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