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