src/gromacs/gmxpreprocess/grompp.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,2017, 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 "grompp.h"
  40
  41 #include <errno.h>
  42 #include <limits.h>
  43 #include <string.h>
  44
  45 #include <cmath>
  46
  47 #include <algorithm>
  48
  49 #include <sys/types.h>
  50
  51 #include "gromacs/commandline/pargs.h"
  52 #include "gromacs/ewald/pme.h"
  53 #include "gromacs/fft/calcgrid.h"
  54 #include "gromacs/fileio/confio.h"
  55 #include "gromacs/fileio/enxio.h"
  56 #include "gromacs/fileio/tpxio.h"
  57 #include "gromacs/fileio/trxio.h"
  58 #include "gromacs/fileio/warninp.h"
  59 #include "gromacs/gmxpreprocess/add_par.h"
  60 #include "gromacs/gmxpreprocess/convparm.h"
  61 #include "gromacs/gmxpreprocess/gen_maxwell_velocities.h"
  62 #include "gromacs/gmxpreprocess/gpp_atomtype.h"
  63 #include "gromacs/gmxpreprocess/grompp-impl.h"
  64 #include "gromacs/gmxpreprocess/notset.h"
  65 #include "gromacs/gmxpreprocess/readir.h"
  66 #include "gromacs/gmxpreprocess/tomorse.h"
  67 #include "gromacs/gmxpreprocess/topio.h"
  68 #include "gromacs/gmxpreprocess/toputil.h"
  69 #include "gromacs/gmxpreprocess/vsite_parm.h"
  70 #include "gromacs/imd/imd.h"
  71 #include "gromacs/math/functions.h"
  72 #include "gromacs/math/invertmatrix.h"
  73 #include "gromacs/math/units.h"
  74 #include "gromacs/math/vec.h"
  75 #include "gromacs/mdlib/calc_verletbuf.h"
  76 #include "gromacs/mdlib/compute_io.h"
  77 #include "gromacs/mdlib/constr.h"
  78 #include "gromacs/mdlib/genborn.h"
  79 #include "gromacs/mdlib/perf_est.h"
  80 #include "gromacs/mdrunutility/mdmodules.h"
  81 #include "gromacs/mdtypes/inputrec.h"
  82 #include "gromacs/mdtypes/md_enums.h"
  83 #include "gromacs/mdtypes/nblist.h"
  84 #include "gromacs/mdtypes/state.h"
  85 #include "gromacs/pbcutil/boxutilities.h"
  86 #include "gromacs/pbcutil/pbc.h"
  87 #include "gromacs/pulling/pull.h"
  88 #include "gromacs/random/seed.h"
  89 #include "gromacs/topology/ifunc.h"
  90 #include "gromacs/topology/mtop_util.h"
  91 #include "gromacs/topology/symtab.h"
  92 #include "gromacs/topology/topology.h"
  93 #include "gromacs/trajectory/trajectoryframe.h"
  94 #include "gromacs/utility/arraysize.h"
  95 #include "gromacs/utility/cstringutil.h"
  96 #include "gromacs/utility/exceptions.h"
  97 #include "gromacs/utility/fatalerror.h"
  98 #include "gromacs/utility/futil.h"
  99 #include "gromacs/utility/gmxassert.h"
 100 #include "gromacs/utility/smalloc.h"
 101 #include "gromacs/utility/snprintf.h"
 102
 103 static int rm_interactions(int ifunc, int nrmols, t_molinfo mols[])
 104 {
 105     int  i, n;
 106
 107     n = 0;
 108     /* For all the molecule types */
 109     for (i = 0; i < nrmols; i++)
 110     {
 111         n += mols[i].plist[ifunc].nr;
 112         mols[i].plist[ifunc].nr = 0;
 113     }
 114     return n;
 115 }
 116
 117 static int check_atom_names(const char *fn1, const char *fn2,
 118                             gmx_mtop_t *mtop, const t_atoms *at)
 119 {
 120     int      mb, m, i, j, nmismatch;
 121     t_atoms *tat;
 122 #define MAXMISMATCH 20
 123
 124     if (mtop->natoms != at->nr)
 125     {
 126         gmx_incons("comparing atom names");
 127     }
 128
 129     nmismatch = 0;
 130     i         = 0;
 131     for (mb = 0; mb < mtop->nmolblock; mb++)
 132     {
 133         tat = &mtop->moltype[mtop->molblock[mb].type].atoms;
 134         for (m = 0; m < mtop->molblock[mb].nmol; m++)
 135         {
 136             for (j = 0; j < tat->nr; j++)
 137             {
 138                 if (strcmp( *(tat->atomname[j]), *(at->atomname[i]) ) != 0)
 139                 {
 140                     if (nmismatch < MAXMISMATCH)
 141                     {
 142                         fprintf(stderr,
 143                                 "Warning: atom name %d in %s and %s does not match (%s - %s)\n",
 144                                 i+1, fn1, fn2, *(tat->atomname[j]), *(at->atomname[i]));
 145                     }
 146                     else if (nmismatch == MAXMISMATCH)
 147                     {
 148                         fprintf(stderr, "(more than %d non-matching atom names)\n", MAXMISMATCH);
 149                     }
 150                     nmismatch++;
 151                 }
 152                 i++;
 153             }
 154         }
 155     }
 156
 157     return nmismatch;
 158 }
 159
 160 static void check_eg_vs_cg(gmx_mtop_t *mtop)
 161 {
 162     int            astart, mb, m, cg, j, firstj;
 163     unsigned char  firsteg, eg;
 164     gmx_moltype_t *molt;
 165
 166     /* Go through all the charge groups and make sure all their
 167      * atoms are in the same energy group.
 168      */
 169
 170     astart = 0;
 171     for (mb = 0; mb < mtop->nmolblock; mb++)
 172     {
 173         molt = &mtop->moltype[mtop->molblock[mb].type];
 174         for (m = 0; m < mtop->molblock[mb].nmol; m++)
 175         {
 176             for (cg = 0; cg < molt->cgs.nr; cg++)
 177             {
 178                 /* Get the energy group of the first atom in this charge group */
 179                 firstj  = astart + molt->cgs.index[cg];
 180                 firsteg = ggrpnr(&mtop->groups, egcENER, firstj);
 181                 for (j = molt->cgs.index[cg]+1; j < molt->cgs.index[cg+1]; j++)
 182                 {
 183                     eg = ggrpnr(&mtop->groups, egcENER, astart+j);
 184                     if (eg != firsteg)
 185                     {
 186                         gmx_fatal(FARGS, "atoms %d and %d in charge group %d of molecule type '%s' are in different energy groups",
 187                                   firstj+1, astart+j+1, cg+1, *molt->name);
 188                     }
 189                 }
 190             }
 191             astart += molt->atoms.nr;
 192         }
 193     }
 194 }
 195
 196 static void check_cg_sizes(const char *topfn, t_block *cgs, warninp_t wi)
 197 {
 198     int  maxsize, cg;
 199     char warn_buf[STRLEN];
 200
 201     maxsize = 0;
 202     for (cg = 0; cg < cgs->nr; cg++)
 203     {
 204         maxsize = std::max(maxsize, cgs->index[cg+1]-cgs->index[cg]);
 205     }
 206
 207     if (maxsize > MAX_CHARGEGROUP_SIZE)
 208     {
 209         gmx_fatal(FARGS, "The largest charge group contains %d atoms. The maximum is %d.", maxsize, MAX_CHARGEGROUP_SIZE);
 210     }
 211     else if (maxsize > 10)
 212     {
 213         set_warning_line(wi, topfn, -1);
 214         sprintf(warn_buf,
 215                 "The largest charge group contains %d atoms.\n"
 216                 "Since atoms only see each other when the centers of geometry of the charge groups they belong to are within the cut-off distance, too large charge groups can lead to serious cut-off artifacts.\n"
 217                 "For efficiency and accuracy, charge group should consist of a few atoms.\n"
 218                 "For all-atom force fields use: CH3, CH2, CH, NH2, NH, OH, CO2, CO, etc.",
 219                 maxsize);
 220         warning_note(wi, warn_buf);
 221     }
 222 }
 223
 224 static void check_bonds_timestep(gmx_mtop_t *mtop, double dt, warninp_t wi)
 225 {
 226     /* This check is not intended to ensure accurate integration,
 227      * rather it is to signal mistakes in the mdp settings.
 228      * A common mistake is to forget to turn on constraints
 229      * for MD after energy minimization with flexible bonds.
 230      * This check can also detect too large time steps for flexible water
 231      * models, but such errors will often be masked by the constraints
 232      * mdp options, which turns flexible water into water with bond constraints,
 233      * but without an angle constraint. Unfortunately such incorrect use
 234      * of water models can not easily be detected without checking
 235      * for specific model names.
 236      *
 237      * The stability limit of leap-frog or velocity verlet is 4.44 steps
 238      * per oscillational period.
 239      * But accurate bonds distributions are lost far before that limit.
 240      * To allow relatively common schemes (although not common with Gromacs)
 241      * of dt=1 fs without constraints and dt=2 fs with only H-bond constraints
 242      * we set the note limit to 10.
 243      */
 244     int            min_steps_warn = 5;
 245     int            min_steps_note = 10;
 246     t_iparams     *ip;
 247     int            molt;
 248     gmx_moltype_t *moltype, *w_moltype;
 249     t_atom        *atom;
 250     t_ilist       *ilist, *ilb, *ilc, *ils;
 251     int            ftype;
 252     int            i, a1, a2, w_a1, w_a2, j;
 253     real           twopi2, limit2, fc, re, m1, m2, period2, w_period2;
 254     gmx_bool       bFound, bWater, bWarn;
 255     char           warn_buf[STRLEN];
 256
 257     ip = mtop->ffparams.iparams;
 258
 259     twopi2 = gmx::square(2*M_PI);
 260
 261     limit2 = gmx::square(min_steps_note*dt);
 262
 263     w_a1      = w_a2 = -1;
 264     w_period2 = -1.0;
 265
 266     w_moltype = nullptr;
 267     for (molt = 0; molt < mtop->nmoltype; molt++)
 268     {
 269         moltype = &mtop->moltype[molt];
 270         atom    = moltype->atoms.atom;
 271         ilist   = moltype->ilist;
 272         ilc     = &ilist[F_CONSTR];
 273         ils     = &ilist[F_SETTLE];
 274         for (ftype = 0; ftype < F_NRE; ftype++)
 275         {
 276             if (!(ftype == F_BONDS || ftype == F_G96BONDS || ftype == F_HARMONIC))
 277             {
 278                 continue;
 279             }
 280
 281             ilb = &ilist[ftype];
 282             for (i = 0; i < ilb->nr; i += 3)
 283             {
 284                 fc = ip[ilb->iatoms[i]].harmonic.krA;
 285                 re = ip[ilb->iatoms[i]].harmonic.rA;
 286                 if (ftype == F_G96BONDS)
 287                 {
 288                     /* Convert squared sqaure fc to harmonic fc */
 289                     fc = 2*fc*re;
 290                 }
 291                 a1 = ilb->iatoms[i+1];
 292                 a2 = ilb->iatoms[i+2];
 293                 m1 = atom[a1].m;
 294                 m2 = atom[a2].m;
 295                 if (fc > 0 && m1 > 0 && m2 > 0)
 296                 {
 297                     period2 = twopi2*m1*m2/((m1 + m2)*fc);
 298                 }
 299                 else
 300                 {
 301                     period2 = GMX_FLOAT_MAX;
 302                 }
 303                 if (debug)
 304                 {
 305                     fprintf(debug, "fc %g m1 %g m2 %g period %g\n",
 306                             fc, m1, m2, std::sqrt(period2));
 307                 }
 308                 if (period2 < limit2)
 309                 {
 310                     bFound = FALSE;
 311                     for (j = 0; j < ilc->nr; j += 3)
 312                     {
 313                         if ((ilc->iatoms[j+1] == a1 && ilc->iatoms[j+2] == a2) ||
 314                             (ilc->iatoms[j+1] == a2 && ilc->iatoms[j+2] == a1))
 315                         {
 316                             bFound = TRUE;
 317                         }
 318                     }
 319                     for (j = 0; j < ils->nr; j += 4)
 320                     {
 321                         if ((a1 == ils->iatoms[j+1] || a1 == ils->iatoms[j+2] || a1 == ils->iatoms[j+3]) &&
 322                             (a2 == ils->iatoms[j+1] || a2 == ils->iatoms[j+2] || a2 == ils->iatoms[j+3]))
 323                         {
 324                             bFound = TRUE;
 325                         }
 326                     }
 327                     if (!bFound &&
 328                         (w_moltype == nullptr || period2 < w_period2))
 329                     {
 330                         w_moltype = moltype;
 331                         w_a1      = a1;
 332                         w_a2      = a2;
 333                         w_period2 = period2;
 334                     }
 335                 }
 336             }
 337         }
 338     }
 339
 340     if (w_moltype != nullptr)
 341     {
 342         bWarn = (w_period2 < gmx::square(min_steps_warn*dt));
 343         /* A check that would recognize most water models */
 344         bWater = ((*w_moltype->atoms.atomname[0])[0] == 'O' &&
 345                   w_moltype->atoms.nr <= 5);
 346         sprintf(warn_buf, "The bond in molecule-type %s between atoms %d %s and %d %s has an estimated oscillational period of %.1e ps, which is less than %d times the time step of %.1e ps.\n"
 347                 "%s",
 348                 *w_moltype->name,
 349                 w_a1+1, *w_moltype->atoms.atomname[w_a1],
 350                 w_a2+1, *w_moltype->atoms.atomname[w_a2],
 351                 std::sqrt(w_period2), bWarn ? min_steps_warn : min_steps_note, dt,
 352                 bWater ?
 353                 "Maybe you asked for fexible water." :
 354                 "Maybe you forgot to change the constraints mdp option.");
 355         if (bWarn)
 356         {
 357             warning(wi, warn_buf);
 358         }
 359         else
 360         {
 361             warning_note(wi, warn_buf);
 362         }
 363     }
 364 }
 365
 366 static void check_vel(gmx_mtop_t *mtop, rvec v[])
 367 {
 368     gmx_mtop_atomloop_all_t aloop;
 369     const t_atom           *atom;
 370     int                     a;
 371
 372     aloop = gmx_mtop_atomloop_all_init(mtop);
 373     while (gmx_mtop_atomloop_all_next(aloop, &a, &atom))
 374     {
 375         if (atom->ptype == eptShell ||
 376             atom->ptype == eptBond  ||
 377             atom->ptype == eptVSite)
 378         {
 379             clear_rvec(v[a]);
 380         }
 381     }
 382 }
 383
 384 static void check_shells_inputrec(gmx_mtop_t *mtop,
 385                                   t_inputrec *ir,
 386                                   warninp_t   wi)
 387 {
 388     gmx_mtop_atomloop_all_t aloop;
 389     const t_atom           *atom;
 390     int                     a, nshells = 0;
 391     char                    warn_buf[STRLEN];
 392
 393     aloop = gmx_mtop_atomloop_all_init(mtop);
 394     while (gmx_mtop_atomloop_all_next(aloop, &a, &atom))
 395     {
 396         if (atom->ptype == eptShell ||
 397             atom->ptype == eptBond)
 398         {
 399             nshells++;
 400         }
 401     }
 402     if ((nshells > 0) && (ir->nstcalcenergy != 1))
 403     {
 404         set_warning_line(wi, "unknown", -1);
 405         snprintf(warn_buf, STRLEN,
 406                  "There are %d shells, changing nstcalcenergy from %d to 1",
 407                  nshells, ir->nstcalcenergy);
 408         ir->nstcalcenergy = 1;
 409         warning(wi, warn_buf);
 410     }
 411 }
 412
 413 /* TODO Decide whether this function can be consolidated with
 414  * gmx_mtop_ftype_count */
 415 static gmx_bool nint_ftype(gmx_mtop_t *mtop, t_molinfo *mi, int ftype)
 416 {
 417     int nint, mb;
 418
 419     nint = 0;
 420     for (mb = 0; mb < mtop->nmolblock; mb++)
 421     {
 422         nint += mtop->molblock[mb].nmol*mi[mtop->molblock[mb].type].plist[ftype].nr;
 423     }
 424
 425     return nint;
 426 }
 427
 428 /* This routine reorders the molecule type array
 429  * in the order of use in the molblocks,
 430  * unused molecule types are deleted.
 431  */
 432 static void renumber_moltypes(gmx_mtop_t *sys,
 433                               int *nmolinfo, t_molinfo **molinfo)
 434 {
 435     int       *order, norder, i;
 436     int        mb, mi;
 437     t_molinfo *minew;
 438
 439     snew(order, *nmolinfo);
 440     norder = 0;
 441     for (mb = 0; mb < sys->nmolblock; mb++)
 442     {
 443         for (i = 0; i < norder; i++)
 444         {
 445             if (order[i] == sys->molblock[mb].type)
 446             {
 447                 break;
 448             }
 449         }
 450         if (i == norder)
 451         {
 452             /* This type did not occur yet, add it */
 453             order[norder] = sys->molblock[mb].type;
 454             /* Renumber the moltype in the topology */
 455             norder++;
 456         }
 457         sys->molblock[mb].type = i;
 458     }
 459
 460     /* We still need to reorder the molinfo structs */
 461     snew(minew, norder);
 462     for (mi = 0; mi < *nmolinfo; mi++)
 463     {
 464         for (i = 0; i < norder; i++)
 465         {
 466             if (order[i] == mi)
 467             {
 468                 break;
 469             }
 470         }
 471         if (i == norder)
 472         {
 473             done_mi(&(*molinfo)[mi]);
 474         }
 475         else
 476         {
 477             minew[i] = (*molinfo)[mi];
 478         }
 479     }
 480     sfree(*molinfo);
 481
 482     *nmolinfo = norder;
 483     *molinfo  = minew;
 484 }
 485
 486 static void molinfo2mtop(int nmi, t_molinfo *mi, gmx_mtop_t *mtop)
 487 {
 488     int            m;
 489     gmx_moltype_t *molt;
 490
 491     mtop->nmoltype = nmi;
 492     snew(mtop->moltype, nmi);
 493     for (m = 0; m < nmi; m++)
 494     {
 495         molt        = &mtop->moltype[m];
 496         molt->name  = mi[m].name;
 497         molt->atoms = mi[m].atoms;
 498         /* ilists are copied later */
 499         molt->cgs   = mi[m].cgs;
 500         molt->excls = mi[m].excls;
 501     }
 502 }
 503
 504 static void
 505 new_status(const char *topfile, const char *topppfile, const char *confin,
 506            t_gromppopts *opts, t_inputrec *ir, gmx_bool bZero,
 507            gmx_bool bGenVel, gmx_bool bVerbose, t_state *state,
 508            gpp_atomtype_t atype, gmx_mtop_t *sys,
 509            int *nmi, t_molinfo **mi, t_molinfo **intermolecular_interactions,
 510            t_params plist[],
 511            int *comb, double *reppow, real *fudgeQQ,
 512            gmx_bool bMorse,
 513            warninp_t wi)
 514 {
 515     t_molinfo      *molinfo = nullptr;
 516     int             nmolblock;
 517     gmx_molblock_t *molblock, *molbs;
 518     int             mb, i, nrmols, nmismatch;
 519     char            buf[STRLEN];
 520     gmx_bool        bGB = FALSE;
 521     char            warn_buf[STRLEN];
 522
 523     init_mtop(sys);
 524
 525     /* Set gmx_boolean for GB */
 526     if (ir->implicit_solvent)
 527     {
 528         bGB = TRUE;
 529     }
 530
 531     /* TOPOLOGY processing */
 532     sys->name = do_top(bVerbose, topfile, topppfile, opts, bZero, &(sys->symtab),
 533                        plist, comb, reppow, fudgeQQ,
 534                        atype, &nrmols, &molinfo, intermolecular_interactions,
 535                        ir,
 536                        &nmolblock, &molblock, bGB,
 537                        wi);
 538
 539     sys->nmolblock = 0;
 540     snew(sys->molblock, nmolblock);
 541
 542     sys->natoms = 0;
 543     for (mb = 0; mb < nmolblock; mb++)
 544     {
 545         if (sys->nmolblock > 0 &&
 546             molblock[mb].type == sys->molblock[sys->nmolblock-1].type)
 547         {
 548             /* Merge consecutive blocks with the same molecule type */
 549             sys->molblock[sys->nmolblock-1].nmol += molblock[mb].nmol;
 550             sys->natoms += molblock[mb].nmol*sys->molblock[sys->nmolblock-1].natoms_mol;
 551         }
 552         else if (molblock[mb].nmol > 0)
 553         {
 554             /* Add a new molblock to the topology */
 555             molbs             = &sys->molblock[sys->nmolblock];
 556             *molbs            = molblock[mb];
 557             molbs->natoms_mol = molinfo[molbs->type].atoms.nr;
 558             molbs->nposres_xA = 0;
 559             molbs->nposres_xB = 0;
 560             sys->natoms      += molbs->nmol*molbs->natoms_mol;
 561             sys->nmolblock++;
 562         }
 563     }
 564     if (sys->nmolblock == 0)
 565     {
 566         gmx_fatal(FARGS, "No molecules were defined in the system");
 567     }
 568
 569     renumber_moltypes(sys, &nrmols, &molinfo);
 570
 571     if (bMorse)
 572     {
 573         convert_harmonics(nrmols, molinfo, atype);
 574     }
 575
 576     if (ir->eDisre == edrNone)
 577     {
 578         i = rm_interactions(F_DISRES, nrmols, molinfo);
 579         if (i > 0)
 580         {
 581             set_warning_line(wi, "unknown", -1);
 582             sprintf(warn_buf, "disre = no, removed %d distance restraints", i);
 583             warning_note(wi, warn_buf);
 584         }
 585     }
 586     if (opts->bOrire == FALSE)
 587     {
 588         i = rm_interactions(F_ORIRES, nrmols, molinfo);
 589         if (i > 0)
 590         {
 591             set_warning_line(wi, "unknown", -1);
 592             sprintf(warn_buf, "orire = no, removed %d orientation restraints", i);
 593             warning_note(wi, warn_buf);
 594         }
 595     }
 596
 597     /* Copy structures from msys to sys */
 598     molinfo2mtop(nrmols, molinfo, sys);
 599
 600     gmx_mtop_finalize(sys);
 601
 602     /* COORDINATE file processing */
 603     if (bVerbose)
 604     {
 605         fprintf(stderr, "processing coordinates...\n");
 606     }
 607
 608     t_topology *conftop;
 609     rvec       *x = nullptr;
 610     rvec       *v = nullptr;
 611     snew(conftop, 1);
 612     read_tps_conf(confin, conftop, nullptr, &x, &v, state->box, FALSE);
 613     state->natoms = conftop->atoms.nr;
 614     if (state->natoms != sys->natoms)
 615     {
 616         gmx_fatal(FARGS, "number of coordinates in coordinate file (%s, %d)\n"
 617                   "             does not match topology (%s, %d)",
 618                   confin, state->natoms, topfile, sys->natoms);
 619     }
 620     /* It would be nice to get rid of the copies below, but we don't know
 621      * a priori if the number of atoms in confin matches what we expect.
 622      */
 623     state->flags |= (1 << estX);
 624     if (v != NULL)
 625     {
 626         state->flags |= (1 << estV);
 627     }
 628     state_change_natoms(state, state->natoms);
 629     for (int i = 0; i < state->natoms; i++)
 630     {
 631         copy_rvec(x[i], state->x[i]);
 632     }
 633     sfree(x);
 634     if (v != nullptr)
 635     {
 636         for (int i = 0; i < state->natoms; i++)
 637         {
 638             copy_rvec(v[i], state->v[i]);
 639         }
 640         sfree(v);
 641     }
 642     /* This call fixes the box shape for runs with pressure scaling */
 643     set_box_rel(ir, state);
 644
 645     nmismatch = check_atom_names(topfile, confin, sys, &conftop->atoms);
 646     done_top(conftop);
 647     sfree(conftop);
 648
 649     if (nmismatch)
 650     {
 651         sprintf(buf, "%d non-matching atom name%s\n"
 652                 "atom names from %s will be used\n"
 653                 "atom names from %s will be ignored\n",
 654                 nmismatch, (nmismatch == 1) ? "" : "s", topfile, confin);
 655         warning(wi, buf);
 656     }
 657
 658     /* Do more checks, mostly related to constraints */
 659     if (bVerbose)
 660     {
 661         fprintf(stderr, "double-checking input for internal consistency...\n");
 662     }
 663     {
 664         int bHasNormalConstraints = 0 < (nint_ftype(sys, molinfo, F_CONSTR) +
 665                                          nint_ftype(sys, molinfo, F_CONSTRNC));
 666         int bHasAnyConstraints = bHasNormalConstraints || 0 < nint_ftype(sys, molinfo, F_SETTLE);
 667         double_check(ir, state->box,
 668                      bHasNormalConstraints,
 669                      bHasAnyConstraints,
 670                      wi);
 671     }
 672
 673     if (bGenVel)
 674     {
 675         real                   *mass;
 676         gmx_mtop_atomloop_all_t aloop;
 677         const t_atom           *atom;
 678
 679         snew(mass, state->natoms);
 680         aloop = gmx_mtop_atomloop_all_init(sys);
 681         while (gmx_mtop_atomloop_all_next(aloop, &i, &atom))
 682         {
 683             mass[i] = atom->m;
 684         }
 685
 686         if (opts->seed == -1)
 687         {
 688             opts->seed = static_cast<int>(gmx::makeRandomSeed());
 689             fprintf(stderr, "Setting gen_seed to %d\n", opts->seed);
 690         }
 691         state->flags |= (1 << estV);
 692         maxwell_speed(opts->tempi, opts->seed, sys, as_rvec_array(state->v.data()));
 693
 694         stop_cm(stdout, state->natoms, mass, as_rvec_array(state->x.data()), as_rvec_array(state->v.data()));
 695         sfree(mass);
 696     }
 697
 698     *nmi = nrmols;
 699     *mi  = molinfo;
 700 }
 701
 702 static void copy_state(const char *slog, t_trxframe *fr,
 703                        gmx_bool bReadVel, t_state *state,
 704                        double *use_time)
 705 {
 706     int i;
 707
 708     if (fr->not_ok & FRAME_NOT_OK)
 709     {
 710         gmx_fatal(FARGS, "Can not start from an incomplete frame");
 711     }
 712     if (!fr->bX)
 713     {
 714         gmx_fatal(FARGS, "Did not find a frame with coordinates in file %s",
 715                   slog);
 716     }
 717
 718     for (i = 0; i < state->natoms; i++)
 719     {
 720         copy_rvec(fr->x[i], state->x[i]);
 721     }
 722     if (bReadVel)
 723     {
 724         if (!fr->bV)
 725         {
 726             gmx_incons("Trajecory frame unexpectedly does not contain velocities");
 727         }
 728         for (i = 0; i < state->natoms; i++)
 729         {
 730             copy_rvec(fr->v[i], state->v[i]);
 731         }
 732     }
 733     if (fr->bBox)
 734     {
 735         copy_mat(fr->box, state->box);
 736     }
 737
 738     *use_time = fr->time;
 739 }
 740
 741 static void cont_status(const char *slog, const char *ener,
 742                         gmx_bool bNeedVel, gmx_bool bGenVel, real fr_time,
 743                         t_inputrec *ir, t_state *state,
 744                         gmx_mtop_t *sys,
 745                         const gmx_output_env_t *oenv)
 746 /* If fr_time == -1 read the last frame available which is complete */
 747 {
 748     gmx_bool     bReadVel;
 749     t_trxframe   fr;
 750     t_trxstatus *fp;
 751     int          i;
 752     double       use_time;
 753
 754     bReadVel = (bNeedVel && !bGenVel);
 755
 756     fprintf(stderr,
 757             "Reading Coordinates%s and Box size from old trajectory\n",
 758             bReadVel ? ", Velocities" : "");
 759     if (fr_time == -1)
 760     {
 761         fprintf(stderr, "Will read whole trajectory\n");
 762     }
 763     else
 764     {
 765         fprintf(stderr, "Will read till time %g\n", fr_time);
 766     }
 767     if (!bReadVel)
 768     {
 769         if (bGenVel)
 770         {
 771             fprintf(stderr, "Velocities generated: "
 772                     "ignoring velocities in input trajectory\n");
 773         }
 774         read_first_frame(oenv, &fp, slog, &fr, TRX_NEED_X);
 775     }
 776     else
 777     {
 778         read_first_frame(oenv, &fp, slog, &fr, TRX_NEED_X | TRX_NEED_V);
 779
 780         if (!fr.bV)
 781         {
 782             fprintf(stderr,
 783                     "\n"
 784                     "WARNING: Did not find a frame with velocities in file %s,\n"
 785                     "         all velocities will be set to zero!\n\n", slog);
 786             for (i = 0; i < sys->natoms; i++)
 787             {
 788                 clear_rvec(state->v[i]);
 789             }
 790             close_trx(fp);
 791             /* Search for a frame without velocities */
 792             bReadVel = FALSE;
 793             read_first_frame(oenv, &fp, slog, &fr, TRX_NEED_X);
 794         }
 795     }
 796
 797     state->natoms = fr.natoms;
 798
 799     if (sys->natoms != state->natoms)
 800     {
 801         gmx_fatal(FARGS, "Number of atoms in Topology "
 802                   "is not the same as in Trajectory");
 803     }
 804     copy_state(slog, &fr, bReadVel, state, &use_time);
 805
 806     /* Find the appropriate frame */
 807     while ((fr_time == -1 || fr.time < fr_time) &&
 808            read_next_frame(oenv, fp, &fr))
 809     {
 810         copy_state(slog, &fr, bReadVel, state, &use_time);
 811     }
 812
 813     close_trx(fp);
 814
 815     /* Set the relative box lengths for preserving the box shape.
 816      * Note that this call can lead to differences in the last bit
 817      * with respect to using gmx convert-tpr to create a [REF].tpx[ref] file.
 818      */
 819     set_box_rel(ir, state);
 820
 821     fprintf(stderr, "Using frame at t = %g ps\n", use_time);
 822     fprintf(stderr, "Starting time for run is %g ps\n", ir->init_t);
 823
 824     if ((ir->epc != epcNO  || ir->etc == etcNOSEHOOVER) && ener)
 825     {
 826         get_enx_state(ener, use_time, &sys->groups, ir, state);
 827         preserve_box_shape(ir, state->box_rel, state->boxv);
 828     }
 829 }
 830
 831 static void read_posres(gmx_mtop_t *mtop, t_molinfo *molinfo, gmx_bool bTopB,
 832                         const char *fn,
 833                         int rc_scaling, int ePBC,
 834                         rvec com,
 835                         warninp_t wi)
 836 {
 837     gmx_bool       *hadAtom;
 838     rvec           *x, *v, *xp;
 839     dvec            sum;
 840     double          totmass;
 841     t_topology     *top;
 842     matrix          box, invbox;
 843     int             natoms, npbcdim = 0;
 844     char            warn_buf[STRLEN];
 845     int             a, i, ai, j, k, mb, nat_molb;
 846     gmx_molblock_t *molb;
 847     t_params       *pr, *prfb;
 848     t_atom         *atom;
 849
 850     snew(top, 1);
 851     read_tps_conf(fn, top, nullptr, &x, &v, box, FALSE);
 852     natoms = top->atoms.nr;
 853     done_top(top);
 854     sfree(top);
 855     if (natoms != mtop->natoms)
 856     {
 857         sprintf(warn_buf, "The number of atoms in %s (%d) does not match the number of atoms in the topology (%d). Will assume that the first %d atoms in the topology and %s match.", fn, natoms, mtop->natoms, std::min(mtop->natoms, natoms), fn);
 858         warning(wi, warn_buf);
 859     }
 860
 861     npbcdim = ePBC2npbcdim(ePBC);
 862     clear_rvec(com);
 863     if (rc_scaling != erscNO)
 864     {
 865         copy_mat(box, invbox);
 866         for (j = npbcdim; j < DIM; j++)
 867         {
 868             clear_rvec(invbox[j]);
 869             invbox[j][j] = 1;
 870         }
 871         gmx::invertBoxMatrix(invbox, invbox);
 872     }
 873
 874     /* Copy the reference coordinates to mtop */
 875     clear_dvec(sum);
 876     totmass = 0;
 877     a       = 0;
 878     snew(hadAtom, natoms);
 879     for (mb = 0; mb < mtop->nmolblock; mb++)
 880     {
 881         molb     = &mtop->molblock[mb];
 882         nat_molb = molb->nmol*mtop->moltype[molb->type].atoms.nr;
 883         pr       = &(molinfo[molb->type].plist[F_POSRES]);
 884         prfb     = &(molinfo[molb->type].plist[F_FBPOSRES]);
 885         if (pr->nr > 0 || prfb->nr > 0)
 886         {
 887             atom = mtop->moltype[molb->type].atoms.atom;
 888             for (i = 0; (i < pr->nr); i++)
 889             {
 890                 ai = pr->param[i].ai();
 891                 if (ai >= natoms)
 892                 {
 893                     gmx_fatal(FARGS, "Position restraint atom index (%d) in moltype '%s' is larger than number of atoms in %s (%d).\n",
 894                               ai+1, *molinfo[molb->type].name, fn, natoms);
 895                 }
 896                 hadAtom[ai] = TRUE;
 897                 if (rc_scaling == erscCOM)
 898                 {
 899                     /* Determine the center of mass of the posres reference coordinates */
 900                     for (j = 0; j < npbcdim; j++)
 901                     {
 902                         sum[j] += atom[ai].m*x[a+ai][j];
 903                     }
 904                     totmass  += atom[ai].m;
 905                 }
 906             }
 907             /* Same for flat-bottomed posres, but do not count an atom twice for COM */
 908             for (i = 0; (i < prfb->nr); i++)
 909             {
 910                 ai = prfb->param[i].ai();
 911                 if (ai >= natoms)
 912                 {
 913                     gmx_fatal(FARGS, "Position restraint atom index (%d) in moltype '%s' is larger than number of atoms in %s (%d).\n",
 914                               ai+1, *molinfo[molb->type].name, fn, natoms);
 915                 }
 916                 if (rc_scaling == erscCOM && hadAtom[ai] == FALSE)
 917                 {
 918                     /* Determine the center of mass of the posres reference coordinates */
 919                     for (j = 0; j < npbcdim; j++)
 920                     {
 921                         sum[j] += atom[ai].m*x[a+ai][j];
 922                     }
 923                     totmass  += atom[ai].m;
 924                 }
 925             }
 926             if (!bTopB)
 927             {
 928                 molb->nposres_xA = nat_molb;
 929                 snew(molb->posres_xA, molb->nposres_xA);
 930                 for (i = 0; i < nat_molb; i++)
 931                 {
 932                     copy_rvec(x[a+i], molb->posres_xA[i]);
 933                 }
 934             }
 935             else
 936             {
 937                 molb->nposres_xB = nat_molb;
 938                 snew(molb->posres_xB, molb->nposres_xB);
 939                 for (i = 0; i < nat_molb; i++)
 940                 {
 941                     copy_rvec(x[a+i], molb->posres_xB[i]);
 942                 }
 943             }
 944         }
 945         a += nat_molb;
 946     }
 947     if (rc_scaling == erscCOM)
 948     {
 949         if (totmass == 0)
 950         {
 951             gmx_fatal(FARGS, "The total mass of the position restraint atoms is 0");
 952         }
 953         for (j = 0; j < npbcdim; j++)
 954         {
 955             com[j] = sum[j]/totmass;
 956         }
 957         fprintf(stderr, "The center of mass of the position restraint coord's is %6.3f %6.3f %6.3f\n", com[XX], com[YY], com[ZZ]);
 958     }
 959
 960     if (rc_scaling != erscNO)
 961     {
 962         GMX_ASSERT(npbcdim <= DIM, "Only DIM dimensions can have PBC");
 963
 964         for (mb = 0; mb < mtop->nmolblock; mb++)
 965         {
 966             molb     = &mtop->molblock[mb];
 967             nat_molb = molb->nmol*mtop->moltype[molb->type].atoms.nr;
 968             if (molb->nposres_xA > 0 || molb->nposres_xB > 0)
 969             {
 970                 xp = (!bTopB ? molb->posres_xA : molb->posres_xB);
 971                 for (i = 0; i < nat_molb; i++)
 972                 {
 973                     for (j = 0; j < npbcdim; j++)
 974                     {
 975                         if (rc_scaling == erscALL)
 976                         {
 977                             /* Convert from Cartesian to crystal coordinates */
 978                             xp[i][j] *= invbox[j][j];
 979                             for (k = j+1; k < npbcdim; k++)
 980                             {
 981                                 xp[i][j] += invbox[k][j]*xp[i][k];
 982                             }
 983                         }
 984                         else if (rc_scaling == erscCOM)
 985                         {
 986                             /* Subtract the center of mass */
 987                             xp[i][j] -= com[j];
 988                         }
 989                     }
 990                 }
 991             }
 992         }
 993
 994         if (rc_scaling == erscCOM)
 995         {
 996             /* Convert the COM from Cartesian to crystal coordinates */
 997             for (j = 0; j < npbcdim; j++)
 998             {
 999                 com[j] *= invbox[j][j];
1000                 for (k = j+1; k < npbcdim; k++)
1001                 {
1002                     com[j] += invbox[k][j]*com[k];
1003                 }
1004             }
1005         }
1006     }
1007
1008     sfree(x);
1009     sfree(v);
1010     sfree(hadAtom);
1011 }
1012
1013 static void gen_posres(gmx_mtop_t *mtop, t_molinfo *mi,
1014                        const char *fnA, const char *fnB,
1015                        int rc_scaling, int ePBC,
1016                        rvec com, rvec comB,
1017                        warninp_t wi)
1018 {
1019     read_posres  (mtop, mi, FALSE, fnA, rc_scaling, ePBC, com, wi);
1020     /* It is safer to simply read the b-state posres rather than trying
1021      * to be smart and copy the positions.
1022      */
1023     read_posres(mtop, mi, TRUE, fnB, rc_scaling, ePBC, comB, wi);
1024 }
1025
1026 static void set_wall_atomtype(gpp_atomtype_t at, t_gromppopts *opts,
1027                               t_inputrec *ir, warninp_t wi)
1028 {
1029     int  i;
1030     char warn_buf[STRLEN];
1031
1032     if (ir->nwall > 0)
1033     {
1034         fprintf(stderr, "Searching the wall atom type(s)\n");
1035     }
1036     for (i = 0; i < ir->nwall; i++)
1037     {
1038         ir->wall_atomtype[i] = get_atomtype_type(opts->wall_atomtype[i], at);
1039         if (ir->wall_atomtype[i] == NOTSET)
1040         {
1041             sprintf(warn_buf, "Specified wall atom type %s is not defined", opts->wall_atomtype[i]);
1042             warning_error(wi, warn_buf);
1043         }
1044     }
1045 }
1046
1047 static int nrdf_internal(t_atoms *atoms)
1048 {
1049     int i, nmass, nrdf;
1050
1051     nmass = 0;
1052     for (i = 0; i < atoms->nr; i++)
1053     {
1054         /* Vsite ptype might not be set here yet, so also check the mass */
1055         if ((atoms->atom[i].ptype == eptAtom ||
1056              atoms->atom[i].ptype == eptNucleus)
1057             && atoms->atom[i].m > 0)
1058         {
1059             nmass++;
1060         }
1061     }
1062     switch (nmass)
1063     {
1064         case 0:  nrdf = 0; break;
1065         case 1:  nrdf = 0; break;
1066         case 2:  nrdf = 1; break;
1067         default: nrdf = nmass*3 - 6; break;
1068     }
1069
1070     return nrdf;
1071 }
1072
1073 static void
1074 spline1d( double        dx,
1075           double *      y,
1076           int           n,
1077           double *      u,
1078           double *      y2 )
1079 {
1080     int    i;
1081     double p, q;
1082
1083     y2[0] = 0.0;
1084     u[0]  = 0.0;
1085
1086     for (i = 1; i < n-1; i++)
1087     {
1088         p     = 0.5*y2[i-1]+2.0;
1089         y2[i] = -0.5/p;
1090         q     = (y[i+1]-2.0*y[i]+y[i-1])/dx;
1091         u[i]  = (3.0*q/dx-0.5*u[i-1])/p;
1092     }
1093
1094     y2[n-1] = 0.0;
1095
1096     for (i = n-2; i >= 0; i--)
1097     {
1098         y2[i] = y2[i]*y2[i+1]+u[i];
1099     }
1100 }
1101
1102
1103 static void
1104 interpolate1d( double     xmin,
1105                double     dx,
1106                double *   ya,
1107                double *   y2a,
1108                double     x,
1109                double *   y,
1110                double *   y1)
1111 {
1112     int    ix;
1113     double a, b;
1114
1115     ix = static_cast<int>((x-xmin)/dx);
1116
1117     a = (xmin+(ix+1)*dx-x)/dx;
1118     b = (x-xmin-ix*dx)/dx;
1119
1120     *y  = a*ya[ix]+b*ya[ix+1]+((a*a*a-a)*y2a[ix]+(b*b*b-b)*y2a[ix+1])*(dx*dx)/6.0;
1121     *y1 = (ya[ix+1]-ya[ix])/dx-(3.0*a*a-1.0)/6.0*dx*y2a[ix]+(3.0*b*b-1.0)/6.0*dx*y2a[ix+1];
1122 }
1123
1124
1125 static void
1126 setup_cmap (int              grid_spacing,
1127             int              nc,
1128             real *           grid,
1129             gmx_cmap_t *     cmap_grid)
1130 {
1131     double *tmp_u, *tmp_u2, *tmp_yy, *tmp_y1, *tmp_t2, *tmp_grid;
1132
1133     int     i, j, k, ii, jj, kk, idx;
1134     int     offset;
1135     double  dx, xmin, v, v1, v2, v12;
1136     double  phi, psi;
1137
1138     snew(tmp_u, 2*grid_spacing);
1139     snew(tmp_u2, 2*grid_spacing);
1140     snew(tmp_yy, 2*grid_spacing);
1141     snew(tmp_y1, 2*grid_spacing);
1142     snew(tmp_t2, 2*grid_spacing*2*grid_spacing);
1143     snew(tmp_grid, 2*grid_spacing*2*grid_spacing);
1144
1145     dx   = 360.0/grid_spacing;
1146     xmin = -180.0-dx*grid_spacing/2;
1147
1148     for (kk = 0; kk < nc; kk++)
1149     {
1150         /* Compute an offset depending on which cmap we are using
1151          * Offset will be the map number multiplied with the
1152          * grid_spacing * grid_spacing * 2
1153          */
1154         offset = kk * grid_spacing * grid_spacing * 2;
1155
1156         for (i = 0; i < 2*grid_spacing; i++)
1157         {
1158             ii = (i+grid_spacing-grid_spacing/2)%grid_spacing;
1159
1160             for (j = 0; j < 2*grid_spacing; j++)
1161             {
1162                 jj = (j+grid_spacing-grid_spacing/2)%grid_spacing;
1163                 tmp_grid[i*grid_spacing*2+j] = grid[offset+ii*grid_spacing+jj];
1164             }
1165         }
1166
1167         for (i = 0; i < 2*grid_spacing; i++)
1168         {
1169             spline1d(dx, &(tmp_grid[2*grid_spacing*i]), 2*grid_spacing, tmp_u, &(tmp_t2[2*grid_spacing*i]));
1170         }
1171
1172         for (i = grid_spacing/2; i < grid_spacing+grid_spacing/2; i++)
1173         {
1174             ii  = i-grid_spacing/2;
1175             phi = ii*dx-180.0;
1176
1177             for (j = grid_spacing/2; j < grid_spacing+grid_spacing/2; j++)
1178             {
1179                 jj  = j-grid_spacing/2;
1180                 psi = jj*dx-180.0;
1181
1182                 for (k = 0; k < 2*grid_spacing; k++)
1183                 {
1184                     interpolate1d(xmin, dx, &(tmp_grid[2*grid_spacing*k]),
1185                                   &(tmp_t2[2*grid_spacing*k]), psi, &tmp_yy[k], &tmp_y1[k]);
1186                 }
1187
1188                 spline1d(dx, tmp_yy, 2*grid_spacing, tmp_u, tmp_u2);
1189                 interpolate1d(xmin, dx, tmp_yy, tmp_u2, phi, &v, &v1);
1190                 spline1d(dx, tmp_y1, 2*grid_spacing, tmp_u, tmp_u2);
1191                 interpolate1d(xmin, dx, tmp_y1, tmp_u2, phi, &v2, &v12);
1192
1193                 idx = ii*grid_spacing+jj;
1194                 cmap_grid->cmapdata[kk].cmap[idx*4]   = grid[offset+ii*grid_spacing+jj];
1195                 cmap_grid->cmapdata[kk].cmap[idx*4+1] = v1;
1196                 cmap_grid->cmapdata[kk].cmap[idx*4+2] = v2;
1197                 cmap_grid->cmapdata[kk].cmap[idx*4+3] = v12;
1198             }
1199         }
1200     }
1201 }
1202
1203 static void init_cmap_grid(gmx_cmap_t *cmap_grid, int ngrid, int grid_spacing)
1204 {
1205     int i, nelem;
1206
1207     cmap_grid->ngrid        = ngrid;
1208     cmap_grid->grid_spacing = grid_spacing;
1209     nelem                   = cmap_grid->grid_spacing*cmap_grid->grid_spacing;
1210
1211     snew(cmap_grid->cmapdata, ngrid);
1212
1213     for (i = 0; i < cmap_grid->ngrid; i++)
1214     {
1215         snew(cmap_grid->cmapdata[i].cmap, 4*nelem);
1216     }
1217 }
1218
1219
1220 static int count_constraints(gmx_mtop_t *mtop, t_molinfo *mi, warninp_t wi)
1221 {
1222     int             count, count_mol, i, mb;
1223     gmx_molblock_t *molb;
1224     t_params       *plist;
1225     char            buf[STRLEN];
1226
1227     count = 0;
1228     for (mb = 0; mb < mtop->nmolblock; mb++)
1229     {
1230         count_mol = 0;
1231         molb      = &mtop->molblock[mb];
1232         plist     = mi[molb->type].plist;
1233
1234         for (i = 0; i < F_NRE; i++)
1235         {
1236             if (i == F_SETTLE)
1237             {
1238                 count_mol += 3*plist[i].nr;
1239             }
1240             else if (interaction_function[i].flags & IF_CONSTRAINT)
1241             {
1242                 count_mol += plist[i].nr;
1243             }
1244         }
1245
1246         if (count_mol > nrdf_internal(&mi[molb->type].atoms))
1247         {
1248             sprintf(buf,
1249                     "Molecule type '%s' has %d constraints.\n"
1250                     "For stability and efficiency there should not be more constraints than internal number of degrees of freedom: %d.\n",
1251                     *mi[molb->type].name, count_mol,
1252                     nrdf_internal(&mi[molb->type].atoms));
1253             warning(wi, buf);
1254         }
1255         count += molb->nmol*count_mol;
1256     }
1257
1258     return count;
1259 }
1260
1261 static void check_gbsa_params_charged(gmx_mtop_t *sys, gpp_atomtype_t atype)
1262 {
1263     int            i, nmiss, natoms, mt;
1264     real           q;
1265     const t_atoms *atoms;
1266
1267     nmiss = 0;
1268     for (mt = 0; mt < sys->nmoltype; mt++)
1269     {
1270         atoms  = &sys->moltype[mt].atoms;
1271         natoms = atoms->nr;
1272
1273         for (i = 0; i < natoms; i++)
1274         {
1275             q = atoms->atom[i].q;
1276             if ((get_atomtype_radius(atoms->atom[i].type, atype)    == 0  ||
1277                  get_atomtype_vol(atoms->atom[i].type, atype)       == 0  ||
1278                  get_atomtype_surftens(atoms->atom[i].type, atype)  == 0  ||
1279                  get_atomtype_gb_radius(atoms->atom[i].type, atype) == 0  ||
1280                  get_atomtype_S_hct(atoms->atom[i].type, atype)     == 0) &&
1281                 q != 0)
1282             {
1283                 fprintf(stderr, "\nGB parameter(s) zero for atom type '%s' while charge is %g\n",
1284                         get_atomtype_name(atoms->atom[i].type, atype), q);
1285                 nmiss++;
1286             }
1287         }
1288     }
1289
1290     if (nmiss > 0)
1291     {
1292         gmx_fatal(FARGS, "Can't do GB electrostatics; the implicit_genborn_params section of the forcefield has parameters with value zero for %d atomtypes that occur as charged atoms.", nmiss);
1293     }
1294 }
1295
1296
1297 static void check_gbsa_params(gpp_atomtype_t atype)
1298 {
1299     int  nmiss, i;
1300
1301     /* If we are doing GBSA, check that we got the parameters we need
1302      * This checking is to see if there are GBSA paratmeters for all
1303      * atoms in the force field. To go around this for testing purposes
1304      * comment out the nerror++ counter temporarily
1305      */
1306     nmiss = 0;
1307     for (i = 0; i < get_atomtype_ntypes(atype); i++)
1308     {
1309         if (get_atomtype_radius(i, atype)    < 0 ||
1310             get_atomtype_vol(i, atype)       < 0 ||
1311             get_atomtype_surftens(i, atype)  < 0 ||
1312             get_atomtype_gb_radius(i, atype) < 0 ||
1313             get_atomtype_S_hct(i, atype)     < 0)
1314         {
1315             fprintf(stderr, "\nGB parameter(s) missing or negative for atom type '%s'\n",
1316                     get_atomtype_name(i, atype));
1317             nmiss++;
1318         }
1319     }
1320
1321     if (nmiss > 0)
1322     {
1323         gmx_fatal(FARGS, "Can't do GB electrostatics; the implicit_genborn_params section of the forcefield is missing parameters for %d atomtypes or they might be negative.", nmiss);
1324     }
1325
1326 }
1327
1328 static real calc_temp(const gmx_mtop_t *mtop,
1329                       const t_inputrec *ir,
1330                       rvec             *v)
1331 {
1332     gmx_mtop_atomloop_all_t aloop;
1333     const t_atom           *atom;
1334     int                     a;
1335
1336     double                  sum_mv2 = 0;
1337     aloop = gmx_mtop_atomloop_all_init(mtop);
1338     while (gmx_mtop_atomloop_all_next(aloop, &a, &atom))
1339     {
1340         sum_mv2 += atom->m*norm2(v[a]);
1341     }
1342
1343     double nrdf = 0;
1344     for (int g = 0; g < ir->opts.ngtc; g++)
1345     {
1346         nrdf += ir->opts.nrdf[g];
1347     }
1348
1349     return sum_mv2/(nrdf*BOLTZ);
1350 }
1351
1352 static real get_max_reference_temp(const t_inputrec *ir,
1353                                    warninp_t         wi)
1354 {
1355     real     ref_t;
1356     int      i;
1357     gmx_bool bNoCoupl;
1358
1359     ref_t    = 0;
1360     bNoCoupl = FALSE;
1361     for (i = 0; i < ir->opts.ngtc; i++)
1362     {
1363         if (ir->opts.tau_t[i] < 0)
1364         {
1365             bNoCoupl = TRUE;
1366         }
1367         else
1368         {
1369             ref_t = std::max(ref_t, ir->opts.ref_t[i]);
1370         }
1371     }
1372
1373     if (bNoCoupl)
1374     {
1375         char buf[STRLEN];
1376
1377         sprintf(buf, "Some temperature coupling groups do not use temperature coupling. We will assume their temperature is not more than %.3f K. If their temperature is higher, the energy error and the Verlet buffer might be underestimated.",
1378                 ref_t);
1379         warning(wi, buf);
1380     }
1381
1382     return ref_t;
1383 }
1384
1385 /* Checks if there are unbound atoms in moleculetype molt.
1386  * Prints a note for each unbound atoms and a warning if any is present.
1387  */
1388 static void checkForUnboundAtoms(const gmx_moltype_t *molt,
1389                                  gmx_bool             bVerbose,
1390                                  warninp_t            wi)
1391 {
1392     const t_atoms *atoms = &molt->atoms;
1393
1394     if (atoms->nr == 1)
1395     {
1396         /* Only one atom, there can't be unbound atoms */
1397         return;
1398     }
1399
1400     std::vector<int> count(atoms->nr, 0);
1401
1402     for (int ftype = 0; ftype < F_NRE; ftype++)
1403     {
1404         if (((interaction_function[ftype].flags & IF_BOND) && ftype != F_CONNBONDS) ||
1405             (interaction_function[ftype].flags & IF_CONSTRAINT) ||
1406             ftype == F_SETTLE)
1407         {
1408             const t_ilist *il   = &molt->ilist[ftype];
1409             int            nral = NRAL(ftype);
1410
1411             for (int i = 0; i < il->nr; i += 1 + nral)
1412             {
1413                 for (int j = 0; j < nral; j++)
1414                 {
1415                     count[il->iatoms[i + 1 + j]]++;
1416                 }
1417             }
1418         }
1419     }
1420
1421     int numDanglingAtoms = 0;
1422     for (int a = 0; a < atoms->nr; a++)
1423     {
1424         if (atoms->atom[a].ptype != eptVSite &&
1425             count[a] == 0)
1426         {
1427             if (bVerbose)
1428             {
1429                 fprintf(stderr, "\nAtom %d '%s' in moleculetype '%s' is not bound by a potential or constraint to any other atom in the same moleculetype.\n",
1430                         a + 1, *atoms->atomname[a], *molt->name);
1431             }
1432             numDanglingAtoms++;
1433         }
1434     }
1435
1436     if (numDanglingAtoms > 0)
1437     {
1438         char buf[STRLEN];
1439         sprintf(buf, "In moleculetype '%s' %d atoms are not bound by a potential or constraint to any other atom in the same moleculetype. Although technically this might not cause issues in a simulation, this often means that the user forgot to add a bond/potential/constraint or put multiple molecules in the same moleculetype definition by mistake. Run with -v to get information for each atom.",
1440                 *molt->name, numDanglingAtoms);
1441         warning_note(wi, buf);
1442     }
1443 }
1444
1445 /* Checks all moleculetypes for unbound atoms */
1446 static void checkForUnboundAtoms(const gmx_mtop_t *mtop,
1447                                  gmx_bool          bVerbose,
1448                                  warninp_t         wi)
1449 {
1450     for (int mt = 0; mt < mtop->nmoltype; mt++)
1451     {
1452         checkForUnboundAtoms(&mtop->moltype[mt], bVerbose, wi);
1453     }
1454 }
1455
1456 /*! \brief Checks if there are decoupled modes in moleculetype \p molt.
1457  *
1458  * The specific decoupled modes this routine check for are angle modes
1459  * where the two bonds are constrained and the atoms a both ends are only
1460  * involved in a single constraint; the mass of the two atoms needs to
1461  * differ by more than \p massFactorThreshold.
1462  */
1463 static bool haveDecoupledModeInMol(const gmx_moltype_t *molt,
1464                                    const t_iparams     *iparams,
1465                                    real                 massFactorThreshold)
1466 {
1467     if (molt->ilist[F_CONSTR].nr == 0 && molt->ilist[F_CONSTRNC].nr == 0)
1468     {
1469         return false;
1470     }
1471
1472     const t_atom * atom = molt->atoms.atom;
1473
1474     int            numFlexibleConstraints;
1475     t_blocka       atomToConstraints = make_at2con(0, molt->atoms.nr,
1476                                                    molt->ilist, iparams,
1477                                                    FALSE,
1478                                                    &numFlexibleConstraints);
1479
1480     bool           haveDecoupledMode = false;
1481     for (int ftype = 0; ftype < F_NRE; ftype++)
1482     {
1483         if (interaction_function[ftype].flags & IF_ATYPE)
1484         {
1485             const int      nral = NRAL(ftype);
1486             const t_ilist *il   = &molt->ilist[ftype];
1487             for (int i = 0; i < il->nr; i += 1 + nral)
1488             {
1489                 /* Here we check for the mass difference between the atoms
1490                  * at both ends of the angle, that the atoms at the ends have
1491                  * 1 contraint and the atom in the middle at least 3; we check
1492                  * that the 3 atoms are linked by constraints below.
1493                  * We check for at least three constraints for the middle atom,
1494                  * since with only the two bonds in the angle, we have 3-atom
1495                  * molecule, which has much more thermal exhange in this single
1496                  * angle mode than molecules with more atoms.
1497                  * Note that this check also catches molecules where more atoms
1498                  * are connected to one or more atoms in the angle, but by
1499                  * bond potentials instead of angles. But such cases will not
1500                  * occur in "normal" molecules and it doesn't hurt running
1501                  * those with higher accuracy settings as well.
1502                  */
1503                 int a0 = il->iatoms[1 + i];
1504                 int a1 = il->iatoms[1 + i + 1];
1505                 int a2 = il->iatoms[1 + i + 2];
1506                 if ((atom[a0].m > atom[a2].m*massFactorThreshold ||
1507                      atom[a2].m > atom[a0].m*massFactorThreshold) &&
1508                     atomToConstraints.index[a0 + 1] - atomToConstraints.index[a0] == 1 &&
1509                     atomToConstraints.index[a2 + 1] - atomToConstraints.index[a2] == 1 &&
1510                     atomToConstraints.index[a1 + 1] - atomToConstraints.index[a1] >= 3)
1511                 {
1512                     int  constraint0 = atomToConstraints.a[atomToConstraints.index[a0]];
1513                     int  constraint2 = atomToConstraints.a[atomToConstraints.index[a2]];
1514
1515                     bool foundAtom0  = false;
1516                     bool foundAtom2  = false;
1517                     for (int conIndex = atomToConstraints.index[a1]; conIndex < atomToConstraints.index[a1 + 1]; conIndex++)
1518                     {
1519                         if (atomToConstraints.a[conIndex] == constraint0)
1520                         {
1521                             foundAtom0 = true;
1522                         }
1523                         if (atomToConstraints.a[conIndex] == constraint2)
1524                         {
1525                             foundAtom2 = true;
1526                         }
1527                     }
1528                     if (foundAtom0 && foundAtom2)
1529                     {
1530                         haveDecoupledMode = true;
1531                     }
1532                 }
1533             }
1534         }
1535     }
1536
1537     done_blocka(&atomToConstraints);
1538
1539     return haveDecoupledMode;
1540 }
1541
1542 /*! \brief Checks if the Verlet buffer and constraint accuracy is sufficient for decoupled dynamic modes.
1543  *
1544  * When decoupled modes are present and the accuray in insufficient,
1545  * this routine issues a warning if the accuracy is insufficient.
1546  */
1547 static void checkDecoupledModeAccuracy(const gmx_mtop_t *mtop,
1548                                        const t_inputrec *ir,
1549                                        warninp_t         wi)
1550 {
1551     /* We only have issues with decoupled modes with normal MD.
1552      * With stochastic dynamics equipartitioning is enforced strongly.
1553      */
1554     if (!EI_MD(ir->eI))
1555     {
1556         return;
1557     }
1558
1559     /* When atoms of very different mass are involved in an angle potential
1560      * and both bonds in the angle are constrained, the dynamic modes in such
1561      * angles have very different periods and significant energy exchange
1562      * takes several nanoseconds. Thus even a small amount of error in
1563      * different algorithms can lead to violation of equipartitioning.
1564      * The parameters below are mainly based on an all-atom chloroform model
1565      * with all bonds constrained. Equipartitioning is off by more than 1%
1566      * (need to run 5-10 ns) when the difference in mass between H and Cl
1567      * is more than a factor 13 and the accuracy is less than the thresholds
1568      * given below. This has been verified on some other molecules.
1569      *
1570      * Note that the buffer and shake parameters have unit length and
1571      * energy/time, respectively, so they will "only" work correctly
1572      * for atomistic force fields using MD units.
1573      */
1574     const real massFactorThreshold      = 13.0;
1575     const real bufferToleranceThreshold = 1e-4;
1576     const int  lincsIterationThreshold  = 2;
1577     const int  lincsOrderThreshold      = 4;
1578     const real shakeToleranceThreshold  = 0.005*ir->delta_t;
1579
1580     bool       lincsWithSufficientTolerance = (ir->eConstrAlg == econtLINCS && ir->nLincsIter >= lincsIterationThreshold && ir->nProjOrder >= lincsOrderThreshold);
1581     bool       shakeWithSufficientTolerance = (ir->eConstrAlg == econtSHAKE && ir->shake_tol <= 1.1*shakeToleranceThreshold);
1582     if (ir->cutoff_scheme == ecutsVERLET &&
1583         ir->verletbuf_tol <= 1.1*bufferToleranceThreshold &&
1584         (lincsWithSufficientTolerance || shakeWithSufficientTolerance))
1585     {
1586         return;
1587     }
1588
1589     bool haveDecoupledMode = false;
1590     for (int mt = 0; mt < mtop->nmoltype; mt++)
1591     {
1592         if (haveDecoupledModeInMol(&mtop->moltype[mt], mtop->ffparams.iparams,
1593                                    massFactorThreshold))
1594         {
1595             haveDecoupledMode = true;
1596         }
1597     }
1598
1599     if (haveDecoupledMode)
1600     {
1601         char modeMessage[STRLEN];
1602         sprintf(modeMessage, "There are atoms at both ends of an angle, connected by constraints and with masses that differ by more than a factor of %g. This means that there are likely dynamic modes that are only very weakly coupled.",
1603                 massFactorThreshold);
1604         char buf[STRLEN];
1605         if (ir->cutoff_scheme == ecutsVERLET)
1606         {
1607             sprintf(buf, "%s To ensure good equipartitioning, you need to either not use constraints on all bonds (but, if possible, only on bonds involving hydrogens) or use integrator = %s or decrease one or more tolerances: verlet-buffer-tolerance <= %g, LINCS iterations >= %d, LINCS order >= %d or SHAKE tolerance <= %g",
1608                     modeMessage,
1609                     ei_names[eiSD1],
1610                     bufferToleranceThreshold,
1611                     lincsIterationThreshold, lincsOrderThreshold,
1612                     shakeToleranceThreshold);
1613         }
1614         else
1615         {
1616             sprintf(buf, "%s To ensure good equipartitioning, we suggest to switch to the %s cutoff-scheme, since that allows for better control over the Verlet buffer size and thus over the energy drift.",
1617                     modeMessage,
1618                     ecutscheme_names[ecutsVERLET]);
1619         }
1620         warning(wi, buf);
1621     }
1622 }
1623
1624 static void set_verlet_buffer(const gmx_mtop_t *mtop,
1625                               t_inputrec       *ir,
1626                               real              buffer_temp,
1627                               matrix            box,
1628                               warninp_t         wi)
1629 {
1630     verletbuf_list_setup_t ls;
1631     real                   rlist_1x1;
1632     int                    n_nonlin_vsite;
1633     char                   warn_buf[STRLEN];
1634
1635     printf("Determining Verlet buffer for a tolerance of %g kJ/mol/ps at %g K\n", ir->verletbuf_tol, buffer_temp);
1636
1637     /* Calculate the buffer size for simple atom vs atoms list */
1638     ls.cluster_size_i = 1;
1639     ls.cluster_size_j = 1;
1640     calc_verlet_buffer_size(mtop, det(box), ir, ir->nstlist, ir->nstlist - 1,
1641                             buffer_temp, &ls, &n_nonlin_vsite, &rlist_1x1);
1642
1643     /* Set the pair-list buffer size in ir */
1644     verletbuf_get_list_setup(FALSE, FALSE, &ls);
1645     calc_verlet_buffer_size(mtop, det(box), ir, ir->nstlist, ir->nstlist - 1,
1646                             buffer_temp, &ls, &n_nonlin_vsite, &ir->rlist);
1647
1648     if (n_nonlin_vsite > 0)
1649     {
1650         sprintf(warn_buf, "There are %d non-linear virtual site constructions. Their contribution to the energy error is approximated. In most cases this does not affect the error significantly.", n_nonlin_vsite);
1651         warning_note(wi, warn_buf);
1652     }
1653
1654     printf("Calculated rlist for %dx%d atom pair-list as %.3f nm, buffer size %.3f nm\n",
1655            1, 1, rlist_1x1, rlist_1x1-std::max(ir->rvdw, ir->rcoulomb));
1656
1657     printf("Set rlist, assuming %dx%d atom pair-list, to %.3f nm, buffer size %.3f nm\n",
1658            ls.cluster_size_i, ls.cluster_size_j,
1659            ir->rlist, ir->rlist-std::max(ir->rvdw, ir->rcoulomb));
1660
1661     printf("Note that mdrun will redetermine rlist based on the actual pair-list setup\n");
1662
1663     if (gmx::square(ir->rlist) >= max_cutoff2(ir->ePBC, box))
1664     {
1665         gmx_fatal(FARGS, "The pair-list cut-off (%g nm) is longer than half the shortest box vector or longer than the smallest box diagonal element (%g nm). Increase the box size or decrease nstlist or increase verlet-buffer-tolerance.", ir->rlist, std::sqrt(max_cutoff2(ir->ePBC, box)));
1666     }
1667 }
1668
1669 int gmx_grompp(int argc, char *argv[])
1670 {
1671     const char        *desc[] = {
1672         "[THISMODULE] (the gromacs preprocessor)",
1673         "reads a molecular topology file, checks the validity of the",
1674         "file, expands the topology from a molecular description to an atomic",
1675         "description. The topology file contains information about",
1676         "molecule types and the number of molecules, the preprocessor",
1677         "copies each molecule as needed. ",
1678         "There is no limitation on the number of molecule types. ",
1679         "Bonds and bond-angles can be converted into constraints, separately",
1680         "for hydrogens and heavy atoms.",
1681         "Then a coordinate file is read and velocities can be generated",
1682         "from a Maxwellian distribution if requested.",
1683         "[THISMODULE] also reads parameters for [gmx-mdrun] ",
1684         "(eg. number of MD steps, time step, cut-off), and others such as",
1685         "NEMD parameters, which are corrected so that the net acceleration",
1686         "is zero.",
1687         "Eventually a binary file is produced that can serve as the sole input",
1688         "file for the MD program.[PAR]",
1689
1690         "[THISMODULE] uses the atom names from the topology file. The atom names",
1691         "in the coordinate file (option [TT]-c[tt]) are only read to generate",
1692         "warnings when they do not match the atom names in the topology.",
1693         "Note that the atom names are irrelevant for the simulation as",
1694         "only the atom types are used for generating interaction parameters.[PAR]",
1695
1696         "[THISMODULE] uses a built-in preprocessor to resolve includes, macros, ",
1697         "etc. The preprocessor supports the following keywords::",
1698         "",
1699         "    #ifdef VARIABLE",
1700         "    #ifndef VARIABLE",
1701         "    #else",
1702         "    #endif",
1703         "    #define VARIABLE",
1704         "    #undef VARIABLE",
1705         "    #include \"filename\"",
1706         "    #include <filename>",
1707         "",
1708         "The functioning of these statements in your topology may be modulated by",
1709         "using the following two flags in your [REF].mdp[ref] file::",
1710         "",
1711         "    define = -DVARIABLE1 -DVARIABLE2",
1712         "    include = -I/home/john/doe",
1713         "",
1714         "For further information a C-programming textbook may help you out.",
1715         "Specifying the [TT]-pp[tt] flag will get the pre-processed",
1716         "topology file written out so that you can verify its contents.[PAR]",
1717
1718         "When using position restraints, a file with restraint coordinates",
1719         "must be supplied with [TT]-r[tt] (can be the same file as supplied",
1720         "for [TT]-c[tt]). For free energy calculations, separate reference",
1721         "coordinates for the B topology can be supplied with [TT]-rb[tt],",
1722         "otherwise they will be equal to those of the A topology.[PAR]",
1723
1724         "Starting coordinates can be read from trajectory with [TT]-t[tt].",
1725         "The last frame with coordinates and velocities will be read,",
1726         "unless the [TT]-time[tt] option is used. Only if this information",
1727         "is absent will the coordinates in the [TT]-c[tt] file be used.",
1728         "Note that these velocities will not be used when [TT]gen_vel = yes[tt]",
1729         "in your [REF].mdp[ref] file. An energy file can be supplied with",
1730         "[TT]-e[tt] to read Nose-Hoover and/or Parrinello-Rahman coupling",
1731         "variables.[PAR]",
1732
1733         "[THISMODULE] can be used to restart simulations (preserving",
1734         "continuity) by supplying just a checkpoint file with [TT]-t[tt].",
1735         "However, for simply changing the number of run steps to extend",
1736         "a run, using [gmx-convert-tpr] is more convenient than [THISMODULE].",
1737         "You then supply the old checkpoint file directly to [gmx-mdrun]",
1738         "with [TT]-cpi[tt]. If you wish to change the ensemble or things",
1739         "like output frequency, then supplying the checkpoint file to",
1740         "[THISMODULE] with [TT]-t[tt] along with a new [REF].mdp[ref] file",
1741         "with [TT]-f[tt] is the recommended procedure. Actually preserving",
1742         "the ensemble (if possible) still requires passing the checkpoint",
1743         "file to [gmx-mdrun] [TT]-cpi[tt].[PAR]",
1744
1745         "By default, all bonded interactions which have constant energy due to",
1746         "virtual site constructions will be removed. If this constant energy is",
1747         "not zero, this will result in a shift in the total energy. All bonded",
1748         "interactions can be kept by turning off [TT]-rmvsbds[tt]. Additionally,",
1749         "all constraints for distances which will be constant anyway because",
1750         "of virtual site constructions will be removed. If any constraints remain",
1751         "which involve virtual sites, a fatal error will result.[PAR]"
1752
1753         "To verify your run input file, please take note of all warnings",
1754         "on the screen, and correct where necessary. Do also look at the contents",
1755         "of the [TT]mdout.mdp[tt] file; this contains comment lines, as well as",
1756         "the input that [THISMODULE] has read. If in doubt, you can start [THISMODULE]",
1757         "with the [TT]-debug[tt] option which will give you more information",
1758         "in a file called [TT]grompp.log[tt] (along with real debug info). You",
1759         "can see the contents of the run input file with the [gmx-dump]",
1760         "program. [gmx-check] can be used to compare the contents of two",
1761         "run input files.[PAR]"
1762
1763         "The [TT]-maxwarn[tt] option can be used to override warnings printed",
1764         "by [THISMODULE] that otherwise halt output. In some cases, warnings are",
1765         "harmless, but usually they are not. The user is advised to carefully",
1766         "interpret the output messages before attempting to bypass them with",
1767         "this option."
1768     };
1769     t_gromppopts      *opts;
1770     gmx_mtop_t        *sys;
1771     int                nmi;
1772     t_molinfo         *mi, *intermolecular_interactions;
1773     gpp_atomtype_t     atype;
1774     int                nvsite, comb, mt;
1775     t_params          *plist;
1776     matrix             box;
1777     real               fudgeQQ;
1778     double             reppow;
1779     const char        *mdparin;
1780     int                ntype;
1781     gmx_bool           bNeedVel, bGenVel;
1782     gmx_bool           have_atomnumber;
1783     gmx_output_env_t  *oenv;
1784     gmx_bool           bVerbose = FALSE;
1785     warninp_t          wi;
1786     char               warn_buf[STRLEN];
1787
1788     t_filenm           fnm[] = {
1789         { efMDP, nullptr,  nullptr,        ffREAD  },
1790         { efMDP, "-po", "mdout",     ffWRITE },
1791         { efSTX, "-c",  nullptr,        ffREAD  },
1792         { efSTX, "-r",  "restraint", ffOPTRD },
1793         { efSTX, "-rb", "restraint", ffOPTRD },
1794         { efNDX, nullptr,  nullptr,        ffOPTRD },
1795         { efTOP, nullptr,  nullptr,        ffREAD  },
1796         { efTOP, "-pp", "processed", ffOPTWR },
1797         { efTPR, "-o",  nullptr,        ffWRITE },
1798         { efTRN, "-t",  nullptr,        ffOPTRD },
1799         { efEDR, "-e",  nullptr,        ffOPTRD },
1800         /* This group is needed by the VMD viewer as the start configuration for IMD sessions: */
1801         { efGRO, "-imd", "imdgroup", ffOPTWR },
1802         { efTRN, "-ref", "rotref",   ffOPTRW }
1803     };
1804 #define NFILE asize(fnm)
1805
1806     /* Command line options */
1807     gmx_bool        bRenum   = TRUE;
1808     gmx_bool        bRmVSBds = TRUE, bZero = FALSE;
1809     int             i, maxwarn = 0;
1810     real            fr_time = -1;
1811     t_pargs         pa[]    = {
1812         { "-v",       FALSE, etBOOL, {&bVerbose},
1813           "Be loud and noisy" },
1814         { "-time",    FALSE, etREAL, {&fr_time},
1815           "Take frame at or first after this time." },
1816         { "-rmvsbds", FALSE, etBOOL, {&bRmVSBds},
1817           "Remove constant bonded interactions with virtual sites" },
1818         { "-maxwarn", FALSE, etINT,  {&maxwarn},
1819           "Number of allowed warnings during input processing. Not for normal use and may generate unstable systems" },
1820         { "-zero",    FALSE, etBOOL, {&bZero},
1821           "Set parameters for bonded interactions without defaults to zero instead of generating an error" },
1822         { "-renum",   FALSE, etBOOL, {&bRenum},
1823           "Renumber atomtypes and minimize number of atomtypes" }
1824     };
1825
1826     /* Parse the command line */
1827     if (!parse_common_args(&argc, argv, 0, NFILE, fnm, asize(pa), pa,
1828                            asize(desc), desc, 0, nullptr, &oenv))
1829     {
1830         return 0;
1831     }
1832
1833     /* Initiate some variables */
1834     gmx::MDModules mdModules;
1835     t_inputrec     irInstance;
1836     t_inputrec    *ir = &irInstance;
1837     snew(opts, 1);
1838     snew(opts->include, STRLEN);
1839     snew(opts->define, STRLEN);
1840
1841     wi = init_warning(TRUE, maxwarn);
1842
1843     /* PARAMETER file processing */
1844     mdparin = opt2fn("-f", NFILE, fnm);
1845     set_warning_line(wi, mdparin, -1);
1846     try
1847     {
1848         get_ir(mdparin, opt2fn("-po", NFILE, fnm), &mdModules, ir, opts, WriteMdpHeader::yes, wi);
1849     }
1850     GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
1851
1852     if (bVerbose)
1853     {
1854         fprintf(stderr, "checking input for internal consistency...\n");
1855     }
1856     check_ir(mdparin, ir, opts, wi);
1857
1858     if (ir->ld_seed == -1)
1859     {
1860         ir->ld_seed = static_cast<int>(gmx::makeRandomSeed());
1861         fprintf(stderr, "Setting the LD random seed to %" GMX_PRId64 "\n", ir->ld_seed);
1862     }
1863
1864     if (ir->expandedvals->lmc_seed == -1)
1865     {
1866         ir->expandedvals->lmc_seed = static_cast<int>(gmx::makeRandomSeed());
1867         fprintf(stderr, "Setting the lambda MC random seed to %d\n", ir->expandedvals->lmc_seed);
1868     }
1869
1870     bNeedVel = EI_STATE_VELOCITY(ir->eI);
1871     bGenVel  = (bNeedVel && opts->bGenVel);
1872     if (bGenVel && ir->bContinuation)
1873     {
1874         sprintf(warn_buf,
1875                 "Generating velocities is inconsistent with attempting "
1876                 "to continue a previous run. Choose only one of "
1877                 "gen-vel = yes and continuation = yes.");
1878         warning_error(wi, warn_buf);
1879     }
1880
1881     snew(plist, F_NRE);
1882     init_plist(plist);
1883     snew(sys, 1);
1884     atype = init_atomtype();
1885     if (debug)
1886     {
1887         pr_symtab(debug, 0, "Just opened", &sys->symtab);
1888     }
1889
1890     const char *fn = ftp2fn(efTOP, NFILE, fnm);
1891     if (!gmx_fexist(fn))
1892     {
1893         gmx_fatal(FARGS, "%s does not exist", fn);
1894     }
1895
1896     t_state state;
1897     new_status(fn, opt2fn_null("-pp", NFILE, fnm), opt2fn("-c", NFILE, fnm),
1898                opts, ir, bZero, bGenVel, bVerbose, &state,
1899                atype, sys, &nmi, &mi, &intermolecular_interactions,
1900                plist, &comb, &reppow, &fudgeQQ,
1901                opts->bMorse,
1902                wi);
1903
1904     if (debug)
1905     {
1906         pr_symtab(debug, 0, "After new_status", &sys->symtab);
1907     }
1908
1909     nvsite = 0;
1910     /* set parameters for virtual site construction (not for vsiten) */
1911     for (mt = 0; mt < sys->nmoltype; mt++)
1912     {
1913         nvsite +=
1914             set_vsites(bVerbose, &sys->moltype[mt].atoms, atype, mi[mt].plist);
1915     }
1916     /* now throw away all obsolete bonds, angles and dihedrals: */
1917     /* note: constraints are ALWAYS removed */
1918     if (nvsite)
1919     {
1920         for (mt = 0; mt < sys->nmoltype; mt++)
1921         {
1922             clean_vsite_bondeds(mi[mt].plist, sys->moltype[mt].atoms.nr, bRmVSBds);
1923         }
1924     }
1925
1926     if (ir->cutoff_scheme == ecutsVERLET)
1927     {
1928         fprintf(stderr, "Removing all charge groups because cutoff-scheme=%s\n",
1929                 ecutscheme_names[ir->cutoff_scheme]);
1930
1931         /* Remove all charge groups */
1932         gmx_mtop_remove_chargegroups(sys);
1933     }
1934
1935     if (count_constraints(sys, mi, wi) && (ir->eConstrAlg == econtSHAKE))
1936     {
1937         if (ir->eI == eiCG || ir->eI == eiLBFGS)
1938         {
1939             sprintf(warn_buf, "Can not do %s with %s, use %s",
1940                     EI(ir->eI), econstr_names[econtSHAKE], econstr_names[econtLINCS]);
1941             warning_error(wi, warn_buf);
1942         }
1943         if (ir->bPeriodicMols)
1944         {
1945             sprintf(warn_buf, "Can not do periodic molecules with %s, use %s",
1946                     econstr_names[econtSHAKE], econstr_names[econtLINCS]);
1947             warning_error(wi, warn_buf);
1948         }
1949     }
1950
1951     if (EI_SD (ir->eI) &&  ir->etc != etcNO)
1952     {
1953         warning_note(wi, "Temperature coupling is ignored with SD integrators.");
1954     }
1955
1956     /* If we are doing QM/MM, check that we got the atom numbers */
1957     have_atomnumber = TRUE;
1958     for (i = 0; i < get_atomtype_ntypes(atype); i++)
1959     {
1960         have_atomnumber = have_atomnumber && (get_atomtype_atomnumber(i, atype) >= 0);
1961     }
1962     if (!have_atomnumber && ir->bQMMM)
1963     {
1964         warning_error(wi,
1965                       "\n"
1966                       "It appears as if you are trying to run a QM/MM calculation, but the force\n"
1967                       "field you are using does not contain atom numbers fields. This is an\n"
1968                       "optional field (introduced in GROMACS 3.3) for general runs, but mandatory\n"
1969                       "for QM/MM. The good news is that it is easy to add - put the atom number as\n"
1970                       "an integer just before the mass column in ffXXXnb.itp.\n"
1971                       "NB: United atoms have the same atom numbers as normal ones.\n\n");
1972     }
1973
1974     /* Check for errors in the input now, since they might cause problems
1975      * during processing further down.
1976      */
1977     check_warning_error(wi, FARGS);
1978
1979     if (nint_ftype(sys, mi, F_POSRES) > 0 ||
1980         nint_ftype(sys, mi, F_FBPOSRES) > 0)
1981     {
1982         const char *fn = opt2fn("-r", NFILE, fnm);
1983         const char *fnB;
1984
1985         if (opt2bSet("-rb", NFILE, fnm))
1986         {
1987             fnB = opt2fn("-rb", NFILE, fnm);
1988         }
1989         else
1990         {
1991             fnB = fn;
1992         }
1993
1994         if (bVerbose)
1995         {
1996             fprintf(stderr, "Reading position restraint coords from %s", fn);
1997             if (strcmp(fn, fnB) == 0)
1998             {
1999                 fprintf(stderr, "\n");
2000             }
2001             else
2002             {
2003                 fprintf(stderr, " and %s\n", fnB);
2004             }
2005         }
2006         gen_posres(sys, mi, fn, fnB,
2007                    ir->refcoord_scaling, ir->ePBC,
2008                    ir->posres_com, ir->posres_comB,
2009                    wi);
2010     }
2011
2012     /* If we are using CMAP, setup the pre-interpolation grid */
2013     if (plist[F_CMAP].ncmap > 0)
2014     {
2015         init_cmap_grid(&sys->ffparams.cmap_grid, plist[F_CMAP].nc, plist[F_CMAP].grid_spacing);
2016         setup_cmap(plist[F_CMAP].grid_spacing, plist[F_CMAP].nc, plist[F_CMAP].cmap, &sys->ffparams.cmap_grid);
2017     }
2018
2019     set_wall_atomtype(atype, opts, ir, wi);
2020     if (bRenum)
2021     {
2022         renum_atype(plist, sys, ir->wall_atomtype, atype, bVerbose);
2023         get_atomtype_ntypes(atype);
2024     }
2025
2026     if (ir->implicit_solvent != eisNO)
2027     {
2028         /* Now we have renumbered the atom types, we can check the GBSA params */
2029         check_gbsa_params(atype);
2030
2031         /* Check that all atoms that have charge and/or LJ-parameters also have
2032          * sensible GB-parameters
2033          */
2034         check_gbsa_params_charged(sys, atype);
2035     }
2036
2037     /* PELA: Copy the atomtype data to the topology atomtype list */
2038     copy_atomtype_atomtypes(atype, &(sys->atomtypes));
2039
2040     if (debug)
2041     {
2042         pr_symtab(debug, 0, "After renum_atype", &sys->symtab);
2043     }
2044
2045     if (bVerbose)
2046     {
2047         fprintf(stderr, "converting bonded parameters...\n");
2048     }
2049
2050     ntype = get_atomtype_ntypes(atype);
2051     convert_params(ntype, plist, mi, intermolecular_interactions,
2052                    comb, reppow, fudgeQQ, sys);
2053
2054     if (debug)
2055     {
2056         pr_symtab(debug, 0, "After convert_params", &sys->symtab);
2057     }
2058
2059     /* set ptype to VSite for virtual sites */
2060     for (mt = 0; mt < sys->nmoltype; mt++)
2061     {
2062         set_vsites_ptype(FALSE, &sys->moltype[mt]);
2063     }
2064     if (debug)
2065     {
2066         pr_symtab(debug, 0, "After virtual sites", &sys->symtab);
2067     }
2068     /* Check velocity for virtual sites and shells */
2069     if (bGenVel)
2070     {
2071         check_vel(sys, as_rvec_array(state.v.data()));
2072     }
2073
2074     /* check for shells and inpurecs */
2075     check_shells_inputrec(sys, ir, wi);
2076
2077     /* check masses */
2078     check_mol(sys, wi);
2079
2080     checkForUnboundAtoms(sys, bVerbose, wi);
2081
2082     for (i = 0; i < sys->nmoltype; i++)
2083     {
2084         check_cg_sizes(ftp2fn(efTOP, NFILE, fnm), &sys->moltype[i].cgs, wi);
2085     }
2086
2087     if (EI_DYNAMICS(ir->eI) && ir->eI != eiBD)
2088     {
2089         check_bonds_timestep(sys, ir->delta_t, wi);
2090     }
2091
2092     checkDecoupledModeAccuracy(sys, ir, wi);
2093
2094     if (EI_ENERGY_MINIMIZATION(ir->eI) && 0 == ir->nsteps)
2095     {
2096         warning_note(wi, "Zero-step energy minimization will alter the coordinates before calculating the energy. If you just want the energy of a single point, try zero-step MD (with unconstrained_start = yes). To do multiple single-point energy evaluations of different configurations of the same topology, use mdrun -rerun.");
2097     }
2098
2099     check_warning_error(wi, FARGS);
2100
2101     if (bVerbose)
2102     {
2103         fprintf(stderr, "initialising group options...\n");
2104     }
2105     do_index(mdparin, ftp2fn_null(efNDX, NFILE, fnm),
2106              sys, bVerbose, ir,
2107              wi);
2108
2109     if (ir->cutoff_scheme == ecutsVERLET && ir->verletbuf_tol > 0)
2110     {
2111         if (EI_DYNAMICS(ir->eI) && inputrec2nboundeddim(ir) == 3)
2112         {
2113             real buffer_temp;
2114
2115             if (EI_MD(ir->eI) && ir->etc == etcNO)
2116             {
2117                 if (bGenVel)
2118                 {
2119                     buffer_temp = opts->tempi;
2120                 }
2121                 else
2122                 {
2123                     buffer_temp = calc_temp(sys, ir, as_rvec_array(state.v.data()));
2124                 }
2125                 if (buffer_temp > 0)
2126                 {
2127                     sprintf(warn_buf, "NVE simulation: will use the initial temperature of %.3f K for determining the Verlet buffer size", buffer_temp);
2128                     warning_note(wi, warn_buf);
2129                 }
2130                 else
2131                 {
2132                     sprintf(warn_buf, "NVE simulation with an initial temperature of zero: will use a Verlet buffer of %d%%. Check your energy drift!",
2133                             (int)(verlet_buffer_ratio_NVE_T0*100 + 0.5));
2134                     warning_note(wi, warn_buf);
2135                 }
2136             }
2137             else
2138             {
2139                 buffer_temp = get_max_reference_temp(ir, wi);
2140             }
2141
2142             if (EI_MD(ir->eI) && ir->etc == etcNO && buffer_temp == 0)
2143             {
2144                 /* NVE with initial T=0: we add a fixed ratio to rlist.
2145                  * Since we don't actually use verletbuf_tol, we set it to -1
2146                  * so it can't be misused later.
2147                  */
2148                 ir->rlist         *= 1.0 + verlet_buffer_ratio_NVE_T0;
2149                 ir->verletbuf_tol  = -1;
2150             }
2151             else
2152             {
2153                 /* We warn for NVE simulations with a drift tolerance that
2154                  * might result in a 1(.1)% drift over the total run-time.
2155                  * Note that we can't warn when nsteps=0, since we don't
2156                  * know how many steps the user intends to run.
2157                  */
2158                 if (EI_MD(ir->eI) && ir->etc == etcNO && ir->nstlist > 1 &&
2159                     ir->nsteps > 0)
2160                 {
2161                     const real driftTolerance = 0.01;
2162                     /* We use 2 DOF per atom = 2kT pot+kin energy,
2163                      * to be on the safe side with constraints.
2164                      */
2165                     const real totalEnergyDriftPerAtomPerPicosecond = 2*BOLTZ*buffer_temp/(ir->nsteps*ir->delta_t);
2166
2167                     if (ir->verletbuf_tol > 1.1*driftTolerance*totalEnergyDriftPerAtomPerPicosecond)
2168                     {
2169                         sprintf(warn_buf, "You are using a Verlet buffer tolerance of %g kJ/mol/ps for an NVE simulation of length %g ps, which can give a final drift of %d%%. For conserving energy to %d%% when using constraints, you might need to set verlet-buffer-tolerance to %.1e.",
2170                                 ir->verletbuf_tol, ir->nsteps*ir->delta_t,
2171                                 (int)(ir->verletbuf_tol/totalEnergyDriftPerAtomPerPicosecond*100 + 0.5),
2172                                 (int)(100*driftTolerance + 0.5),
2173                                 driftTolerance*totalEnergyDriftPerAtomPerPicosecond);
2174                         warning_note(wi, warn_buf);
2175                     }
2176                 }
2177
2178                 set_verlet_buffer(sys, ir, buffer_temp, state.box, wi);
2179             }
2180         }
2181     }
2182
2183     /* Init the temperature coupling state */
2184     init_gtc_state(&state, ir->opts.ngtc, 0, ir->opts.nhchainlength); /* need to add nnhpres here? */
2185
2186     if (bVerbose)
2187     {
2188         fprintf(stderr, "Checking consistency between energy and charge groups...\n");
2189     }
2190     check_eg_vs_cg(sys);
2191
2192     if (debug)
2193     {
2194         pr_symtab(debug, 0, "After index", &sys->symtab);
2195     }
2196
2197     triple_check(mdparin, ir, sys, wi);
2198     close_symtab(&sys->symtab);
2199     if (debug)
2200     {
2201         pr_symtab(debug, 0, "After close", &sys->symtab);
2202     }
2203
2204     /* make exclusions between QM atoms */
2205     if (ir->bQMMM)
2206     {
2207         if (ir->QMMMscheme == eQMMMschemenormal && ir->ns_type == ensSIMPLE)
2208         {
2209             gmx_fatal(FARGS, "electrostatic embedding only works with grid neighboursearching, use ns-type=grid instead\n");
2210         }
2211         else
2212         {
2213             generate_qmexcl(sys, ir, wi);
2214         }
2215     }
2216
2217     if (ftp2bSet(efTRN, NFILE, fnm))
2218     {
2219         if (bVerbose)
2220         {
2221             fprintf(stderr, "getting data from old trajectory ...\n");
2222         }
2223         cont_status(ftp2fn(efTRN, NFILE, fnm), ftp2fn_null(efEDR, NFILE, fnm),
2224                     bNeedVel, bGenVel, fr_time, ir, &state, sys, oenv);
2225     }
2226
2227     if (ir->ePBC == epbcXY && ir->nwall != 2)
2228     {
2229         clear_rvec(state.box[ZZ]);
2230     }
2231
2232     if (ir->cutoff_scheme != ecutsVERLET && ir->rlist > 0)
2233     {
2234         set_warning_line(wi, mdparin, -1);
2235         check_chargegroup_radii(sys, ir, as_rvec_array(state.x.data()), wi);
2236     }
2237
2238     if (EEL_FULL(ir->coulombtype) || EVDW_PME(ir->vdwtype))
2239     {
2240         /* Calculate the optimal grid dimensions */
2241         copy_mat(state.box, box);
2242         if (ir->ePBC == epbcXY && ir->nwall == 2)
2243         {
2244             svmul(ir->wall_ewald_zfac, box[ZZ], box[ZZ]);
2245         }
2246         if (ir->nkx > 0 && ir->nky > 0 && ir->nkz > 0)
2247         {
2248             /* Mark fourier_spacing as not used */
2249             ir->fourier_spacing = 0;
2250         }
2251         else if (ir->nkx != 0 && ir->nky != 0 && ir->nkz != 0)
2252         {
2253             set_warning_line(wi, mdparin, -1);
2254             warning_error(wi, "Some of the Fourier grid sizes are set, but all of them need to be set.");
2255         }
2256         const int minGridSize = minimalPmeGridSize(ir->pme_order);
2257         calcFftGrid(stdout, box, ir->fourier_spacing, minGridSize,
2258                     &(ir->nkx), &(ir->nky), &(ir->nkz));
2259         if (ir->nkx < minGridSize ||
2260             ir->nky < minGridSize ||
2261             ir->nkz < minGridSize)
2262         {
2263             warning_error(wi, "The PME grid size should be >= 2*(pme-order - 1); either manually increase the grid size or decrease pme-order");
2264         }
2265     }
2266
2267     /* MRS: eventually figure out better logic for initializing the fep
2268        values that makes declaring the lambda and declaring the state not
2269        potentially conflict if not handled correctly. */
2270     if (ir->efep != efepNO)
2271     {
2272         state.fep_state = ir->fepvals->init_fep_state;
2273         for (i = 0; i < efptNR; i++)
2274         {
2275             /* init_lambda trumps state definitions*/
2276             if (ir->fepvals->init_lambda >= 0)
2277             {
2278                 state.lambda[i] = ir->fepvals->init_lambda;
2279             }
2280             else
2281             {
2282                 if (ir->fepvals->all_lambda[i] == nullptr)
2283                 {
2284                     gmx_fatal(FARGS, "Values of lambda not set for a free energy calculation!");
2285                 }
2286                 else
2287                 {
2288                     state.lambda[i] = ir->fepvals->all_lambda[i][state.fep_state];
2289                 }
2290             }
2291         }
2292     }
2293
2294     struct pull_t *pull = nullptr;
2295
2296     if (ir->bPull)
2297     {
2298         pull = set_pull_init(ir, sys, as_rvec_array(state.x.data()), state.box, state.lambda[efptMASS], oenv);
2299     }
2300
2301     /* Modules that supply external potential for pull coordinates
2302      * should register those potentials here. finish_pull() will check
2303      * that providers have been registerd for all external potentials.
2304      */
2305
2306     if (ir->bPull)
2307     {
2308         finish_pull(pull);
2309     }
2310
2311     if (ir->bRot)
2312     {
2313         set_reference_positions(ir->rot, as_rvec_array(state.x.data()), state.box,
2314                                 opt2fn("-ref", NFILE, fnm), opt2bSet("-ref", NFILE, fnm),
2315                                 wi);
2316     }
2317
2318     /*  reset_multinr(sys); */
2319
2320     if (EEL_PME(ir->coulombtype))
2321     {
2322         float ratio = pme_load_estimate(sys, ir, state.box);
2323         fprintf(stderr, "Estimate for the relative computational load of the PME mesh part: %.2f\n", ratio);
2324         /* With free energy we might need to do PME both for the A and B state
2325          * charges. This will double the cost, but the optimal performance will
2326          * then probably be at a slightly larger cut-off and grid spacing.
2327          */
2328         if ((ir->efep == efepNO && ratio > 1.0/2.0) ||
2329             (ir->efep != efepNO && ratio > 2.0/3.0))
2330         {
2331             warning_note(wi,
2332                          "The optimal PME mesh load for parallel simulations is below 0.5\n"
2333                          "and for highly parallel simulations between 0.25 and 0.33,\n"
2334                          "for higher performance, increase the cut-off and the PME grid spacing.\n");
2335             if (ir->efep != efepNO)
2336             {
2337                 warning_note(wi,
2338                              "For free energy simulations, the optimal load limit increases from 0.5 to 0.667\n");
2339             }
2340         }
2341     }
2342
2343     {
2344         char   warn_buf[STRLEN];
2345         double cio = compute_io(ir, sys->natoms, &sys->groups, F_NRE, 1);
2346         sprintf(warn_buf, "This run will generate roughly %.0f Mb of data", cio);
2347         if (cio > 2000)
2348         {
2349             set_warning_line(wi, mdparin, -1);
2350             warning_note(wi, warn_buf);
2351         }
2352         else
2353         {
2354             printf("%s\n", warn_buf);
2355         }
2356     }
2357
2358     if (bVerbose)
2359     {
2360         fprintf(stderr, "writing run input file...\n");
2361     }
2362
2363     done_warning(wi, FARGS);
2364     write_tpx_state(ftp2fn(efTPR, NFILE, fnm), ir, &state, sys);
2365
2366     /* Output IMD group, if bIMD is TRUE */
2367     write_IMDgroup_to_file(ir->bIMD, ir, &state, sys, NFILE, fnm);
2368
2369     done_atomtype(atype);
2370     done_mtop(sys);
2371     done_inputrec_strings();
2372
2373     return 0;
2374 }