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