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