src/gromacs/mdlib/constr.cpp

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