src/gromacs/mdlib/constr.c

   1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
   2  *
   3  *
   4  *                This source code is part of
   5  *
   6  *                 G   R   O   M   A   C   S
   7  *
   8  *          GROningen MAchine for Chemical Simulations
   9  *
  10  *                        VERSION 3.2.0
  11  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  12  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  13  * Copyright (c) 2001-2004, The GROMACS development team,
  14  * check out http://www.gromacs.org for more information.
  15
  16  * This program is free software; you can redistribute it and/or
  17  * modify it under the terms of the GNU General Public License
  18  * as published by the Free Software Foundation; either version 2
  19  * of the License, or (at your option) any later version.
  20  *
  21  * If you want to redistribute modifications, please consider that
  22  * scientific software is very special. Version control is crucial -
  23  * bugs must be traceable. We will be happy to consider code for
  24  * inclusion in the official distribution, but derived work must not
  25  * be called official GROMACS. Details are found in the README & COPYING
  26  * files - if they are missing, get the official version at www.gromacs.org.
  27  *
  28  * To help us fund GROMACS development, we humbly ask that you cite
  29  * the papers on the package - you can find them in the top README file.
  30  *
  31  * For more info, check our website at http://www.gromacs.org
  32  *
  33  * And Hey:
  34  * GROwing Monsters And Cloning Shrimps
  35  */
  36 #ifdef HAVE_CONFIG_H
  37 #include <config.h>
  38 #endif
  39
  40 #include "confio.h"
  41 #include "constr.h"
  42 #include "copyrite.h"
  43 #include "invblock.h"
  44 #include "main.h"
  45 #include "mdrun.h"
  46 #include "nrnb.h"
  47 #include "smalloc.h"
  48 #include "vec.h"
  49 #include "physics.h"
  50 #include "names.h"
  51 #include "txtdump.h"
  52 #include "domdec.h"
  53 #include "pdbio.h"
  54 #include "partdec.h"
  55 #include "splitter.h"
  56 #include "mtop_util.h"
  57 #include "gmxfio.h"
  58 #include "macros.h"
  59
  60 typedef struct gmx_constr {
  61   int              ncon_tot;     /* The total number of constraints    */
  62   int              nflexcon;     /* The number of flexible constraints */
  63   int              n_at2con_mt;  /* The size of at2con = #moltypes     */
  64   t_blocka         *at2con_mt;   /* A list of atoms to constraints     */
  65   gmx_lincsdata_t  lincsd;       /* LINCS data                         */
  66   gmx_shakedata_t  shaked;       /* SHAKE data                         */
  67   gmx_settledata_t settled;      /* SETTLE data                        */
  68   int              nblocks;      /* The number of SHAKE blocks         */
  69   int              *sblock;      /* The SHAKE blocks                   */
  70   int              sblock_nalloc;/* The allocation size of sblock      */
  71   real             *lagr;        /* Lagrange multipliers for SHAKE     */
  72   int              lagr_nalloc;  /* The allocation size of lagr        */
  73   int              maxwarn;      /* The maximum number of warnings     */
  74   int              warncount_lincs;
  75   int              warncount_settle;
  76   gmx_edsam_t      ed;           /* The essential dynamics data        */
  77
  78   gmx_mtop_t       *warn_mtop;   /* Only used for printing warnings    */
  79 } t_gmx_constr;
  80
  81 typedef struct {
  82   atom_id iatom[3];
  83   atom_id blocknr;
  84 } t_sortblock;
  85
  86 static void *init_vetavars(t_vetavars *vars,
  87                            gmx_bool constr_deriv,
  88                            real veta,real vetanew, t_inputrec *ir, gmx_ekindata_t *ekind, gmx_bool bPscal)
  89 {
  90     double g;
  91     int i;
  92
  93     /* first, set the alpha integrator variable */
  94     if ((ir->opts.nrdf[0] > 0) && bPscal)
  95     {
  96         vars->alpha = 1.0 + DIM/((double)ir->opts.nrdf[0]);
  97     } else {
  98         vars->alpha = 1.0;
  99     }
 100     g = 0.5*veta*ir->delta_t;
 101     vars->rscale = exp(g)*series_sinhx(g);
 102     g = -0.25*vars->alpha*veta*ir->delta_t;
 103     vars->vscale = exp(g)*series_sinhx(g);
 104     vars->rvscale = vars->vscale*vars->rscale;
 105     vars->veta = vetanew;
 106
 107     if (constr_deriv)
 108     {
 109         snew(vars->vscale_nhc,ir->opts.ngtc);
 110         if ((ekind==NULL) || (!bPscal))
 111         {
 112             for (i=0;i<ir->opts.ngtc;i++)
 113             {
 114                 vars->vscale_nhc[i] = 1;
 115             }
 116         }
 117         else
 118         {
 119             for (i=0;i<ir->opts.ngtc;i++)
 120             {
 121                 vars->vscale_nhc[i] = ekind->tcstat[i].vscale_nhc;
 122             }
 123         }
 124     }
 125     else
 126     {
 127         vars->vscale_nhc = NULL;
 128     }
 129
 130     return vars;
 131 }
 132
 133 static void free_vetavars(t_vetavars *vars)
 134 {
 135     if (vars->vscale_nhc != NULL)
 136     {
 137         sfree(vars->vscale_nhc);
 138     }
 139 }
 140
 141 static int pcomp(const void *p1, const void *p2)
 142 {
 143   int     db;
 144   atom_id min1,min2,max1,max2;
 145   t_sortblock *a1=(t_sortblock *)p1;
 146   t_sortblock *a2=(t_sortblock *)p2;
 147
 148   db=a1->blocknr-a2->blocknr;
 149
 150   if (db != 0)
 151     return db;
 152
 153   min1=min(a1->iatom[1],a1->iatom[2]);
 154   max1=max(a1->iatom[1],a1->iatom[2]);
 155   min2=min(a2->iatom[1],a2->iatom[2]);
 156   max2=max(a2->iatom[1],a2->iatom[2]);
 157
 158   if (min1 == min2)
 159     return max1-max2;
 160   else
 161     return min1-min2;
 162 }
 163
 164 static int icomp(const void *p1, const void *p2)
 165 {
 166   atom_id *a1=(atom_id *)p1;
 167   atom_id *a2=(atom_id *)p2;
 168
 169   return (*a1)-(*a2);
 170 }
 171
 172 int n_flexible_constraints(struct gmx_constr *constr)
 173 {
 174   int nflexcon;
 175
 176   if (constr)
 177     nflexcon = constr->nflexcon;
 178   else
 179     nflexcon = 0;
 180
 181   return nflexcon;
 182 }
 183
 184 void too_many_constraint_warnings(int eConstrAlg,int warncount)
 185 {
 186   const char *abort="- aborting to avoid logfile runaway.\n"
 187     "This normally happens when your system is not sufficiently equilibrated,"
 188     "or if you are changing lambda too fast in free energy simulations.\n";
 189
 190   gmx_fatal(FARGS,
 191             "Too many %s warnings (%d)\n"
 192             "If you know what you are doing you can %s"
 193             "set the environment variable GMX_MAXCONSTRWARN to -1,\n"
 194             "but normally it is better to fix the problem",
 195             (eConstrAlg == econtLINCS) ? "LINCS" : "SETTLE",warncount,
 196             (eConstrAlg == econtLINCS) ?
 197             "adjust the lincs warning threshold in your mdp file\nor " : "\n");
 198 }
 199
 200 static void write_constr_pdb(const char *fn,const char *title,
 201                              gmx_mtop_t *mtop,
 202                              int start,int homenr,t_commrec *cr,
 203                              rvec x[],matrix box)
 204 {
 205     char fname[STRLEN],format[STRLEN];
 206     FILE *out;
 207     int  dd_ac0=0,dd_ac1=0,i,ii,resnr;
 208     gmx_domdec_t *dd;
 209     char *anm,*resnm;
 210
 211     dd = NULL;
 212     if (PAR(cr))
 213     {
 214         sprintf(fname,"%s_n%d.pdb",fn,cr->sim_nodeid);
 215         if (DOMAINDECOMP(cr))
 216         {
 217             dd = cr->dd;
 218             dd_get_constraint_range(dd,&dd_ac0,&dd_ac1);
 219             start = 0;
 220             homenr = dd_ac1;
 221         }
 222     }
 223     else
 224     {
 225         sprintf(fname,"%s.pdb",fn);
 226     }
 227     sprintf(format,"%s\n",get_pdbformat());
 228
 229     out = gmx_fio_fopen(fname,"w");
 230
 231     fprintf(out,"TITLE     %s\n",title);
 232     gmx_write_pdb_box(out,-1,box);
 233     for(i=start; i<start+homenr; i++)
 234     {
 235         if (dd != NULL)
 236         {
 237             if (i >= dd->nat_home && i < dd_ac0)
 238             {
 239                 continue;
 240             }
 241             ii = dd->gatindex[i];
 242         }
 243         else
 244         {
 245             ii = i;
 246         }
 247         gmx_mtop_atominfo_global(mtop,ii,&anm,&resnr,&resnm);
 248         fprintf(out,format,"ATOM",(ii+1)%100000,
 249                 anm,resnm,' ',resnr%10000,' ',
 250                 10*x[i][XX],10*x[i][YY],10*x[i][ZZ]);
 251     }
 252     fprintf(out,"TER\n");
 253
 254     gmx_fio_fclose(out);
 255 }
 256
 257 static void dump_confs(FILE *fplog,gmx_large_int_t step,gmx_mtop_t *mtop,
 258                        int start,int homenr,t_commrec *cr,
 259                        rvec x[],rvec xprime[],matrix box)
 260 {
 261   char buf[256],buf2[22];
 262
 263   char *env=getenv("GMX_SUPPRESS_DUMP");
 264   if (env)
 265       return;
 266
 267   sprintf(buf,"step%sb",gmx_step_str(step,buf2));
 268   write_constr_pdb(buf,"initial coordinates",
 269                    mtop,start,homenr,cr,x,box);
 270   sprintf(buf,"step%sc",gmx_step_str(step,buf2));
 271   write_constr_pdb(buf,"coordinates after constraining",
 272                    mtop,start,homenr,cr,xprime,box);
 273   if (fplog)
 274   {
 275       fprintf(fplog,"Wrote pdb files with previous and current coordinates\n");
 276   }
 277   fprintf(stderr,"Wrote pdb files with previous and current coordinates\n");
 278 }
 279
 280 static void pr_sortblock(FILE *fp,const char *title,int nsb,t_sortblock sb[])
 281 {
 282   int i;
 283
 284   fprintf(fp,"%s\n",title);
 285   for(i=0; (i<nsb); i++)
 286     fprintf(fp,"i: %5d, iatom: (%5d %5d %5d), blocknr: %5d\n",
 287             i,sb[i].iatom[0],sb[i].iatom[1],sb[i].iatom[2],
 288             sb[i].blocknr);
 289 }
 290
 291 gmx_bool constrain(FILE *fplog,gmx_bool bLog,gmx_bool bEner,
 292                struct gmx_constr *constr,
 293                t_idef *idef,t_inputrec *ir,gmx_ekindata_t *ekind,
 294                t_commrec *cr,
 295                gmx_large_int_t step,int delta_step,
 296                t_mdatoms *md,
 297                rvec *x,rvec *xprime,rvec *min_proj,matrix box,
 298                real lambda,real *dvdlambda,
 299                rvec *v,tensor *vir,
 300                t_nrnb *nrnb,int econq,gmx_bool bPscal,real veta, real vetanew)
 301 {
 302     gmx_bool    bOK,bDump;
 303     int     start,homenr,nrend;
 304     int     i,j,d;
 305     int     ncons,error;
 306     tensor  rmdr;
 307     rvec    *vstor;
 308     real    invdt,vir_fac,t;
 309     t_ilist *settle;
 310     int     nsettle;
 311     t_pbc   pbc;
 312     char    buf[22];
 313     t_vetavars vetavar;
 314
 315     if (econq == econqForceDispl && !EI_ENERGY_MINIMIZATION(ir->eI))
 316     {
 317         gmx_incons("constrain called for forces displacements while not doing energy minimization, can not do this while the LINCS and SETTLE constraint connection matrices are mass weighted");
 318     }
 319
 320     bOK   = TRUE;
 321     bDump = FALSE;
 322
 323     start  = md->start;
 324     homenr = md->homenr;
 325     nrend = start+homenr;
 326
 327     /* set constants for pressure control integration */
 328     init_vetavars(&vetavar,econq!=econqCoord,
 329                   veta,vetanew,ir,ekind,bPscal);
 330
 331     if (ir->delta_t == 0)
 332     {
 333         invdt = 0;
 334     }
 335     else
 336     {
 337         invdt  = 1/ir->delta_t;
 338     }
 339
 340     if (ir->efep != efepNO && EI_DYNAMICS(ir->eI))
 341     {
 342         /* Set the constraint lengths for the step at which this configuration
 343          * is meant to be. The invmasses should not be changed.
 344          */
 345         lambda += delta_step*ir->fepvals->delta_lambda;
 346     }
 347
 348     if (vir != NULL)
 349     {
 350         clear_mat(rmdr);
 351     }
 352
 353     where();
 354     if (constr->lincsd)
 355     {
 356         bOK = constrain_lincs(fplog,bLog,bEner,ir,step,constr->lincsd,md,cr,
 357                               x,xprime,min_proj,box,lambda,dvdlambda,
 358                               invdt,v,vir!=NULL,rmdr,
 359                               econq,nrnb,
 360                               constr->maxwarn,&constr->warncount_lincs);
 361         if (!bOK && constr->maxwarn >= 0)
 362         {
 363             if (fplog != NULL)
 364             {
 365                 fprintf(fplog,"Constraint error in algorithm %s at step %s\n",
 366                         econstr_names[econtLINCS],gmx_step_str(step,buf));
 367             }
 368             bDump = TRUE;
 369         }
 370     }
 371
 372     if (constr->nblocks > 0)
 373     {
 374         switch (econq) {
 375         case (econqCoord):
 376             bOK = bshakef(fplog,constr->shaked,
 377                           homenr,md->invmass,constr->nblocks,constr->sblock,
 378                           idef,ir,box,x,xprime,nrnb,
 379                           constr->lagr,lambda,dvdlambda,
 380                           invdt,v,vir!=NULL,rmdr,constr->maxwarn>=0,econq,
 381                           &vetavar);
 382             break;
 383         case (econqVeloc):
 384             bOK = bshakef(fplog,constr->shaked,
 385                           homenr,md->invmass,constr->nblocks,constr->sblock,
 386                           idef,ir,box,x,min_proj,nrnb,
 387                           constr->lagr,lambda,dvdlambda,
 388                           invdt,NULL,vir!=NULL,rmdr,constr->maxwarn>=0,econq,
 389                           &vetavar);
 390             break;
 391         default:
 392             gmx_fatal(FARGS,"Internal error, SHAKE called for constraining something else than coordinates");
 393             break;
 394         }
 395
 396         if (!bOK && constr->maxwarn >= 0)
 397         {
 398             if (fplog != NULL)
 399             {
 400                 fprintf(fplog,"Constraint error in algorithm %s at step %s\n",
 401                         econstr_names[econtSHAKE],gmx_step_str(step,buf));
 402             }
 403             bDump = TRUE;
 404         }
 405     }
 406
 407     settle  = &idef->il[F_SETTLE];
 408     if (settle->nr > 0)
 409     {
 410         nsettle = settle->nr/4;
 411
 412         switch (econq)
 413         {
 414         case econqCoord:
 415             csettle(constr->settled,
 416                     nsettle,settle->iatoms,x[0],xprime[0],
 417                     invdt,v[0],vir!=NULL,rmdr,&error,&vetavar);
 418             inc_nrnb(nrnb,eNR_SETTLE,nsettle);
 419             if (v != NULL)
 420             {
 421                 inc_nrnb(nrnb,eNR_CONSTR_V,nsettle*3);
 422             }
 423             if (vir != NULL)
 424             {
 425                 inc_nrnb(nrnb,eNR_CONSTR_VIR,nsettle*3);
 426             }
 427
 428             bOK = (error < 0);
 429             if (!bOK && constr->maxwarn >= 0)
 430             {
 431                 char buf[256];
 432                 sprintf(buf,
 433                         "\nstep " gmx_large_int_pfmt ": Water molecule starting at atom %d can not be "
 434                         "settled.\nCheck for bad contacts and/or reduce the timestep if appropriate.\n",
 435                         step,ddglatnr(cr->dd,settle->iatoms[error*4+1]));
 436                 if (fplog)
 437                 {
 438                     fprintf(fplog,"%s",buf);
 439                 }
 440                 fprintf(stderr,"%s",buf);
 441                 constr->warncount_settle++;
 442                 if (constr->warncount_settle > constr->maxwarn)
 443                 {
 444                     too_many_constraint_warnings(-1,constr->warncount_settle);
 445                 }
 446                 bDump = TRUE;
 447                 break;
 448             case econqVeloc:
 449             case econqDeriv:
 450             case econqForce:
 451             case econqForceDispl:
 452                 settle_proj(fplog,constr->settled,econq,
 453                             nsettle,settle->iatoms,x,
 454                             xprime,min_proj,vir!=NULL,rmdr,&vetavar);
 455                 /* This is an overestimate */
 456                 inc_nrnb(nrnb,eNR_SETTLE,nsettle);
 457                 break;
 458             case econqDeriv_FlexCon:
 459                 /* Nothing to do, since the are no flexible constraints in settles */
 460                 break;
 461             default:
 462                 gmx_incons("Unknown constraint quantity for settle");
 463             }
 464         }
 465     }
 466
 467     free_vetavars(&vetavar);
 468
 469     if (vir != NULL)
 470     {
 471         switch (econq)
 472         {
 473         case econqCoord:
 474             vir_fac = 0.5/(ir->delta_t*ir->delta_t);
 475             break;
 476         case econqVeloc:
 477             vir_fac = 0.5/ir->delta_t;
 478             break;
 479         case econqForce:
 480         case econqForceDispl:
 481             vir_fac = 0.5;
 482             break;
 483         default:
 484             vir_fac = 0;
 485             gmx_incons("Unsupported constraint quantity for virial");
 486         }
 487
 488         if (EI_VV(ir->eI))
 489         {
 490             vir_fac *= 2;  /* only constraining over half the distance here */
 491         }
 492         for(i=0; i<DIM; i++)
 493         {
 494             for(j=0; j<DIM; j++)
 495             {
 496                 (*vir)[i][j] = vir_fac*rmdr[i][j];
 497             }
 498         }
 499     }
 500
 501     if (bDump)
 502     {
 503         dump_confs(fplog,step,constr->warn_mtop,start,homenr,cr,x,xprime,box);
 504     }
 505
 506     if (econq == econqCoord)
 507     {
 508         if (ir->ePull == epullCONSTRAINT)
 509         {
 510             if (EI_DYNAMICS(ir->eI))
 511             {
 512                 t = ir->init_t + (step + delta_step)*ir->delta_t;
 513             }
 514             else
 515             {
 516                 t = ir->init_t;
 517             }
 518             set_pbc(&pbc,ir->ePBC,box);
 519             pull_constraint(ir->pull,md,&pbc,cr,ir->delta_t,t,x,xprime,v,*vir);
 520         }
 521         if (constr->ed && delta_step > 0)
 522         {
 523             /* apply the essential dynamcs constraints here */
 524             do_edsam(ir,step,md,cr,xprime,v,box,constr->ed);
 525         }
 526     }
 527
 528     return bOK;
 529 }
 530
 531 real *constr_rmsd_data(struct gmx_constr *constr)
 532 {
 533   if (constr->lincsd)
 534     return lincs_rmsd_data(constr->lincsd);
 535   else
 536     return NULL;
 537 }
 538
 539 real constr_rmsd(struct gmx_constr *constr,gmx_bool bSD2)
 540 {
 541   if (constr->lincsd)
 542     return lincs_rmsd(constr->lincsd,bSD2);
 543   else
 544     return 0;
 545 }
 546
 547 static void make_shake_sblock_pd(struct gmx_constr *constr,
 548                                  t_idef *idef,t_mdatoms *md)
 549 {
 550   int  i,j,m,ncons;
 551   int  bstart,bnr;
 552   t_blocka    sblocks;
 553   t_sortblock *sb;
 554   t_iatom     *iatom;
 555   atom_id     *inv_sblock;
 556
 557   /* Since we are processing the local topology,
 558    * the F_CONSTRNC ilist has been concatenated to the F_CONSTR ilist.
 559    */
 560   ncons = idef->il[F_CONSTR].nr/3;
 561
 562   init_blocka(&sblocks);
 563   gen_sblocks(NULL,md->start,md->start+md->homenr,idef,&sblocks,FALSE);
 564
 565   /*
 566     bstart=(idef->nodeid > 0) ? blocks->multinr[idef->nodeid-1] : 0;
 567     nblocks=blocks->multinr[idef->nodeid] - bstart;
 568   */
 569   bstart  = 0;
 570   constr->nblocks = sblocks.nr;
 571   if (debug)
 572     fprintf(debug,"ncons: %d, bstart: %d, nblocks: %d\n",
 573             ncons,bstart,constr->nblocks);
 574
 575   /* Calculate block number for each atom */
 576   inv_sblock = make_invblocka(&sblocks,md->nr);
 577
 578   done_blocka(&sblocks);
 579
 580   /* Store the block number in temp array and
 581    * sort the constraints in order of the sblock number
 582    * and the atom numbers, really sorting a segment of the array!
 583    */
 584 #ifdef DEBUGIDEF
 585   pr_idef(fplog,0,"Before Sort",idef);
 586 #endif
 587   iatom=idef->il[F_CONSTR].iatoms;
 588   snew(sb,ncons);
 589   for(i=0; (i<ncons); i++,iatom+=3) {
 590     for(m=0; (m<3); m++)
 591       sb[i].iatom[m] = iatom[m];
 592     sb[i].blocknr = inv_sblock[iatom[1]];
 593   }
 594
 595   /* Now sort the blocks */
 596   if (debug) {
 597     pr_sortblock(debug,"Before sorting",ncons,sb);
 598     fprintf(debug,"Going to sort constraints\n");
 599   }
 600
 601   qsort(sb,ncons,(size_t)sizeof(*sb),pcomp);
 602
 603   if (debug) {
 604     pr_sortblock(debug,"After sorting",ncons,sb);
 605   }
 606
 607   iatom=idef->il[F_CONSTR].iatoms;
 608   for(i=0; (i<ncons); i++,iatom+=3)
 609     for(m=0; (m<3); m++)
 610       iatom[m]=sb[i].iatom[m];
 611 #ifdef DEBUGIDEF
 612   pr_idef(fplog,0,"After Sort",idef);
 613 #endif
 614
 615   j=0;
 616   snew(constr->sblock,constr->nblocks+1);
 617   bnr=-2;
 618   for(i=0; (i<ncons); i++) {
 619     if (sb[i].blocknr != bnr) {
 620       bnr=sb[i].blocknr;
 621       constr->sblock[j++]=3*i;
 622     }
 623   }
 624   /* Last block... */
 625   constr->sblock[j++] = 3*ncons;
 626
 627   if (j != (constr->nblocks+1)) {
 628     fprintf(stderr,"bstart: %d\n",bstart);
 629     fprintf(stderr,"j: %d, nblocks: %d, ncons: %d\n",
 630             j,constr->nblocks,ncons);
 631     for(i=0; (i<ncons); i++)
 632       fprintf(stderr,"i: %5d  sb[i].blocknr: %5u\n",i,sb[i].blocknr);
 633     for(j=0; (j<=constr->nblocks); j++)
 634       fprintf(stderr,"sblock[%3d]=%5d\n",j,(int)constr->sblock[j]);
 635     gmx_fatal(FARGS,"DEATH HORROR: "
 636               "sblocks does not match idef->il[F_CONSTR]");
 637   }
 638   sfree(sb);
 639   sfree(inv_sblock);
 640 }
 641
 642 static void make_shake_sblock_dd(struct gmx_constr *constr,
 643                                  t_ilist *ilcon,t_block *cgs,
 644                                  gmx_domdec_t *dd)
 645 {
 646   int ncons,c,cg;
 647   t_iatom *iatom;
 648
 649   if (dd->ncg_home+1 > constr->sblock_nalloc) {
 650     constr->sblock_nalloc = over_alloc_dd(dd->ncg_home+1);
 651     srenew(constr->sblock,constr->sblock_nalloc);
 652   }
 653
 654   ncons = ilcon->nr/3;
 655   iatom = ilcon->iatoms;
 656   constr->nblocks = 0;
 657   cg = 0;
 658   for(c=0; c<ncons; c++) {
 659     if (c == 0 || iatom[1] >= cgs->index[cg+1]) {
 660       constr->sblock[constr->nblocks++] = 3*c;
 661       while (iatom[1] >= cgs->index[cg+1])
 662         cg++;
 663     }
 664     iatom += 3;
 665   }
 666   constr->sblock[constr->nblocks] = 3*ncons;
 667 }
 668
 669 t_blocka make_at2con(int start,int natoms,
 670                      t_ilist *ilist,t_iparams *iparams,
 671                      gmx_bool bDynamics,int *nflexiblecons)
 672 {
 673   int *count,ncon,con,con_tot,nflexcon,ftype,i,a;
 674   t_iatom  *ia;
 675   t_blocka at2con;
 676   gmx_bool bFlexCon;
 677
 678   snew(count,natoms);
 679   nflexcon = 0;
 680   for(ftype=F_CONSTR; ftype<=F_CONSTRNC; ftype++) {
 681     ncon = ilist[ftype].nr/3;
 682     ia   = ilist[ftype].iatoms;
 683     for(con=0; con<ncon; con++) {
 684       bFlexCon = (iparams[ia[0]].constr.dA == 0 &&
 685                   iparams[ia[0]].constr.dB == 0);
 686       if (bFlexCon)
 687         nflexcon++;
 688       if (bDynamics || !bFlexCon) {
 689         for(i=1; i<3; i++) {
 690           a = ia[i] - start;
 691           count[a]++;
 692         }
 693       }
 694       ia += 3;
 695     }
 696   }
 697   *nflexiblecons = nflexcon;
 698
 699   at2con.nr = natoms;
 700   at2con.nalloc_index = at2con.nr+1;
 701   snew(at2con.index,at2con.nalloc_index);
 702   at2con.index[0] = 0;
 703   for(a=0; a<natoms; a++) {
 704     at2con.index[a+1] = at2con.index[a] + count[a];
 705     count[a] = 0;
 706   }
 707   at2con.nra = at2con.index[natoms];
 708   at2con.nalloc_a = at2con.nra;
 709   snew(at2con.a,at2con.nalloc_a);
 710
 711   /* The F_CONSTRNC constraints have constraint numbers
 712    * that continue after the last F_CONSTR constraint.
 713    */
 714   con_tot = 0;
 715   for(ftype=F_CONSTR; ftype<=F_CONSTRNC; ftype++) {
 716     ncon = ilist[ftype].nr/3;
 717     ia   = ilist[ftype].iatoms;
 718     for(con=0; con<ncon; con++) {
 719       bFlexCon = (iparams[ia[0]].constr.dA == 0 &&
 720                   iparams[ia[0]].constr.dB == 0);
 721       if (bDynamics || !bFlexCon) {
 722         for(i=1; i<3; i++) {
 723           a = ia[i] - start;
 724           at2con.a[at2con.index[a]+count[a]++] = con_tot;
 725         }
 726       }
 727       con_tot++;
 728       ia += 3;
 729     }
 730   }
 731
 732   sfree(count);
 733
 734   return at2con;
 735 }
 736
 737 void set_constraints(struct gmx_constr *constr,
 738                      gmx_localtop_t *top,t_inputrec *ir,
 739                      t_mdatoms *md,t_commrec *cr)
 740 {
 741     t_idef *idef;
 742     int    ncons;
 743     t_ilist *settle;
 744     int    iO,iH;
 745
 746     idef = &top->idef;
 747
 748     if (constr->ncon_tot > 0)
 749     {
 750         /* We are using the local topology,
 751          * so there are only F_CONSTR constraints.
 752          */
 753         ncons = idef->il[F_CONSTR].nr/3;
 754
 755         /* With DD we might also need to call LINCS with ncons=0 for
 756          * communicating coordinates to other nodes that do have constraints.
 757          */
 758         if (ir->eConstrAlg == econtLINCS)
 759         {
 760             set_lincs(idef,md,EI_DYNAMICS(ir->eI),cr,constr->lincsd);
 761         }
 762         if (ir->eConstrAlg == econtSHAKE)
 763         {
 764             if (cr->dd)
 765             {
 766                 make_shake_sblock_dd(constr,&idef->il[F_CONSTR],&top->cgs,cr->dd);
 767             }
 768             else
 769             {
 770                 make_shake_sblock_pd(constr,idef,md);
 771             }
 772             if (ncons > constr->lagr_nalloc)
 773             {
 774                 constr->lagr_nalloc = over_alloc_dd(ncons);
 775                 srenew(constr->lagr,constr->lagr_nalloc);
 776             }
 777         }
 778     }
 779
 780     if (idef->il[F_SETTLE].nr > 0 && constr->settled == NULL)
 781     {
 782         settle = &idef->il[F_SETTLE];
 783         iO = settle->iatoms[1];
 784         iH = settle->iatoms[2];
 785         constr->settled =
 786             settle_init(md->massT[iO],md->massT[iH],
 787                         md->invmass[iO],md->invmass[iH],
 788                         idef->iparams[settle->iatoms[0]].settle.doh,
 789                         idef->iparams[settle->iatoms[0]].settle.dhh);
 790     }
 791
 792     /* Make a selection of the local atoms for essential dynamics */
 793     if (constr->ed && cr->dd)
 794     {
 795         dd_make_local_ed_indices(cr->dd,constr->ed);
 796     }
 797 }
 798
 799 static void constr_recur(t_blocka *at2con,
 800                          t_ilist *ilist,t_iparams *iparams,gmx_bool bTopB,
 801                          int at,int depth,int nc,int *path,
 802                          real r0,real r1,real *r2max,
 803                          int *count)
 804 {
 805   int  ncon1;
 806   t_iatom *ia1,*ia2;
 807   int  c,con,a1;
 808   gmx_bool bUse;
 809   t_iatom *ia;
 810   real len,rn0,rn1;
 811
 812   (*count)++;
 813
 814   ncon1 = ilist[F_CONSTR].nr/3;
 815   ia1   = ilist[F_CONSTR].iatoms;
 816   ia2   = ilist[F_CONSTRNC].iatoms;
 817
 818   /* Loop over all constraints connected to this atom */
 819   for(c=at2con->index[at]; c<at2con->index[at+1]; c++) {
 820     con = at2con->a[c];
 821     /* Do not walk over already used constraints */
 822     bUse = TRUE;
 823     for(a1=0; a1<depth; a1++) {
 824       if (con == path[a1])
 825         bUse = FALSE;
 826     }
 827     if (bUse) {
 828       ia = constr_iatomptr(ncon1,ia1,ia2,con);
 829       /* Flexible constraints currently have length 0, which is incorrect */
 830       if (!bTopB)
 831         len = iparams[ia[0]].constr.dA;
 832       else
 833         len = iparams[ia[0]].constr.dB;
 834       /* In the worst case the bond directions alternate */
 835       if (nc % 2 == 0) {
 836         rn0 = r0 + len;
 837         rn1 = r1;
 838       } else {
 839         rn0 = r0;
 840         rn1 = r1 + len;
 841       }
 842       /* Assume angles of 120 degrees between all bonds */
 843       if (rn0*rn0 + rn1*rn1 + rn0*rn1 > *r2max) {
 844         *r2max = rn0*rn0 + rn1*rn1 + r0*rn1;
 845         if (debug) {
 846           fprintf(debug,"Found longer constraint distance: r0 %5.3f r1 %5.3f rmax %5.3f\n", rn0,rn1,sqrt(*r2max));
 847           for(a1=0; a1<depth; a1++)
 848             fprintf(debug," %d %5.3f",
 849                     path[a1],
 850                     iparams[constr_iatomptr(ncon1,ia1,ia2,con)[0]].constr.dA);
 851           fprintf(debug," %d %5.3f\n",con,len);
 852         }
 853       }
 854       /* Limit the number of recursions to 1000*nc,
 855        * so a call does not take more than a second,
 856        * even for highly connected systems.
 857        */
 858       if (depth + 1 < nc && *count < 1000*nc) {
 859         if (ia[1] == at)
 860           a1 = ia[2];
 861         else
 862           a1 = ia[1];
 863         /* Recursion */
 864         path[depth] = con;
 865         constr_recur(at2con,ilist,iparams,
 866                      bTopB,a1,depth+1,nc,path,rn0,rn1,r2max,count);
 867         path[depth] = -1;
 868       }
 869     }
 870   }
 871 }
 872
 873 static real constr_r_max_moltype(FILE *fplog,
 874                                  gmx_moltype_t *molt,t_iparams *iparams,
 875                                  t_inputrec *ir)
 876 {
 877   int natoms,nflexcon,*path,at,count;
 878
 879   t_blocka at2con;
 880   real r0,r1,r2maxA,r2maxB,rmax,lam0,lam1;
 881
 882   if (molt->ilist[F_CONSTR].nr   == 0 &&
 883       molt->ilist[F_CONSTRNC].nr == 0) {
 884     return 0;
 885   }
 886
 887   natoms = molt->atoms.nr;
 888
 889   at2con = make_at2con(0,natoms,molt->ilist,iparams,
 890                        EI_DYNAMICS(ir->eI),&nflexcon);
 891   snew(path,1+ir->nProjOrder);
 892   for(at=0; at<1+ir->nProjOrder; at++)
 893     path[at] = -1;
 894
 895   r2maxA = 0;
 896   for(at=0; at<natoms; at++) {
 897     r0 = 0;
 898     r1 = 0;
 899
 900     count = 0;
 901     constr_recur(&at2con,molt->ilist,iparams,
 902                  FALSE,at,0,1+ir->nProjOrder,path,r0,r1,&r2maxA,&count);
 903   }
 904   if (ir->efep == efepNO) {
 905     rmax = sqrt(r2maxA);
 906   } else {
 907     r2maxB = 0;
 908     for(at=0; at<natoms; at++) {
 909       r0 = 0;
 910       r1 = 0;
 911       count = 0;
 912       constr_recur(&at2con,molt->ilist,iparams,
 913                    TRUE,at,0,1+ir->nProjOrder,path,r0,r1,&r2maxB,&count);
 914     }
 915     lam0 = ir->fepvals->init_lambda;
 916     if (EI_DYNAMICS(ir->eI))
 917       lam0 += ir->init_step*ir->fepvals->delta_lambda;
 918     rmax = (1 - lam0)*sqrt(r2maxA) + lam0*sqrt(r2maxB);
 919     if (EI_DYNAMICS(ir->eI)) {
 920       lam1 = ir->fepvals->init_lambda + (ir->init_step + ir->nsteps)*ir->fepvals->delta_lambda;
 921       rmax = max(rmax,(1 - lam1)*sqrt(r2maxA) + lam1*sqrt(r2maxB));
 922     }
 923   }
 924
 925   done_blocka(&at2con);
 926   sfree(path);
 927
 928   return rmax;
 929 }
 930
 931 real constr_r_max(FILE *fplog,gmx_mtop_t *mtop,t_inputrec *ir)
 932 {
 933   int mt;
 934   real rmax;
 935
 936   rmax = 0;
 937   for(mt=0; mt<mtop->nmoltype; mt++) {
 938     rmax = max(rmax,
 939                constr_r_max_moltype(fplog,&mtop->moltype[mt],
 940                                     mtop->ffparams.iparams,ir));
 941   }
 942
 943   if (fplog)
 944     fprintf(fplog,"Maximum distance for %d constraints, at 120 deg. angles, all-trans: %.3f nm\n",1+ir->nProjOrder,rmax);
 945
 946   return rmax;
 947 }
 948
 949 gmx_constr_t init_constraints(FILE *fplog,
 950                               gmx_mtop_t *mtop,t_inputrec *ir,
 951                               gmx_edsam_t ed,t_state *state,
 952                               t_commrec *cr)
 953 {
 954     int  ncon,nset,nmol,settle_type,i,natoms,mt,nflexcon;
 955     struct gmx_constr *constr;
 956     char *env;
 957     t_ilist *ilist;
 958     gmx_mtop_ilistloop_t iloop;
 959
 960     ncon =
 961         gmx_mtop_ftype_count(mtop,F_CONSTR) +
 962         gmx_mtop_ftype_count(mtop,F_CONSTRNC);
 963     nset = gmx_mtop_ftype_count(mtop,F_SETTLE);
 964
 965     if (ncon+nset == 0 && ir->ePull != epullCONSTRAINT && ed == NULL)
 966     {
 967         return NULL;
 968     }
 969
 970     snew(constr,1);
 971
 972     constr->ncon_tot = ncon;
 973     constr->nflexcon = 0;
 974     if (ncon > 0)
 975     {
 976         constr->n_at2con_mt = mtop->nmoltype;
 977         snew(constr->at2con_mt,constr->n_at2con_mt);
 978         for(mt=0; mt<mtop->nmoltype; mt++)
 979         {
 980             constr->at2con_mt[mt] = make_at2con(0,mtop->moltype[mt].atoms.nr,
 981                                                 mtop->moltype[mt].ilist,
 982                                                 mtop->ffparams.iparams,
 983                                                 EI_DYNAMICS(ir->eI),&nflexcon);
 984             for(i=0; i<mtop->nmolblock; i++)
 985             {
 986                 if (mtop->molblock[i].type == mt)
 987                 {
 988                     constr->nflexcon += mtop->molblock[i].nmol*nflexcon;
 989                 }
 990             }
 991         }
 992
 993         if (constr->nflexcon > 0)
 994         {
 995             if (fplog)
 996             {
 997                 fprintf(fplog,"There are %d flexible constraints\n",
 998                         constr->nflexcon);
 999                 if (ir->fc_stepsize == 0)
1000                 {
1001                     fprintf(fplog,"\n"
1002                             "WARNING: step size for flexible constraining = 0\n"
1003                             "         All flexible constraints will be rigid.\n"
1004                             "         Will try to keep all flexible constraints at their original length,\n"
1005                             "         but the lengths may exhibit some drift.\n\n");
1006                     constr->nflexcon = 0;
1007                 }
1008             }
1009             if (constr->nflexcon > 0)
1010             {
1011                 please_cite(fplog,"Hess2002");
1012             }
1013         }
1014
1015         if (ir->eConstrAlg == econtLINCS)
1016         {
1017             constr->lincsd = init_lincs(fplog,mtop,
1018                                         constr->nflexcon,constr->at2con_mt,
1019                                         DOMAINDECOMP(cr) && cr->dd->bInterCGcons,
1020                                         ir->nLincsIter,ir->nProjOrder);
1021         }
1022
1023         if (ir->eConstrAlg == econtSHAKE) {
1024             if (DOMAINDECOMP(cr) && cr->dd->bInterCGcons)
1025             {
1026                 gmx_fatal(FARGS,"SHAKE is not supported with domain decomposition and constraint that cross charge group boundaries, use LINCS");
1027             }
1028             if (constr->nflexcon)
1029             {
1030                 gmx_fatal(FARGS,"For this system also velocities and/or forces need to be constrained, this can not be done with SHAKE, you should select LINCS");
1031             }
1032             please_cite(fplog,"Ryckaert77a");
1033             if (ir->bShakeSOR)
1034             {
1035                 please_cite(fplog,"Barth95a");
1036             }
1037
1038             constr->shaked = shake_init();
1039         }
1040     }
1041
1042     if (nset > 0) {
1043         please_cite(fplog,"Miyamoto92a");
1044
1045         /* Check that we have only one settle type */
1046         settle_type = -1;
1047         iloop = gmx_mtop_ilistloop_init(mtop);
1048         while (gmx_mtop_ilistloop_next(iloop,&ilist,&nmol))
1049         {
1050             for (i=0; i<ilist[F_SETTLE].nr; i+=4)
1051             {
1052                 if (settle_type == -1)
1053                 {
1054                     settle_type = ilist[F_SETTLE].iatoms[i];
1055                 }
1056                 else if (ilist[F_SETTLE].iatoms[i] != settle_type)
1057                 {
1058                     gmx_fatal(FARGS,
1059                               "The [molecules] section of your topology specifies more than one block of\n"
1060                               "a [moleculetype] with a [settles] block. Only one such is allowed. If you\n"
1061                               "are trying to partition your solvent into different *groups* (e.g. for\n"
1062                               "freezing, T-coupling, etc.) then you are using the wrong approach. Index\n"
1063                               "files specify groups. Otherwise, you may wish to change the least-used\n"
1064                               "block of molecules with SETTLE constraints into 3 normal constraints.");
1065                 }
1066             }
1067         }
1068     }
1069
1070     constr->maxwarn = 999;
1071     env = getenv("GMX_MAXCONSTRWARN");
1072     if (env)
1073     {
1074         constr->maxwarn = 0;
1075         sscanf(env,"%d",&constr->maxwarn);
1076         if (fplog)
1077         {
1078             fprintf(fplog,
1079                     "Setting the maximum number of constraint warnings to %d\n",
1080                     constr->maxwarn);
1081         }
1082         if (MASTER(cr))
1083         {
1084             fprintf(stderr,
1085                     "Setting the maximum number of constraint warnings to %d\n",
1086                     constr->maxwarn);
1087         }
1088     }
1089     if (constr->maxwarn < 0 && fplog)
1090     {
1091         fprintf(fplog,"maxwarn < 0, will not stop on constraint errors\n");
1092     }
1093     constr->warncount_lincs  = 0;
1094     constr->warncount_settle = 0;
1095
1096     /* Initialize the essential dynamics sampling.
1097      * Put the pointer to the ED struct in constr */
1098     constr->ed = ed;
1099     if (ed != NULL)
1100     {
1101         init_edsam(mtop,ir,cr,ed,state->x,state->box);
1102     }
1103
1104     constr->warn_mtop = mtop;
1105
1106     return constr;
1107 }
1108
1109 t_blocka *atom2constraints_moltype(gmx_constr_t constr)
1110 {
1111   return constr->at2con_mt;
1112 }
1113
1114
1115 gmx_bool inter_charge_group_constraints(gmx_mtop_t *mtop)
1116 {
1117   const gmx_moltype_t *molt;
1118   const t_block *cgs;
1119   const t_ilist *il;
1120   int  mb;
1121   int  nat,*at2cg,cg,a,ftype,i;
1122   gmx_bool bInterCG;
1123
1124   bInterCG = FALSE;
1125   for(mb=0; mb<mtop->nmolblock && !bInterCG; mb++) {
1126     molt = &mtop->moltype[mtop->molblock[mb].type];
1127
1128     if (molt->ilist[F_CONSTR].nr   > 0 ||
1129         molt->ilist[F_CONSTRNC].nr > 0) {
1130       cgs  = &molt->cgs;
1131       snew(at2cg,molt->atoms.nr);
1132       for(cg=0; cg<cgs->nr; cg++) {
1133         for(a=cgs->index[cg]; a<cgs->index[cg+1]; a++)
1134           at2cg[a] = cg;
1135       }
1136
1137       for(ftype=F_CONSTR; ftype<=F_CONSTRNC; ftype++) {
1138         il = &molt->ilist[ftype];
1139         for(i=0; i<il->nr && !bInterCG; i+=3) {
1140           if (at2cg[il->iatoms[i+1]] != at2cg[il->iatoms[i+2]])
1141             bInterCG = TRUE;
1142         }
1143       }
1144       sfree(at2cg);
1145     }
1146   }
1147
1148   return bInterCG;
1149 }
1150
1151 /* helper functions for andersen temperature control, because the
1152  * gmx_constr construct is only defined in constr.c. Return the list
1153  * of blocks (get_sblock) and the number of blocks (get_nblocks).  */
1154
1155 extern int *get_sblock(struct gmx_constr *constr)
1156 {
1157     return constr->sblock;
1158 }
1159
1160 extern int get_nblocks(struct gmx_constr *constr)
1161 {
1162     return constr->nblocks;
1163 }