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