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