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