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