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