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