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