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