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