src/gromacs/mdlib/constr.cpp

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