src/gromacs/mdlib/constr.c

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