src/gromacs/mdlib/constr.c

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