src/gromacs/mdlib/constr.cpp

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