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