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