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