src/gromacs/mdlib/constr.c

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