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