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