src/gromacs/mdlib/constr.c

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