src/gromacs/mdlib/constr.c

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