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           *vir_r_m_dr_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  vir_r_m_dr;
 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(vir_r_m_dr);
 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->vir_r_m_dr_th == NULL)
 369     {
 370         snew(constr->vir_r_m_dr_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,vir_r_m_dr,
 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,vir_r_m_dr,
 435                           constr->maxwarn>=0,econq,&vetavar);
 436             break;
 437         case (econqVeloc):
 438             bOK = bshakef(fplog,constr->shaked,
 439                           homenr,md->invmass,constr->nblocks,constr->sblock,
 440                           idef,ir,x,min_proj,nrnb,
 441                           constr->lagr,lambda,dvdlambda,
 442                           invdt,NULL,vir!=NULL,vir_r_m_dr,
 443                           constr->maxwarn>=0,econq,&vetavar);
 444             break;
 445         default:
 446             gmx_fatal(FARGS,"Internal error, SHAKE called for constraining something else than coordinates");
 447             break;
 448         }
 449
 450         if (!bOK && constr->maxwarn >= 0)
 451         {
 452             if (fplog != NULL)
 453             {
 454                 fprintf(fplog,"Constraint error in algorithm %s at step %s\n",
 455                         econstr_names[econtSHAKE],gmx_step_str(step,buf));
 456             }
 457             bDump = TRUE;
 458         }
 459     }
 460
 461     if (nsettle > 0)
 462     {
 463         int calcvir_atom_end;
 464
 465         if (vir == NULL)
 466         {
 467             calcvir_atom_end = 0;
 468         }
 469         else
 470         {
 471             calcvir_atom_end = md->start + md->homenr;
 472         }
 473
 474         switch (econq)
 475         {
 476         case econqCoord:
 477 #pragma omp parallel for num_threads(nth) schedule(static)
 478             for(th=0; th<nth; th++)
 479             {
 480                 int start_th,end_th;
 481
 482                 if (th > 0)
 483                 {
 484                     clear_mat(constr->vir_r_m_dr_th[th]);
 485                 }
 486
 487                 start_th = (nsettle* th   )/nth;
 488                 end_th   = (nsettle*(th+1))/nth;
 489                 if (start_th >= 0 && end_th - start_th > 0)
 490                 {
 491                     csettle(constr->settled,
 492                             end_th-start_th,
 493                             settle->iatoms+start_th*(1+NRAL(F_SETTLE)),
 494                             pbc_null,
 495                             x[0],xprime[0],
 496                             invdt,v?v[0]:NULL,calcvir_atom_end,
 497                             th == 0 ? vir_r_m_dr : constr->vir_r_m_dr_th[th],
 498                             th == 0 ? &settle_error : &constr->settle_error[th],
 499                             &vetavar);
 500                 }
 501             }
 502             inc_nrnb(nrnb,eNR_SETTLE,nsettle);
 503             if (v != NULL)
 504             {
 505                 inc_nrnb(nrnb,eNR_CONSTR_V,nsettle*3);
 506             }
 507             if (vir != NULL)
 508             {
 509                 inc_nrnb(nrnb,eNR_CONSTR_VIR,nsettle*3);
 510             }
 511             break;
 512         case econqVeloc:
 513         case econqDeriv:
 514         case econqForce:
 515         case econqForceDispl:
 516 #pragma omp parallel for num_threads(nth) schedule(static)
 517             for(th=0; th<nth; th++)
 518             {
 519                 int start_th,end_th;
 520
 521                 if (th > 0)
 522                 {
 523                     clear_mat(constr->vir_r_m_dr_th[th]);
 524                 }
 525
 526                 start_th = (nsettle* th   )/nth;
 527                 end_th   = (nsettle*(th+1))/nth;
 528
 529                 if (start_th >= 0 && end_th - start_th > 0)
 530                 {
 531                     settle_proj(fplog,constr->settled,econq,
 532                                 end_th-start_th,
 533                                 settle->iatoms+start_th*(1+NRAL(F_SETTLE)),
 534                                 pbc_null,
 535                                 x,
 536                                 xprime,min_proj,calcvir_atom_end,
 537                                 th == 0 ? vir_r_m_dr : constr->vir_r_m_dr_th[th],
 538                                 &vetavar);
 539                 }
 540             }
 541             /* This is an overestimate */
 542             inc_nrnb(nrnb,eNR_SETTLE,nsettle);
 543             break;
 544         case econqDeriv_FlexCon:
 545             /* Nothing to do, since the are no flexible constraints in settles */
 546             break;
 547         default:
 548             gmx_incons("Unknown constraint quantity for settle");
 549         }
 550     }
 551
 552     if (settle->nr > 0)
 553     {
 554         /* Combine virial and error info of the other threads */
 555         for(i=1; i<nth; i++)
 556         {
 557             m_add(vir_r_m_dr,constr->vir_r_m_dr_th[i],vir_r_m_dr);
 558             settle_error = constr->settle_error[i];
 559         }
 560
 561         if (econq == econqCoord && settle_error >= 0)
 562         {
 563             bOK = FALSE;
 564             if (constr->maxwarn >= 0)
 565             {
 566                 char buf[256];
 567                 sprintf(buf,
 568                         "\nstep " gmx_large_int_pfmt ": Water molecule starting at atom %d can not be "
 569                         "settled.\nCheck for bad contacts and/or reduce the timestep if appropriate.\n",
 570                         step,ddglatnr(cr->dd,settle->iatoms[settle_error*(1+NRAL(F_SETTLE))+1]));
 571                 if (fplog)
 572                 {
 573                     fprintf(fplog,"%s",buf);
 574                 }
 575                 fprintf(stderr,"%s",buf);
 576                 constr->warncount_settle++;
 577                 if (constr->warncount_settle > constr->maxwarn)
 578                 {
 579                     too_many_constraint_warnings(-1,constr->warncount_settle);
 580                 }
 581                 bDump = TRUE;
 582             }
 583         }
 584     }
 585
 586     free_vetavars(&vetavar);
 587
 588     if (vir != NULL)
 589     {
 590         switch (econq)
 591         {
 592         case econqCoord:
 593             vir_fac = 0.5/(ir->delta_t*ir->delta_t);
 594             break;
 595         case econqVeloc:
 596             vir_fac = 0.5/ir->delta_t;
 597             break;
 598         case econqForce:
 599         case econqForceDispl:
 600             vir_fac = 0.5;
 601             break;
 602         default:
 603             vir_fac = 0;
 604             gmx_incons("Unsupported constraint quantity for virial");
 605         }
 606
 607         if (EI_VV(ir->eI))
 608         {
 609             vir_fac *= 2;  /* only constraining over half the distance here */
 610         }
 611         for(i=0; i<DIM; i++)
 612         {
 613             for(j=0; j<DIM; j++)
 614             {
 615                 (*vir)[i][j] = vir_fac*vir_r_m_dr[i][j];
 616             }
 617         }
 618     }
 619
 620     if (bDump)
 621     {
 622         dump_confs(fplog,step,constr->warn_mtop,start,homenr,cr,x,xprime,box);
 623     }
 624
 625     if (econq == econqCoord)
 626     {
 627         if (ir->ePull == epullCONSTRAINT)
 628         {
 629             if (EI_DYNAMICS(ir->eI))
 630             {
 631                 t = ir->init_t + (step + delta_step)*ir->delta_t;
 632             }
 633             else
 634             {
 635                 t = ir->init_t;
 636             }
 637             set_pbc(&pbc,ir->ePBC,box);
 638             pull_constraint(ir->pull,md,&pbc,cr,ir->delta_t,t,x,xprime,v,*vir);
 639         }
 640         if (constr->ed && delta_step > 0)
 641         {
 642             /* apply the essential dynamcs constraints here */
 643             do_edsam(ir,step,md,cr,xprime,v,box,constr->ed);
 644         }
 645     }
 646
 647     return bOK;
 648 }
 649
 650 real *constr_rmsd_data(struct gmx_constr *constr)
 651 {
 652   if (constr->lincsd)
 653     return lincs_rmsd_data(constr->lincsd);
 654   else
 655     return NULL;
 656 }
 657
 658 real constr_rmsd(struct gmx_constr *constr,gmx_bool bSD2)
 659 {
 660   if (constr->lincsd)
 661     return lincs_rmsd(constr->lincsd,bSD2);
 662   else
 663     return 0;
 664 }
 665
 666 static void make_shake_sblock_pd(struct gmx_constr *constr,
 667                                  t_idef *idef,t_mdatoms *md)
 668 {
 669   int  i,j,m,ncons;
 670   int  bstart,bnr;
 671   t_blocka    sblocks;
 672   t_sortblock *sb;
 673   t_iatom     *iatom;
 674   atom_id     *inv_sblock;
 675
 676   /* Since we are processing the local topology,
 677    * the F_CONSTRNC ilist has been concatenated to the F_CONSTR ilist.
 678    */
 679   ncons = idef->il[F_CONSTR].nr/3;
 680
 681   init_blocka(&sblocks);
 682   gen_sblocks(NULL,md->start,md->start+md->homenr,idef,&sblocks,FALSE);
 683
 684   /*
 685     bstart=(idef->nodeid > 0) ? blocks->multinr[idef->nodeid-1] : 0;
 686     nblocks=blocks->multinr[idef->nodeid] - bstart;
 687   */
 688   bstart  = 0;
 689   constr->nblocks = sblocks.nr;
 690   if (debug)
 691     fprintf(debug,"ncons: %d, bstart: %d, nblocks: %d\n",
 692             ncons,bstart,constr->nblocks);
 693
 694   /* Calculate block number for each atom */
 695   inv_sblock = make_invblocka(&sblocks,md->nr);
 696
 697   done_blocka(&sblocks);
 698
 699   /* Store the block number in temp array and
 700    * sort the constraints in order of the sblock number
 701    * and the atom numbers, really sorting a segment of the array!
 702    */
 703 #ifdef DEBUGIDEF
 704   pr_idef(fplog,0,"Before Sort",idef);
 705 #endif
 706   iatom=idef->il[F_CONSTR].iatoms;
 707   snew(sb,ncons);
 708   for(i=0; (i<ncons); i++,iatom+=3) {
 709     for(m=0; (m<3); m++)
 710       sb[i].iatom[m] = iatom[m];
 711     sb[i].blocknr = inv_sblock[iatom[1]];
 712   }
 713
 714   /* Now sort the blocks */
 715   if (debug) {
 716     pr_sortblock(debug,"Before sorting",ncons,sb);
 717     fprintf(debug,"Going to sort constraints\n");
 718   }
 719
 720   qsort(sb,ncons,(size_t)sizeof(*sb),pcomp);
 721
 722   if (debug) {
 723     pr_sortblock(debug,"After sorting",ncons,sb);
 724   }
 725
 726   iatom=idef->il[F_CONSTR].iatoms;
 727   for(i=0; (i<ncons); i++,iatom+=3)
 728     for(m=0; (m<3); m++)
 729       iatom[m]=sb[i].iatom[m];
 730 #ifdef DEBUGIDEF
 731   pr_idef(fplog,0,"After Sort",idef);
 732 #endif
 733
 734   j=0;
 735   snew(constr->sblock,constr->nblocks+1);
 736   bnr=-2;
 737   for(i=0; (i<ncons); i++) {
 738     if (sb[i].blocknr != bnr) {
 739       bnr=sb[i].blocknr;
 740       constr->sblock[j++]=3*i;
 741     }
 742   }
 743   /* Last block... */
 744   constr->sblock[j++] = 3*ncons;
 745
 746   if (j != (constr->nblocks+1)) {
 747     fprintf(stderr,"bstart: %d\n",bstart);
 748     fprintf(stderr,"j: %d, nblocks: %d, ncons: %d\n",
 749             j,constr->nblocks,ncons);
 750     for(i=0; (i<ncons); i++)
 751       fprintf(stderr,"i: %5d  sb[i].blocknr: %5u\n",i,sb[i].blocknr);
 752     for(j=0; (j<=constr->nblocks); j++)
 753       fprintf(stderr,"sblock[%3d]=%5d\n",j,(int)constr->sblock[j]);
 754     gmx_fatal(FARGS,"DEATH HORROR: "
 755               "sblocks does not match idef->il[F_CONSTR]");
 756   }
 757   sfree(sb);
 758   sfree(inv_sblock);
 759 }
 760
 761 static void make_shake_sblock_dd(struct gmx_constr *constr,
 762                                  t_ilist *ilcon,t_block *cgs,
 763                                  gmx_domdec_t *dd)
 764 {
 765   int ncons,c,cg;
 766   t_iatom *iatom;
 767
 768   if (dd->ncg_home+1 > constr->sblock_nalloc) {
 769     constr->sblock_nalloc = over_alloc_dd(dd->ncg_home+1);
 770     srenew(constr->sblock,constr->sblock_nalloc);
 771   }
 772
 773   ncons = ilcon->nr/3;
 774   iatom = ilcon->iatoms;
 775   constr->nblocks = 0;
 776   cg = 0;
 777   for(c=0; c<ncons; c++) {
 778     if (c == 0 || iatom[1] >= cgs->index[cg+1]) {
 779       constr->sblock[constr->nblocks++] = 3*c;
 780       while (iatom[1] >= cgs->index[cg+1])
 781         cg++;
 782     }
 783     iatom += 3;
 784   }
 785   constr->sblock[constr->nblocks] = 3*ncons;
 786 }
 787
 788 t_blocka make_at2con(int start,int natoms,
 789                      t_ilist *ilist,t_iparams *iparams,
 790                      gmx_bool bDynamics,int *nflexiblecons)
 791 {
 792   int *count,ncon,con,con_tot,nflexcon,ftype,i,a;
 793   t_iatom  *ia;
 794   t_blocka at2con;
 795   gmx_bool bFlexCon;
 796
 797   snew(count,natoms);
 798   nflexcon = 0;
 799   for(ftype=F_CONSTR; ftype<=F_CONSTRNC; ftype++) {
 800     ncon = ilist[ftype].nr/3;
 801     ia   = ilist[ftype].iatoms;
 802     for(con=0; con<ncon; con++) {
 803       bFlexCon = (iparams[ia[0]].constr.dA == 0 &&
 804                   iparams[ia[0]].constr.dB == 0);
 805       if (bFlexCon)
 806         nflexcon++;
 807       if (bDynamics || !bFlexCon) {
 808         for(i=1; i<3; i++) {
 809           a = ia[i] - start;
 810           count[a]++;
 811         }
 812       }
 813       ia += 3;
 814     }
 815   }
 816   *nflexiblecons = nflexcon;
 817
 818   at2con.nr = natoms;
 819   at2con.nalloc_index = at2con.nr+1;
 820   snew(at2con.index,at2con.nalloc_index);
 821   at2con.index[0] = 0;
 822   for(a=0; a<natoms; a++) {
 823     at2con.index[a+1] = at2con.index[a] + count[a];
 824     count[a] = 0;
 825   }
 826   at2con.nra = at2con.index[natoms];
 827   at2con.nalloc_a = at2con.nra;
 828   snew(at2con.a,at2con.nalloc_a);
 829
 830   /* The F_CONSTRNC constraints have constraint numbers
 831    * that continue after the last F_CONSTR constraint.
 832    */
 833   con_tot = 0;
 834   for(ftype=F_CONSTR; ftype<=F_CONSTRNC; ftype++) {
 835     ncon = ilist[ftype].nr/3;
 836     ia   = ilist[ftype].iatoms;
 837     for(con=0; con<ncon; con++) {
 838       bFlexCon = (iparams[ia[0]].constr.dA == 0 &&
 839                   iparams[ia[0]].constr.dB == 0);
 840       if (bDynamics || !bFlexCon) {
 841         for(i=1; i<3; i++) {
 842           a = ia[i] - start;
 843           at2con.a[at2con.index[a]+count[a]++] = con_tot;
 844         }
 845       }
 846       con_tot++;
 847       ia += 3;
 848     }
 849   }
 850
 851   sfree(count);
 852
 853   return at2con;
 854 }
 855
 856 static int *make_at2settle(int natoms,const t_ilist *ilist)
 857 {
 858     int *at2s;
 859     int a,stride,s;
 860
 861     snew(at2s,natoms);
 862     /* Set all to no settle */
 863     for(a=0; a<natoms; a++)
 864     {
 865         at2s[a] = -1;
 866     }
 867
 868     stride = 1 + NRAL(F_SETTLE);
 869
 870     for(s=0; s<ilist->nr; s+=stride)
 871     {
 872         at2s[ilist->iatoms[s+1]] = s/stride;
 873         at2s[ilist->iatoms[s+2]] = s/stride;
 874         at2s[ilist->iatoms[s+3]] = s/stride;
 875     }
 876
 877     return at2s;
 878 }
 879
 880 void set_constraints(struct gmx_constr *constr,
 881                      gmx_localtop_t *top,t_inputrec *ir,
 882                      t_mdatoms *md,t_commrec *cr)
 883 {
 884     t_idef *idef;
 885     int    ncons;
 886     t_ilist *settle;
 887     int    iO,iH;
 888
 889     idef = &top->idef;
 890
 891     if (constr->ncon_tot > 0)
 892     {
 893         /* We are using the local topology,
 894          * so there are only F_CONSTR constraints.
 895          */
 896         ncons = idef->il[F_CONSTR].nr/3;
 897
 898         /* With DD we might also need to call LINCS with ncons=0 for
 899          * communicating coordinates to other nodes that do have constraints.
 900          */
 901         if (ir->eConstrAlg == econtLINCS)
 902         {
 903             set_lincs(idef,md,EI_DYNAMICS(ir->eI),cr,constr->lincsd);
 904         }
 905         if (ir->eConstrAlg == econtSHAKE)
 906         {
 907             if (cr->dd)
 908             {
 909                 make_shake_sblock_dd(constr,&idef->il[F_CONSTR],&top->cgs,cr->dd);
 910             }
 911             else
 912             {
 913                 make_shake_sblock_pd(constr,idef,md);
 914             }
 915             if (ncons > constr->lagr_nalloc)
 916             {
 917                 constr->lagr_nalloc = over_alloc_dd(ncons);
 918                 srenew(constr->lagr,constr->lagr_nalloc);
 919             }
 920         }
 921     }
 922
 923     if (idef->il[F_SETTLE].nr > 0 && constr->settled == NULL)
 924     {
 925         settle = &idef->il[F_SETTLE];
 926         iO = settle->iatoms[1];
 927         iH = settle->iatoms[2];
 928         constr->settled =
 929             settle_init(md->massT[iO],md->massT[iH],
 930                         md->invmass[iO],md->invmass[iH],
 931                         idef->iparams[settle->iatoms[0]].settle.doh,
 932                         idef->iparams[settle->iatoms[0]].settle.dhh);
 933     }
 934
 935     /* Make a selection of the local atoms for essential dynamics */
 936     if (constr->ed && cr->dd)
 937     {
 938         dd_make_local_ed_indices(cr->dd,constr->ed);
 939     }
 940 }
 941
 942 static void constr_recur(t_blocka *at2con,
 943                          t_ilist *ilist,t_iparams *iparams,gmx_bool bTopB,
 944                          int at,int depth,int nc,int *path,
 945                          real r0,real r1,real *r2max,
 946                          int *count)
 947 {
 948   int  ncon1;
 949   t_iatom *ia1,*ia2;
 950   int  c,con,a1;
 951   gmx_bool bUse;
 952   t_iatom *ia;
 953   real len,rn0,rn1;
 954
 955   (*count)++;
 956
 957   ncon1 = ilist[F_CONSTR].nr/3;
 958   ia1   = ilist[F_CONSTR].iatoms;
 959   ia2   = ilist[F_CONSTRNC].iatoms;
 960
 961   /* Loop over all constraints connected to this atom */
 962   for(c=at2con->index[at]; c<at2con->index[at+1]; c++) {
 963     con = at2con->a[c];
 964     /* Do not walk over already used constraints */
 965     bUse = TRUE;
 966     for(a1=0; a1<depth; a1++) {
 967       if (con == path[a1])
 968         bUse = FALSE;
 969     }
 970     if (bUse) {
 971       ia = constr_iatomptr(ncon1,ia1,ia2,con);
 972       /* Flexible constraints currently have length 0, which is incorrect */
 973       if (!bTopB)
 974         len = iparams[ia[0]].constr.dA;
 975       else
 976         len = iparams[ia[0]].constr.dB;
 977       /* In the worst case the bond directions alternate */
 978       if (nc % 2 == 0) {
 979         rn0 = r0 + len;
 980         rn1 = r1;
 981       } else {
 982         rn0 = r0;
 983         rn1 = r1 + len;
 984       }
 985       /* Assume angles of 120 degrees between all bonds */
 986       if (rn0*rn0 + rn1*rn1 + rn0*rn1 > *r2max) {
 987         *r2max = rn0*rn0 + rn1*rn1 + r0*rn1;
 988         if (debug) {
 989           fprintf(debug,"Found longer constraint distance: r0 %5.3f r1 %5.3f rmax %5.3f\n", rn0,rn1,sqrt(*r2max));
 990           for(a1=0; a1<depth; a1++)
 991             fprintf(debug," %d %5.3f",
 992                     path[a1],
 993                     iparams[constr_iatomptr(ncon1,ia1,ia2,con)[0]].constr.dA);
 994           fprintf(debug," %d %5.3f\n",con,len);
 995         }
 996       }
 997       /* Limit the number of recursions to 1000*nc,
 998        * so a call does not take more than a second,
 999        * even for highly connected systems.
1000        */
1001       if (depth + 1 < nc && *count < 1000*nc) {
1002         if (ia[1] == at)
1003           a1 = ia[2];
1004         else
1005           a1 = ia[1];
1006         /* Recursion */
1007         path[depth] = con;
1008         constr_recur(at2con,ilist,iparams,
1009                      bTopB,a1,depth+1,nc,path,rn0,rn1,r2max,count);
1010         path[depth] = -1;
1011       }
1012     }
1013   }
1014 }
1015
1016 static real constr_r_max_moltype(FILE *fplog,
1017                                  gmx_moltype_t *molt,t_iparams *iparams,
1018                                  t_inputrec *ir)
1019 {
1020   int natoms,nflexcon,*path,at,count;
1021
1022   t_blocka at2con;
1023   real r0,r1,r2maxA,r2maxB,rmax,lam0,lam1;
1024
1025   if (molt->ilist[F_CONSTR].nr   == 0 &&
1026       molt->ilist[F_CONSTRNC].nr == 0) {
1027     return 0;
1028   }
1029
1030   natoms = molt->atoms.nr;
1031
1032   at2con = make_at2con(0,natoms,molt->ilist,iparams,
1033                        EI_DYNAMICS(ir->eI),&nflexcon);
1034   snew(path,1+ir->nProjOrder);
1035   for(at=0; at<1+ir->nProjOrder; at++)
1036     path[at] = -1;
1037
1038   r2maxA = 0;
1039   for(at=0; at<natoms; at++) {
1040     r0 = 0;
1041     r1 = 0;
1042
1043     count = 0;
1044     constr_recur(&at2con,molt->ilist,iparams,
1045                  FALSE,at,0,1+ir->nProjOrder,path,r0,r1,&r2maxA,&count);
1046   }
1047   if (ir->efep == efepNO) {
1048     rmax = sqrt(r2maxA);
1049   } else {
1050     r2maxB = 0;
1051     for(at=0; at<natoms; at++) {
1052       r0 = 0;
1053       r1 = 0;
1054       count = 0;
1055       constr_recur(&at2con,molt->ilist,iparams,
1056                    TRUE,at,0,1+ir->nProjOrder,path,r0,r1,&r2maxB,&count);
1057     }
1058     lam0 = ir->fepvals->init_lambda;
1059     if (EI_DYNAMICS(ir->eI))
1060       lam0 += ir->init_step*ir->fepvals->delta_lambda;
1061     rmax = (1 - lam0)*sqrt(r2maxA) + lam0*sqrt(r2maxB);
1062     if (EI_DYNAMICS(ir->eI)) {
1063       lam1 = ir->fepvals->init_lambda + (ir->init_step + ir->nsteps)*ir->fepvals->delta_lambda;
1064       rmax = max(rmax,(1 - lam1)*sqrt(r2maxA) + lam1*sqrt(r2maxB));
1065     }
1066   }
1067
1068   done_blocka(&at2con);
1069   sfree(path);
1070
1071   return rmax;
1072 }
1073
1074 real constr_r_max(FILE *fplog,gmx_mtop_t *mtop,t_inputrec *ir)
1075 {
1076   int mt;
1077   real rmax;
1078
1079   rmax = 0;
1080   for(mt=0; mt<mtop->nmoltype; mt++) {
1081     rmax = max(rmax,
1082                constr_r_max_moltype(fplog,&mtop->moltype[mt],
1083                                     mtop->ffparams.iparams,ir));
1084   }
1085
1086   if (fplog)
1087     fprintf(fplog,"Maximum distance for %d constraints, at 120 deg. angles, all-trans: %.3f nm\n",1+ir->nProjOrder,rmax);
1088
1089   return rmax;
1090 }
1091
1092 gmx_constr_t init_constraints(FILE *fplog,
1093                               gmx_mtop_t *mtop,t_inputrec *ir,
1094                               gmx_edsam_t ed,t_state *state,
1095                               t_commrec *cr)
1096 {
1097     int  ncon,nset,nmol,settle_type,i,natoms,mt,nflexcon;
1098     struct gmx_constr *constr;
1099     char *env;
1100     t_ilist *ilist;
1101     gmx_mtop_ilistloop_t iloop;
1102
1103     ncon =
1104         gmx_mtop_ftype_count(mtop,F_CONSTR) +
1105         gmx_mtop_ftype_count(mtop,F_CONSTRNC);
1106     nset = gmx_mtop_ftype_count(mtop,F_SETTLE);
1107
1108     if (ncon+nset == 0 && ir->ePull != epullCONSTRAINT && ed == NULL)
1109     {
1110         return NULL;
1111     }
1112
1113     snew(constr,1);
1114
1115     constr->ncon_tot = ncon;
1116     constr->nflexcon = 0;
1117     if (ncon > 0)
1118     {
1119         constr->n_at2con_mt = mtop->nmoltype;
1120         snew(constr->at2con_mt,constr->n_at2con_mt);
1121         for(mt=0; mt<mtop->nmoltype; mt++)
1122         {
1123             constr->at2con_mt[mt] = make_at2con(0,mtop->moltype[mt].atoms.nr,
1124                                                 mtop->moltype[mt].ilist,
1125                                                 mtop->ffparams.iparams,
1126                                                 EI_DYNAMICS(ir->eI),&nflexcon);
1127             for(i=0; i<mtop->nmolblock; i++)
1128             {
1129                 if (mtop->molblock[i].type == mt)
1130                 {
1131                     constr->nflexcon += mtop->molblock[i].nmol*nflexcon;
1132                 }
1133             }
1134         }
1135
1136         if (constr->nflexcon > 0)
1137         {
1138             if (fplog)
1139             {
1140                 fprintf(fplog,"There are %d flexible constraints\n",
1141                         constr->nflexcon);
1142                 if (ir->fc_stepsize == 0)
1143                 {
1144                     fprintf(fplog,"\n"
1145                             "WARNING: step size for flexible constraining = 0\n"
1146                             "         All flexible constraints will be rigid.\n"
1147                             "         Will try to keep all flexible constraints at their original length,\n"
1148                             "         but the lengths may exhibit some drift.\n\n");
1149                     constr->nflexcon = 0;
1150                 }
1151             }
1152             if (constr->nflexcon > 0)
1153             {
1154                 please_cite(fplog,"Hess2002");
1155             }
1156         }
1157
1158         if (ir->eConstrAlg == econtLINCS)
1159         {
1160             constr->lincsd = init_lincs(fplog,mtop,
1161                                         constr->nflexcon,constr->at2con_mt,
1162                                         DOMAINDECOMP(cr) && cr->dd->bInterCGcons,
1163                                         ir->nLincsIter,ir->nProjOrder);
1164         }
1165
1166         if (ir->eConstrAlg == econtSHAKE) {
1167             if (DOMAINDECOMP(cr) && cr->dd->bInterCGcons)
1168             {
1169                 gmx_fatal(FARGS,"SHAKE is not supported with domain decomposition and constraint that cross charge group boundaries, use LINCS");
1170             }
1171             if (constr->nflexcon)
1172             {
1173                 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");
1174             }
1175             please_cite(fplog,"Ryckaert77a");
1176             if (ir->bShakeSOR)
1177             {
1178                 please_cite(fplog,"Barth95a");
1179             }
1180
1181             constr->shaked = shake_init();
1182         }
1183     }
1184
1185     if (nset > 0) {
1186         please_cite(fplog,"Miyamoto92a");
1187
1188         constr->bInterCGsettles = inter_charge_group_settles(mtop);
1189
1190         /* Check that we have only one settle type */
1191         settle_type = -1;
1192         iloop = gmx_mtop_ilistloop_init(mtop);
1193         while (gmx_mtop_ilistloop_next(iloop,&ilist,&nmol))
1194         {
1195             for (i=0; i<ilist[F_SETTLE].nr; i+=4)
1196             {
1197                 if (settle_type == -1)
1198                 {
1199                     settle_type = ilist[F_SETTLE].iatoms[i];
1200                 }
1201                 else if (ilist[F_SETTLE].iatoms[i] != settle_type)
1202                 {
1203                     gmx_fatal(FARGS,
1204                               "The [molecules] section of your topology specifies more than one block of\n"
1205                               "a [moleculetype] with a [settles] block. Only one such is allowed. If you\n"
1206                               "are trying to partition your solvent into different *groups* (e.g. for\n"
1207                               "freezing, T-coupling, etc.) then you are using the wrong approach. Index\n"
1208                               "files specify groups. Otherwise, you may wish to change the least-used\n"
1209                               "block of molecules with SETTLE constraints into 3 normal constraints.");
1210                 }
1211             }
1212         }
1213
1214         constr->n_at2settle_mt = mtop->nmoltype;
1215         snew(constr->at2settle_mt,constr->n_at2settle_mt);
1216         for(mt=0; mt<mtop->nmoltype; mt++)
1217         {
1218             constr->at2settle_mt[mt] =
1219                 make_at2settle(mtop->moltype[mt].atoms.nr,
1220                                &mtop->moltype[mt].ilist[F_SETTLE]);
1221         }
1222     }
1223
1224     constr->maxwarn = 999;
1225     env = getenv("GMX_MAXCONSTRWARN");
1226     if (env)
1227     {
1228         constr->maxwarn = 0;
1229         sscanf(env,"%d",&constr->maxwarn);
1230         if (fplog)
1231         {
1232             fprintf(fplog,
1233                     "Setting the maximum number of constraint warnings to %d\n",
1234                     constr->maxwarn);
1235         }
1236         if (MASTER(cr))
1237         {
1238             fprintf(stderr,
1239                     "Setting the maximum number of constraint warnings to %d\n",
1240                     constr->maxwarn);
1241         }
1242     }
1243     if (constr->maxwarn < 0 && fplog)
1244     {
1245         fprintf(fplog,"maxwarn < 0, will not stop on constraint errors\n");
1246     }
1247     constr->warncount_lincs  = 0;
1248     constr->warncount_settle = 0;
1249
1250     /* Initialize the essential dynamics sampling.
1251      * Put the pointer to the ED struct in constr */
1252     constr->ed = ed;
1253     if (ed != NULL)
1254     {
1255         init_edsam(mtop,ir,cr,ed,state->x,state->box);
1256     }
1257
1258     constr->warn_mtop = mtop;
1259
1260     return constr;
1261 }
1262
1263 const t_blocka *atom2constraints_moltype(gmx_constr_t constr)
1264 {
1265   return constr->at2con_mt;
1266 }
1267
1268 const int **atom2settle_moltype(gmx_constr_t constr)
1269 {
1270     return (const int **)constr->at2settle_mt;
1271 }
1272
1273
1274 gmx_bool inter_charge_group_constraints(const gmx_mtop_t *mtop)
1275 {
1276     const gmx_moltype_t *molt;
1277     const t_block *cgs;
1278     const t_ilist *il;
1279     int  mb;
1280     int  nat,*at2cg,cg,a,ftype,i;
1281     gmx_bool bInterCG;
1282
1283     bInterCG = FALSE;
1284     for(mb=0; mb<mtop->nmolblock && !bInterCG; mb++)
1285     {
1286         molt = &mtop->moltype[mtop->molblock[mb].type];
1287
1288         if (molt->ilist[F_CONSTR].nr   > 0 ||
1289             molt->ilist[F_CONSTRNC].nr > 0 ||
1290             molt->ilist[F_SETTLE].nr > 0)
1291         {
1292             cgs  = &molt->cgs;
1293             snew(at2cg,molt->atoms.nr);
1294             for(cg=0; cg<cgs->nr; cg++)
1295             {
1296                 for(a=cgs->index[cg]; a<cgs->index[cg+1]; a++)
1297                     at2cg[a] = cg;
1298             }
1299
1300             for(ftype=F_CONSTR; ftype<=F_CONSTRNC; ftype++)
1301             {
1302                 il = &molt->ilist[ftype];
1303                 for(i=0; i<il->nr && !bInterCG; i+=1+NRAL(ftype))
1304                 {
1305                     if (at2cg[il->iatoms[i+1]] != at2cg[il->iatoms[i+2]])
1306                     {
1307                         bInterCG = TRUE;
1308                     }
1309                 }
1310             }
1311
1312             sfree(at2cg);
1313         }
1314     }
1315
1316     return bInterCG;
1317 }
1318
1319 gmx_bool inter_charge_group_settles(const gmx_mtop_t *mtop)
1320 {
1321     const gmx_moltype_t *molt;
1322     const t_block *cgs;
1323     const t_ilist *il;
1324     int  mb;
1325     int  nat,*at2cg,cg,a,ftype,i;
1326     gmx_bool bInterCG;
1327
1328     bInterCG = FALSE;
1329     for(mb=0; mb<mtop->nmolblock && !bInterCG; mb++)
1330     {
1331         molt = &mtop->moltype[mtop->molblock[mb].type];
1332
1333         if (molt->ilist[F_SETTLE].nr > 0)
1334         {
1335             cgs  = &molt->cgs;
1336             snew(at2cg,molt->atoms.nr);
1337             for(cg=0; cg<cgs->nr; cg++)
1338             {
1339                 for(a=cgs->index[cg]; a<cgs->index[cg+1]; a++)
1340                     at2cg[a] = cg;
1341             }
1342
1343             for(ftype=F_SETTLE; ftype<=F_SETTLE; ftype++)
1344             {
1345                 il = &molt->ilist[ftype];
1346                 for(i=0; i<il->nr && !bInterCG; i+=1+NRAL(F_SETTLE))
1347                 {
1348                     if (at2cg[il->iatoms[i+1]] != at2cg[il->iatoms[i+2]] ||
1349                         at2cg[il->iatoms[i+1]] != at2cg[il->iatoms[i+3]])
1350                     {
1351                         bInterCG = TRUE;
1352                     }
1353                 }
1354             }
1355
1356             sfree(at2cg);
1357         }
1358     }
1359
1360     return bInterCG;
1361 }
1362
1363 /* helper functions for andersen temperature control, because the
1364  * gmx_constr construct is only defined in constr.c. Return the list
1365  * of blocks (get_sblock) and the number of blocks (get_nblocks).  */
1366
1367 extern int *get_sblock(struct gmx_constr *constr)
1368 {
1369     return constr->sblock;
1370 }
1371
1372 extern int get_nblocks(struct gmx_constr *constr)
1373 {
1374     return constr->nblocks;
1375 }