src/mdlib/constr.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2004, The GROMACS development team,
   6  * check out http://www.gromacs.org for more information.
   7  * Copyright (c) 2012, by the GROMACS development team, led by
   8  * David van der Spoel, Berk Hess, Erik Lindahl, and including many
   9  * others, as listed in the AUTHORS file in the top-level source
  10  * directory and at http://www.gromacs.org.
  11  *
  12  * GROMACS is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU Lesser General Public License
  14  * as published by the Free Software Foundation; either version 2.1
  15  * of the License, or (at your option) any later version.
  16  *
  17  * GROMACS is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20  * Lesser General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Lesser General Public
  23  * License along with GROMACS; if not, see
  24  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  25  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  26  *
  27  * If you want to redistribute modifications to GROMACS, please
  28  * consider that scientific software is very special. Version
  29  * control is crucial - bugs must be traceable. We will be happy to
  30  * consider code for inclusion in the official distribution, but
  31  * derived work must not be called official GROMACS. Details are found
  32  * in the README & COPYING files - if they are missing, get the
  33  * official version at http://www.gromacs.org.
  34  *
  35  * To help us fund GROMACS development, we humbly ask that you cite
  36  * the research papers on the package. Check out http://www.gromacs.org.
  37  */
  38 #ifdef HAVE_CONFIG_H
  39 #include <config.h>
  40 #endif
  41
  42 #include "confio.h"
  43 #include "constr.h"
  44 #include "copyrite.h"
  45 #include "invblock.h"
  46 #include "main.h"
  47 #include "mdrun.h"
  48 #include "nrnb.h"
  49 #include "smalloc.h"
  50 #include "vec.h"
  51 #include "physics.h"
  52 #include "names.h"
  53 #include "txtdump.h"
  54 #include "domdec.h"
  55 #include "pdbio.h"
  56 #include "partdec.h"
  57 #include "splitter.h"
  58 #include "mtop_util.h"
  59 #include "gmxfio.h"
  60 #include "gmx_omp_nthreads.h"
  61
  62 typedef struct gmx_constr {
  63   int              ncon_tot;     /* The total number of constraints    */
  64   int              nflexcon;     /* The number of flexible constraints */
  65   int              n_at2con_mt;  /* The size of at2con = #moltypes     */
  66   t_blocka         *at2con_mt;   /* A list of atoms to constraints     */
  67   int              n_at2settle_mt; /* The size of at2settle = #moltypes  */
  68   int              **at2settle_mt; /* A list of atoms to settles         */
  69   gmx_bool         bInterCGsettles;
  70   gmx_lincsdata_t  lincsd;       /* LINCS data                         */
  71   gmx_shakedata_t  shaked;       /* SHAKE data                         */
  72   gmx_settledata_t settled;      /* SETTLE data                        */
  73   int              nblocks;      /* The number of SHAKE blocks         */
  74   int              *sblock;      /* The SHAKE blocks                   */
  75   int              sblock_nalloc;/* The allocation size of sblock      */
  76   real             *lagr;        /* Lagrange multipliers for SHAKE     */
  77   int              lagr_nalloc;  /* The allocation size of lagr        */
  78   int              maxwarn;      /* The maximum number of warnings     */
  79   int              warncount_lincs;
  80   int              warncount_settle;
  81   gmx_edsam_t      ed;           /* The essential dynamics data        */
  82
  83     tensor           *vir_r_m_dr_th;/* Thread local working data          */
  84     int              *settle_error; /* Thread local working data          */
  85
  86   gmx_mtop_t       *warn_mtop;   /* Only used for printing warnings    */
  87 } t_gmx_constr;
  88
  89 typedef struct {
  90   atom_id iatom[3];
  91   atom_id blocknr;
  92 } t_sortblock;
  93
  94 static void *init_vetavars(t_vetavars *vars,
  95                            gmx_bool constr_deriv,
  96                            real veta,real vetanew, t_inputrec *ir, gmx_ekindata_t *ekind, gmx_bool bPscal)
  97 {
  98     double g;
  99     int i;
 100
 101     /* first, set the alpha integrator variable */
 102     if ((ir->opts.nrdf[0] > 0) && bPscal)
 103     {
 104         vars->alpha = 1.0 + DIM/((double)ir->opts.nrdf[0]);
 105     } else {
 106         vars->alpha = 1.0;
 107     }
 108     g = 0.5*veta*ir->delta_t;
 109     vars->rscale = exp(g)*series_sinhx(g);
 110     g = -0.25*vars->alpha*veta*ir->delta_t;
 111     vars->vscale = exp(g)*series_sinhx(g);
 112     vars->rvscale = vars->vscale*vars->rscale;
 113     vars->veta = vetanew;
 114
 115     if (constr_deriv)
 116     {
 117         snew(vars->vscale_nhc,ir->opts.ngtc);
 118         if ((ekind==NULL) || (!bPscal))
 119         {
 120             for (i=0;i<ir->opts.ngtc;i++)
 121             {
 122                 vars->vscale_nhc[i] = 1;
 123             }
 124         }
 125         else
 126         {
 127             for (i=0;i<ir->opts.ngtc;i++)
 128             {
 129                 vars->vscale_nhc[i] = ekind->tcstat[i].vscale_nhc;
 130             }
 131         }
 132     }
 133     else
 134     {
 135         vars->vscale_nhc = NULL;
 136     }
 137
 138     return vars;
 139 }
 140
 141 static void free_vetavars(t_vetavars *vars)
 142 {
 143     if (vars->vscale_nhc != NULL)
 144     {
 145         sfree(vars->vscale_nhc);
 146     }
 147 }
 148
 149 static int pcomp(const void *p1, const void *p2)
 150 {
 151   int     db;
 152   atom_id min1,min2,max1,max2;
 153   t_sortblock *a1=(t_sortblock *)p1;
 154   t_sortblock *a2=(t_sortblock *)p2;
 155
 156   db=a1->blocknr-a2->blocknr;
 157
 158   if (db != 0)
 159     return db;
 160
 161   min1=min(a1->iatom[1],a1->iatom[2]);
 162   max1=max(a1->iatom[1],a1->iatom[2]);
 163   min2=min(a2->iatom[1],a2->iatom[2]);
 164   max2=max(a2->iatom[1],a2->iatom[2]);
 165
 166   if (min1 == min2)
 167     return max1-max2;
 168   else
 169     return min1-min2;
 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",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,
 298                    gmx_bool bMolPBC,matrix box,
 299                    real lambda,real *dvdlambda,
 300                    rvec *v,tensor *vir,
 301                    t_nrnb *nrnb,int econq,gmx_bool bPscal,
 302                    real veta, real vetanew)
 303 {
 304     gmx_bool    bOK,bDump;
 305     int     start,homenr,nrend;
 306     int     i,j,d;
 307     int     ncons,settle_error;
 308     tensor  vir_r_m_dr;
 309     rvec    *vstor;
 310     real    invdt,vir_fac,t;
 311     t_ilist *settle;
 312     int     nsettle;
 313     t_pbc   pbc,*pbc_null;
 314     char    buf[22];
 315     t_vetavars vetavar;
 316     int     nth,th;
 317
 318     if (econq == econqForceDispl && !EI_ENERGY_MINIMIZATION(ir->eI))
 319     {
 320         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");
 321     }
 322
 323     bOK   = TRUE;
 324     bDump = FALSE;
 325
 326     start  = md->start;
 327     homenr = md->homenr;
 328     nrend = start+homenr;
 329
 330     /* set constants for pressure control integration */
 331     init_vetavars(&vetavar,econq!=econqCoord,
 332                   veta,vetanew,ir,ekind,bPscal);
 333
 334     if (ir->delta_t == 0)
 335     {
 336         invdt = 0;
 337     }
 338     else
 339     {
 340         invdt  = 1/ir->delta_t;
 341     }
 342
 343     if (ir->efep != efepNO && EI_DYNAMICS(ir->eI))
 344     {
 345         /* Set the constraint lengths for the step at which this configuration
 346          * is meant to be. The invmasses should not be changed.
 347          */
 348         lambda += delta_step*ir->fepvals->delta_lambda;
 349     }
 350
 351     if (vir != NULL)
 352     {
 353         clear_mat(vir_r_m_dr);
 354     }
 355
 356     where();
 357
 358     settle  = &idef->il[F_SETTLE];
 359     nsettle = settle->nr/(1+NRAL(F_SETTLE));
 360
 361     if (nsettle > 0)
 362     {
 363         nth = gmx_omp_nthreads_get(emntSETTLE);
 364     }
 365     else
 366     {
 367         nth = 1;
 368     }
 369
 370     if (nth > 1 && constr->vir_r_m_dr_th == NULL)
 371     {
 372         snew(constr->vir_r_m_dr_th,nth);
 373         snew(constr->settle_error,nth);
 374     }
 375
 376     settle_error = -1;
 377
 378     /* We do not need full pbc when constraints do not cross charge groups,
 379      * i.e. when dd->constraint_comm==NULL.
 380      * Note that PBC for constraints is different from PBC for bondeds.
 381      * For constraints there is both forward and backward communication.
 382      */
 383     if (ir->ePBC != epbcNONE &&
 384         (cr->dd || bMolPBC) && !(cr->dd && cr->dd->constraint_comm==NULL))
 385     {
 386         /* With pbc=screw the screw has been changed to a shift
 387          * by the constraint coordinate communication routine,
 388          * so that here we can use normal pbc.
 389          */
 390         pbc_null = set_pbc_dd(&pbc,ir->ePBC,cr->dd,FALSE,box);
 391     }
 392     else
 393     {
 394         pbc_null = NULL;
 395     }
 396
 397     /* Communicate the coordinates required for the non-local constraints
 398      * for LINCS and/or SETTLE.
 399      */
 400     if (cr->dd)
 401     {
 402         dd_move_x_constraints(cr->dd,box,x,xprime);
 403     }
 404         else if (PARTDECOMP(cr))
 405         {
 406                 pd_move_x_constraints(cr,x,xprime);
 407         }
 408
 409     if (constr->lincsd != NULL)
 410     {
 411         bOK = constrain_lincs(fplog,bLog,bEner,ir,step,constr->lincsd,md,cr,
 412                               x,xprime,min_proj,
 413                               box,pbc_null,lambda,dvdlambda,
 414                               invdt,v,vir!=NULL,vir_r_m_dr,
 415                               econq,nrnb,
 416                               constr->maxwarn,&constr->warncount_lincs);
 417         if (!bOK && constr->maxwarn >= 0)
 418         {
 419             if (fplog != NULL)
 420             {
 421                 fprintf(fplog,"Constraint error in algorithm %s at step %s\n",
 422                         econstr_names[econtLINCS],gmx_step_str(step,buf));
 423             }
 424             bDump = TRUE;
 425         }
 426     }
 427
 428     if (constr->nblocks > 0)
 429     {
 430         switch (econq) {
 431         case (econqCoord):
 432             bOK = bshakef(fplog,constr->shaked,
 433                           homenr,md->invmass,constr->nblocks,constr->sblock,
 434                           idef,ir,x,xprime,nrnb,
 435                           constr->lagr,lambda,dvdlambda,
 436                           invdt,v,vir!=NULL,vir_r_m_dr,
 437                           constr->maxwarn>=0,econq,&vetavar);
 438             break;
 439         case (econqVeloc):
 440             bOK = bshakef(fplog,constr->shaked,
 441                           homenr,md->invmass,constr->nblocks,constr->sblock,
 442                           idef,ir,x,min_proj,nrnb,
 443                           constr->lagr,lambda,dvdlambda,
 444                           invdt,NULL,vir!=NULL,vir_r_m_dr,
 445                           constr->maxwarn>=0,econq,&vetavar);
 446             break;
 447         default:
 448             gmx_fatal(FARGS,"Internal error, SHAKE called for constraining something else than coordinates");
 449             break;
 450         }
 451
 452         if (!bOK && constr->maxwarn >= 0)
 453         {
 454             if (fplog != NULL)
 455             {
 456                 fprintf(fplog,"Constraint error in algorithm %s at step %s\n",
 457                         econstr_names[econtSHAKE],gmx_step_str(step,buf));
 458             }
 459             bDump = TRUE;
 460         }
 461     }
 462
 463     if (nsettle > 0)
 464     {
 465         int calcvir_atom_end;
 466
 467         if (vir == NULL)
 468         {
 469             calcvir_atom_end = 0;
 470         }
 471         else
 472         {
 473             calcvir_atom_end = md->start + md->homenr;
 474         }
 475
 476         switch (econq)
 477         {
 478         case econqCoord:
 479 #pragma omp parallel for num_threads(nth) schedule(static)
 480             for(th=0; th<nth; th++)
 481             {
 482                 int start_th,end_th;
 483
 484                 if (th > 0)
 485                 {
 486                     clear_mat(constr->vir_r_m_dr_th[th]);
 487                 }
 488
 489                 start_th = (nsettle* th   )/nth;
 490                 end_th   = (nsettle*(th+1))/nth;
 491                 if (start_th >= 0 && end_th - start_th > 0)
 492                 {
 493                     csettle(constr->settled,
 494                             end_th-start_th,
 495                             settle->iatoms+start_th*(1+NRAL(F_SETTLE)),
 496                             pbc_null,
 497                             x[0],xprime[0],
 498                             invdt,v?v[0]:NULL,calcvir_atom_end,
 499                             th == 0 ? vir_r_m_dr : constr->vir_r_m_dr_th[th],
 500                             th == 0 ? &settle_error : &constr->settle_error[th],
 501                             &vetavar);
 502                 }
 503             }
 504             inc_nrnb(nrnb,eNR_SETTLE,nsettle);
 505             if (v != NULL)
 506             {
 507                 inc_nrnb(nrnb,eNR_CONSTR_V,nsettle*3);
 508             }
 509             if (vir != NULL)
 510             {
 511                 inc_nrnb(nrnb,eNR_CONSTR_VIR,nsettle*3);
 512             }
 513             break;
 514         case econqVeloc:
 515         case econqDeriv:
 516         case econqForce:
 517         case econqForceDispl:
 518 #pragma omp parallel for num_threads(nth) schedule(static)
 519             for(th=0; th<nth; th++)
 520             {
 521                 int start_th,end_th;
 522
 523                 if (th > 0)
 524                 {
 525                     clear_mat(constr->vir_r_m_dr_th[th]);
 526                 }
 527
 528                 start_th = (nsettle* th   )/nth;
 529                 end_th   = (nsettle*(th+1))/nth;
 530
 531                 if (start_th >= 0 && end_th - start_th > 0)
 532                 {
 533                     settle_proj(fplog,constr->settled,econq,
 534                                 end_th-start_th,
 535                                 settle->iatoms+start_th*(1+NRAL(F_SETTLE)),
 536                                 pbc_null,
 537                                 x,
 538                                 xprime,min_proj,calcvir_atom_end,
 539                                 th == 0 ? vir_r_m_dr : constr->vir_r_m_dr_th[th],
 540                                 &vetavar);
 541                 }
 542             }
 543             /* This is an overestimate */
 544             inc_nrnb(nrnb,eNR_SETTLE,nsettle);
 545             break;
 546         case econqDeriv_FlexCon:
 547             /* Nothing to do, since the are no flexible constraints in settles */
 548             break;
 549         default:
 550             gmx_incons("Unknown constraint quantity for settle");
 551         }
 552     }
 553
 554     if (settle->nr > 0)
 555     {
 556         /* Combine virial and error info of the other threads */
 557         for(i=1; i<nth; i++)
 558         {
 559             m_add(vir_r_m_dr,constr->vir_r_m_dr_th[i],vir_r_m_dr);
 560             settle_error = constr->settle_error[i];
 561         }
 562
 563         if (econq == econqCoord && settle_error >= 0)
 564         {
 565             bOK = FALSE;
 566             if (constr->maxwarn >= 0)
 567             {
 568                 char buf[256];
 569                 sprintf(buf,
 570                         "\nstep " gmx_large_int_pfmt ": Water molecule starting at atom %d can not be "
 571                         "settled.\nCheck for bad contacts and/or reduce the timestep if appropriate.\n",
 572                         step,ddglatnr(cr->dd,settle->iatoms[settle_error*(1+NRAL(F_SETTLE))+1]));
 573                 if (fplog)
 574                 {
 575                     fprintf(fplog,"%s",buf);
 576                 }
 577                 fprintf(stderr,"%s",buf);
 578                 constr->warncount_settle++;
 579                 if (constr->warncount_settle > constr->maxwarn)
 580                 {
 581                     too_many_constraint_warnings(-1,constr->warncount_settle);
 582                 }
 583                 bDump = TRUE;
 584             }
 585         }
 586     }
 587
 588     free_vetavars(&vetavar);
 589
 590     if (vir != NULL)
 591     {
 592         switch (econq)
 593         {
 594         case econqCoord:
 595             vir_fac = 0.5/(ir->delta_t*ir->delta_t);
 596             break;
 597         case econqVeloc:
 598             vir_fac = 0.5/ir->delta_t;
 599             break;
 600         case econqForce:
 601         case econqForceDispl:
 602             vir_fac = 0.5;
 603             break;
 604         default:
 605             vir_fac = 0;
 606             gmx_incons("Unsupported constraint quantity for virial");
 607         }
 608
 609         if (EI_VV(ir->eI))
 610         {
 611             vir_fac *= 2;  /* only constraining over half the distance here */
 612         }
 613         for(i=0; i<DIM; i++)
 614         {
 615             for(j=0; j<DIM; j++)
 616             {
 617                 (*vir)[i][j] = vir_fac*vir_r_m_dr[i][j];
 618             }
 619         }
 620     }
 621
 622     if (bDump)
 623     {
 624         dump_confs(fplog,step,constr->warn_mtop,start,homenr,cr,x,xprime,box);
 625     }
 626
 627     if (econq == econqCoord)
 628     {
 629         if (ir->ePull == epullCONSTRAINT)
 630         {
 631             if (EI_DYNAMICS(ir->eI))
 632             {
 633                 t = ir->init_t + (step + delta_step)*ir->delta_t;
 634             }
 635             else
 636             {
 637                 t = ir->init_t;
 638             }
 639             set_pbc(&pbc,ir->ePBC,box);
 640             pull_constraint(ir->pull,md,&pbc,cr,ir->delta_t,t,x,xprime,v,*vir);
 641         }
 642         if (constr->ed && delta_step > 0)
 643         {
 644             /* apply the essential dynamcs constraints here */
 645             do_edsam(ir,step,cr,xprime,v,box,constr->ed);
 646         }
 647     }
 648
 649     return bOK;
 650 }
 651
 652 real *constr_rmsd_data(struct gmx_constr *constr)
 653 {
 654   if (constr->lincsd)
 655     return lincs_rmsd_data(constr->lincsd);
 656   else
 657     return NULL;
 658 }
 659
 660 real constr_rmsd(struct gmx_constr *constr,gmx_bool bSD2)
 661 {
 662   if (constr->lincsd)
 663     return lincs_rmsd(constr->lincsd,bSD2);
 664   else
 665     return 0;
 666 }
 667
 668 static void make_shake_sblock_pd(struct gmx_constr *constr,
 669                                  t_idef *idef,t_mdatoms *md)
 670 {
 671   int  i,j,m,ncons;
 672   int  bstart,bnr;
 673   t_blocka    sblocks;
 674   t_sortblock *sb;
 675   t_iatom     *iatom;
 676   atom_id     *inv_sblock;
 677
 678   /* Since we are processing the local topology,
 679    * the F_CONSTRNC ilist has been concatenated to the F_CONSTR ilist.
 680    */
 681   ncons = idef->il[F_CONSTR].nr/3;
 682
 683   init_blocka(&sblocks);
 684   gen_sblocks(NULL,md->start,md->start+md->homenr,idef,&sblocks,FALSE);
 685
 686   /*
 687     bstart=(idef->nodeid > 0) ? blocks->multinr[idef->nodeid-1] : 0;
 688     nblocks=blocks->multinr[idef->nodeid] - bstart;
 689   */
 690   bstart  = 0;
 691   constr->nblocks = sblocks.nr;
 692   if (debug)
 693     fprintf(debug,"ncons: %d, bstart: %d, nblocks: %d\n",
 694             ncons,bstart,constr->nblocks);
 695
 696   /* Calculate block number for each atom */
 697   inv_sblock = make_invblocka(&sblocks,md->nr);
 698
 699   done_blocka(&sblocks);
 700
 701   /* Store the block number in temp array and
 702    * sort the constraints in order of the sblock number
 703    * and the atom numbers, really sorting a segment of the array!
 704    */
 705 #ifdef DEBUGIDEF
 706   pr_idef(fplog,0,"Before Sort",idef);
 707 #endif
 708   iatom=idef->il[F_CONSTR].iatoms;
 709   snew(sb,ncons);
 710   for(i=0; (i<ncons); i++,iatom+=3) {
 711     for(m=0; (m<3); m++)
 712       sb[i].iatom[m] = iatom[m];
 713     sb[i].blocknr = inv_sblock[iatom[1]];
 714   }
 715
 716   /* Now sort the blocks */
 717   if (debug) {
 718     pr_sortblock(debug,"Before sorting",ncons,sb);
 719     fprintf(debug,"Going to sort constraints\n");
 720   }
 721
 722   qsort(sb,ncons,(size_t)sizeof(*sb),pcomp);
 723
 724   if (debug) {
 725     pr_sortblock(debug,"After sorting",ncons,sb);
 726   }
 727
 728   iatom=idef->il[F_CONSTR].iatoms;
 729   for(i=0; (i<ncons); i++,iatom+=3)
 730     for(m=0; (m<3); m++)
 731       iatom[m]=sb[i].iatom[m];
 732 #ifdef DEBUGIDEF
 733   pr_idef(fplog,0,"After Sort",idef);
 734 #endif
 735
 736   j=0;
 737   snew(constr->sblock,constr->nblocks+1);
 738   bnr=-2;
 739   for(i=0; (i<ncons); i++) {
 740     if (sb[i].blocknr != bnr) {
 741       bnr=sb[i].blocknr;
 742       constr->sblock[j++]=3*i;
 743     }
 744   }
 745   /* Last block... */
 746   constr->sblock[j++] = 3*ncons;
 747
 748   if (j != (constr->nblocks+1)) {
 749     fprintf(stderr,"bstart: %d\n",bstart);
 750     fprintf(stderr,"j: %d, nblocks: %d, ncons: %d\n",
 751             j,constr->nblocks,ncons);
 752     for(i=0; (i<ncons); i++)
 753       fprintf(stderr,"i: %5d  sb[i].blocknr: %5u\n",i,sb[i].blocknr);
 754     for(j=0; (j<=constr->nblocks); j++)
 755       fprintf(stderr,"sblock[%3d]=%5d\n",j,(int)constr->sblock[j]);
 756     gmx_fatal(FARGS,"DEATH HORROR: "
 757               "sblocks does not match idef->il[F_CONSTR]");
 758   }
 759   sfree(sb);
 760   sfree(inv_sblock);
 761 }
 762
 763 static void make_shake_sblock_dd(struct gmx_constr *constr,
 764                                  t_ilist *ilcon,t_block *cgs,
 765                                  gmx_domdec_t *dd)
 766 {
 767   int ncons,c,cg;
 768   t_iatom *iatom;
 769
 770   if (dd->ncg_home+1 > constr->sblock_nalloc) {
 771     constr->sblock_nalloc = over_alloc_dd(dd->ncg_home+1);
 772     srenew(constr->sblock,constr->sblock_nalloc);
 773   }
 774
 775   ncons = ilcon->nr/3;
 776   iatom = ilcon->iatoms;
 777   constr->nblocks = 0;
 778   cg = 0;
 779   for(c=0; c<ncons; c++) {
 780     if (c == 0 || iatom[1] >= cgs->index[cg+1]) {
 781       constr->sblock[constr->nblocks++] = 3*c;
 782       while (iatom[1] >= cgs->index[cg+1])
 783         cg++;
 784     }
 785     iatom += 3;
 786   }
 787   constr->sblock[constr->nblocks] = 3*ncons;
 788 }
 789
 790 t_blocka make_at2con(int start,int natoms,
 791                      t_ilist *ilist,t_iparams *iparams,
 792                      gmx_bool bDynamics,int *nflexiblecons)
 793 {
 794   int *count,ncon,con,con_tot,nflexcon,ftype,i,a;
 795   t_iatom  *ia;
 796   t_blocka at2con;
 797   gmx_bool bFlexCon;
 798
 799   snew(count,natoms);
 800   nflexcon = 0;
 801   for(ftype=F_CONSTR; ftype<=F_CONSTRNC; ftype++) {
 802     ncon = ilist[ftype].nr/3;
 803     ia   = ilist[ftype].iatoms;
 804     for(con=0; con<ncon; con++) {
 805       bFlexCon = (iparams[ia[0]].constr.dA == 0 &&
 806                   iparams[ia[0]].constr.dB == 0);
 807       if (bFlexCon)
 808         nflexcon++;
 809       if (bDynamics || !bFlexCon) {
 810         for(i=1; i<3; i++) {
 811           a = ia[i] - start;
 812           count[a]++;
 813         }
 814       }
 815       ia += 3;
 816     }
 817   }
 818   *nflexiblecons = nflexcon;
 819
 820   at2con.nr = natoms;
 821   at2con.nalloc_index = at2con.nr+1;
 822   snew(at2con.index,at2con.nalloc_index);
 823   at2con.index[0] = 0;
 824   for(a=0; a<natoms; a++) {
 825     at2con.index[a+1] = at2con.index[a] + count[a];
 826     count[a] = 0;
 827   }
 828   at2con.nra = at2con.index[natoms];
 829   at2con.nalloc_a = at2con.nra;
 830   snew(at2con.a,at2con.nalloc_a);
 831
 832   /* The F_CONSTRNC constraints have constraint numbers
 833    * that continue after the last F_CONSTR constraint.
 834    */
 835   con_tot = 0;
 836   for(ftype=F_CONSTR; ftype<=F_CONSTRNC; ftype++) {
 837     ncon = ilist[ftype].nr/3;
 838     ia   = ilist[ftype].iatoms;
 839     for(con=0; con<ncon; con++) {
 840       bFlexCon = (iparams[ia[0]].constr.dA == 0 &&
 841                   iparams[ia[0]].constr.dB == 0);
 842       if (bDynamics || !bFlexCon) {
 843         for(i=1; i<3; i++) {
 844           a = ia[i] - start;
 845           at2con.a[at2con.index[a]+count[a]++] = con_tot;
 846         }
 847       }
 848       con_tot++;
 849       ia += 3;
 850     }
 851   }
 852
 853   sfree(count);
 854
 855   return at2con;
 856 }
 857
 858 static int *make_at2settle(int natoms,const t_ilist *ilist)
 859 {
 860     int *at2s;
 861     int a,stride,s;
 862
 863     snew(at2s,natoms);
 864     /* Set all to no settle */
 865     for(a=0; a<natoms; a++)
 866     {
 867         at2s[a] = -1;
 868     }
 869
 870     stride = 1 + NRAL(F_SETTLE);
 871
 872     for(s=0; s<ilist->nr; s+=stride)
 873     {
 874         at2s[ilist->iatoms[s+1]] = s/stride;
 875         at2s[ilist->iatoms[s+2]] = s/stride;
 876         at2s[ilist->iatoms[s+3]] = s/stride;
 877     }
 878
 879     return at2s;
 880 }
 881
 882 void set_constraints(struct gmx_constr *constr,
 883                      gmx_localtop_t *top,t_inputrec *ir,
 884                      t_mdatoms *md,t_commrec *cr)
 885 {
 886     t_idef *idef;
 887     int    ncons;
 888     t_ilist *settle;
 889     int    iO,iH;
 890
 891     idef = &top->idef;
 892
 893     if (constr->ncon_tot > 0)
 894     {
 895         /* We are using the local topology,
 896          * so there are only F_CONSTR constraints.
 897          */
 898         ncons = idef->il[F_CONSTR].nr/3;
 899
 900         /* With DD we might also need to call LINCS with ncons=0 for
 901          * communicating coordinates to other nodes that do have constraints.
 902          */
 903         if (ir->eConstrAlg == econtLINCS)
 904         {
 905             set_lincs(idef,md,EI_DYNAMICS(ir->eI),cr,constr->lincsd);
 906         }
 907         if (ir->eConstrAlg == econtSHAKE)
 908         {
 909             if (cr->dd)
 910             {
 911                 make_shake_sblock_dd(constr,&idef->il[F_CONSTR],&top->cgs,cr->dd);
 912             }
 913             else
 914             {
 915                 make_shake_sblock_pd(constr,idef,md);
 916             }
 917             if (ncons > constr->lagr_nalloc)
 918             {
 919                 constr->lagr_nalloc = over_alloc_dd(ncons);
 920                 srenew(constr->lagr,constr->lagr_nalloc);
 921             }
 922         }
 923     }
 924
 925     if (idef->il[F_SETTLE].nr > 0 && constr->settled == NULL)
 926     {
 927         settle = &idef->il[F_SETTLE];
 928         iO = settle->iatoms[1];
 929         iH = settle->iatoms[2];
 930         constr->settled =
 931             settle_init(md->massT[iO],md->massT[iH],
 932                         md->invmass[iO],md->invmass[iH],
 933                         idef->iparams[settle->iatoms[0]].settle.doh,
 934                         idef->iparams[settle->iatoms[0]].settle.dhh);
 935     }
 936
 937     /* Make a selection of the local atoms for essential dynamics */
 938     if (constr->ed && cr->dd)
 939     {
 940         dd_make_local_ed_indices(cr->dd,constr->ed);
 941     }
 942 }
 943
 944 static void constr_recur(t_blocka *at2con,
 945                          t_ilist *ilist,t_iparams *iparams,gmx_bool bTopB,
 946                          int at,int depth,int nc,int *path,
 947                          real r0,real r1,real *r2max,
 948                          int *count)
 949 {
 950   int  ncon1;
 951   t_iatom *ia1,*ia2;
 952   int  c,con,a1;
 953   gmx_bool bUse;
 954   t_iatom *ia;
 955   real len,rn0,rn1;
 956
 957   (*count)++;
 958
 959   ncon1 = ilist[F_CONSTR].nr/3;
 960   ia1   = ilist[F_CONSTR].iatoms;
 961   ia2   = ilist[F_CONSTRNC].iatoms;
 962
 963   /* Loop over all constraints connected to this atom */
 964   for(c=at2con->index[at]; c<at2con->index[at+1]; c++) {
 965     con = at2con->a[c];
 966     /* Do not walk over already used constraints */
 967     bUse = TRUE;
 968     for(a1=0; a1<depth; a1++) {
 969       if (con == path[a1])
 970         bUse = FALSE;
 971     }
 972     if (bUse) {
 973       ia = constr_iatomptr(ncon1,ia1,ia2,con);
 974       /* Flexible constraints currently have length 0, which is incorrect */
 975       if (!bTopB)
 976         len = iparams[ia[0]].constr.dA;
 977       else
 978         len = iparams[ia[0]].constr.dB;
 979       /* In the worst case the bond directions alternate */
 980       if (nc % 2 == 0) {
 981         rn0 = r0 + len;
 982         rn1 = r1;
 983       } else {
 984         rn0 = r0;
 985         rn1 = r1 + len;
 986       }
 987       /* Assume angles of 120 degrees between all bonds */
 988       if (rn0*rn0 + rn1*rn1 + rn0*rn1 > *r2max) {
 989         *r2max = rn0*rn0 + rn1*rn1 + r0*rn1;
 990         if (debug) {
 991           fprintf(debug,"Found longer constraint distance: r0 %5.3f r1 %5.3f rmax %5.3f\n", rn0,rn1,sqrt(*r2max));
 992           for(a1=0; a1<depth; a1++)
 993             fprintf(debug," %d %5.3f",
 994                     path[a1],
 995                     iparams[constr_iatomptr(ncon1,ia1,ia2,con)[0]].constr.dA);
 996           fprintf(debug," %d %5.3f\n",con,len);
 997         }
 998       }
 999       /* Limit the number of recursions to 1000*nc,
1000        * so a call does not take more than a second,
1001        * even for highly connected systems.
1002        */
1003       if (depth + 1 < nc && *count < 1000*nc) {
1004         if (ia[1] == at)
1005           a1 = ia[2];
1006         else
1007           a1 = ia[1];
1008         /* Recursion */
1009         path[depth] = con;
1010         constr_recur(at2con,ilist,iparams,
1011                      bTopB,a1,depth+1,nc,path,rn0,rn1,r2max,count);
1012         path[depth] = -1;
1013       }
1014     }
1015   }
1016 }
1017
1018 static real constr_r_max_moltype(FILE *fplog,
1019                                  gmx_moltype_t *molt,t_iparams *iparams,
1020                                  t_inputrec *ir)
1021 {
1022   int natoms,nflexcon,*path,at,count;
1023
1024   t_blocka at2con;
1025   real r0,r1,r2maxA,r2maxB,rmax,lam0,lam1;
1026
1027   if (molt->ilist[F_CONSTR].nr   == 0 &&
1028       molt->ilist[F_CONSTRNC].nr == 0) {
1029     return 0;
1030   }
1031
1032   natoms = molt->atoms.nr;
1033
1034   at2con = make_at2con(0,natoms,molt->ilist,iparams,
1035                        EI_DYNAMICS(ir->eI),&nflexcon);
1036   snew(path,1+ir->nProjOrder);
1037   for(at=0; at<1+ir->nProjOrder; at++)
1038     path[at] = -1;
1039
1040   r2maxA = 0;
1041   for(at=0; at<natoms; at++) {
1042     r0 = 0;
1043     r1 = 0;
1044
1045     count = 0;
1046     constr_recur(&at2con,molt->ilist,iparams,
1047                  FALSE,at,0,1+ir->nProjOrder,path,r0,r1,&r2maxA,&count);
1048   }
1049   if (ir->efep == efepNO) {
1050     rmax = sqrt(r2maxA);
1051   } else {
1052     r2maxB = 0;
1053     for(at=0; at<natoms; at++) {
1054       r0 = 0;
1055       r1 = 0;
1056       count = 0;
1057       constr_recur(&at2con,molt->ilist,iparams,
1058                    TRUE,at,0,1+ir->nProjOrder,path,r0,r1,&r2maxB,&count);
1059     }
1060     lam0 = ir->fepvals->init_lambda;
1061     if (EI_DYNAMICS(ir->eI))
1062       lam0 += ir->init_step*ir->fepvals->delta_lambda;
1063     rmax = (1 - lam0)*sqrt(r2maxA) + lam0*sqrt(r2maxB);
1064     if (EI_DYNAMICS(ir->eI)) {
1065       lam1 = ir->fepvals->init_lambda + (ir->init_step + ir->nsteps)*ir->fepvals->delta_lambda;
1066       rmax = max(rmax,(1 - lam1)*sqrt(r2maxA) + lam1*sqrt(r2maxB));
1067     }
1068   }
1069
1070   done_blocka(&at2con);
1071   sfree(path);
1072
1073   return rmax;
1074 }
1075
1076 real constr_r_max(FILE *fplog,gmx_mtop_t *mtop,t_inputrec *ir)
1077 {
1078   int mt;
1079   real rmax;
1080
1081   rmax = 0;
1082   for(mt=0; mt<mtop->nmoltype; mt++) {
1083     rmax = max(rmax,
1084                constr_r_max_moltype(fplog,&mtop->moltype[mt],
1085                                     mtop->ffparams.iparams,ir));
1086   }
1087
1088   if (fplog)
1089     fprintf(fplog,"Maximum distance for %d constraints, at 120 deg. angles, all-trans: %.3f nm\n",1+ir->nProjOrder,rmax);
1090
1091   return rmax;
1092 }
1093
1094 gmx_constr_t init_constraints(FILE *fplog,
1095                               gmx_mtop_t *mtop,t_inputrec *ir,
1096                               gmx_edsam_t ed,t_state *state,
1097                               t_commrec *cr)
1098 {
1099     int  ncon,nset,nmol,settle_type,i,natoms,mt,nflexcon;
1100     struct gmx_constr *constr;
1101     char *env;
1102     t_ilist *ilist;
1103     gmx_mtop_ilistloop_t iloop;
1104
1105     ncon =
1106         gmx_mtop_ftype_count(mtop,F_CONSTR) +
1107         gmx_mtop_ftype_count(mtop,F_CONSTRNC);
1108     nset = gmx_mtop_ftype_count(mtop,F_SETTLE);
1109
1110     if (ncon+nset == 0 && ir->ePull != epullCONSTRAINT && ed == NULL)
1111     {
1112         return NULL;
1113     }
1114
1115     snew(constr,1);
1116
1117     constr->ncon_tot = ncon;
1118     constr->nflexcon = 0;
1119     if (ncon > 0)
1120     {
1121         constr->n_at2con_mt = mtop->nmoltype;
1122         snew(constr->at2con_mt,constr->n_at2con_mt);
1123         for(mt=0; mt<mtop->nmoltype; mt++)
1124         {
1125             constr->at2con_mt[mt] = make_at2con(0,mtop->moltype[mt].atoms.nr,
1126                                                 mtop->moltype[mt].ilist,
1127                                                 mtop->ffparams.iparams,
1128                                                 EI_DYNAMICS(ir->eI),&nflexcon);
1129             for(i=0; i<mtop->nmolblock; i++)
1130             {
1131                 if (mtop->molblock[i].type == mt)
1132                 {
1133                     constr->nflexcon += mtop->molblock[i].nmol*nflexcon;
1134                 }
1135             }
1136         }
1137
1138         if (constr->nflexcon > 0)
1139         {
1140             if (fplog)
1141             {
1142                 fprintf(fplog,"There are %d flexible constraints\n",
1143                         constr->nflexcon);
1144                 if (ir->fc_stepsize == 0)
1145                 {
1146                     fprintf(fplog,"\n"
1147                             "WARNING: step size for flexible constraining = 0\n"
1148                             "         All flexible constraints will be rigid.\n"
1149                             "         Will try to keep all flexible constraints at their original length,\n"
1150                             "         but the lengths may exhibit some drift.\n\n");
1151                     constr->nflexcon = 0;
1152                 }
1153             }
1154             if (constr->nflexcon > 0)
1155             {
1156                 please_cite(fplog,"Hess2002");
1157             }
1158         }
1159
1160         if (ir->eConstrAlg == econtLINCS)
1161         {
1162             constr->lincsd = init_lincs(fplog,mtop,
1163                                         constr->nflexcon,constr->at2con_mt,
1164                                         DOMAINDECOMP(cr) && cr->dd->bInterCGcons,
1165                                         ir->nLincsIter,ir->nProjOrder);
1166         }
1167
1168         if (ir->eConstrAlg == econtSHAKE) {
1169             if (DOMAINDECOMP(cr) && cr->dd->bInterCGcons)
1170             {
1171                 gmx_fatal(FARGS,"SHAKE is not supported with domain decomposition and constraint that cross charge group boundaries, use LINCS");
1172             }
1173             if (constr->nflexcon)
1174             {
1175                 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");
1176             }
1177             please_cite(fplog,"Ryckaert77a");
1178             if (ir->bShakeSOR)
1179             {
1180                 please_cite(fplog,"Barth95a");
1181             }
1182
1183             constr->shaked = shake_init();
1184         }
1185     }
1186
1187     if (nset > 0) {
1188         please_cite(fplog,"Miyamoto92a");
1189
1190         constr->bInterCGsettles = inter_charge_group_settles(mtop);
1191
1192         /* Check that we have only one settle type */
1193         settle_type = -1;
1194         iloop = gmx_mtop_ilistloop_init(mtop);
1195         while (gmx_mtop_ilistloop_next(iloop,&ilist,&nmol))
1196         {
1197             for (i=0; i<ilist[F_SETTLE].nr; i+=4)
1198             {
1199                 if (settle_type == -1)
1200                 {
1201                     settle_type = ilist[F_SETTLE].iatoms[i];
1202                 }
1203                 else if (ilist[F_SETTLE].iatoms[i] != settle_type)
1204                 {
1205                     gmx_fatal(FARGS,
1206                               "The [molecules] section of your topology specifies more than one block of\n"
1207                               "a [moleculetype] with a [settles] block. Only one such is allowed. If you\n"
1208                               "are trying to partition your solvent into different *groups* (e.g. for\n"
1209                               "freezing, T-coupling, etc.) then you are using the wrong approach. Index\n"
1210                               "files specify groups. Otherwise, you may wish to change the least-used\n"
1211                               "block of molecules with SETTLE constraints into 3 normal constraints.");
1212                 }
1213             }
1214         }
1215
1216         constr->n_at2settle_mt = mtop->nmoltype;
1217         snew(constr->at2settle_mt,constr->n_at2settle_mt);
1218         for(mt=0; mt<mtop->nmoltype; mt++)
1219         {
1220             constr->at2settle_mt[mt] =
1221                 make_at2settle(mtop->moltype[mt].atoms.nr,
1222                                &mtop->moltype[mt].ilist[F_SETTLE]);
1223         }
1224     }
1225
1226     constr->maxwarn = 999;
1227     env = getenv("GMX_MAXCONSTRWARN");
1228     if (env)
1229     {
1230         constr->maxwarn = 0;
1231         sscanf(env,"%d",&constr->maxwarn);
1232         if (fplog)
1233         {
1234             fprintf(fplog,
1235                     "Setting the maximum number of constraint warnings to %d\n",
1236                     constr->maxwarn);
1237         }
1238         if (MASTER(cr))
1239         {
1240             fprintf(stderr,
1241                     "Setting the maximum number of constraint warnings to %d\n",
1242                     constr->maxwarn);
1243         }
1244     }
1245     if (constr->maxwarn < 0 && fplog)
1246     {
1247         fprintf(fplog,"maxwarn < 0, will not stop on constraint errors\n");
1248     }
1249     constr->warncount_lincs  = 0;
1250     constr->warncount_settle = 0;
1251
1252     /* Initialize the essential dynamics sampling.
1253      * Put the pointer to the ED struct in constr */
1254     constr->ed = ed;
1255     if (ed != NULL || state->edsamstate.nED > 0)
1256     {
1257         init_edsam(mtop,ir,cr,ed,state->x,state->box,&state->edsamstate);
1258     }
1259
1260     constr->warn_mtop = mtop;
1261
1262     return constr;
1263 }
1264
1265 const t_blocka *atom2constraints_moltype(gmx_constr_t constr)
1266 {
1267   return constr->at2con_mt;
1268 }
1269
1270 const int **atom2settle_moltype(gmx_constr_t constr)
1271 {
1272     return (const int **)constr->at2settle_mt;
1273 }
1274
1275
1276 gmx_bool inter_charge_group_constraints(const gmx_mtop_t *mtop)
1277 {
1278     const gmx_moltype_t *molt;
1279     const t_block *cgs;
1280     const t_ilist *il;
1281     int  mb;
1282     int  nat,*at2cg,cg,a,ftype,i;
1283     gmx_bool bInterCG;
1284
1285     bInterCG = FALSE;
1286     for(mb=0; mb<mtop->nmolblock && !bInterCG; mb++)
1287     {
1288         molt = &mtop->moltype[mtop->molblock[mb].type];
1289
1290         if (molt->ilist[F_CONSTR].nr   > 0 ||
1291             molt->ilist[F_CONSTRNC].nr > 0 ||
1292             molt->ilist[F_SETTLE].nr > 0)
1293         {
1294             cgs  = &molt->cgs;
1295             snew(at2cg,molt->atoms.nr);
1296             for(cg=0; cg<cgs->nr; cg++)
1297             {
1298                 for(a=cgs->index[cg]; a<cgs->index[cg+1]; a++)
1299                     at2cg[a] = cg;
1300             }
1301
1302             for(ftype=F_CONSTR; ftype<=F_CONSTRNC; ftype++)
1303             {
1304                 il = &molt->ilist[ftype];
1305                 for(i=0; i<il->nr && !bInterCG; i+=1+NRAL(ftype))
1306                 {
1307                     if (at2cg[il->iatoms[i+1]] != at2cg[il->iatoms[i+2]])
1308                     {
1309                         bInterCG = TRUE;
1310                     }
1311                 }
1312             }
1313
1314             sfree(at2cg);
1315         }
1316     }
1317
1318     return bInterCG;
1319 }
1320
1321 gmx_bool inter_charge_group_settles(const gmx_mtop_t *mtop)
1322 {
1323     const gmx_moltype_t *molt;
1324     const t_block *cgs;
1325     const t_ilist *il;
1326     int  mb;
1327     int  nat,*at2cg,cg,a,ftype,i;
1328     gmx_bool bInterCG;
1329
1330     bInterCG = FALSE;
1331     for(mb=0; mb<mtop->nmolblock && !bInterCG; mb++)
1332     {
1333         molt = &mtop->moltype[mtop->molblock[mb].type];
1334
1335         if (molt->ilist[F_SETTLE].nr > 0)
1336         {
1337             cgs  = &molt->cgs;
1338             snew(at2cg,molt->atoms.nr);
1339             for(cg=0; cg<cgs->nr; cg++)
1340             {
1341                 for(a=cgs->index[cg]; a<cgs->index[cg+1]; a++)
1342                     at2cg[a] = cg;
1343             }
1344
1345             for(ftype=F_SETTLE; ftype<=F_SETTLE; ftype++)
1346             {
1347                 il = &molt->ilist[ftype];
1348                 for(i=0; i<il->nr && !bInterCG; i+=1+NRAL(F_SETTLE))
1349                 {
1350                     if (at2cg[il->iatoms[i+1]] != at2cg[il->iatoms[i+2]] ||
1351                         at2cg[il->iatoms[i+1]] != at2cg[il->iatoms[i+3]])
1352                     {
1353                         bInterCG = TRUE;
1354                     }
1355                 }
1356             }
1357
1358             sfree(at2cg);
1359         }
1360     }
1361
1362     return bInterCG;
1363 }
1364
1365 /* helper functions for andersen temperature control, because the
1366  * gmx_constr construct is only defined in constr.c. Return the list
1367  * of blocks (get_sblock) and the number of blocks (get_nblocks).  */
1368
1369 extern int *get_sblock(struct gmx_constr *constr)
1370 {
1371     return constr->sblock;
1372 }
1373
1374 extern int get_nblocks(struct gmx_constr *constr)
1375 {
1376     return constr->nblocks;
1377 }