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