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