src/gromacs/gmxpreprocess/grompp.c

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