--- /dev/null
- CTYPE ("e.g.: -DI_Want_Cookies -DMe_Too");
+/* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
+ *
+ *
+ * This source code is part of
+ *
+ * G R O M A C S
+ *
+ * GROningen MAchine for Chemical Simulations
+ *
+ * VERSION 3.2.0
+ * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
+ * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
+ * Copyright (c) 2001-2004, The GROMACS development team,
+ * check out http://www.gromacs.org for more information.
+
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version.
+ *
+ * If you want to redistribute modifications, please consider that
+ * scientific software is very special. Version control is crucial -
+ * bugs must be traceable. We will be happy to consider code for
+ * inclusion in the official distribution, but derived work must not
+ * be called official GROMACS. Details are found in the README & COPYING
+ * files - if they are missing, get the official version at www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the papers on the package - you can find them in the top README file.
+ *
+ * For more info, check our website at http://www.gromacs.org
+ *
+ * And Hey:
+ * Gallium Rubidium Oxygen Manganese Argon Carbon Silicon
+ */
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <ctype.h>
+#include <stdlib.h>
+#include <limits.h>
+#include "sysstuff.h"
+#include "smalloc.h"
+#include "typedefs.h"
+#include "physics.h"
+#include "names.h"
+#include "gmx_fatal.h"
+#include "macros.h"
+#include "index.h"
+#include "symtab.h"
+#include "string2.h"
+#include "readinp.h"
+#include "warninp.h"
+#include "readir.h"
+#include "toputil.h"
+#include "index.h"
+#include "network.h"
+#include "vec.h"
+#include "pbc.h"
+#include "mtop_util.h"
+#include "chargegroup.h"
+#include "inputrec.h"
+
+#define MAXPTR 254
+#define NOGID 255
+
+/* Resource parameters
+ * Do not change any of these until you read the instruction
+ * in readinp.h. Some cpp's do not take spaces after the backslash
+ * (like the c-shell), which will give you a very weird compiler
+ * message.
+ */
+
+static char tcgrps[STRLEN],tau_t[STRLEN],ref_t[STRLEN],
+ acc[STRLEN],accgrps[STRLEN],freeze[STRLEN],frdim[STRLEN],
+ energy[STRLEN],user1[STRLEN],user2[STRLEN],vcm[STRLEN],xtc_grps[STRLEN],
+ couple_moltype[STRLEN],orirefitgrp[STRLEN],egptable[STRLEN],egpexcl[STRLEN],
+ wall_atomtype[STRLEN],wall_density[STRLEN],deform[STRLEN],QMMM[STRLEN];
+static char foreign_lambda[STRLEN];
+static char **pull_grp;
+static char **rot_grp;
+static char anneal[STRLEN],anneal_npoints[STRLEN],
+ anneal_time[STRLEN],anneal_temp[STRLEN];
+static char QMmethod[STRLEN],QMbasis[STRLEN],QMcharge[STRLEN],QMmult[STRLEN],
+ bSH[STRLEN],CASorbitals[STRLEN], CASelectrons[STRLEN],SAon[STRLEN],
+ SAoff[STRLEN],SAsteps[STRLEN],bTS[STRLEN],bOPT[STRLEN];
+static char efield_x[STRLEN],efield_xt[STRLEN],efield_y[STRLEN],
+ efield_yt[STRLEN],efield_z[STRLEN],efield_zt[STRLEN];
+
+enum {
+ egrptpALL, /* All particles have to be a member of a group. */
+ egrptpALL_GENREST, /* A rest group with name is generated for particles *
+ * that are not part of any group. */
+ egrptpPART, /* As egrptpALL_GENREST, but no name is generated *
+ * for the rest group. */
+ egrptpONE /* Merge all selected groups into one group, *
+ * make a rest group for the remaining particles. */
+};
+
+
+void init_ir(t_inputrec *ir, t_gromppopts *opts)
+{
+ snew(opts->include,STRLEN);
+ snew(opts->define,STRLEN);
+}
+
+static void _low_check(gmx_bool b,char *s,warninp_t wi)
+{
+ if (b)
+ {
+ warning_error(wi,s);
+ }
+}
+
+static void check_nst(const char *desc_nst,int nst,
+ const char *desc_p,int *p,
+ warninp_t wi)
+{
+ char buf[STRLEN];
+
+ if (*p > 0 && *p % nst != 0)
+ {
+ /* Round up to the next multiple of nst */
+ *p = ((*p)/nst + 1)*nst;
+ sprintf(buf,"%s should be a multiple of %s, changing %s to %d\n",
+ desc_p,desc_nst,desc_p,*p);
+ warning(wi,buf);
+ }
+}
+
+static gmx_bool ir_NVE(const t_inputrec *ir)
+{
+ return ((ir->eI == eiMD || EI_VV(ir->eI)) && ir->etc == etcNO);
+}
+
+static int lcd(int n1,int n2)
+{
+ int d,i;
+
+ d = 1;
+ for(i=2; (i<=n1 && i<=n2); i++)
+ {
+ if (n1 % i == 0 && n2 % i == 0)
+ {
+ d = i;
+ }
+ }
+
+ return d;
+}
+
+void check_ir(const char *mdparin,t_inputrec *ir, t_gromppopts *opts,
+ warninp_t wi)
+/* Check internal consistency */
+{
+ /* Strange macro: first one fills the err_buf, and then one can check
+ * the condition, which will print the message and increase the error
+ * counter.
+ */
+#define CHECK(b) _low_check(b,err_buf,wi)
+ char err_buf[256],warn_buf[STRLEN];
+ int ns_type=0;
+ real dt_pcoupl;
+
+ set_warning_line(wi,mdparin,-1);
+
+ /* BASIC CUT-OFF STUFF */
+ if (ir->rlist == 0 ||
+ !((EEL_MIGHT_BE_ZERO_AT_CUTOFF(ir->coulombtype) && ir->rcoulomb > ir->rlist) ||
+ (EVDW_MIGHT_BE_ZERO_AT_CUTOFF(ir->vdwtype) && ir->rvdw > ir->rlist))) {
+ /* No switched potential and/or no twin-range:
+ * we can set the long-range cut-off to the maximum of the other cut-offs.
+ */
+ ir->rlistlong = max_cutoff(ir->rlist,max_cutoff(ir->rvdw,ir->rcoulomb));
+ } else if (ir->rlistlong < 0) {
+ ir->rlistlong = max_cutoff(ir->rlist,max_cutoff(ir->rvdw,ir->rcoulomb));
+ sprintf(warn_buf,"rlistlong was not set, setting it to %g (no buffer)",
+ ir->rlistlong);
+ warning(wi,warn_buf);
+ }
+ if (ir->rlistlong == 0 && ir->ePBC != epbcNONE) {
+ warning_error(wi,"Can not have an infinite cut-off with PBC");
+ }
+ if (ir->rlistlong > 0 && (ir->rlist == 0 || ir->rlistlong < ir->rlist)) {
+ warning_error(wi,"rlistlong can not be shorter than rlist");
+ }
+ if (IR_TWINRANGE(*ir) && ir->nstlist <= 0) {
+ warning_error(wi,"Can not have nstlist<=0 with twin-range interactions");
+ }
+
+ /* GENERAL INTEGRATOR STUFF */
+ if (!(ir->eI == eiMD || EI_VV(ir->eI)))
+ {
+ ir->etc = etcNO;
+ }
+ if (!EI_DYNAMICS(ir->eI))
+ {
+ ir->epc = epcNO;
+ }
+ if (EI_DYNAMICS(ir->eI))
+ {
+ if (ir->nstcalcenergy < 0)
+ {
+ ir->nstcalcenergy = ir_optimal_nstcalcenergy(ir);
+ if (ir->nstenergy != 0 && ir->nstenergy < ir->nstcalcenergy)
+ {
+ /* nstcalcenergy larger than nstener does not make sense.
+ * We ideally want nstcalcenergy=nstener.
+ */
+ if (ir->nstlist > 0)
+ {
+ ir->nstcalcenergy = lcd(ir->nstenergy,ir->nstlist);
+ }
+ else
+ {
+ ir->nstcalcenergy = ir->nstenergy;
+ }
+ }
+ }
+ if (ir->epc != epcNO)
+ {
+ if (ir->nstpcouple < 0)
+ {
+ ir->nstpcouple = ir_optimal_nstpcouple(ir);
+ }
+ }
+ if (IR_TWINRANGE(*ir))
+ {
+ check_nst("nstlist",ir->nstlist,
+ "nstcalcenergy",&ir->nstcalcenergy,wi);
+ if (ir->epc != epcNO)
+ {
+ check_nst("nstlist",ir->nstlist,
+ "nstpcouple",&ir->nstpcouple,wi);
+ }
+ }
+
+ if (ir->nstcalcenergy > 1)
+ {
+ /* for storing exact averages nstenergy should be
+ * a multiple of nstcalcenergy
+ */
+ check_nst("nstcalcenergy",ir->nstcalcenergy,
+ "nstenergy",&ir->nstenergy,wi);
+ if (ir->efep != efepNO)
+ {
+ /* nstdhdl should be a multiple of nstcalcenergy */
+ check_nst("nstcalcenergy",ir->nstcalcenergy,
+ "nstdhdl",&ir->nstdhdl,wi);
+ }
+ }
+ }
+
+ /* LD STUFF */
+ if ((EI_SD(ir->eI) || ir->eI == eiBD) &&
+ ir->bContinuation && ir->ld_seed != -1) {
+ warning_note(wi,"You are doing a continuation with SD or BD, make sure that ld_seed is different from the previous run (using ld_seed=-1 will ensure this)");
+ }
+
+ /* TPI STUFF */
+ if (EI_TPI(ir->eI)) {
+ sprintf(err_buf,"TPI only works with pbc = %s",epbc_names[epbcXYZ]);
+ CHECK(ir->ePBC != epbcXYZ);
+ sprintf(err_buf,"TPI only works with ns = %s",ens_names[ensGRID]);
+ CHECK(ir->ns_type != ensGRID);
+ sprintf(err_buf,"with TPI nstlist should be larger than zero");
+ CHECK(ir->nstlist <= 0);
+ sprintf(err_buf,"TPI does not work with full electrostatics other than PME");
+ CHECK(EEL_FULL(ir->coulombtype) && !EEL_PME(ir->coulombtype));
+ }
+
+ /* SHAKE / LINCS */
+ if ( (opts->nshake > 0) && (opts->bMorse) ) {
+ sprintf(warn_buf,
+ "Using morse bond-potentials while constraining bonds is useless");
+ warning(wi,warn_buf);
+ }
+
+ sprintf(err_buf,"shake_tol must be > 0 instead of %g while using shake",
+ ir->shake_tol);
+ CHECK(((ir->shake_tol <= 0.0) && (opts->nshake>0) &&
+ (ir->eConstrAlg == econtSHAKE)));
+
+ /* PBC/WALLS */
+ sprintf(err_buf,"walls only work with pbc=%s",epbc_names[epbcXY]);
+ CHECK(ir->nwall && ir->ePBC!=epbcXY);
+
+ /* VACUUM STUFF */
+ if (ir->ePBC != epbcXYZ && ir->nwall != 2) {
+ if (ir->ePBC == epbcNONE) {
+ if (ir->epc != epcNO) {
+ warning(wi,"Turning off pressure coupling for vacuum system");
+ ir->epc = epcNO;
+ }
+ } else {
+ sprintf(err_buf,"Can not have pressure coupling with pbc=%s",
+ epbc_names[ir->ePBC]);
+ CHECK(ir->epc != epcNO);
+ }
+ sprintf(err_buf,"Can not have Ewald with pbc=%s",epbc_names[ir->ePBC]);
+ CHECK(EEL_FULL(ir->coulombtype));
+
+ sprintf(err_buf,"Can not have dispersion correction with pbc=%s",
+ epbc_names[ir->ePBC]);
+ CHECK(ir->eDispCorr != edispcNO);
+ }
+
+ if (ir->rlist == 0.0) {
+ sprintf(err_buf,"can only have neighborlist cut-off zero (=infinite)\n"
+ "with coulombtype = %s or coulombtype = %s\n"
+ "without periodic boundary conditions (pbc = %s) and\n"
+ "rcoulomb and rvdw set to zero",
+ eel_names[eelCUT],eel_names[eelUSER],epbc_names[epbcNONE]);
+ CHECK(((ir->coulombtype != eelCUT) && (ir->coulombtype != eelUSER)) ||
+ (ir->ePBC != epbcNONE) ||
+ (ir->rcoulomb != 0.0) || (ir->rvdw != 0.0));
+
+ if (ir->nstlist < 0) {
+ warning_error(wi,"Can not have heuristic neighborlist updates without cut-off");
+ }
+ if (ir->nstlist > 0) {
+ warning_note(wi,"Simulating without cut-offs is usually (slightly) faster with nstlist=0, nstype=simple and particle decomposition");
+ }
+ }
+
+ /* COMM STUFF */
+ if (ir->nstcomm == 0) {
+ ir->comm_mode = ecmNO;
+ }
+ if (ir->comm_mode != ecmNO) {
+ if (ir->nstcomm < 0) {
+ warning(wi,"If you want to remove the rotation around the center of mass, you should set comm_mode = Angular instead of setting nstcomm < 0. nstcomm is modified to its absolute value");
+ ir->nstcomm = abs(ir->nstcomm);
+ }
+
+ if (ir->nstcalcenergy > 0 && ir->nstcomm < ir->nstcalcenergy) {
+ warning_note(wi,"nstcomm < nstcalcenergy defeats the purpose of nstcalcenergy, setting nstcomm to nstcalcenergy");
+ ir->nstcomm = ir->nstcalcenergy;
+ }
+
+ if (ir->comm_mode == ecmANGULAR) {
+ sprintf(err_buf,"Can not remove the rotation around the center of mass with periodic molecules");
+ CHECK(ir->bPeriodicMols);
+ if (ir->ePBC != epbcNONE)
+ warning(wi,"Removing the rotation around the center of mass in a periodic system (this is not a problem when you have only one molecule).");
+ }
+ }
+
+ if (EI_STATE_VELOCITY(ir->eI) && ir->ePBC == epbcNONE && ir->comm_mode != ecmANGULAR) {
+ warning_note(wi,"Tumbling and or flying ice-cubes: We are not removing rotation around center of mass in a non-periodic system. You should probably set comm_mode = ANGULAR.");
+ }
+
+ sprintf(err_buf,"Free-energy not implemented for Ewald and PPPM");
+ CHECK((ir->coulombtype==eelEWALD || ir->coulombtype==eelPPPM)
+ && (ir->efep!=efepNO));
+
+ sprintf(err_buf,"Twin-range neighbour searching (NS) with simple NS"
+ " algorithm not implemented");
+ CHECK(((ir->rcoulomb > ir->rlist) || (ir->rvdw > ir->rlist))
+ && (ir->ns_type == ensSIMPLE));
+
+ /* TEMPERATURE COUPLING */
+ if (ir->etc == etcYES)
+ {
+ ir->etc = etcBERENDSEN;
+ warning_note(wi,"Old option for temperature coupling given: "
+ "changing \"yes\" to \"Berendsen\"\n");
+ }
+
+ if (ir->etc == etcNOSEHOOVER)
+ {
+ if (ir->opts.nhchainlength < 1)
+ {
+ sprintf(warn_buf,"number of Nose-Hoover chains (currently %d) cannot be less than 1,reset to 1\n",ir->opts.nhchainlength);
+ ir->opts.nhchainlength =1;
+ warning(wi,warn_buf);
+ }
+
+ if (ir->etc==etcNOSEHOOVER && !EI_VV(ir->eI) && ir->opts.nhchainlength > 1)
+ {
+ warning_note(wi,"leapfrog does not yet support Nose-Hoover chains, nhchainlength reset to 1");
+ ir->opts.nhchainlength = 1;
+ }
+ }
+ else
+ {
+ ir->opts.nhchainlength = 0;
+ }
+
+ if (ir->etc == etcBERENDSEN)
+ {
+ sprintf(warn_buf,"The %s thermostat does not generate the correct kinetic energy distribution. You might want to consider using the %s thermostat.",
+ ETCOUPLTYPE(ir->etc),ETCOUPLTYPE(etcVRESCALE));
+ warning_note(wi,warn_buf);
+ }
+
+ if ((ir->etc==etcNOSEHOOVER || ir->etc==etcANDERSEN || ir->etc==etcANDERSENINTERVAL)
+ && ir->epc==epcBERENDSEN)
+ {
+ sprintf(warn_buf,"Using Berendsen pressure coupling invalidates the "
+ "true ensemble for the thermostat");
+ warning(wi,warn_buf);
+ }
+
+ /* PRESSURE COUPLING */
+ if (ir->epc == epcISOTROPIC)
+ {
+ ir->epc = epcBERENDSEN;
+ warning_note(wi,"Old option for pressure coupling given: "
+ "changing \"Isotropic\" to \"Berendsen\"\n");
+ }
+
+ if (ir->epc != epcNO)
+ {
+ dt_pcoupl = ir->nstpcouple*ir->delta_t;
+
+ sprintf(err_buf,"tau_p must be > 0 instead of %g\n",ir->tau_p);
+ CHECK(ir->tau_p <= 0);
+
+ if (ir->tau_p/dt_pcoupl < pcouple_min_integration_steps(ir->epc))
+ {
+ sprintf(warn_buf,"For proper integration of the %s barostat, tau_p (%g) should be at least %d times larger than nstpcouple*dt (%g)",
+ EPCOUPLTYPE(ir->epc),ir->tau_p,pcouple_min_integration_steps(ir->epc),dt_pcoupl);
+ warning(wi,warn_buf);
+ }
+
+ sprintf(err_buf,"compressibility must be > 0 when using pressure"
+ " coupling %s\n",EPCOUPLTYPE(ir->epc));
+ CHECK(ir->compress[XX][XX] < 0 || ir->compress[YY][YY] < 0 ||
+ ir->compress[ZZ][ZZ] < 0 ||
+ (trace(ir->compress) == 0 && ir->compress[YY][XX] <= 0 &&
+ ir->compress[ZZ][XX] <= 0 && ir->compress[ZZ][YY] <= 0));
+
+ sprintf(err_buf,"pressure coupling with PPPM not implemented, use PME");
+ CHECK(ir->coulombtype == eelPPPM);
+
+ }
+ else if (ir->coulombtype == eelPPPM)
+ {
+ sprintf(warn_buf,"The pressure with PPPM is incorrect, if you need the pressure use PME");
+ warning(wi,warn_buf);
+ }
+
+ if (EI_VV(ir->eI))
+ {
+ if (ir->epc > epcNO)
+ {
+ if (ir->epc!=epcMTTK)
+ {
+ warning_error(wi,"NPT only defined for vv using Martyna-Tuckerman-Tobias-Klein equations");
+ }
+ }
+ }
+
+ /* ELECTROSTATICS */
+ /* More checks are in triple check (grompp.c) */
+ if (ir->coulombtype == eelPPPM)
+ {
+ warning_error(wi,"PPPM is not functional in the current version, we plan to implement PPPM through a small modification of the PME code");
+ }
+
+ if (ir->coulombtype == eelSWITCH) {
+ sprintf(warn_buf,"coulombtype = %s is only for testing purposes and can lead to serious artifacts, advice: use coulombtype = %s",
+ eel_names[ir->coulombtype],
+ eel_names[eelRF_ZERO]);
+ warning(wi,warn_buf);
+ }
+
+ if (ir->epsilon_r!=1 && ir->implicit_solvent==eisGBSA) {
+ sprintf(warn_buf,"epsilon_r = %g with GB implicit solvent, will use this value for inner dielectric",ir->epsilon_r);
+ warning_note(wi,warn_buf);
+ }
+
+ if (EEL_RF(ir->coulombtype) && ir->epsilon_rf==1 && ir->epsilon_r!=1) {
+ sprintf(warn_buf,"epsilon_r = %g and epsilon_rf = 1 with reaction field, assuming old format and exchanging epsilon_r and epsilon_rf",ir->epsilon_r);
+ warning(wi,warn_buf);
+ ir->epsilon_rf = ir->epsilon_r;
+ ir->epsilon_r = 1.0;
+ }
+
+ if (getenv("GALACTIC_DYNAMICS") == NULL) {
+ sprintf(err_buf,"epsilon_r must be >= 0 instead of %g\n",ir->epsilon_r);
+ CHECK(ir->epsilon_r < 0);
+ }
+
+ if (EEL_RF(ir->coulombtype)) {
+ /* reaction field (at the cut-off) */
+
+ if (ir->coulombtype == eelRF_ZERO) {
+ sprintf(err_buf,"With coulombtype = %s, epsilon_rf must be 0",
+ eel_names[ir->coulombtype]);
+ CHECK(ir->epsilon_rf != 0);
+ }
+
+ sprintf(err_buf,"epsilon_rf must be >= epsilon_r");
+ CHECK((ir->epsilon_rf < ir->epsilon_r && ir->epsilon_rf != 0) ||
+ (ir->epsilon_r == 0));
+ if (ir->epsilon_rf == ir->epsilon_r) {
+ sprintf(warn_buf,"Using epsilon_rf = epsilon_r with %s does not make sense",
+ eel_names[ir->coulombtype]);
+ warning(wi,warn_buf);
+ }
+ }
+ /* Allow rlist>rcoulomb for tabulated long range stuff. This just
+ * means the interaction is zero outside rcoulomb, but it helps to
+ * provide accurate energy conservation.
+ */
+ if (EEL_MIGHT_BE_ZERO_AT_CUTOFF(ir->coulombtype)) {
+ if (EEL_SWITCHED(ir->coulombtype)) {
+ sprintf(err_buf,
+ "With coulombtype = %s rcoulomb_switch must be < rcoulomb",
+ eel_names[ir->coulombtype]);
+ CHECK(ir->rcoulomb_switch >= ir->rcoulomb);
+ }
+ } else if (ir->coulombtype == eelCUT || EEL_RF(ir->coulombtype)) {
+ sprintf(err_buf,"With coulombtype = %s, rcoulomb must be >= rlist",
+ eel_names[ir->coulombtype]);
+ CHECK(ir->rlist > ir->rcoulomb);
+ }
+
+ if (EEL_FULL(ir->coulombtype)) {
+ if (ir->coulombtype==eelPMESWITCH || ir->coulombtype==eelPMEUSER ||
+ ir->coulombtype==eelPMEUSERSWITCH) {
+ sprintf(err_buf,"With coulombtype = %s, rcoulomb must be <= rlist",
+ eel_names[ir->coulombtype]);
+ CHECK(ir->rcoulomb > ir->rlist);
+ } else {
+ if (ir->coulombtype == eelPME) {
+ sprintf(err_buf,
+ "With coulombtype = %s, rcoulomb must be equal to rlist\n"
+ "If you want optimal energy conservation or exact integration use %s",
+ eel_names[ir->coulombtype],eel_names[eelPMESWITCH]);
+ } else {
+ sprintf(err_buf,
+ "With coulombtype = %s, rcoulomb must be equal to rlist",
+ eel_names[ir->coulombtype]);
+ }
+ CHECK(ir->rcoulomb != ir->rlist);
+ }
+ }
+
+ if (EEL_PME(ir->coulombtype)) {
+ if (ir->pme_order < 3) {
+ warning_error(wi,"pme_order can not be smaller than 3");
+ }
+ }
+
+ if (ir->nwall==2 && EEL_FULL(ir->coulombtype)) {
+ if (ir->ewald_geometry == eewg3D) {
+ sprintf(warn_buf,"With pbc=%s you should use ewald_geometry=%s",
+ epbc_names[ir->ePBC],eewg_names[eewg3DC]);
+ warning(wi,warn_buf);
+ }
+ /* This check avoids extra pbc coding for exclusion corrections */
+ sprintf(err_buf,"wall_ewald_zfac should be >= 2");
+ CHECK(ir->wall_ewald_zfac < 2);
+ }
+
+ if (EVDW_SWITCHED(ir->vdwtype)) {
+ sprintf(err_buf,"With vdwtype = %s rvdw_switch must be < rvdw",
+ evdw_names[ir->vdwtype]);
+ CHECK(ir->rvdw_switch >= ir->rvdw);
+ } else if (ir->vdwtype == evdwCUT) {
+ sprintf(err_buf,"With vdwtype = %s, rvdw must be >= rlist",evdw_names[ir->vdwtype]);
+ CHECK(ir->rlist > ir->rvdw);
+ }
+ if (EEL_IS_ZERO_AT_CUTOFF(ir->coulombtype)
+ && (ir->rlistlong <= ir->rcoulomb)) {
+ sprintf(warn_buf,"For energy conservation with switch/shift potentials, %s should be 0.1 to 0.3 nm larger than rcoulomb.",
+ IR_TWINRANGE(*ir) ? "rlistlong" : "rlist");
+ warning_note(wi,warn_buf);
+ }
+ if (EVDW_SWITCHED(ir->vdwtype) && (ir->rlistlong <= ir->rvdw)) {
+ sprintf(warn_buf,"For energy conservation with switch/shift potentials, %s should be 0.1 to 0.3 nm larger than rvdw.",
+ IR_TWINRANGE(*ir) ? "rlistlong" : "rlist");
+ warning_note(wi,warn_buf);
+ }
+
+ if (ir->vdwtype == evdwUSER && ir->eDispCorr != edispcNO) {
+ warning_note(wi,"You have selected user tables with dispersion correction, the dispersion will be corrected to -C6/r^6 beyond rvdw_switch (the tabulated interaction between rvdw_switch and rvdw will not be double counted). Make sure that you really want dispersion correction to -C6/r^6.");
+ }
+
+ if (ir->nstlist == -1) {
+ sprintf(err_buf,
+ "nstlist=-1 only works with switched or shifted potentials,\n"
+ "suggestion: use vdw-type=%s and coulomb-type=%s",
+ evdw_names[evdwSHIFT],eel_names[eelPMESWITCH]);
+ CHECK(!(EEL_MIGHT_BE_ZERO_AT_CUTOFF(ir->coulombtype) &&
+ EVDW_MIGHT_BE_ZERO_AT_CUTOFF(ir->vdwtype)));
+
+ sprintf(err_buf,"With nstlist=-1 rvdw and rcoulomb should be smaller than rlist to account for diffusion and possibly charge-group radii");
+ CHECK(ir->rvdw >= ir->rlist || ir->rcoulomb >= ir->rlist);
+ }
+ sprintf(err_buf,"nstlist can not be smaller than -1");
+ CHECK(ir->nstlist < -1);
+
+ if (ir->eI == eiLBFGS && (ir->coulombtype==eelCUT || ir->vdwtype==evdwCUT)
+ && ir->rvdw != 0) {
+ warning(wi,"For efficient BFGS minimization, use switch/shift/pme instead of cut-off.");
+ }
+
+ if (ir->eI == eiLBFGS && ir->nbfgscorr <= 0) {
+ warning(wi,"Using L-BFGS with nbfgscorr<=0 just gets you steepest descent.");
+ }
+
+ /* FREE ENERGY */
+ if (ir->efep != efepNO) {
+ sprintf(err_buf,"The soft-core power is %d and can only be 1 or 2",
+ ir->sc_power);
+ CHECK(ir->sc_alpha!=0 && ir->sc_power!=1 && ir->sc_power!=2);
+ }
+
+ /* ENERGY CONSERVATION */
+ if (ir_NVE(ir))
+ {
+ if (!EVDW_MIGHT_BE_ZERO_AT_CUTOFF(ir->vdwtype) && ir->rvdw > 0)
+ {
+ sprintf(warn_buf,"You are using a cut-off for VdW interactions with NVE, for good energy conservation use vdwtype = %s (possibly with DispCorr)",
+ evdw_names[evdwSHIFT]);
+ warning_note(wi,warn_buf);
+ }
+ if (!EEL_MIGHT_BE_ZERO_AT_CUTOFF(ir->coulombtype) && ir->rcoulomb > 0)
+ {
+ sprintf(warn_buf,"You are using a cut-off for electrostatics with NVE, for good energy conservation use coulombtype = %s or %s",
+ eel_names[eelPMESWITCH],eel_names[eelRF_ZERO]);
+ warning_note(wi,warn_buf);
+ }
+ }
+
+ /* IMPLICIT SOLVENT */
+ if(ir->coulombtype==eelGB_NOTUSED)
+ {
+ ir->coulombtype=eelCUT;
+ ir->implicit_solvent=eisGBSA;
+ fprintf(stderr,"Note: Old option for generalized born electrostatics given:\n"
+ "Changing coulombtype from \"generalized-born\" to \"cut-off\" and instead\n"
+ "setting implicit_solvent value to \"GBSA\" in input section.\n");
+ }
+
+ if(ir->sa_algorithm==esaSTILL)
+ {
+ sprintf(err_buf,"Still SA algorithm not available yet, use %s or %s instead\n",esa_names[esaAPPROX],esa_names[esaNO]);
+ CHECK(ir->sa_algorithm == esaSTILL);
+ }
+
+ if(ir->implicit_solvent==eisGBSA)
+ {
+ sprintf(err_buf,"With GBSA implicit solvent, rgbradii must be equal to rlist.");
+ CHECK(ir->rgbradii != ir->rlist);
+
+ if(ir->coulombtype!=eelCUT)
+ {
+ sprintf(err_buf,"With GBSA, coulombtype must be equal to %s\n",eel_names[eelCUT]);
+ CHECK(ir->coulombtype!=eelCUT);
+ }
+ if(ir->vdwtype!=evdwCUT)
+ {
+ sprintf(err_buf,"With GBSA, vdw-type must be equal to %s\n",evdw_names[evdwCUT]);
+ CHECK(ir->vdwtype!=evdwCUT);
+ }
+ if(ir->nstgbradii<1)
+ {
+ sprintf(warn_buf,"Using GBSA with nstgbradii<1, setting nstgbradii=1");
+ warning_note(wi,warn_buf);
+ ir->nstgbradii=1;
+ }
+ if(ir->sa_algorithm==esaNO)
+ {
+ sprintf(warn_buf,"No SA (non-polar) calculation requested together with GB. Are you sure this is what you want?\n");
+ warning_note(wi,warn_buf);
+ }
+ if(ir->sa_surface_tension<0 && ir->sa_algorithm!=esaNO)
+ {
+ sprintf(warn_buf,"Value of sa_surface_tension is < 0. Changing it to 2.05016 or 2.25936 kJ/nm^2/mol for Still and HCT/OBC respectively\n");
+ warning_note(wi,warn_buf);
+
+ if(ir->gb_algorithm==egbSTILL)
+ {
+ ir->sa_surface_tension = 0.0049 * CAL2JOULE * 100;
+ }
+ else
+ {
+ ir->sa_surface_tension = 0.0054 * CAL2JOULE * 100;
+ }
+ }
+ if(ir->sa_surface_tension==0 && ir->sa_algorithm!=esaNO)
+ {
+ sprintf(err_buf, "Surface tension set to 0 while SA-calculation requested\n");
+ CHECK(ir->sa_surface_tension==0 && ir->sa_algorithm!=esaNO);
+ }
+
+ }
+}
+
+static int str_nelem(const char *str,int maxptr,char *ptr[])
+{
+ int np=0;
+ char *copy0,*copy;
+
+ copy0=strdup(str);
+ copy=copy0;
+ ltrim(copy);
+ while (*copy != '\0') {
+ if (np >= maxptr)
+ gmx_fatal(FARGS,"Too many groups on line: '%s' (max is %d)",
+ str,maxptr);
+ if (ptr)
+ ptr[np]=copy;
+ np++;
+ while ((*copy != '\0') && !isspace(*copy))
+ copy++;
+ if (*copy != '\0') {
+ *copy='\0';
+ copy++;
+ }
+ ltrim(copy);
+ }
+ if (ptr == NULL)
+ sfree(copy0);
+
+ return np;
+}
+
+static void parse_n_double(char *str,int *n,double **r)
+{
+ char *ptr[MAXPTR];
+ int i;
+
+ *n = str_nelem(str,MAXPTR,ptr);
+
+ snew(*r,*n);
+ for(i=0; i<*n; i++) {
+ (*r)[i] = strtod(ptr[i],NULL);
+ }
+}
+
+static void do_wall_params(t_inputrec *ir,
+ char *wall_atomtype, char *wall_density,
+ t_gromppopts *opts)
+{
+ int nstr,i;
+ char *names[MAXPTR];
+ double dbl;
+
+ opts->wall_atomtype[0] = NULL;
+ opts->wall_atomtype[1] = NULL;
+
+ ir->wall_atomtype[0] = -1;
+ ir->wall_atomtype[1] = -1;
+ ir->wall_density[0] = 0;
+ ir->wall_density[1] = 0;
+
+ if (ir->nwall > 0)
+ {
+ nstr = str_nelem(wall_atomtype,MAXPTR,names);
+ if (nstr != ir->nwall)
+ {
+ gmx_fatal(FARGS,"Expected %d elements for wall_atomtype, found %d",
+ ir->nwall,nstr);
+ }
+ for(i=0; i<ir->nwall; i++)
+ {
+ opts->wall_atomtype[i] = strdup(names[i]);
+ }
+
+ if (ir->wall_type == ewt93 || ir->wall_type == ewt104) {
+ nstr = str_nelem(wall_density,MAXPTR,names);
+ if (nstr != ir->nwall)
+ {
+ gmx_fatal(FARGS,"Expected %d elements for wall_density, found %d",ir->nwall,nstr);
+ }
+ for(i=0; i<ir->nwall; i++)
+ {
+ sscanf(names[i],"%lf",&dbl);
+ if (dbl <= 0)
+ {
+ gmx_fatal(FARGS,"wall_density[%d] = %f\n",i,dbl);
+ }
+ ir->wall_density[i] = dbl;
+ }
+ }
+ }
+}
+
+static void add_wall_energrps(gmx_groups_t *groups,int nwall,t_symtab *symtab)
+{
+ int i;
+ t_grps *grps;
+ char str[STRLEN];
+
+ if (nwall > 0) {
+ srenew(groups->grpname,groups->ngrpname+nwall);
+ grps = &(groups->grps[egcENER]);
+ srenew(grps->nm_ind,grps->nr+nwall);
+ for(i=0; i<nwall; i++) {
+ sprintf(str,"wall%d",i);
+ groups->grpname[groups->ngrpname] = put_symtab(symtab,str);
+ grps->nm_ind[grps->nr++] = groups->ngrpname++;
+ }
+ }
+}
+
+void get_ir(const char *mdparin,const char *mdparout,
+ t_inputrec *ir,t_gromppopts *opts,
+ warninp_t wi)
+{
+ char *dumstr[2];
+ double dumdub[2][6];
+ t_inpfile *inp;
+ const char *tmp;
+ int i,j,m,ninp;
+ char warn_buf[STRLEN];
+
+ inp = read_inpfile(mdparin, &ninp, NULL, wi);
+
+ snew(dumstr[0],STRLEN);
+ snew(dumstr[1],STRLEN);
+
+ REM_TYPE("title");
+ REM_TYPE("cpp");
+ REM_TYPE("domain-decomposition");
+ REPL_TYPE("unconstrained-start","continuation");
+ REM_TYPE("dihre-tau");
+ REM_TYPE("nstdihreout");
+ REM_TYPE("nstcheckpoint");
+
+ CCTYPE ("VARIOUS PREPROCESSING OPTIONS");
+ CTYPE ("Preprocessor information: use cpp syntax.");
+ CTYPE ("e.g.: -I/home/joe/doe -I/home/mary/roe");
+ STYPE ("include", opts->include, NULL);
++ CTYPE ("e.g.: -DPOSRES -DFLEXIBLE (note these variable names are case sensitive)");
+ STYPE ("define", opts->define, NULL);
+
+ CCTYPE ("RUN CONTROL PARAMETERS");
+ EETYPE("integrator", ir->eI, ei_names);
+ CTYPE ("Start time and timestep in ps");
+ RTYPE ("tinit", ir->init_t, 0.0);
+ RTYPE ("dt", ir->delta_t, 0.001);
+ STEPTYPE ("nsteps", ir->nsteps, 0);
+ CTYPE ("For exact run continuation or redoing part of a run");
+ STEPTYPE ("init_step",ir->init_step, 0);
+ CTYPE ("Part index is updated automatically on checkpointing (keeps files separate)");
+ ITYPE ("simulation_part", ir->simulation_part, 1);
+ CTYPE ("mode for center of mass motion removal");
+ EETYPE("comm-mode", ir->comm_mode, ecm_names);
+ CTYPE ("number of steps for center of mass motion removal");
+ ITYPE ("nstcomm", ir->nstcomm, 10);
+ CTYPE ("group(s) for center of mass motion removal");
+ STYPE ("comm-grps", vcm, NULL);
+
+ CCTYPE ("LANGEVIN DYNAMICS OPTIONS");
+ CTYPE ("Friction coefficient (amu/ps) and random seed");
+ RTYPE ("bd-fric", ir->bd_fric, 0.0);
+ ITYPE ("ld-seed", ir->ld_seed, 1993);
+
+ /* Em stuff */
+ CCTYPE ("ENERGY MINIMIZATION OPTIONS");
+ CTYPE ("Force tolerance and initial step-size");
+ RTYPE ("emtol", ir->em_tol, 10.0);
+ RTYPE ("emstep", ir->em_stepsize,0.01);
+ CTYPE ("Max number of iterations in relax_shells");
+ ITYPE ("niter", ir->niter, 20);
+ CTYPE ("Step size (ps^2) for minimization of flexible constraints");
+ RTYPE ("fcstep", ir->fc_stepsize, 0);
+ CTYPE ("Frequency of steepest descents steps when doing CG");
+ ITYPE ("nstcgsteep", ir->nstcgsteep, 1000);
+ ITYPE ("nbfgscorr", ir->nbfgscorr, 10);
+
+ CCTYPE ("TEST PARTICLE INSERTION OPTIONS");
+ RTYPE ("rtpi", ir->rtpi, 0.05);
+
+ /* Output options */
+ CCTYPE ("OUTPUT CONTROL OPTIONS");
+ CTYPE ("Output frequency for coords (x), velocities (v) and forces (f)");
+ ITYPE ("nstxout", ir->nstxout, 100);
+ ITYPE ("nstvout", ir->nstvout, 100);
+ ITYPE ("nstfout", ir->nstfout, 0);
+ ir->nstcheckpoint = 1000;
+ CTYPE ("Output frequency for energies to log file and energy file");
+ ITYPE ("nstlog", ir->nstlog, 100);
+ ITYPE ("nstcalcenergy",ir->nstcalcenergy, -1);
+ ITYPE ("nstenergy", ir->nstenergy, 100);
+ CTYPE ("Output frequency and precision for xtc file");
+ ITYPE ("nstxtcout", ir->nstxtcout, 0);
+ RTYPE ("xtc-precision",ir->xtcprec, 1000.0);
+ CTYPE ("This selects the subset of atoms for the xtc file. You can");
+ CTYPE ("select multiple groups. By default all atoms will be written.");
+ STYPE ("xtc-grps", xtc_grps, NULL);
+ CTYPE ("Selection of energy groups");
+ STYPE ("energygrps", energy, NULL);
+
+ /* Neighbor searching */
+ CCTYPE ("NEIGHBORSEARCHING PARAMETERS");
+ CTYPE ("nblist update frequency");
+ ITYPE ("nstlist", ir->nstlist, 10);
+ CTYPE ("ns algorithm (simple or grid)");
+ EETYPE("ns-type", ir->ns_type, ens_names);
+ /* set ndelta to the optimal value of 2 */
+ ir->ndelta = 2;
+ CTYPE ("Periodic boundary conditions: xyz, no, xy");
+ EETYPE("pbc", ir->ePBC, epbc_names);
+ EETYPE("periodic_molecules", ir->bPeriodicMols, yesno_names);
+ CTYPE ("nblist cut-off");
+ RTYPE ("rlist", ir->rlist, 1.0);
+ CTYPE ("long-range cut-off for switched potentials");
+ RTYPE ("rlistlong", ir->rlistlong, -1);
+
+ /* Electrostatics */
+ CCTYPE ("OPTIONS FOR ELECTROSTATICS AND VDW");
+ CTYPE ("Method for doing electrostatics");
+ EETYPE("coulombtype", ir->coulombtype, eel_names);
+ CTYPE ("cut-off lengths");
+ RTYPE ("rcoulomb-switch", ir->rcoulomb_switch, 0.0);
+ RTYPE ("rcoulomb", ir->rcoulomb, 1.0);
+ CTYPE ("Relative dielectric constant for the medium and the reaction field");
+ RTYPE ("epsilon_r", ir->epsilon_r, 1.0);
+ RTYPE ("epsilon_rf", ir->epsilon_rf, 1.0);
+ CTYPE ("Method for doing Van der Waals");
+ EETYPE("vdw-type", ir->vdwtype, evdw_names);
+ CTYPE ("cut-off lengths");
+ RTYPE ("rvdw-switch", ir->rvdw_switch, 0.0);
+ RTYPE ("rvdw", ir->rvdw, 1.0);
+ CTYPE ("Apply long range dispersion corrections for Energy and Pressure");
+ EETYPE("DispCorr", ir->eDispCorr, edispc_names);
+ CTYPE ("Extension of the potential lookup tables beyond the cut-off");
+ RTYPE ("table-extension", ir->tabext, 1.0);
+ CTYPE ("Seperate tables between energy group pairs");
+ STYPE ("energygrp_table", egptable, NULL);
+ CTYPE ("Spacing for the PME/PPPM FFT grid");
+ RTYPE ("fourierspacing", opts->fourierspacing,0.12);
+ CTYPE ("FFT grid size, when a value is 0 fourierspacing will be used");
+ ITYPE ("fourier_nx", ir->nkx, 0);
+ ITYPE ("fourier_ny", ir->nky, 0);
+ ITYPE ("fourier_nz", ir->nkz, 0);
+ CTYPE ("EWALD/PME/PPPM parameters");
+ ITYPE ("pme_order", ir->pme_order, 4);
+ RTYPE ("ewald_rtol", ir->ewald_rtol, 0.00001);
+ EETYPE("ewald_geometry", ir->ewald_geometry, eewg_names);
+ RTYPE ("epsilon_surface", ir->epsilon_surface, 0.0);
+ EETYPE("optimize_fft",ir->bOptFFT, yesno_names);
+
+ CCTYPE("IMPLICIT SOLVENT ALGORITHM");
+ EETYPE("implicit_solvent", ir->implicit_solvent, eis_names);
+
+ CCTYPE ("GENERALIZED BORN ELECTROSTATICS");
+ CTYPE ("Algorithm for calculating Born radii");
+ EETYPE("gb_algorithm", ir->gb_algorithm, egb_names);
+ CTYPE ("Frequency of calculating the Born radii inside rlist");
+ ITYPE ("nstgbradii", ir->nstgbradii, 1);
+ CTYPE ("Cutoff for Born radii calculation; the contribution from atoms");
+ CTYPE ("between rlist and rgbradii is updated every nstlist steps");
+ RTYPE ("rgbradii", ir->rgbradii, 1.0);
+ CTYPE ("Dielectric coefficient of the implicit solvent");
+ RTYPE ("gb_epsilon_solvent",ir->gb_epsilon_solvent, 80.0);
+ CTYPE ("Salt concentration in M for Generalized Born models");
+ RTYPE ("gb_saltconc", ir->gb_saltconc, 0.0);
+ CTYPE ("Scaling factors used in the OBC GB model. Default values are OBC(II)");
+ RTYPE ("gb_obc_alpha", ir->gb_obc_alpha, 1.0);
+ RTYPE ("gb_obc_beta", ir->gb_obc_beta, 0.8);
+ RTYPE ("gb_obc_gamma", ir->gb_obc_gamma, 4.85);
+ RTYPE ("gb_dielectric_offset", ir->gb_dielectric_offset, 0.009);
+ EETYPE("sa_algorithm", ir->sa_algorithm, esa_names);
+ CTYPE ("Surface tension (kJ/mol/nm^2) for the SA (nonpolar surface) part of GBSA");
+ CTYPE ("The value -1 will set default value for Still/HCT/OBC GB-models.");
+ RTYPE ("sa_surface_tension", ir->sa_surface_tension, -1);
+
+ /* Coupling stuff */
+ CCTYPE ("OPTIONS FOR WEAK COUPLING ALGORITHMS");
+ CTYPE ("Temperature coupling");
+ EETYPE("tcoupl", ir->etc, etcoupl_names);
+ ITYPE ("nsttcouple", ir->nsttcouple, -1);
+ ITYPE("nh-chain-length", ir->opts.nhchainlength, NHCHAINLENGTH);
+ CTYPE ("Groups to couple separately");
+ STYPE ("tc-grps", tcgrps, NULL);
+ CTYPE ("Time constant (ps) and reference temperature (K)");
+ STYPE ("tau-t", tau_t, NULL);
+ STYPE ("ref-t", ref_t, NULL);
+ CTYPE ("Pressure coupling");
+ EETYPE("Pcoupl", ir->epc, epcoupl_names);
+ EETYPE("Pcoupltype", ir->epct, epcoupltype_names);
+ ITYPE ("nstpcouple", ir->nstpcouple, -1);
+ CTYPE ("Time constant (ps), compressibility (1/bar) and reference P (bar)");
+ RTYPE ("tau-p", ir->tau_p, 1.0);
+ STYPE ("compressibility", dumstr[0], NULL);
+ STYPE ("ref-p", dumstr[1], NULL);
+ CTYPE ("Scaling of reference coordinates, No, All or COM");
+ EETYPE ("refcoord_scaling",ir->refcoord_scaling,erefscaling_names);
+
+ CTYPE ("Random seed for Andersen thermostat");
+ ITYPE ("andersen_seed", ir->andersen_seed, 815131);
+
+ /* QMMM */
+ CCTYPE ("OPTIONS FOR QMMM calculations");
+ EETYPE("QMMM", ir->bQMMM, yesno_names);
+ CTYPE ("Groups treated Quantum Mechanically");
+ STYPE ("QMMM-grps", QMMM, NULL);
+ CTYPE ("QM method");
+ STYPE("QMmethod", QMmethod, NULL);
+ CTYPE ("QMMM scheme");
+ EETYPE("QMMMscheme", ir->QMMMscheme, eQMMMscheme_names);
+ CTYPE ("QM basisset");
+ STYPE("QMbasis", QMbasis, NULL);
+ CTYPE ("QM charge");
+ STYPE ("QMcharge", QMcharge,NULL);
+ CTYPE ("QM multiplicity");
+ STYPE ("QMmult", QMmult,NULL);
+ CTYPE ("Surface Hopping");
+ STYPE ("SH", bSH, NULL);
+ CTYPE ("CAS space options");
+ STYPE ("CASorbitals", CASorbitals, NULL);
+ STYPE ("CASelectrons", CASelectrons, NULL);
+ STYPE ("SAon", SAon, NULL);
+ STYPE ("SAoff",SAoff,NULL);
+ STYPE ("SAsteps", SAsteps, NULL);
+ CTYPE ("Scale factor for MM charges");
+ RTYPE ("MMChargeScaleFactor", ir->scalefactor, 1.0);
+ CTYPE ("Optimization of QM subsystem");
+ STYPE ("bOPT", bOPT, NULL);
+ STYPE ("bTS", bTS, NULL);
+
+ /* Simulated annealing */
+ CCTYPE("SIMULATED ANNEALING");
+ CTYPE ("Type of annealing for each temperature group (no/single/periodic)");
+ STYPE ("annealing", anneal, NULL);
+ CTYPE ("Number of time points to use for specifying annealing in each group");
+ STYPE ("annealing_npoints", anneal_npoints, NULL);
+ CTYPE ("List of times at the annealing points for each group");
+ STYPE ("annealing_time", anneal_time, NULL);
+ CTYPE ("Temp. at each annealing point, for each group.");
+ STYPE ("annealing_temp", anneal_temp, NULL);
+
+ /* Startup run */
+ CCTYPE ("GENERATE VELOCITIES FOR STARTUP RUN");
+ EETYPE("gen-vel", opts->bGenVel, yesno_names);
+ RTYPE ("gen-temp", opts->tempi, 300.0);
+ ITYPE ("gen-seed", opts->seed, 173529);
+
+ /* Shake stuff */
+ CCTYPE ("OPTIONS FOR BONDS");
+ EETYPE("constraints", opts->nshake, constraints);
+ CTYPE ("Type of constraint algorithm");
+ EETYPE("constraint-algorithm", ir->eConstrAlg, econstr_names);
+ CTYPE ("Do not constrain the start configuration");
+ EETYPE("continuation", ir->bContinuation, yesno_names);
+ CTYPE ("Use successive overrelaxation to reduce the number of shake iterations");
+ EETYPE("Shake-SOR", ir->bShakeSOR, yesno_names);
+ CTYPE ("Relative tolerance of shake");
+ RTYPE ("shake-tol", ir->shake_tol, 0.0001);
+ CTYPE ("Highest order in the expansion of the constraint coupling matrix");
+ ITYPE ("lincs-order", ir->nProjOrder, 4);
+ CTYPE ("Number of iterations in the final step of LINCS. 1 is fine for");
+ CTYPE ("normal simulations, but use 2 to conserve energy in NVE runs.");
+ CTYPE ("For energy minimization with constraints it should be 4 to 8.");
+ ITYPE ("lincs-iter", ir->nLincsIter, 1);
+ CTYPE ("Lincs will write a warning to the stderr if in one step a bond");
+ CTYPE ("rotates over more degrees than");
+ RTYPE ("lincs-warnangle", ir->LincsWarnAngle, 30.0);
+ CTYPE ("Convert harmonic bonds to morse potentials");
+ EETYPE("morse", opts->bMorse,yesno_names);
+
+ /* Energy group exclusions */
+ CCTYPE ("ENERGY GROUP EXCLUSIONS");
+ CTYPE ("Pairs of energy groups for which all non-bonded interactions are excluded");
+ STYPE ("energygrp_excl", egpexcl, NULL);
+
+ /* Walls */
+ CCTYPE ("WALLS");
+ CTYPE ("Number of walls, type, atom types, densities and box-z scale factor for Ewald");
+ ITYPE ("nwall", ir->nwall, 0);
+ EETYPE("wall_type", ir->wall_type, ewt_names);
+ RTYPE ("wall_r_linpot", ir->wall_r_linpot, -1);
+ STYPE ("wall_atomtype", wall_atomtype, NULL);
+ STYPE ("wall_density", wall_density, NULL);
+ RTYPE ("wall_ewald_zfac", ir->wall_ewald_zfac, 3);
+
+ /* COM pulling */
+ CCTYPE("COM PULLING");
+ CTYPE("Pull type: no, umbrella, constraint or constant_force");
+ EETYPE("pull", ir->ePull, epull_names);
+ if (ir->ePull != epullNO) {
+ snew(ir->pull,1);
+ pull_grp = read_pullparams(&ninp,&inp,ir->pull,&opts->pull_start,wi);
+ }
+
+ /* Enforced rotation */
+ CCTYPE("ENFORCED ROTATION");
+ CTYPE("Enforced rotation: No or Yes");
+ EETYPE("rotation", ir->bRot, yesno_names);
+ if (ir->bRot) {
+ snew(ir->rot,1);
+ rot_grp = read_rotparams(&ninp,&inp,ir->rot,wi);
+ }
+
+ /* Refinement */
+ CCTYPE("NMR refinement stuff");
+ CTYPE ("Distance restraints type: No, Simple or Ensemble");
+ EETYPE("disre", ir->eDisre, edisre_names);
+ CTYPE ("Force weighting of pairs in one distance restraint: Conservative or Equal");
+ EETYPE("disre-weighting", ir->eDisreWeighting, edisreweighting_names);
+ CTYPE ("Use sqrt of the time averaged times the instantaneous violation");
+ EETYPE("disre-mixed", ir->bDisreMixed, yesno_names);
+ RTYPE ("disre-fc", ir->dr_fc, 1000.0);
+ RTYPE ("disre-tau", ir->dr_tau, 0.0);
+ CTYPE ("Output frequency for pair distances to energy file");
+ ITYPE ("nstdisreout", ir->nstdisreout, 100);
+ CTYPE ("Orientation restraints: No or Yes");
+ EETYPE("orire", opts->bOrire, yesno_names);
+ CTYPE ("Orientation restraints force constant and tau for time averaging");
+ RTYPE ("orire-fc", ir->orires_fc, 0.0);
+ RTYPE ("orire-tau", ir->orires_tau, 0.0);
+ STYPE ("orire-fitgrp",orirefitgrp, NULL);
+ CTYPE ("Output frequency for trace(SD) and S to energy file");
+ ITYPE ("nstorireout", ir->nstorireout, 100);
+ CTYPE ("Dihedral angle restraints: No or Yes");
+ EETYPE("dihre", opts->bDihre, yesno_names);
+ RTYPE ("dihre-fc", ir->dihre_fc, 1000.0);
+
+ /* Free energy stuff */
+ CCTYPE ("Free energy control stuff");
+ EETYPE("free-energy", ir->efep, efep_names);
+ RTYPE ("init-lambda", ir->init_lambda,0.0);
+ RTYPE ("delta-lambda",ir->delta_lambda,0.0);
+ STYPE ("foreign_lambda", foreign_lambda, NULL);
+ RTYPE ("sc-alpha",ir->sc_alpha,0.0);
+ ITYPE ("sc-power",ir->sc_power,0);
+ RTYPE ("sc-sigma",ir->sc_sigma,0.3);
+ ITYPE ("nstdhdl", ir->nstdhdl, 10);
+ EETYPE("separate-dhdl-file", ir->separate_dhdl_file,
+ separate_dhdl_file_names);
+ EETYPE("dhdl-derivatives", ir->dhdl_derivatives, dhdl_derivatives_names);
+ ITYPE ("dh_hist_size", ir->dh_hist_size, 0);
+ RTYPE ("dh_hist_spacing", ir->dh_hist_spacing, 0.1);
+ STYPE ("couple-moltype", couple_moltype, NULL);
+ EETYPE("couple-lambda0", opts->couple_lam0, couple_lam);
+ EETYPE("couple-lambda1", opts->couple_lam1, couple_lam);
+ EETYPE("couple-intramol", opts->bCoupleIntra, yesno_names);
+
+ /* Non-equilibrium MD stuff */
+ CCTYPE("Non-equilibrium MD stuff");
+ STYPE ("acc-grps", accgrps, NULL);
+ STYPE ("accelerate", acc, NULL);
+ STYPE ("freezegrps", freeze, NULL);
+ STYPE ("freezedim", frdim, NULL);
+ RTYPE ("cos-acceleration", ir->cos_accel, 0);
+ STYPE ("deform", deform, NULL);
+
+ /* Electric fields */
+ CCTYPE("Electric fields");
+ CTYPE ("Format is number of terms (int) and for all terms an amplitude (real)");
+ CTYPE ("and a phase angle (real)");
+ STYPE ("E-x", efield_x, NULL);
+ STYPE ("E-xt", efield_xt, NULL);
+ STYPE ("E-y", efield_y, NULL);
+ STYPE ("E-yt", efield_yt, NULL);
+ STYPE ("E-z", efield_z, NULL);
+ STYPE ("E-zt", efield_zt, NULL);
+
+ /* User defined thingies */
+ CCTYPE ("User defined thingies");
+ STYPE ("user1-grps", user1, NULL);
+ STYPE ("user2-grps", user2, NULL);
+ ITYPE ("userint1", ir->userint1, 0);
+ ITYPE ("userint2", ir->userint2, 0);
+ ITYPE ("userint3", ir->userint3, 0);
+ ITYPE ("userint4", ir->userint4, 0);
+ RTYPE ("userreal1", ir->userreal1, 0);
+ RTYPE ("userreal2", ir->userreal2, 0);
+ RTYPE ("userreal3", ir->userreal3, 0);
+ RTYPE ("userreal4", ir->userreal4, 0);
+#undef CTYPE
+
+ write_inpfile(mdparout,ninp,inp,FALSE,wi);
+ for (i=0; (i<ninp); i++) {
+ sfree(inp[i].name);
+ sfree(inp[i].value);
+ }
+ sfree(inp);
+
+ /* Process options if necessary */
+ for(m=0; m<2; m++) {
+ for(i=0; i<2*DIM; i++)
+ dumdub[m][i]=0.0;
+ if(ir->epc) {
+ switch (ir->epct) {
+ case epctISOTROPIC:
+ if (sscanf(dumstr[m],"%lf",&(dumdub[m][XX]))!=1) {
+ warning_error(wi,"Pressure coupling not enough values (I need 1)");
+ }
+ dumdub[m][YY]=dumdub[m][ZZ]=dumdub[m][XX];
+ break;
+ case epctSEMIISOTROPIC:
+ case epctSURFACETENSION:
+ if (sscanf(dumstr[m],"%lf%lf",
+ &(dumdub[m][XX]),&(dumdub[m][ZZ]))!=2) {
+ warning_error(wi,"Pressure coupling not enough values (I need 2)");
+ }
+ dumdub[m][YY]=dumdub[m][XX];
+ break;
+ case epctANISOTROPIC:
+ if (sscanf(dumstr[m],"%lf%lf%lf%lf%lf%lf",
+ &(dumdub[m][XX]),&(dumdub[m][YY]),&(dumdub[m][ZZ]),
+ &(dumdub[m][3]),&(dumdub[m][4]),&(dumdub[m][5]))!=6) {
+ warning_error(wi,"Pressure coupling not enough values (I need 6)");
+ }
+ break;
+ default:
+ gmx_fatal(FARGS,"Pressure coupling type %s not implemented yet",
+ epcoupltype_names[ir->epct]);
+ }
+ }
+ }
+ clear_mat(ir->ref_p);
+ clear_mat(ir->compress);
+ for(i=0; i<DIM; i++) {
+ ir->ref_p[i][i] = dumdub[1][i];
+ ir->compress[i][i] = dumdub[0][i];
+ }
+ if (ir->epct == epctANISOTROPIC) {
+ ir->ref_p[XX][YY] = dumdub[1][3];
+ ir->ref_p[XX][ZZ] = dumdub[1][4];
+ ir->ref_p[YY][ZZ] = dumdub[1][5];
+ if (ir->ref_p[XX][YY]!=0 && ir->ref_p[XX][ZZ]!=0 && ir->ref_p[YY][ZZ]!=0) {
+ warning(wi,"All off-diagonal reference pressures are non-zero. Are you sure you want to apply a threefold shear stress?\n");
+ }
+ ir->compress[XX][YY] = dumdub[0][3];
+ ir->compress[XX][ZZ] = dumdub[0][4];
+ ir->compress[YY][ZZ] = dumdub[0][5];
+ for(i=0; i<DIM; i++) {
+ for(m=0; m<i; m++) {
+ ir->ref_p[i][m] = ir->ref_p[m][i];
+ ir->compress[i][m] = ir->compress[m][i];
+ }
+ }
+ }
+
+ if (ir->comm_mode == ecmNO)
+ ir->nstcomm = 0;
+
+ opts->couple_moltype = NULL;
+ if (strlen(couple_moltype) > 0) {
+ if (ir->efep != efepNO) {
+ opts->couple_moltype = strdup(couple_moltype);
+ if (opts->couple_lam0 == opts->couple_lam1)
+ warning(wi,"The lambda=0 and lambda=1 states for coupling are identical");
+ if (ir->eI == eiMD && (opts->couple_lam0 == ecouplamNONE ||
+ opts->couple_lam1 == ecouplamNONE)) {
+ warning(wi,"For proper sampling of the (nearly) decoupled state, stochastic dynamics should be used");
+ }
+ } else {
+ warning(wi,"Can not couple a molecule with free_energy = no");
+ }
+ }
+
+ do_wall_params(ir,wall_atomtype,wall_density,opts);
+
+ if (opts->bOrire && str_nelem(orirefitgrp,MAXPTR,NULL)!=1) {
+ warning_error(wi,"ERROR: Need one orientation restraint fit group\n");
+ }
+
+ clear_mat(ir->deform);
+ for(i=0; i<6; i++)
+ dumdub[0][i] = 0;
+ m = sscanf(deform,"%lf %lf %lf %lf %lf %lf",
+ &(dumdub[0][0]),&(dumdub[0][1]),&(dumdub[0][2]),
+ &(dumdub[0][3]),&(dumdub[0][4]),&(dumdub[0][5]));
+ for(i=0; i<3; i++)
+ ir->deform[i][i] = dumdub[0][i];
+ ir->deform[YY][XX] = dumdub[0][3];
+ ir->deform[ZZ][XX] = dumdub[0][4];
+ ir->deform[ZZ][YY] = dumdub[0][5];
+ if (ir->epc != epcNO) {
+ for(i=0; i<3; i++)
+ for(j=0; j<=i; j++)
+ if (ir->deform[i][j]!=0 && ir->compress[i][j]!=0) {
+ warning_error(wi,"A box element has deform set and compressibility > 0");
+ }
+ for(i=0; i<3; i++)
+ for(j=0; j<i; j++)
+ if (ir->deform[i][j]!=0) {
+ for(m=j; m<DIM; m++)
+ if (ir->compress[m][j]!=0) {
+ sprintf(warn_buf,"An off-diagonal box element has deform set while compressibility > 0 for the same component of another box vector, this might lead to spurious periodicity effects.");
+ warning(wi,warn_buf);
+ }
+ }
+ }
+
+ if (ir->efep != efepNO) {
+ parse_n_double(foreign_lambda,&ir->n_flambda,&ir->flambda);
+ if (ir->n_flambda > 0 && ir->rlist < max(ir->rvdw,ir->rcoulomb)) {
+ warning_note(wi,"For foreign lambda free energy differences it is assumed that the soft-core interactions have no effect beyond the neighborlist cut-off");
+ }
+ } else {
+ ir->n_flambda = 0;
+ }
+
+ sfree(dumstr[0]);
+ sfree(dumstr[1]);
+}
+
+static int search_QMstring(char *s,int ng,const char *gn[])
+{
+ /* same as normal search_string, but this one searches QM strings */
+ int i;
+
+ for(i=0; (i<ng); i++)
+ if (gmx_strcasecmp(s,gn[i]) == 0)
+ return i;
+
+ gmx_fatal(FARGS,"this QM method or basisset (%s) is not implemented\n!",s);
+
+ return -1;
+
+} /* search_QMstring */
+
+
+int search_string(char *s,int ng,char *gn[])
+{
+ int i;
+
+ for(i=0; (i<ng); i++)
+ if (gmx_strcasecmp(s,gn[i]) == 0)
+ return i;
+
+ gmx_fatal(FARGS,"Group %s not found in indexfile.\nMaybe you have non-default goups in your mdp file, while not using the '-n' option of grompp.\nIn that case use the '-n' option.\n",s);
+
+ return -1;
+}
+
+static gmx_bool do_numbering(int natoms,gmx_groups_t *groups,int ng,char *ptrs[],
+ t_blocka *block,char *gnames[],
+ int gtype,int restnm,
+ int grptp,gmx_bool bVerbose,
+ warninp_t wi)
+{
+ unsigned short *cbuf;
+ t_grps *grps=&(groups->grps[gtype]);
+ int i,j,gid,aj,ognr,ntot=0;
+ const char *title;
+ gmx_bool bRest;
+ char warn_buf[STRLEN];
+
+ if (debug)
+ {
+ fprintf(debug,"Starting numbering %d groups of type %d\n",ng,gtype);
+ }
+
+ title = gtypes[gtype];
+
+ snew(cbuf,natoms);
+ /* Mark all id's as not set */
+ for(i=0; (i<natoms); i++)
+ {
+ cbuf[i] = NOGID;
+ }
+
+ snew(grps->nm_ind,ng+1); /* +1 for possible rest group */
+ for(i=0; (i<ng); i++)
+ {
+ /* Lookup the group name in the block structure */
+ gid = search_string(ptrs[i],block->nr,gnames);
+ if ((grptp != egrptpONE) || (i == 0))
+ {
+ grps->nm_ind[grps->nr++]=gid;
+ }
+ if (debug)
+ {
+ fprintf(debug,"Found gid %d for group %s\n",gid,ptrs[i]);
+ }
+
+ /* Now go over the atoms in the group */
+ for(j=block->index[gid]; (j<block->index[gid+1]); j++)
+ {
+
+ aj=block->a[j];
+
+ /* Range checking */
+ if ((aj < 0) || (aj >= natoms))
+ {
+ gmx_fatal(FARGS,"Invalid atom number %d in indexfile",aj);
+ }
+ /* Lookup up the old group number */
+ ognr = cbuf[aj];
+ if (ognr != NOGID)
+ {
+ gmx_fatal(FARGS,"Atom %d in multiple %s groups (%d and %d)",
+ aj+1,title,ognr+1,i+1);
+ }
+ else
+ {
+ /* Store the group number in buffer */
+ if (grptp == egrptpONE)
+ {
+ cbuf[aj] = 0;
+ }
+ else
+ {
+ cbuf[aj] = i;
+ }
+ ntot++;
+ }
+ }
+ }
+
+ /* Now check whether we have done all atoms */
+ bRest = FALSE;
+ if (ntot != natoms)
+ {
+ if (grptp == egrptpALL)
+ {
+ gmx_fatal(FARGS,"%d atoms are not part of any of the %s groups",
+ natoms-ntot,title);
+ }
+ else if (grptp == egrptpPART)
+ {
+ sprintf(warn_buf,"%d atoms are not part of any of the %s groups",
+ natoms-ntot,title);
+ warning_note(wi,warn_buf);
+ }
+ /* Assign all atoms currently unassigned to a rest group */
+ for(j=0; (j<natoms); j++)
+ {
+ if (cbuf[j] == NOGID)
+ {
+ cbuf[j] = grps->nr;
+ bRest = TRUE;
+ }
+ }
+ if (grptp != egrptpPART)
+ {
+ if (bVerbose)
+ {
+ fprintf(stderr,
+ "Making dummy/rest group for %s containing %d elements\n",
+ title,natoms-ntot);
+ }
+ /* Add group name "rest" */
+ grps->nm_ind[grps->nr] = restnm;
+
+ /* Assign the rest name to all atoms not currently assigned to a group */
+ for(j=0; (j<natoms); j++)
+ {
+ if (cbuf[j] == NOGID)
+ {
+ cbuf[j] = grps->nr;
+ }
+ }
+ grps->nr++;
+ }
+ }
+
+ if (grps->nr == 1)
+ {
+ groups->ngrpnr[gtype] = 0;
+ groups->grpnr[gtype] = NULL;
+ }
+ else
+ {
+ groups->ngrpnr[gtype] = natoms;
+ snew(groups->grpnr[gtype],natoms);
+ for(j=0; (j<natoms); j++)
+ {
+ groups->grpnr[gtype][j] = cbuf[j];
+ }
+ }
+
+ sfree(cbuf);
+
+ return (bRest && grptp == egrptpPART);
+}
+
+static void calc_nrdf(gmx_mtop_t *mtop,t_inputrec *ir,char **gnames)
+{
+ t_grpopts *opts;
+ gmx_groups_t *groups;
+ t_pull *pull;
+ int natoms,ai,aj,i,j,d,g,imin,jmin,nc;
+ t_iatom *ia;
+ int *nrdf2,*na_vcm,na_tot;
+ double *nrdf_tc,*nrdf_vcm,nrdf_uc,n_sub=0;
+ gmx_mtop_atomloop_all_t aloop;
+ t_atom *atom;
+ int mb,mol,ftype,as;
+ gmx_molblock_t *molb;
+ gmx_moltype_t *molt;
+
+ /* Calculate nrdf.
+ * First calc 3xnr-atoms for each group
+ * then subtract half a degree of freedom for each constraint
+ *
+ * Only atoms and nuclei contribute to the degrees of freedom...
+ */
+
+ opts = &ir->opts;
+
+ groups = &mtop->groups;
+ natoms = mtop->natoms;
+
+ /* Allocate one more for a possible rest group */
+ /* We need to sum degrees of freedom into doubles,
+ * since floats give too low nrdf's above 3 million atoms.
+ */
+ snew(nrdf_tc,groups->grps[egcTC].nr+1);
+ snew(nrdf_vcm,groups->grps[egcVCM].nr+1);
+ snew(na_vcm,groups->grps[egcVCM].nr+1);
+
+ for(i=0; i<groups->grps[egcTC].nr; i++)
+ nrdf_tc[i] = 0;
+ for(i=0; i<groups->grps[egcVCM].nr+1; i++)
+ nrdf_vcm[i] = 0;
+
+ snew(nrdf2,natoms);
+ aloop = gmx_mtop_atomloop_all_init(mtop);
+ while (gmx_mtop_atomloop_all_next(aloop,&i,&atom)) {
+ nrdf2[i] = 0;
+ if (atom->ptype == eptAtom || atom->ptype == eptNucleus) {
+ g = ggrpnr(groups,egcFREEZE,i);
+ /* Double count nrdf for particle i */
+ for(d=0; d<DIM; d++) {
+ if (opts->nFreeze[g][d] == 0) {
+ nrdf2[i] += 2;
+ }
+ }
+ nrdf_tc [ggrpnr(groups,egcTC ,i)] += 0.5*nrdf2[i];
+ nrdf_vcm[ggrpnr(groups,egcVCM,i)] += 0.5*nrdf2[i];
+ }
+ }
+
+ as = 0;
+ for(mb=0; mb<mtop->nmolblock; mb++) {
+ molb = &mtop->molblock[mb];
+ molt = &mtop->moltype[molb->type];
+ atom = molt->atoms.atom;
+ for(mol=0; mol<molb->nmol; mol++) {
+ for (ftype=F_CONSTR; ftype<=F_CONSTRNC; ftype++) {
+ ia = molt->ilist[ftype].iatoms;
+ for(i=0; i<molt->ilist[ftype].nr; ) {
+ /* Subtract degrees of freedom for the constraints,
+ * if the particles still have degrees of freedom left.
+ * If one of the particles is a vsite or a shell, then all
+ * constraint motion will go there, but since they do not
+ * contribute to the constraints the degrees of freedom do not
+ * change.
+ */
+ ai = as + ia[1];
+ aj = as + ia[2];
+ if (((atom[ia[1]].ptype == eptNucleus) ||
+ (atom[ia[1]].ptype == eptAtom)) &&
+ ((atom[ia[2]].ptype == eptNucleus) ||
+ (atom[ia[2]].ptype == eptAtom))) {
+ if (nrdf2[ai] > 0)
+ jmin = 1;
+ else
+ jmin = 2;
+ if (nrdf2[aj] > 0)
+ imin = 1;
+ else
+ imin = 2;
+ imin = min(imin,nrdf2[ai]);
+ jmin = min(jmin,nrdf2[aj]);
+ nrdf2[ai] -= imin;
+ nrdf2[aj] -= jmin;
+ nrdf_tc [ggrpnr(groups,egcTC ,ai)] -= 0.5*imin;
+ nrdf_tc [ggrpnr(groups,egcTC ,aj)] -= 0.5*jmin;
+ nrdf_vcm[ggrpnr(groups,egcVCM,ai)] -= 0.5*imin;
+ nrdf_vcm[ggrpnr(groups,egcVCM,aj)] -= 0.5*jmin;
+ }
+ ia += interaction_function[ftype].nratoms+1;
+ i += interaction_function[ftype].nratoms+1;
+ }
+ }
+ ia = molt->ilist[F_SETTLE].iatoms;
+ for(i=0; i<molt->ilist[F_SETTLE].nr; ) {
+ /* Subtract 1 dof from every atom in the SETTLE */
+ for(ai=as+ia[1]; ai<as+ia[1]+3; ai++) {
+ imin = min(2,nrdf2[ai]);
+ nrdf2[ai] -= imin;
+ nrdf_tc [ggrpnr(groups,egcTC ,ai)] -= 0.5*imin;
+ nrdf_vcm[ggrpnr(groups,egcVCM,ai)] -= 0.5*imin;
+ }
+ ia += 2;
+ i += 2;
+ }
+ as += molt->atoms.nr;
+ }
+ }
+
+ if (ir->ePull == epullCONSTRAINT) {
+ /* Correct nrdf for the COM constraints.
+ * We correct using the TC and VCM group of the first atom
+ * in the reference and pull group. If atoms in one pull group
+ * belong to different TC or VCM groups it is anyhow difficult
+ * to determine the optimal nrdf assignment.
+ */
+ pull = ir->pull;
+ if (pull->eGeom == epullgPOS) {
+ nc = 0;
+ for(i=0; i<DIM; i++) {
+ if (pull->dim[i])
+ nc++;
+ }
+ } else {
+ nc = 1;
+ }
+ for(i=0; i<pull->ngrp; i++) {
+ imin = 2*nc;
+ if (pull->grp[0].nat > 0) {
+ /* Subtract 1/2 dof from the reference group */
+ ai = pull->grp[0].ind[0];
+ if (nrdf_tc[ggrpnr(groups,egcTC,ai)] > 1) {
+ nrdf_tc [ggrpnr(groups,egcTC ,ai)] -= 0.5;
+ nrdf_vcm[ggrpnr(groups,egcVCM,ai)] -= 0.5;
+ imin--;
+ }
+ }
+ /* Subtract 1/2 dof from the pulled group */
+ ai = pull->grp[1+i].ind[0];
+ nrdf_tc [ggrpnr(groups,egcTC ,ai)] -= 0.5*imin;
+ nrdf_vcm[ggrpnr(groups,egcVCM,ai)] -= 0.5*imin;
+ if (nrdf_tc[ggrpnr(groups,egcTC,ai)] < 0)
+ gmx_fatal(FARGS,"Center of mass pulling constraints caused the number of degrees of freedom for temperature coupling group %s to be negative",gnames[groups->grps[egcTC].nm_ind[ggrpnr(groups,egcTC,ai)]]);
+ }
+ }
+
+ if (ir->nstcomm != 0) {
+ /* Subtract 3 from the number of degrees of freedom in each vcm group
+ * when com translation is removed and 6 when rotation is removed
+ * as well.
+ */
+ switch (ir->comm_mode) {
+ case ecmLINEAR:
+ n_sub = ndof_com(ir);
+ break;
+ case ecmANGULAR:
+ n_sub = 6;
+ break;
+ default:
+ n_sub = 0;
+ gmx_incons("Checking comm_mode");
+ }
+
+ for(i=0; i<groups->grps[egcTC].nr; i++) {
+ /* Count the number of atoms of TC group i for every VCM group */
+ for(j=0; j<groups->grps[egcVCM].nr+1; j++)
+ na_vcm[j] = 0;
+ na_tot = 0;
+ for(ai=0; ai<natoms; ai++)
+ if (ggrpnr(groups,egcTC,ai) == i) {
+ na_vcm[ggrpnr(groups,egcVCM,ai)]++;
+ na_tot++;
+ }
+ /* Correct for VCM removal according to the fraction of each VCM
+ * group present in this TC group.
+ */
+ nrdf_uc = nrdf_tc[i];
+ if (debug) {
+ fprintf(debug,"T-group[%d] nrdf_uc = %g, n_sub = %g\n",
+ i,nrdf_uc,n_sub);
+ }
+ nrdf_tc[i] = 0;
+ for(j=0; j<groups->grps[egcVCM].nr+1; j++) {
+ if (nrdf_vcm[j] > n_sub) {
+ nrdf_tc[i] += nrdf_uc*((double)na_vcm[j]/(double)na_tot)*
+ (nrdf_vcm[j] - n_sub)/nrdf_vcm[j];
+ }
+ if (debug) {
+ fprintf(debug," nrdf_vcm[%d] = %g, nrdf = %g\n",
+ j,nrdf_vcm[j],nrdf_tc[i]);
+ }
+ }
+ }
+ }
+ for(i=0; (i<groups->grps[egcTC].nr); i++) {
+ opts->nrdf[i] = nrdf_tc[i];
+ if (opts->nrdf[i] < 0)
+ opts->nrdf[i] = 0;
+ fprintf(stderr,
+ "Number of degrees of freedom in T-Coupling group %s is %.2f\n",
+ gnames[groups->grps[egcTC].nm_ind[i]],opts->nrdf[i]);
+ }
+
+ sfree(nrdf2);
+ sfree(nrdf_tc);
+ sfree(nrdf_vcm);
+ sfree(na_vcm);
+}
+
+static void decode_cos(char *s,t_cosines *cosine,gmx_bool bTime)
+{
+ char *t;
+ char format[STRLEN],f1[STRLEN];
+ double a,phi;
+ int i;
+
+ t=strdup(s);
+ trim(t);
+
+ cosine->n=0;
+ cosine->a=NULL;
+ cosine->phi=NULL;
+ if (strlen(t)) {
+ sscanf(t,"%d",&(cosine->n));
+ if (cosine->n <= 0) {
+ cosine->n=0;
+ } else {
+ snew(cosine->a,cosine->n);
+ snew(cosine->phi,cosine->n);
+
+ sprintf(format,"%%*d");
+ for(i=0; (i<cosine->n); i++) {
+ strcpy(f1,format);
+ strcat(f1,"%lf%lf");
+ if (sscanf(t,f1,&a,&phi) < 2)
+ gmx_fatal(FARGS,"Invalid input for electric field shift: '%s'",t);
+ cosine->a[i]=a;
+ cosine->phi[i]=phi;
+ strcat(format,"%*lf%*lf");
+ }
+ }
+ }
+ sfree(t);
+}
+
+static gmx_bool do_egp_flag(t_inputrec *ir,gmx_groups_t *groups,
+ const char *option,const char *val,int flag)
+{
+ /* The maximum number of energy group pairs would be MAXPTR*(MAXPTR+1)/2.
+ * But since this is much larger than STRLEN, such a line can not be parsed.
+ * The real maximum is the number of names that fit in a string: STRLEN/2.
+ */
+#define EGP_MAX (STRLEN/2)
+ int nelem,i,j,k,nr;
+ char *names[EGP_MAX];
+ char ***gnames;
+ gmx_bool bSet;
+
+ gnames = groups->grpname;
+
+ nelem = str_nelem(val,EGP_MAX,names);
+ if (nelem % 2 != 0)
+ gmx_fatal(FARGS,"The number of groups for %s is odd",option);
+ nr = groups->grps[egcENER].nr;
+ bSet = FALSE;
+ for(i=0; i<nelem/2; i++) {
+ j = 0;
+ while ((j < nr) &&
+ gmx_strcasecmp(names[2*i],*(gnames[groups->grps[egcENER].nm_ind[j]])))
+ j++;
+ if (j == nr)
+ gmx_fatal(FARGS,"%s in %s is not an energy group\n",
+ names[2*i],option);
+ k = 0;
+ while ((k < nr) &&
+ gmx_strcasecmp(names[2*i+1],*(gnames[groups->grps[egcENER].nm_ind[k]])))
+ k++;
+ if (k==nr)
+ gmx_fatal(FARGS,"%s in %s is not an energy group\n",
+ names[2*i+1],option);
+ if ((j < nr) && (k < nr)) {
+ ir->opts.egp_flags[nr*j+k] |= flag;
+ ir->opts.egp_flags[nr*k+j] |= flag;
+ bSet = TRUE;
+ }
+ }
+
+ return bSet;
+}
+
+void do_index(const char* mdparin, const char *ndx,
+ gmx_mtop_t *mtop,
+ gmx_bool bVerbose,
+ t_inputrec *ir,rvec *v,
+ warninp_t wi)
+{
+ t_blocka *grps;
+ gmx_groups_t *groups;
+ int natoms;
+ t_symtab *symtab;
+ t_atoms atoms_all;
+ char warnbuf[STRLEN],**gnames;
+ int nr,ntcg,ntau_t,nref_t,nacc,nofg,nSA,nSA_points,nSA_time,nSA_temp;
+ real tau_min;
+ int nstcmin;
+ int nacg,nfreeze,nfrdim,nenergy,nvcm,nuser;
+ char *ptr1[MAXPTR],*ptr2[MAXPTR],*ptr3[MAXPTR];
+ int i,j,k,restnm;
+ real SAtime;
+ gmx_bool bExcl,bTable,bSetTCpar,bAnneal,bRest;
+ int nQMmethod,nQMbasis,nQMcharge,nQMmult,nbSH,nCASorb,nCASelec,
+ nSAon,nSAoff,nSAsteps,nQMg,nbOPT,nbTS;
+ char warn_buf[STRLEN];
+
+ if (bVerbose)
+ fprintf(stderr,"processing index file...\n");
+ debug_gmx();
+ if (ndx == NULL) {
+ snew(grps,1);
+ snew(grps->index,1);
+ snew(gnames,1);
+ atoms_all = gmx_mtop_global_atoms(mtop);
+ analyse(&atoms_all,grps,&gnames,FALSE,TRUE);
+ free_t_atoms(&atoms_all,FALSE);
+ } else {
+ grps = init_index(ndx,&gnames);
+ }
+
+ groups = &mtop->groups;
+ natoms = mtop->natoms;
+ symtab = &mtop->symtab;
+
+ snew(groups->grpname,grps->nr+1);
+
+ for(i=0; (i<grps->nr); i++) {
+ groups->grpname[i] = put_symtab(symtab,gnames[i]);
+ }
+ groups->grpname[i] = put_symtab(symtab,"rest");
+ restnm=i;
+ srenew(gnames,grps->nr+1);
+ gnames[restnm] = *(groups->grpname[i]);
+ groups->ngrpname = grps->nr+1;
+
+ set_warning_line(wi,mdparin,-1);
+
+ ntau_t = str_nelem(tau_t,MAXPTR,ptr1);
+ nref_t = str_nelem(ref_t,MAXPTR,ptr2);
+ ntcg = str_nelem(tcgrps,MAXPTR,ptr3);
+ if ((ntau_t != ntcg) || (nref_t != ntcg)) {
+ gmx_fatal(FARGS,"Invalid T coupling input: %d groups, %d ref_t values and "
+ "%d tau_t values",ntcg,nref_t,ntau_t);
+ }
+
+ bSetTCpar = (ir->etc || EI_SD(ir->eI) || ir->eI==eiBD || EI_TPI(ir->eI));
+ do_numbering(natoms,groups,ntcg,ptr3,grps,gnames,egcTC,
+ restnm,bSetTCpar ? egrptpALL : egrptpALL_GENREST,bVerbose,wi);
+ nr = groups->grps[egcTC].nr;
+ ir->opts.ngtc = nr;
+ snew(ir->opts.nrdf,nr);
+ snew(ir->opts.tau_t,nr);
+ snew(ir->opts.ref_t,nr);
+ if (ir->eI==eiBD && ir->bd_fric==0) {
+ fprintf(stderr,"bd_fric=0, so tau_t will be used as the inverse friction constant(s)\n");
+ }
+
+ if (bSetTCpar)
+ {
+ if (nr != nref_t)
+ {
+ gmx_fatal(FARGS,"Not enough ref_t and tau_t values!");
+ }
+
+ tau_min = 1e20;
+ for(i=0; (i<nr); i++)
+ {
+ ir->opts.tau_t[i] = strtod(ptr1[i],NULL);
+ if ((ir->eI == eiBD || ir->eI == eiSD2) && ir->opts.tau_t[i] <= 0)
+ {
+ sprintf(warn_buf,"With integrator %s tau_t should be larger than 0",ei_names[ir->eI]);
+ warning_error(wi,warn_buf);
+ }
+ if ((ir->etc == etcVRESCALE && ir->opts.tau_t[i] >= 0) ||
+ (ir->etc != etcVRESCALE && ir->opts.tau_t[i] > 0))
+ {
+ tau_min = min(tau_min,ir->opts.tau_t[i]);
+ }
+ }
+ if (ir->etc != etcNO && ir->nsttcouple == -1)
+ {
+ ir->nsttcouple = ir_optimal_nsttcouple(ir);
+ }
+ if (EI_VV(ir->eI))
+ {
+ if ((ir->epc==epcMTTK) && (ir->etc>etcNO))
+ {
+ int mincouple;
+ mincouple = ir->nsttcouple;
+ if (ir->nstpcouple < mincouple)
+ {
+ mincouple = ir->nstpcouple;
+ }
+ ir->nstpcouple = mincouple;
+ ir->nsttcouple = mincouple;
+ warning_note(wi,"for current Trotter decomposition methods with vv, nsttcouple and nstpcouple must be equal. Both have been reset to min(nsttcouple,nstpcouple)");
+ }
+ }
+ nstcmin = tcouple_min_integration_steps(ir->etc);
+ if (nstcmin > 1)
+ {
+ if (tau_min/(ir->delta_t*ir->nsttcouple) < nstcmin)
+ {
+ sprintf(warn_buf,"For proper integration of the %s thermostat, tau_t (%g) should be at least %d times larger than nsttcouple*dt (%g)",
+ ETCOUPLTYPE(ir->etc),
+ tau_min,nstcmin,
+ ir->nsttcouple*ir->delta_t);
+ warning(wi,warn_buf);
+ }
+ }
+ for(i=0; (i<nr); i++)
+ {
+ ir->opts.ref_t[i] = strtod(ptr2[i],NULL);
+ if (ir->opts.ref_t[i] < 0)
+ {
+ gmx_fatal(FARGS,"ref_t for group %d negative",i);
+ }
+ }
+ }
+
+ /* Simulated annealing for each group. There are nr groups */
+ nSA = str_nelem(anneal,MAXPTR,ptr1);
+ if (nSA == 1 && (ptr1[0][0]=='n' || ptr1[0][0]=='N'))
+ nSA = 0;
+ if(nSA>0 && nSA != nr)
+ gmx_fatal(FARGS,"Not enough annealing values: %d (for %d groups)\n",nSA,nr);
+ else {
+ snew(ir->opts.annealing,nr);
+ snew(ir->opts.anneal_npoints,nr);
+ snew(ir->opts.anneal_time,nr);
+ snew(ir->opts.anneal_temp,nr);
+ for(i=0;i<nr;i++) {
+ ir->opts.annealing[i]=eannNO;
+ ir->opts.anneal_npoints[i]=0;
+ ir->opts.anneal_time[i]=NULL;
+ ir->opts.anneal_temp[i]=NULL;
+ }
+ if (nSA > 0) {
+ bAnneal=FALSE;
+ for(i=0;i<nr;i++) {
+ if(ptr1[i][0]=='n' || ptr1[i][0]=='N') {
+ ir->opts.annealing[i]=eannNO;
+ } else if(ptr1[i][0]=='s'|| ptr1[i][0]=='S') {
+ ir->opts.annealing[i]=eannSINGLE;
+ bAnneal=TRUE;
+ } else if(ptr1[i][0]=='p'|| ptr1[i][0]=='P') {
+ ir->opts.annealing[i]=eannPERIODIC;
+ bAnneal=TRUE;
+ }
+ }
+ if(bAnneal) {
+ /* Read the other fields too */
+ nSA_points = str_nelem(anneal_npoints,MAXPTR,ptr1);
+ if(nSA_points!=nSA)
+ gmx_fatal(FARGS,"Found %d annealing_npoints values for %d groups\n",nSA_points,nSA);
+ for(k=0,i=0;i<nr;i++) {
+ ir->opts.anneal_npoints[i]=strtol(ptr1[i],NULL,10);
+ if(ir->opts.anneal_npoints[i]==1)
+ gmx_fatal(FARGS,"Please specify at least a start and an end point for annealing\n");
+ snew(ir->opts.anneal_time[i],ir->opts.anneal_npoints[i]);
+ snew(ir->opts.anneal_temp[i],ir->opts.anneal_npoints[i]);
+ k += ir->opts.anneal_npoints[i];
+ }
+
+ nSA_time = str_nelem(anneal_time,MAXPTR,ptr1);
+ if(nSA_time!=k)
+ gmx_fatal(FARGS,"Found %d annealing_time values, wanter %d\n",nSA_time,k);
+ nSA_temp = str_nelem(anneal_temp,MAXPTR,ptr2);
+ if(nSA_temp!=k)
+ gmx_fatal(FARGS,"Found %d annealing_temp values, wanted %d\n",nSA_temp,k);
+
+ for(i=0,k=0;i<nr;i++) {
+
+ for(j=0;j<ir->opts.anneal_npoints[i];j++) {
+ ir->opts.anneal_time[i][j]=strtod(ptr1[k],NULL);
+ ir->opts.anneal_temp[i][j]=strtod(ptr2[k],NULL);
+ if(j==0) {
+ if(ir->opts.anneal_time[i][0] > (ir->init_t+GMX_REAL_EPS))
+ gmx_fatal(FARGS,"First time point for annealing > init_t.\n");
+ } else {
+ /* j>0 */
+ if(ir->opts.anneal_time[i][j]<ir->opts.anneal_time[i][j-1])
+ gmx_fatal(FARGS,"Annealing timepoints out of order: t=%f comes after t=%f\n",
+ ir->opts.anneal_time[i][j],ir->opts.anneal_time[i][j-1]);
+ }
+ if(ir->opts.anneal_temp[i][j]<0)
+ gmx_fatal(FARGS,"Found negative temperature in annealing: %f\n",ir->opts.anneal_temp[i][j]);
+ k++;
+ }
+ }
+ /* Print out some summary information, to make sure we got it right */
+ for(i=0,k=0;i<nr;i++) {
+ if(ir->opts.annealing[i]!=eannNO) {
+ j = groups->grps[egcTC].nm_ind[i];
+ fprintf(stderr,"Simulated annealing for group %s: %s, %d timepoints\n",
+ *(groups->grpname[j]),eann_names[ir->opts.annealing[i]],
+ ir->opts.anneal_npoints[i]);
+ fprintf(stderr,"Time (ps) Temperature (K)\n");
+ /* All terms except the last one */
+ for(j=0;j<(ir->opts.anneal_npoints[i]-1);j++)
+ fprintf(stderr,"%9.1f %5.1f\n",ir->opts.anneal_time[i][j],ir->opts.anneal_temp[i][j]);
+
+ /* Finally the last one */
+ j = ir->opts.anneal_npoints[i]-1;
+ if(ir->opts.annealing[i]==eannSINGLE)
+ fprintf(stderr,"%9.1f- %5.1f\n",ir->opts.anneal_time[i][j],ir->opts.anneal_temp[i][j]);
+ else {
+ fprintf(stderr,"%9.1f %5.1f\n",ir->opts.anneal_time[i][j],ir->opts.anneal_temp[i][j]);
+ if(fabs(ir->opts.anneal_temp[i][j]-ir->opts.anneal_temp[i][0])>GMX_REAL_EPS)
+ warning_note(wi,"There is a temperature jump when your annealing loops back.\n");
+ }
+ }
+ }
+ }
+ }
+ }
+
+ if (ir->ePull != epullNO) {
+ make_pull_groups(ir->pull,pull_grp,grps,gnames);
+ }
+
+ if (ir->bRot) {
+ make_rotation_groups(ir->rot,rot_grp,grps,gnames);
+ }
+
+ nacc = str_nelem(acc,MAXPTR,ptr1);
+ nacg = str_nelem(accgrps,MAXPTR,ptr2);
+ if (nacg*DIM != nacc)
+ gmx_fatal(FARGS,"Invalid Acceleration input: %d groups and %d acc. values",
+ nacg,nacc);
+ do_numbering(natoms,groups,nacg,ptr2,grps,gnames,egcACC,
+ restnm,egrptpALL_GENREST,bVerbose,wi);
+ nr = groups->grps[egcACC].nr;
+ snew(ir->opts.acc,nr);
+ ir->opts.ngacc=nr;
+
+ for(i=k=0; (i<nacg); i++)
+ for(j=0; (j<DIM); j++,k++)
+ ir->opts.acc[i][j]=strtod(ptr1[k],NULL);
+ for( ;(i<nr); i++)
+ for(j=0; (j<DIM); j++)
+ ir->opts.acc[i][j]=0;
+
+ nfrdim = str_nelem(frdim,MAXPTR,ptr1);
+ nfreeze = str_nelem(freeze,MAXPTR,ptr2);
+ if (nfrdim != DIM*nfreeze)
+ gmx_fatal(FARGS,"Invalid Freezing input: %d groups and %d freeze values",
+ nfreeze,nfrdim);
+ do_numbering(natoms,groups,nfreeze,ptr2,grps,gnames,egcFREEZE,
+ restnm,egrptpALL_GENREST,bVerbose,wi);
+ nr = groups->grps[egcFREEZE].nr;
+ ir->opts.ngfrz=nr;
+ snew(ir->opts.nFreeze,nr);
+ for(i=k=0; (i<nfreeze); i++)
+ for(j=0; (j<DIM); j++,k++) {
+ ir->opts.nFreeze[i][j]=(gmx_strncasecmp(ptr1[k],"Y",1)==0);
+ if (!ir->opts.nFreeze[i][j]) {
+ if (gmx_strncasecmp(ptr1[k],"N",1) != 0) {
+ sprintf(warnbuf,"Please use Y(ES) or N(O) for freezedim only "
+ "(not %s)", ptr1[k]);
+ warning(wi,warn_buf);
+ }
+ }
+ }
+ for( ; (i<nr); i++)
+ for(j=0; (j<DIM); j++)
+ ir->opts.nFreeze[i][j]=0;
+
+ nenergy=str_nelem(energy,MAXPTR,ptr1);
+ do_numbering(natoms,groups,nenergy,ptr1,grps,gnames,egcENER,
+ restnm,egrptpALL_GENREST,bVerbose,wi);
+ add_wall_energrps(groups,ir->nwall,symtab);
+ ir->opts.ngener = groups->grps[egcENER].nr;
+ nvcm=str_nelem(vcm,MAXPTR,ptr1);
+ bRest =
+ do_numbering(natoms,groups,nvcm,ptr1,grps,gnames,egcVCM,
+ restnm,nvcm==0 ? egrptpALL_GENREST : egrptpPART,bVerbose,wi);
+ if (bRest) {
+ warning(wi,"Some atoms are not part of any center of mass motion removal group.\n"
+ "This may lead to artifacts.\n"
+ "In most cases one should use one group for the whole system.");
+ }
+
+ /* Now we have filled the freeze struct, so we can calculate NRDF */
+ calc_nrdf(mtop,ir,gnames);
+
+ if (v && NULL) {
+ real fac,ntot=0;
+
+ /* Must check per group! */
+ for(i=0; (i<ir->opts.ngtc); i++)
+ ntot += ir->opts.nrdf[i];
+ if (ntot != (DIM*natoms)) {
+ fac = sqrt(ntot/(DIM*natoms));
+ if (bVerbose)
+ fprintf(stderr,"Scaling velocities by a factor of %.3f to account for constraints\n"
+ "and removal of center of mass motion\n",fac);
+ for(i=0; (i<natoms); i++)
+ svmul(fac,v[i],v[i]);
+ }
+ }
+
+ nuser=str_nelem(user1,MAXPTR,ptr1);
+ do_numbering(natoms,groups,nuser,ptr1,grps,gnames,egcUser1,
+ restnm,egrptpALL_GENREST,bVerbose,wi);
+ nuser=str_nelem(user2,MAXPTR,ptr1);
+ do_numbering(natoms,groups,nuser,ptr1,grps,gnames,egcUser2,
+ restnm,egrptpALL_GENREST,bVerbose,wi);
+ nuser=str_nelem(xtc_grps,MAXPTR,ptr1);
+ do_numbering(natoms,groups,nuser,ptr1,grps,gnames,egcXTC,
+ restnm,egrptpONE,bVerbose,wi);
+ nofg = str_nelem(orirefitgrp,MAXPTR,ptr1);
+ do_numbering(natoms,groups,nofg,ptr1,grps,gnames,egcORFIT,
+ restnm,egrptpALL_GENREST,bVerbose,wi);
+
+ /* QMMM input processing */
+ nQMg = str_nelem(QMMM,MAXPTR,ptr1);
+ nQMmethod = str_nelem(QMmethod,MAXPTR,ptr2);
+ nQMbasis = str_nelem(QMbasis,MAXPTR,ptr3);
+ if((nQMmethod != nQMg)||(nQMbasis != nQMg)){
+ gmx_fatal(FARGS,"Invalid QMMM input: %d groups %d basissets"
+ " and %d methods\n",nQMg,nQMbasis,nQMmethod);
+ }
+ /* group rest, if any, is always MM! */
+ do_numbering(natoms,groups,nQMg,ptr1,grps,gnames,egcQMMM,
+ restnm,egrptpALL_GENREST,bVerbose,wi);
+ nr = nQMg; /*atoms->grps[egcQMMM].nr;*/
+ ir->opts.ngQM = nQMg;
+ snew(ir->opts.QMmethod,nr);
+ snew(ir->opts.QMbasis,nr);
+ for(i=0;i<nr;i++){
+ /* input consists of strings: RHF CASSCF PM3 .. These need to be
+ * converted to the corresponding enum in names.c
+ */
+ ir->opts.QMmethod[i] = search_QMstring(ptr2[i],eQMmethodNR,
+ eQMmethod_names);
+ ir->opts.QMbasis[i] = search_QMstring(ptr3[i],eQMbasisNR,
+ eQMbasis_names);
+
+ }
+ nQMmult = str_nelem(QMmult,MAXPTR,ptr1);
+ nQMcharge = str_nelem(QMcharge,MAXPTR,ptr2);
+ nbSH = str_nelem(bSH,MAXPTR,ptr3);
+ snew(ir->opts.QMmult,nr);
+ snew(ir->opts.QMcharge,nr);
+ snew(ir->opts.bSH,nr);
+
+ for(i=0;i<nr;i++){
+ ir->opts.QMmult[i] = strtol(ptr1[i],NULL,10);
+ ir->opts.QMcharge[i] = strtol(ptr2[i],NULL,10);
+ ir->opts.bSH[i] = (gmx_strncasecmp(ptr3[i],"Y",1)==0);
+ }
+
+ nCASelec = str_nelem(CASelectrons,MAXPTR,ptr1);
+ nCASorb = str_nelem(CASorbitals,MAXPTR,ptr2);
+ snew(ir->opts.CASelectrons,nr);
+ snew(ir->opts.CASorbitals,nr);
+ for(i=0;i<nr;i++){
+ ir->opts.CASelectrons[i]= strtol(ptr1[i],NULL,10);
+ ir->opts.CASorbitals[i] = strtol(ptr2[i],NULL,10);
+ }
+ /* special optimization options */
+
+ nbOPT = str_nelem(bOPT,MAXPTR,ptr1);
+ nbTS = str_nelem(bTS,MAXPTR,ptr2);
+ snew(ir->opts.bOPT,nr);
+ snew(ir->opts.bTS,nr);
+ for(i=0;i<nr;i++){
+ ir->opts.bOPT[i] = (gmx_strncasecmp(ptr1[i],"Y",1)==0);
+ ir->opts.bTS[i] = (gmx_strncasecmp(ptr2[i],"Y",1)==0);
+ }
+ nSAon = str_nelem(SAon,MAXPTR,ptr1);
+ nSAoff = str_nelem(SAoff,MAXPTR,ptr2);
+ nSAsteps = str_nelem(SAsteps,MAXPTR,ptr3);
+ snew(ir->opts.SAon,nr);
+ snew(ir->opts.SAoff,nr);
+ snew(ir->opts.SAsteps,nr);
+
+ for(i=0;i<nr;i++){
+ ir->opts.SAon[i] = strtod(ptr1[i],NULL);
+ ir->opts.SAoff[i] = strtod(ptr2[i],NULL);
+ ir->opts.SAsteps[i] = strtol(ptr3[i],NULL,10);
+ }
+ /* end of QMMM input */
+
+ if (bVerbose)
+ for(i=0; (i<egcNR); i++) {
+ fprintf(stderr,"%-16s has %d element(s):",gtypes[i],groups->grps[i].nr);
+ for(j=0; (j<groups->grps[i].nr); j++)
+ fprintf(stderr," %s",*(groups->grpname[groups->grps[i].nm_ind[j]]));
+ fprintf(stderr,"\n");
+ }
+
+ nr = groups->grps[egcENER].nr;
+ snew(ir->opts.egp_flags,nr*nr);
+
+ bExcl = do_egp_flag(ir,groups,"energygrp_excl",egpexcl,EGP_EXCL);
+ if (bExcl && EEL_FULL(ir->coulombtype))
+ warning(wi,"Can not exclude the lattice Coulomb energy between energy groups");
+
+ bTable = do_egp_flag(ir,groups,"energygrp_table",egptable,EGP_TABLE);
+ if (bTable && !(ir->vdwtype == evdwUSER) &&
+ !(ir->coulombtype == eelUSER) && !(ir->coulombtype == eelPMEUSER) &&
+ !(ir->coulombtype == eelPMEUSERSWITCH))
+ gmx_fatal(FARGS,"Can only have energy group pair tables in combination with user tables for VdW and/or Coulomb");
+
+ decode_cos(efield_x,&(ir->ex[XX]),FALSE);
+ decode_cos(efield_xt,&(ir->et[XX]),TRUE);
+ decode_cos(efield_y,&(ir->ex[YY]),FALSE);
+ decode_cos(efield_yt,&(ir->et[YY]),TRUE);
+ decode_cos(efield_z,&(ir->ex[ZZ]),FALSE);
+ decode_cos(efield_zt,&(ir->et[ZZ]),TRUE);
+
+ for(i=0; (i<grps->nr); i++)
+ sfree(gnames[i]);
+ sfree(gnames);
+ done_blocka(grps);
+ sfree(grps);
+
+}
+
+
+
+static void check_disre(gmx_mtop_t *mtop)
+{
+ gmx_ffparams_t *ffparams;
+ t_functype *functype;
+ t_iparams *ip;
+ int i,ndouble,ftype;
+ int label,old_label;
+
+ if (gmx_mtop_ftype_count(mtop,F_DISRES) > 0) {
+ ffparams = &mtop->ffparams;
+ functype = ffparams->functype;
+ ip = ffparams->iparams;
+ ndouble = 0;
+ old_label = -1;
+ for(i=0; i<ffparams->ntypes; i++) {
+ ftype = functype[i];
+ if (ftype == F_DISRES) {
+ label = ip[i].disres.label;
+ if (label == old_label) {
+ fprintf(stderr,"Distance restraint index %d occurs twice\n",label);
+ ndouble++;
+ }
+ old_label = label;
+ }
+ }
+ if (ndouble>0)
+ gmx_fatal(FARGS,"Found %d double distance restraint indices,\n"
+ "probably the parameters for multiple pairs in one restraint "
+ "are not identical\n",ndouble);
+ }
+}
+
+static gmx_bool absolute_reference(t_inputrec *ir,gmx_mtop_t *sys,ivec AbsRef)
+{
+ int d,g,i;
+ gmx_mtop_ilistloop_t iloop;
+ t_ilist *ilist;
+ int nmol;
+ t_iparams *pr;
+
+ /* Check the COM */
+ for(d=0; d<DIM; d++) {
+ AbsRef[d] = (d < ndof_com(ir) ? 0 : 1);
+ }
+ /* Check for freeze groups */
+ for(g=0; g<ir->opts.ngfrz; g++) {
+ for(d=0; d<DIM; d++) {
+ if (ir->opts.nFreeze[g][d] != 0) {
+ AbsRef[d] = 1;
+ }
+ }
+ }
+ /* Check for position restraints */
+ iloop = gmx_mtop_ilistloop_init(sys);
+ while (gmx_mtop_ilistloop_next(iloop,&ilist,&nmol)) {
+ if (nmol > 0) {
+ for(i=0; i<ilist[F_POSRES].nr; i+=2) {
+ pr = &sys->ffparams.iparams[ilist[F_POSRES].iatoms[i]];
+ for(d=0; d<DIM; d++) {
+ if (pr->posres.fcA[d] != 0) {
+ AbsRef[d] = 1;
+ }
+ }
+ }
+ }
+ }
+
+ return (AbsRef[XX] != 0 && AbsRef[YY] != 0 && AbsRef[ZZ] != 0);
+}
+
+void triple_check(const char *mdparin,t_inputrec *ir,gmx_mtop_t *sys,
+ warninp_t wi)
+{
+ char err_buf[256];
+ int i,m,g,nmol,npct;
+ gmx_bool bCharge,bAcc;
+ real gdt_max,*mgrp,mt;
+ rvec acc;
+ gmx_mtop_atomloop_block_t aloopb;
+ gmx_mtop_atomloop_all_t aloop;
+ t_atom *atom;
+ ivec AbsRef;
+ char warn_buf[STRLEN];
+
+ set_warning_line(wi,mdparin,-1);
+
+ if (EI_DYNAMICS(ir->eI) && !EI_SD(ir->eI) && ir->eI != eiBD &&
+ ir->comm_mode == ecmNO &&
+ !(absolute_reference(ir,sys,AbsRef) || ir->nsteps <= 10)) {
+ warning(wi,"You are not using center of mass motion removal (mdp option comm-mode), numerical rounding errors can lead to build up of kinetic energy of the center of mass");
+ }
+
+ bCharge = FALSE;
+ aloopb = gmx_mtop_atomloop_block_init(sys);
+ while (gmx_mtop_atomloop_block_next(aloopb,&atom,&nmol)) {
+ if (atom->q != 0 || atom->qB != 0) {
+ bCharge = TRUE;
+ }
+ }
+
+ if (!bCharge) {
+ if (EEL_FULL(ir->coulombtype)) {
+ sprintf(err_buf,
+ "You are using full electrostatics treatment %s for a system without charges.\n"
+ "This costs a lot of performance for just processing zeros, consider using %s instead.\n",
+ EELTYPE(ir->coulombtype),EELTYPE(eelCUT));
+ warning(wi,err_buf);
+ }
+ } else {
+ if (ir->coulombtype == eelCUT && ir->rcoulomb > 0 && !ir->implicit_solvent) {
+ sprintf(err_buf,
+ "You are using a plain Coulomb cut-off, which might produce artifacts.\n"
+ "You might want to consider using %s electrostatics.\n",
+ EELTYPE(eelPME));
+ warning_note(wi,err_buf);
+ }
+ }
+
+ /* Generalized reaction field */
+ if (ir->opts.ngtc == 0) {
+ sprintf(err_buf,"No temperature coupling while using coulombtype %s",
+ eel_names[eelGRF]);
+ CHECK(ir->coulombtype == eelGRF);
+ }
+ else {
+ sprintf(err_buf,"When using coulombtype = %s"
+ " ref_t for temperature coupling should be > 0",
+ eel_names[eelGRF]);
+ CHECK((ir->coulombtype == eelGRF) && (ir->opts.ref_t[0] <= 0));
+ }
+
+ if (ir->eI == eiSD1) {
+ gdt_max = 0;
+ for(i=0; (i<ir->opts.ngtc); i++)
+ gdt_max = max(gdt_max,ir->delta_t/ir->opts.tau_t[i]);
+ if (0.5*gdt_max > 0.0015) {
+ sprintf(warn_buf,"The relative error with integrator %s is 0.5*delta_t/tau_t = %g, you might want to switch to integrator %s\n",
+ ei_names[ir->eI],0.5*gdt_max,ei_names[eiSD2]);
+ warning_note(wi,warn_buf);
+ }
+ }
+
+ bAcc = FALSE;
+ for(i=0; (i<sys->groups.grps[egcACC].nr); i++) {
+ for(m=0; (m<DIM); m++) {
+ if (fabs(ir->opts.acc[i][m]) > 1e-6) {
+ bAcc = TRUE;
+ }
+ }
+ }
+ if (bAcc) {
+ clear_rvec(acc);
+ snew(mgrp,sys->groups.grps[egcACC].nr);
+ aloop = gmx_mtop_atomloop_all_init(sys);
+ while (gmx_mtop_atomloop_all_next(aloop,&i,&atom)) {
+ mgrp[ggrpnr(&sys->groups,egcACC,i)] += atom->m;
+ }
+ mt = 0.0;
+ for(i=0; (i<sys->groups.grps[egcACC].nr); i++) {
+ for(m=0; (m<DIM); m++)
+ acc[m] += ir->opts.acc[i][m]*mgrp[i];
+ mt += mgrp[i];
+ }
+ for(m=0; (m<DIM); m++) {
+ if (fabs(acc[m]) > 1e-6) {
+ const char *dim[DIM] = { "X", "Y", "Z" };
+ fprintf(stderr,
+ "Net Acceleration in %s direction, will %s be corrected\n",
+ dim[m],ir->nstcomm != 0 ? "" : "not");
+ if (ir->nstcomm != 0 && m < ndof_com(ir)) {
+ acc[m] /= mt;
+ for (i=0; (i<sys->groups.grps[egcACC].nr); i++)
+ ir->opts.acc[i][m] -= acc[m];
+ }
+ }
+ }
+ sfree(mgrp);
+ }
+
+ if (ir->efep != efepNO && ir->sc_alpha != 0 &&
+ !gmx_within_tol(sys->ffparams.reppow,12.0,10*GMX_DOUBLE_EPS)) {
+ gmx_fatal(FARGS,"Soft-core interactions are only supported with VdW repulsion power 12");
+ }
+
+ if (ir->ePull != epullNO) {
+ if (ir->pull->grp[0].nat == 0) {
+ absolute_reference(ir,sys,AbsRef);
+ for(m=0; m<DIM; m++) {
+ if (ir->pull->dim[m] && !AbsRef[m]) {
+ warning(wi,"You are using an absolute reference for pulling, but the rest of the system does not have an absolute reference. This will lead to artifacts.");
+ break;
+ }
+ }
+ }
+
+ if (ir->pull->eGeom == epullgDIRPBC) {
+ for(i=0; i<3; i++) {
+ for(m=0; m<=i; m++) {
+ if ((ir->epc != epcNO && ir->compress[i][m] != 0) ||
+ ir->deform[i][m] != 0) {
+ for(g=1; g<ir->pull->ngrp; g++) {
+ if (ir->pull->grp[g].vec[m] != 0) {
+ gmx_fatal(FARGS,"Can not have dynamic box while using pull geometry '%s' (dim %c)",EPULLGEOM(ir->pull->eGeom),'x'+m);
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+
+ check_disre(sys);
+}
+
+void double_check(t_inputrec *ir,matrix box,gmx_bool bConstr,warninp_t wi)
+{
+ real min_size;
+ gmx_bool bTWIN;
+ char warn_buf[STRLEN];
+ const char *ptr;
+
+ ptr = check_box(ir->ePBC,box);
+ if (ptr) {
+ warning_error(wi,ptr);
+ }
+
+ if (bConstr && ir->eConstrAlg == econtSHAKE) {
+ if (ir->shake_tol <= 0.0) {
+ sprintf(warn_buf,"ERROR: shake_tol must be > 0 instead of %g\n",
+ ir->shake_tol);
+ warning_error(wi,warn_buf);
+ }
+
+ if (IR_TWINRANGE(*ir) && ir->nstlist > 1) {
+ sprintf(warn_buf,"With twin-range cut-off's and SHAKE the virial and the pressure are incorrect.");
+ if (ir->epc == epcNO) {
+ warning(wi,warn_buf);
+ } else {
+ warning_error(wi,warn_buf);
+ }
+ }
+ }
+
+ if( (ir->eConstrAlg == econtLINCS) && bConstr) {
+ /* If we have Lincs constraints: */
+ if(ir->eI==eiMD && ir->etc==etcNO &&
+ ir->eConstrAlg==econtLINCS && ir->nLincsIter==1) {
+ sprintf(warn_buf,"For energy conservation with LINCS, lincs_iter should be 2 or larger.\n");
+ warning_note(wi,warn_buf);
+ }
+
+ if ((ir->eI == eiCG || ir->eI == eiLBFGS) && (ir->nProjOrder<8)) {
+ sprintf(warn_buf,"For accurate %s with LINCS constraints, lincs_order should be 8 or more.",ei_names[ir->eI]);
+ warning_note(wi,warn_buf);
+ }
+ if (ir->epc==epcMTTK) {
+ warning_error(wi,"MTTK not compatible with lincs -- use shake instead.");
+ }
+ }
+
+ if (ir->LincsWarnAngle > 90.0) {
+ sprintf(warn_buf,"lincs-warnangle can not be larger than 90 degrees, setting it to 90.\n");
+ warning(wi,warn_buf);
+ ir->LincsWarnAngle = 90.0;
+ }
+
+ if (ir->ePBC != epbcNONE) {
+ if (ir->nstlist == 0) {
+ warning(wi,"With nstlist=0 atoms are only put into the box at step 0, therefore drifting atoms might cause the simulation to crash.");
+ }
+ bTWIN = (ir->rlistlong > ir->rlist);
+ if (ir->ns_type == ensGRID) {
+ if (sqr(ir->rlistlong) >= max_cutoff2(ir->ePBC,box)) {
+ sprintf(warn_buf,"ERROR: The cut-off length is longer than half the shortest box vector or longer than the smallest box diagonal element. Increase the box size or decrease %s.\n",
+ bTWIN ? (ir->rcoulomb==ir->rlistlong ? "rcoulomb" : "rvdw"):"rlist");
+ warning_error(wi,warn_buf);
+ }
+ } else {
+ min_size = min(box[XX][XX],min(box[YY][YY],box[ZZ][ZZ]));
+ if (2*ir->rlistlong >= min_size) {
+ sprintf(warn_buf,"ERROR: One of the box lengths is smaller than twice the cut-off length. Increase the box size or decrease rlist.");
+ warning_error(wi,warn_buf);
+ if (TRICLINIC(box))
+ fprintf(stderr,"Grid search might allow larger cut-off's than simple search with triclinic boxes.");
+ }
+ }
+ }
+}
+
+void check_chargegroup_radii(const gmx_mtop_t *mtop,const t_inputrec *ir,
+ rvec *x,
+ warninp_t wi)
+{
+ real rvdw1,rvdw2,rcoul1,rcoul2;
+ char warn_buf[STRLEN];
+
+ calc_chargegroup_radii(mtop,x,&rvdw1,&rvdw2,&rcoul1,&rcoul2);
+
+ if (rvdw1 > 0)
+ {
+ printf("Largest charge group radii for Van der Waals: %5.3f, %5.3f nm\n",
+ rvdw1,rvdw2);
+ }
+ if (rcoul1 > 0)
+ {
+ printf("Largest charge group radii for Coulomb: %5.3f, %5.3f nm\n",
+ rcoul1,rcoul2);
+ }
+
+ if (ir->rlist > 0)
+ {
+ if (rvdw1 + rvdw2 > ir->rlist ||
+ rcoul1 + rcoul2 > ir->rlist)
+ {
+ sprintf(warn_buf,"The sum of the two largest charge group radii (%f) is larger than rlist (%f)\n",max(rvdw1+rvdw2,rcoul1+rcoul2),ir->rlist);
+ warning(wi,warn_buf);
+ }
+ else
+ {
+ /* Here we do not use the zero at cut-off macro,
+ * since user defined interactions might purposely
+ * not be zero at the cut-off.
+ */
+ if (EVDW_IS_ZERO_AT_CUTOFF(ir->vdwtype) &&
+ rvdw1 + rvdw2 > ir->rlist - ir->rvdw)
+ {
+ sprintf(warn_buf,"The sum of the two largest charge group radii (%f) is larger than rlist (%f) - rvdw (%f)\n",
+ rvdw1+rvdw2,
+ ir->rlist,ir->rvdw);
+ if (ir_NVE(ir))
+ {
+ warning(wi,warn_buf);
+ }
+ else
+ {
+ warning_note(wi,warn_buf);
+ }
+ }
+ if (EEL_IS_ZERO_AT_CUTOFF(ir->coulombtype) &&
+ rcoul1 + rcoul2 > ir->rlistlong - ir->rcoulomb)
+ {
+ sprintf(warn_buf,"The sum of the two largest charge group radii (%f) is larger than %s (%f) - rcoulomb (%f)\n",
+ rcoul1+rcoul2,
+ ir->rlistlong > ir->rlist ? "rlistlong" : "rlist",
+ ir->rlistlong,ir->rcoulomb);
+ if (ir_NVE(ir))
+ {
+ warning(wi,warn_buf);
+ }
+ else
+ {
+ warning_note(wi,warn_buf);
+ }
+ }
+ }
+ }
+}
--- /dev/null
- "program."
+/* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
+ *
+ *
+ * This source code is part of
+ *
+ * G R O M A C S
+ *
+ * GROningen MAchine for Chemical Simulations
+ *
+ * VERSION 3.2.03
+ * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
+ * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
+ * Copyright (c) 2001-2004, The GROMACS development team,
+ * check out http://www.gromacs.org for more information.
+
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version.
+ *
+ * If you want to redistribute modifications, please consider that
+ * scientific software is very special. Version control is crucial -
+ * bugs must be traceable. We will be happy to consider code for
+ * inclusion in the official distribution, but derived work must not
+ * be called official GROMACS. Details are found in the README & COPYING
+ * files - if they are missing, get the official version at www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the papers on the package - you can find them in the top README file.
+ *
+ * For more info, check our website at http://www.gromacs.org
+ *
+ * And Hey:
+ * Gallium Rubidium Oxygen Manganese Argon Carbon Silicon
+ */
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <sys/types.h>
+#include <math.h>
+#include <string.h>
+#include <errno.h>
+#include <limits.h>
+
+#include "sysstuff.h"
+#include "smalloc.h"
+#include "macros.h"
+#include "string2.h"
+#include "readir.h"
+#include "toputil.h"
+#include "topio.h"
+#include "confio.h"
+#include "copyrite.h"
+#include "readir.h"
+#include "symtab.h"
+#include "names.h"
+#include "grompp.h"
+#include "random.h"
+#include "vec.h"
+#include "futil.h"
+#include "statutil.h"
+#include "splitter.h"
+#include "sortwater.h"
+#include "convparm.h"
+#include "gmx_fatal.h"
+#include "warninp.h"
+#include "index.h"
+#include "gmxfio.h"
+#include "trnio.h"
+#include "tpxio.h"
+#include "vsite_parm.h"
+#include "txtdump.h"
+#include "calcgrid.h"
+#include "add_par.h"
+#include "enxio.h"
+#include "perf_est.h"
+#include "compute_io.h"
+#include "gpp_atomtype.h"
+#include "gpp_tomorse.h"
+#include "mtop_util.h"
+#include "genborn.h"
+
+static int rm_interactions(int ifunc,int nrmols,t_molinfo mols[])
+{
+ int i,n;
+
+ n=0;
+ /* For all the molecule types */
+ for(i=0; i<nrmols; i++) {
+ n += mols[i].plist[ifunc].nr;
+ mols[i].plist[ifunc].nr=0;
+ }
+ return n;
+}
+
+static int check_atom_names(const char *fn1, const char *fn2,
+ gmx_mtop_t *mtop, t_atoms *at)
+{
+ int mb,m,i,j,nmismatch;
+ t_atoms *tat;
+#define MAXMISMATCH 20
+
+ if (mtop->natoms != at->nr)
+ gmx_incons("comparing atom names");
+
+ nmismatch=0;
+ i = 0;
+ for(mb=0; mb<mtop->nmolblock; mb++) {
+ tat = &mtop->moltype[mtop->molblock[mb].type].atoms;
+ for(m=0; m<mtop->molblock[mb].nmol; m++) {
+ for(j=0; j < tat->nr; j++) {
+ if (strcmp( *(tat->atomname[j]) , *(at->atomname[i]) ) != 0) {
+ if (nmismatch < MAXMISMATCH) {
+ fprintf(stderr,
+ "Warning: atom name %d in %s and %s does not match (%s - %s)\n",
+ i+1, fn1, fn2, *(tat->atomname[j]), *(at->atomname[i]));
+ } else if (nmismatch == MAXMISMATCH) {
+ fprintf(stderr,"(more than %d non-matching atom names)\n",MAXMISMATCH);
+ }
+ nmismatch++;
+ }
+ i++;
+ }
+ }
+ }
+
+ return nmismatch;
+}
+
+static void check_eg_vs_cg(gmx_mtop_t *mtop)
+{
+ int astart,mb,m,cg,j,firstj;
+ unsigned char firsteg,eg;
+ gmx_moltype_t *molt;
+
+ /* Go through all the charge groups and make sure all their
+ * atoms are in the same energy group.
+ */
+
+ astart = 0;
+ for(mb=0; mb<mtop->nmolblock; mb++) {
+ molt = &mtop->moltype[mtop->molblock[mb].type];
+ for(m=0; m<mtop->molblock[mb].nmol; m++) {
+ for(cg=0; cg<molt->cgs.nr;cg++) {
+ /* Get the energy group of the first atom in this charge group */
+ firstj = astart + molt->cgs.index[cg];
+ firsteg = ggrpnr(&mtop->groups,egcENER,firstj);
+ for(j=molt->cgs.index[cg]+1;j<molt->cgs.index[cg+1];j++) {
+ eg = ggrpnr(&mtop->groups,egcENER,astart+j);
+ if(eg != firsteg) {
+ gmx_fatal(FARGS,"atoms %d and %d in charge group %d of molecule type '%s' are in different energy groups",
+ firstj+1,astart+j+1,cg+1,*molt->name);
+ }
+ }
+ }
+ astart += molt->atoms.nr;
+ }
+ }
+}
+
+static void check_cg_sizes(const char *topfn,t_block *cgs,warninp_t wi)
+{
+ int maxsize,cg;
+ char warn_buf[STRLEN];
+
+ maxsize = 0;
+ for(cg=0; cg<cgs->nr; cg++)
+ {
+ maxsize = max(maxsize,cgs->index[cg+1]-cgs->index[cg]);
+ }
+
+ if (maxsize > MAX_CHARGEGROUP_SIZE)
+ {
+ gmx_fatal(FARGS,"The largest charge group contains %d atoms. The maximum is %d.",maxsize,MAX_CHARGEGROUP_SIZE);
+ }
+ else if (maxsize > 10)
+ {
+ set_warning_line(wi,topfn,-1);
+ sprintf(warn_buf,
+ "The largest charge group contains %d atoms.\n"
+ "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"
+ "For efficiency and accuracy, charge group should consist of a few atoms.\n"
+ "For all-atom force fields use: CH3, CH2, CH, NH2, NH, OH, CO2, CO, etc.",
+ maxsize);
+ warning_note(wi,warn_buf);
+ }
+}
+
+static void check_bonds_timestep(gmx_mtop_t *mtop,double dt,warninp_t wi)
+{
+ /* This check is not intended to ensure accurate integration,
+ * rather it is to signal mistakes in the mdp settings.
+ * A common mistake is to forget to turn on constraints
+ * for MD after energy minimization with flexible bonds.
+ * This check can also detect too large time steps for flexible water
+ * models, but such errors will often be masked by the constraints
+ * mdp options, which turns flexible water into water with bond constraints,
+ * but without an angle constraint. Unfortunately such incorrect use
+ * of water models can not easily be detected without checking
+ * for specific model names.
+ *
+ * The stability limit of leap-frog or velocity verlet is 4.44 steps
+ * per oscillational period.
+ * But accurate bonds distributions are lost far before that limit.
+ * To allow relatively common schemes (although not common with Gromacs)
+ * of dt=1 fs without constraints and dt=2 fs with only H-bond constraints
+ * we set the note limit to 10.
+ */
+ int min_steps_warn=5;
+ int min_steps_note=10;
+ t_iparams *ip;
+ int molt;
+ gmx_moltype_t *moltype,*w_moltype;
+ t_atom *atom;
+ t_ilist *ilist,*ilb,*ilc,*ils;
+ int ftype;
+ int i,a1,a2,w_a1,w_a2,j;
+ real twopi2,limit2,fc,re,m1,m2,period2,w_period2;
+ gmx_bool bFound,bWater,bWarn;
+ char warn_buf[STRLEN];
+
+ ip = mtop->ffparams.iparams;
+
+ twopi2 = sqr(2*M_PI);
+
+ limit2 = sqr(min_steps_note*dt);
+
+ w_a1 = w_a2 = -1;
+ w_period2 = -1.0;
+
+ w_moltype = NULL;
+ for(molt=0; molt<mtop->nmoltype; molt++)
+ {
+ moltype = &mtop->moltype[molt];
+ atom = moltype->atoms.atom;
+ ilist = moltype->ilist;
+ ilc = &ilist[F_CONSTR];
+ ils = &ilist[F_SETTLE];
+ for(ftype=0; ftype<F_NRE; ftype++)
+ {
+ if (!(ftype == F_BONDS || ftype == F_G96BONDS || ftype == F_HARMONIC))
+ {
+ continue;
+ }
+
+ ilb = &ilist[ftype];
+ for(i=0; i<ilb->nr; i+=3)
+ {
+ fc = ip[ilb->iatoms[i]].harmonic.krA;
+ re = ip[ilb->iatoms[i]].harmonic.rA;
+ if (ftype == F_G96BONDS)
+ {
+ /* Convert squared sqaure fc to harmonic fc */
+ fc = 2*fc*re;
+ }
+ a1 = ilb->iatoms[i+1];
+ a2 = ilb->iatoms[i+2];
+ m1 = atom[a1].m;
+ m2 = atom[a2].m;
+ if (fc > 0 && m1 > 0 && m2 > 0)
+ {
+ period2 = twopi2*m1*m2/((m1 + m2)*fc);
+ }
+ else
+ {
+ period2 = GMX_FLOAT_MAX;
+ }
+ if (debug)
+ {
+ fprintf(debug,"fc %g m1 %g m2 %g period %g\n",
+ fc,m1,m2,sqrt(period2));
+ }
+ if (period2 < limit2)
+ {
+ bFound = FALSE;
+ for(j=0; j<ilc->nr; j+=3)
+ {
+ if ((ilc->iatoms[j+1] == a1 && ilc->iatoms[j+2] == a2) ||
+ (ilc->iatoms[j+1] == a2 && ilc->iatoms[j+2] == a1))
+ {
+ bFound = TRUE;
+ }
+ }
+ for(j=0; j<ils->nr; j+=2)
+ {
+ if ((a1 >= ils->iatoms[j+1] && a1 < ils->iatoms[j+1]+3) &&
+ (a2 >= ils->iatoms[j+1] && a2 < ils->iatoms[j+1]+3))
+ {
+ bFound = TRUE;
+ }
+ }
+ if (!bFound &&
+ (w_moltype == NULL || period2 < w_period2))
+ {
+ w_moltype = moltype;
+ w_a1 = a1;
+ w_a2 = a2;
+ w_period2 = period2;
+ }
+ }
+ }
+ }
+ }
+
+ if (w_moltype != NULL)
+ {
+ bWarn = (w_period2 < sqr(min_steps_warn*dt));
+ /* A check that would recognize most water models */
+ bWater = ((*w_moltype->atoms.atomname[0])[0] == 'O' &&
+ w_moltype->atoms.nr <= 5);
+ 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"
+ "%s",
+ *w_moltype->name,
+ w_a1+1,*w_moltype->atoms.atomname[w_a1],
+ w_a2+1,*w_moltype->atoms.atomname[w_a2],
+ sqrt(w_period2),bWarn ? min_steps_warn : min_steps_note,dt,
+ bWater ?
+ "Maybe you asked for fexible water." :
+ "Maybe you forgot to change the constraints mdp option.");
+ if (bWarn)
+ {
+ warning(wi,warn_buf);
+ }
+ else
+ {
+ warning_note(wi,warn_buf);
+ }
+ }
+}
+
+static void check_vel(gmx_mtop_t *mtop,rvec v[])
+{
+ gmx_mtop_atomloop_all_t aloop;
+ t_atom *atom;
+ int a;
+
+ aloop = gmx_mtop_atomloop_all_init(mtop);
+ while (gmx_mtop_atomloop_all_next(aloop,&a,&atom)) {
+ if (atom->ptype == eptShell ||
+ atom->ptype == eptBond ||
+ atom->ptype == eptVSite) {
+ clear_rvec(v[a]);
+ }
+ }
+}
+
+static gmx_bool nint_ftype(gmx_mtop_t *mtop,t_molinfo *mi,int ftype)
+{
+ int nint,mb;
+
+ nint = 0;
+ for(mb=0; mb<mtop->nmolblock; mb++) {
+ nint += mtop->molblock[mb].nmol*mi[mtop->molblock[mb].type].plist[ftype].nr;
+ }
+
+ return nint;
+}
+
+/* This routine reorders the molecule type array
+ * in the order of use in the molblocks,
+ * unused molecule types are deleted.
+ */
+static void renumber_moltypes(gmx_mtop_t *sys,
+ int *nmolinfo,t_molinfo **molinfo)
+{
+ int *order,norder,i;
+ int mb,mi;
+ t_molinfo *minew;
+
+ snew(order,*nmolinfo);
+ norder = 0;
+ for(mb=0; mb<sys->nmolblock; mb++) {
+ for(i=0; i<norder; i++) {
+ if (order[i] == sys->molblock[mb].type) {
+ break;
+ }
+ }
+ if (i == norder) {
+ /* This type did not occur yet, add it */
+ order[norder] = sys->molblock[mb].type;
+ /* Renumber the moltype in the topology */
+ norder++;
+ }
+ sys->molblock[mb].type = i;
+ }
+
+ /* We still need to reorder the molinfo structs */
+ snew(minew,norder);
+ for(mi=0; mi<*nmolinfo; mi++) {
+ for(i=0; i<norder; i++) {
+ if (order[i] == mi) {
+ break;
+ }
+ }
+ if (i == norder) {
+ done_mi(&(*molinfo)[mi]);
+ } else {
+ minew[i] = (*molinfo)[mi];
+ }
+ }
+ sfree(*molinfo);
+
+ *nmolinfo = norder;
+ *molinfo = minew;
+}
+
+static void molinfo2mtop(int nmi,t_molinfo *mi,gmx_mtop_t *mtop)
+{
+ int m;
+ gmx_moltype_t *molt;
+
+ mtop->nmoltype = nmi;
+ snew(mtop->moltype,nmi);
+ for(m=0; m<nmi; m++) {
+ molt = &mtop->moltype[m];
+ molt->name = mi[m].name;
+ molt->atoms = mi[m].atoms;
+ /* ilists are copied later */
+ molt->cgs = mi[m].cgs;
+ molt->excls = mi[m].excls;
+ }
+}
+
+static void
+new_status(const char *topfile,const char *topppfile,const char *confin,
+ t_gromppopts *opts,t_inputrec *ir,gmx_bool bZero,
+ gmx_bool bGenVel,gmx_bool bVerbose,t_state *state,
+ gpp_atomtype_t atype,gmx_mtop_t *sys,
+ int *nmi,t_molinfo **mi,t_params plist[],
+ int *comb,double *reppow,real *fudgeQQ,
+ gmx_bool bMorse,
+ warninp_t wi)
+{
+ t_molinfo *molinfo=NULL;
+ int nmolblock;
+ gmx_molblock_t *molblock,*molbs;
+ t_atoms *confat;
+ int mb,i,nrmols,nmismatch;
+ char buf[STRLEN];
+ gmx_bool bGB=FALSE;
+ char warn_buf[STRLEN];
+
+ init_mtop(sys);
+
+ /* Set gmx_boolean for GB */
+ if(ir->implicit_solvent)
+ bGB=TRUE;
+
+ /* TOPOLOGY processing */
+ sys->name = do_top(bVerbose,topfile,topppfile,opts,bZero,&(sys->symtab),
+ plist,comb,reppow,fudgeQQ,
+ atype,&nrmols,&molinfo,ir,
+ &nmolblock,&molblock,bGB,
+ wi);
+
+ sys->nmolblock = 0;
+ snew(sys->molblock,nmolblock);
+
+ sys->natoms = 0;
+ for(mb=0; mb<nmolblock; mb++) {
+ if (sys->nmolblock > 0 &&
+ molblock[mb].type == sys->molblock[sys->nmolblock-1].type) {
+ /* Merge consecutive blocks with the same molecule type */
+ sys->molblock[sys->nmolblock-1].nmol += molblock[mb].nmol;
+ sys->natoms += molblock[mb].nmol*sys->molblock[sys->nmolblock-1].natoms_mol;
+ } else if (molblock[mb].nmol > 0) {
+ /* Add a new molblock to the topology */
+ molbs = &sys->molblock[sys->nmolblock];
+ *molbs = molblock[mb];
+ molbs->natoms_mol = molinfo[molbs->type].atoms.nr;
+ molbs->nposres_xA = 0;
+ molbs->nposres_xB = 0;
+ sys->natoms += molbs->nmol*molbs->natoms_mol;
+ sys->nmolblock++;
+ }
+ }
+ if (sys->nmolblock == 0) {
+ gmx_fatal(FARGS,"No molecules were defined in the system");
+ }
+
+ renumber_moltypes(sys,&nrmols,&molinfo);
+
+ if (bMorse)
+ convert_harmonics(nrmols,molinfo,atype);
+
+ if (ir->eDisre == edrNone) {
+ i = rm_interactions(F_DISRES,nrmols,molinfo);
+ if (i > 0) {
+ set_warning_line(wi,"unknown",-1);
+ sprintf(warn_buf,"disre = no, removed %d distance restraints",i);
+ warning_note(wi,warn_buf);
+ }
+ }
+ if (opts->bOrire == FALSE) {
+ i = rm_interactions(F_ORIRES,nrmols,molinfo);
+ if (i > 0) {
+ set_warning_line(wi,"unknown",-1);
+ sprintf(warn_buf,"orire = no, removed %d orientation restraints",i);
+ warning_note(wi,warn_buf);
+ }
+ }
+ if (opts->bDihre == FALSE) {
+ i = rm_interactions(F_DIHRES,nrmols,molinfo);
+ if (i > 0) {
+ set_warning_line(wi,"unknown",-1);
+ sprintf(warn_buf,"dihre = no, removed %d dihedral restraints",i);
+ warning_note(wi,warn_buf);
+ }
+ }
+
+ /* Copy structures from msys to sys */
+ molinfo2mtop(nrmols,molinfo,sys);
+
+ gmx_mtop_finalize(sys);
+
+ /* COORDINATE file processing */
+ if (bVerbose)
+ fprintf(stderr,"processing coordinates...\n");
+
+ get_stx_coordnum(confin,&state->natoms);
+ if (state->natoms != sys->natoms)
+ gmx_fatal(FARGS,"number of coordinates in coordinate file (%s, %d)\n"
+ " does not match topology (%s, %d)",
+ confin,state->natoms,topfile,sys->natoms);
+ else {
+ /* make space for coordinates and velocities */
+ char title[STRLEN];
+ snew(confat,1);
+ init_t_atoms(confat,state->natoms,FALSE);
+ init_state(state,state->natoms,0,0,0);
+ read_stx_conf(confin,title,confat,state->x,state->v,NULL,state->box);
+ /* This call fixes the box shape for runs with pressure scaling */
+ set_box_rel(ir,state);
+
+ nmismatch = check_atom_names(topfile, confin, sys, confat);
+ free_t_atoms(confat,TRUE);
+ sfree(confat);
+
+ if (nmismatch) {
+ sprintf(buf,"%d non-matching atom name%s\n"
+ "atom names from %s will be used\n"
+ "atom names from %s will be ignored\n",
+ nmismatch,(nmismatch == 1) ? "" : "s",topfile,confin);
+ warning(wi,buf);
+ }
+ if (bVerbose)
+ fprintf(stderr,"double-checking input for internal consistency...\n");
+ double_check(ir,state->box,nint_ftype(sys,molinfo,F_CONSTR),wi);
+ }
+
+ if (bGenVel) {
+ real *mass;
+ gmx_mtop_atomloop_all_t aloop;
+ t_atom *atom;
+
+ snew(mass,state->natoms);
+ aloop = gmx_mtop_atomloop_all_init(sys);
+ while (gmx_mtop_atomloop_all_next(aloop,&i,&atom)) {
+ mass[i] = atom->m;
+ }
+
+ if (opts->seed == -1) {
+ opts->seed = make_seed();
+ fprintf(stderr,"Setting gen_seed to %d\n",opts->seed);
+ }
+ maxwell_speed(opts->tempi,opts->seed,sys,state->v);
+
+ stop_cm(stdout,state->natoms,mass,state->x,state->v);
+ sfree(mass);
+ }
+
+ *nmi = nrmols;
+ *mi = molinfo;
+}
+
+static void copy_state(const char *slog,t_trxframe *fr,
+ gmx_bool bReadVel,t_state *state,
+ double *use_time)
+{
+ int i;
+
+ if (fr->not_ok & FRAME_NOT_OK)
+ {
+ gmx_fatal(FARGS,"Can not start from an incomplete frame");
+ }
+ if (!fr->bX)
+ {
+ gmx_fatal(FARGS,"Did not find a frame with coordinates in file %s",
+ slog);
+ }
+
+ for(i=0; i<state->natoms; i++)
+ {
+ copy_rvec(fr->x[i],state->x[i]);
+ }
+ if (bReadVel)
+ {
+ if (!fr->bV)
+ {
+ gmx_incons("Trajecory frame unexpectedly does not contain velocities");
+ }
+ for(i=0; i<state->natoms; i++)
+ {
+ copy_rvec(fr->v[i],state->v[i]);
+ }
+ }
+ if (fr->bBox)
+ {
+ copy_mat(fr->box,state->box);
+ }
+
+ *use_time = fr->time;
+}
+
+static void cont_status(const char *slog,const char *ener,
+ gmx_bool bNeedVel,gmx_bool bGenVel, real fr_time,
+ t_inputrec *ir,t_state *state,
+ gmx_mtop_t *sys,
+ const output_env_t oenv)
+ /* If fr_time == -1 read the last frame available which is complete */
+{
+ gmx_bool bReadVel;
+ t_trxframe fr;
+ t_trxstatus *fp;
+ int i;
+ double use_time;
+
+ bReadVel = (bNeedVel && !bGenVel);
+
+ fprintf(stderr,
+ "Reading Coordinates%s and Box size from old trajectory\n",
+ bReadVel ? ", Velocities" : "");
+ if (fr_time == -1)
+ {
+ fprintf(stderr,"Will read whole trajectory\n");
+ }
+ else
+ {
+ fprintf(stderr,"Will read till time %g\n",fr_time);
+ }
+ if (!bReadVel)
+ {
+ if (bGenVel)
+ {
+ fprintf(stderr,"Velocities generated: "
+ "ignoring velocities in input trajectory\n");
+ }
+ read_first_frame(oenv,&fp,slog,&fr,TRX_NEED_X);
+ }
+ else
+ {
+ read_first_frame(oenv,&fp,slog,&fr,TRX_NEED_X | TRX_NEED_V);
+
+ if (!fr.bV)
+ {
+ fprintf(stderr,
+ "\n"
+ "WARNING: Did not find a frame with velocities in file %s,\n"
+ " all velocities will be set to zero!\n\n",slog);
+ for(i=0; i<sys->natoms; i++)
+ {
+ clear_rvec(state->v[i]);
+ }
+ close_trj(fp);
+ /* Search for a frame without velocities */
+ bReadVel = FALSE;
+ read_first_frame(oenv,&fp,slog,&fr,TRX_NEED_X);
+ }
+ }
+
+ state->natoms = fr.natoms;
+
+ if (sys->natoms != state->natoms)
+ {
+ gmx_fatal(FARGS,"Number of atoms in Topology "
+ "is not the same as in Trajectory");
+ }
+ copy_state(slog,&fr,bReadVel,state,&use_time);
+
+ /* Find the appropriate frame */
+ while ((fr_time == -1 || fr.time < fr_time) &&
+ read_next_frame(oenv,fp,&fr))
+ {
+ copy_state(slog,&fr,bReadVel,state,&use_time);
+ }
+
+ close_trj(fp);
+
+ /* Set the relative box lengths for preserving the box shape.
+ * Note that this call can lead to differences in the last bit
+ * with respect to using tpbconv to create a tpx file.
+ */
+ set_box_rel(ir,state);
+
+ fprintf(stderr,"Using frame at t = %g ps\n",use_time);
+ fprintf(stderr,"Starting time for run is %g ps\n",ir->init_t);
+
+ if ((ir->epc != epcNO || ir->etc ==etcNOSEHOOVER) && ener)
+ {
+ get_enx_state(ener,use_time,&sys->groups,ir,state);
+ preserve_box_shape(ir,state->box_rel,state->boxv);
+ }
+}
+
+static void read_posres(gmx_mtop_t *mtop,t_molinfo *molinfo,gmx_bool bTopB,
+ char *fn,
+ int rc_scaling, int ePBC,
+ rvec com,
+ warninp_t wi)
+{
+ gmx_bool bFirst = TRUE;
+ rvec *x,*v,*xp;
+ dvec sum;
+ double totmass;
+ t_atoms dumat;
+ matrix box,invbox;
+ int natoms,npbcdim=0;
+ char warn_buf[STRLEN],title[STRLEN];
+ int a,i,ai,j,k,mb,nat_molb;
+ gmx_molblock_t *molb;
+ t_params *pr;
+ t_atom *atom;
+
+ get_stx_coordnum(fn,&natoms);
+ if (natoms != mtop->natoms) {
+ 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);
+ warning(wi,warn_buf);
+ }
+ snew(x,natoms);
+ snew(v,natoms);
+ init_t_atoms(&dumat,natoms,FALSE);
+ read_stx_conf(fn,title,&dumat,x,v,NULL,box);
+
+ npbcdim = ePBC2npbcdim(ePBC);
+ clear_rvec(com);
+ if (rc_scaling != erscNO) {
+ copy_mat(box,invbox);
+ for(j=npbcdim; j<DIM; j++) {
+ clear_rvec(invbox[j]);
+ invbox[j][j] = 1;
+ }
+ m_inv_ur0(invbox,invbox);
+ }
+
+ /* Copy the reference coordinates to mtop */
+ clear_dvec(sum);
+ totmass = 0;
+ a = 0;
+ for(mb=0; mb<mtop->nmolblock; mb++) {
+ molb = &mtop->molblock[mb];
+ nat_molb = molb->nmol*mtop->moltype[molb->type].atoms.nr;
+ pr = &(molinfo[molb->type].plist[F_POSRES]);
+ if (pr->nr > 0) {
+ atom = mtop->moltype[molb->type].atoms.atom;
+ for(i=0; (i<pr->nr); i++) {
+ ai=pr->param[i].AI;
+ if (ai >= natoms) {
+ gmx_fatal(FARGS,"Position restraint atom index (%d) in moltype '%s' is larger than number of atoms in %s (%d).\n",
+ ai+1,*molinfo[molb->type].name,fn,natoms);
+ }
+ if (rc_scaling == erscCOM) {
+ /* Determine the center of mass of the posres reference coordinates */
+ for(j=0; j<npbcdim; j++) {
+ sum[j] += atom[ai].m*x[a+ai][j];
+ }
+ totmass += atom[ai].m;
+ }
+ }
+ if (!bTopB) {
+ molb->nposres_xA = nat_molb;
+ snew(molb->posres_xA,molb->nposres_xA);
+ for(i=0; i<nat_molb; i++) {
+ copy_rvec(x[a+i],molb->posres_xA[i]);
+ }
+ } else {
+ molb->nposres_xB = nat_molb;
+ snew(molb->posres_xB,molb->nposres_xB);
+ for(i=0; i<nat_molb; i++) {
+ copy_rvec(x[a+i],molb->posres_xB[i]);
+ }
+ }
+ }
+ a += nat_molb;
+ }
+ if (rc_scaling == erscCOM) {
+ if (totmass == 0)
+ gmx_fatal(FARGS,"The total mass of the position restraint atoms is 0");
+ for(j=0; j<npbcdim; j++)
+ com[j] = sum[j]/totmass;
+ 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]);
+ }
+
+ if (rc_scaling != erscNO) {
+ for(mb=0; mb<mtop->nmolblock; mb++) {
+ molb = &mtop->molblock[mb];
+ nat_molb = molb->nmol*mtop->moltype[molb->type].atoms.nr;
+ if (molb->nposres_xA > 0 || molb->nposres_xB > 0) {
+ xp = (!bTopB ? molb->posres_xA : molb->posres_xB);
+ for(i=0; i<nat_molb; i++) {
+ for(j=0; j<npbcdim; j++) {
+ if (rc_scaling == erscALL) {
+ /* Convert from Cartesian to crystal coordinates */
+ xp[i][j] *= invbox[j][j];
+ for(k=j+1; k<npbcdim; k++) {
+ xp[i][j] += invbox[k][j]*xp[i][k];
+ }
+ } else if (rc_scaling == erscCOM) {
+ /* Subtract the center of mass */
+ xp[i][j] -= com[j];
+ }
+ }
+ }
+ }
+ }
+
+ if (rc_scaling == erscCOM) {
+ /* Convert the COM from Cartesian to crystal coordinates */
+ for(j=0; j<npbcdim; j++) {
+ com[j] *= invbox[j][j];
+ for(k=j+1; k<npbcdim; k++) {
+ com[j] += invbox[k][j]*com[k];
+ }
+ }
+ }
+ }
+
+ free_t_atoms(&dumat,TRUE);
+ sfree(x);
+ sfree(v);
+}
+
+static void gen_posres(gmx_mtop_t *mtop,t_molinfo *mi,
+ char *fnA, char *fnB,
+ int rc_scaling, int ePBC,
+ rvec com, rvec comB,
+ warninp_t wi)
+{
+ int i,j;
+
+ read_posres (mtop,mi,FALSE,fnA,rc_scaling,ePBC,com,wi);
+ if (strcmp(fnA,fnB) != 0) {
+ read_posres(mtop,mi,TRUE ,fnB,rc_scaling,ePBC,comB,wi);
+ }
+}
+
+static void set_wall_atomtype(gpp_atomtype_t at,t_gromppopts *opts,
+ t_inputrec *ir)
+{
+ int i;
+
+ if (ir->nwall > 0)
+ fprintf(stderr,"Searching the wall atom type(s)\n");
+ for(i=0; i<ir->nwall; i++)
+ ir->wall_atomtype[i] = get_atomtype_type(opts->wall_atomtype[i],at);
+}
+
+static int nrdf_internal(t_atoms *atoms)
+{
+ int i,nmass,nrdf;
+
+ nmass = 0;
+ for(i=0; i<atoms->nr; i++) {
+ /* Vsite ptype might not be set here yet, so also check the mass */
+ if ((atoms->atom[i].ptype == eptAtom ||
+ atoms->atom[i].ptype == eptNucleus)
+ && atoms->atom[i].m > 0) {
+ nmass++;
+ }
+ }
+ switch (nmass) {
+ case 0: nrdf = 0; break;
+ case 1: nrdf = 0; break;
+ case 2: nrdf = 1; break;
+ default: nrdf = nmass*3 - 6; break;
+ }
+
+ return nrdf;
+}
+
+void
+spline1d( double dx,
+ double * y,
+ int n,
+ double * u,
+ double * y2 )
+{
+ int i;
+ double p,q;
+
+ y2[0] = 0.0;
+ u[0] = 0.0;
+
+ for(i=1;i<n-1;i++)
+ {
+ p = 0.5*y2[i-1]+2.0;
+ y2[i] = -0.5/p;
+ q = (y[i+1]-2.0*y[i]+y[i-1])/dx;
+ u[i] = (3.0*q/dx-0.5*u[i-1])/p;
+ }
+
+ y2[n-1] = 0.0;
+
+ for(i=n-2;i>=0;i--)
+ {
+ y2[i] = y2[i]*y2[i+1]+u[i];
+ }
+}
+
+
+void
+interpolate1d( double xmin,
+ double dx,
+ double * ya,
+ double * y2a,
+ double x,
+ double * y,
+ double * y1)
+{
+ int ix;
+ double a,b;
+
+ ix = (x-xmin)/dx;
+
+ a = (xmin+(ix+1)*dx-x)/dx;
+ b = (x-xmin-ix*dx)/dx;
+
+ *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;
+ *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];
+}
+
+
+void
+setup_cmap (int grid_spacing,
+ int nc,
+ real * grid ,
+ gmx_cmap_t * cmap_grid)
+{
+ double *tmp_u,*tmp_u2,*tmp_yy,*tmp_y1,*tmp_t2,*tmp_grid;
+
+ int i,j,k,ii,jj,kk,idx;
+ int offset;
+ double dx,xmin,v,v1,v2,v12;
+ double phi,psi;
+
+ snew(tmp_u,2*grid_spacing);
+ snew(tmp_u2,2*grid_spacing);
+ snew(tmp_yy,2*grid_spacing);
+ snew(tmp_y1,2*grid_spacing);
+ snew(tmp_t2,2*grid_spacing*2*grid_spacing);
+ snew(tmp_grid,2*grid_spacing*2*grid_spacing);
+
+ dx = 360.0/grid_spacing;
+ xmin = -180.0-dx*grid_spacing/2;
+
+ for(kk=0;kk<nc;kk++)
+ {
+ /* Compute an offset depending on which cmap we are using
+ * Offset will be the map number multiplied with the grid_spacing * grid_spacing * 2
+ */
+ offset = kk * grid_spacing * grid_spacing * 2;
+
+ for(i=0;i<2*grid_spacing;i++)
+ {
+ ii=(i+grid_spacing-grid_spacing/2)%grid_spacing;
+
+ for(j=0;j<2*grid_spacing;j++)
+ {
+ jj=(j+grid_spacing-grid_spacing/2)%grid_spacing;
+ tmp_grid[i*grid_spacing*2+j] = grid[offset+ii*grid_spacing+jj];
+ }
+ }
+
+ for(i=0;i<2*grid_spacing;i++)
+ {
+ spline1d(dx,&(tmp_grid[2*grid_spacing*i]),2*grid_spacing,tmp_u,&(tmp_t2[2*grid_spacing*i]));
+ }
+
+ for(i=grid_spacing/2;i<grid_spacing+grid_spacing/2;i++)
+ {
+ ii = i-grid_spacing/2;
+ phi = ii*dx-180.0;
+
+ for(j=grid_spacing/2;j<grid_spacing+grid_spacing/2;j++)
+ {
+ jj = j-grid_spacing/2;
+ psi = jj*dx-180.0;
+
+ for(k=0;k<2*grid_spacing;k++)
+ {
+ interpolate1d(xmin,dx,&(tmp_grid[2*grid_spacing*k]),
+ &(tmp_t2[2*grid_spacing*k]),psi,&tmp_yy[k],&tmp_y1[k]);
+ }
+
+ spline1d(dx,tmp_yy,2*grid_spacing,tmp_u,tmp_u2);
+ interpolate1d(xmin,dx,tmp_yy,tmp_u2,phi,&v,&v1);
+ spline1d(dx,tmp_y1,2*grid_spacing,tmp_u,tmp_u2);
+ interpolate1d(xmin,dx,tmp_y1,tmp_u2,phi,&v2,&v12);
+
+ idx = ii*grid_spacing+jj;
+ cmap_grid->cmapdata[kk].cmap[idx*4] = grid[offset+ii*grid_spacing+jj];
+ cmap_grid->cmapdata[kk].cmap[idx*4+1] = v1;
+ cmap_grid->cmapdata[kk].cmap[idx*4+2] = v2;
+ cmap_grid->cmapdata[kk].cmap[idx*4+3] = v12;
+ }
+ }
+ }
+}
+
+void init_cmap_grid(gmx_cmap_t *cmap_grid, int ngrid, int grid_spacing)
+{
+ int i,k,nelem;
+
+ cmap_grid->ngrid = ngrid;
+ cmap_grid->grid_spacing = grid_spacing;
+ nelem = cmap_grid->grid_spacing*cmap_grid->grid_spacing;
+
+ snew(cmap_grid->cmapdata,ngrid);
+
+ for(i=0;i<cmap_grid->ngrid;i++)
+ {
+ snew(cmap_grid->cmapdata[i].cmap,4*nelem);
+ }
+}
+
+
+static int count_constraints(gmx_mtop_t *mtop,t_molinfo *mi,warninp_t wi)
+{
+ int count,count_mol,i,mb;
+ gmx_molblock_t *molb;
+ t_params *plist;
+ char buf[STRLEN];
+
+ count = 0;
+ for(mb=0; mb<mtop->nmolblock; mb++) {
+ count_mol = 0;
+ molb = &mtop->molblock[mb];
+ plist = mi[molb->type].plist;
+
+ for(i=0; i<F_NRE; i++) {
+ if (i == F_SETTLE)
+ count_mol += 3*plist[i].nr;
+ else if (interaction_function[i].flags & IF_CONSTRAINT)
+ count_mol += plist[i].nr;
+ }
+
+ if (count_mol > nrdf_internal(&mi[molb->type].atoms)) {
+ sprintf(buf,
+ "Molecule type '%s' has %d constraints.\n"
+ "For stability and efficiency there should not be more constraints than internal number of degrees of freedom: %d.\n",
+ *mi[molb->type].name,count_mol,
+ nrdf_internal(&mi[molb->type].atoms));
+ warning(wi,buf);
+ }
+ count += molb->nmol*count_mol;
+ }
+
+ return count;
+}
+
+static void check_gbsa_params_charged(gmx_mtop_t *sys, gpp_atomtype_t atype)
+{
+ int i,nmiss,natoms,mt;
+ real q;
+ const t_atoms *atoms;
+
+ nmiss = 0;
+ for(mt=0;mt<sys->nmoltype;mt++)
+ {
+ atoms = &sys->moltype[mt].atoms;
+ natoms = atoms->nr;
+
+ for(i=0;i<natoms;i++)
+ {
+ q = atoms->atom[i].q;
+ if ((get_atomtype_radius(atoms->atom[i].type,atype) == 0 ||
+ get_atomtype_vol(atoms->atom[i].type,atype) == 0 ||
+ get_atomtype_surftens(atoms->atom[i].type,atype) == 0 ||
+ get_atomtype_gb_radius(atoms->atom[i].type,atype) == 0 ||
+ get_atomtype_S_hct(atoms->atom[i].type,atype) == 0) &&
+ q != 0)
+ {
+ fprintf(stderr,"\nGB parameter(s) zero for atom type '%s' while charge is %g\n",
+ get_atomtype_name(atoms->atom[i].type,atype),q);
+ nmiss++;
+ }
+ }
+ }
+
+ if (nmiss > 0)
+ {
+ 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);
+ }
+}
+
+
+static void check_gbsa_params(t_inputrec *ir,gpp_atomtype_t atype)
+{
+ int nmiss,i;
+
+ /* If we are doing GBSA, check that we got the parameters we need
+ * This checking is to see if there are GBSA paratmeters for all
+ * atoms in the force field. To go around this for testing purposes
+ * comment out the nerror++ counter temporarily
+ */
+ nmiss = 0;
+ for(i=0;i<get_atomtype_ntypes(atype);i++)
+ {
+ if (get_atomtype_radius(i,atype) < 0 ||
+ get_atomtype_vol(i,atype) < 0 ||
+ get_atomtype_surftens(i,atype) < 0 ||
+ get_atomtype_gb_radius(i,atype) < 0 ||
+ get_atomtype_S_hct(i,atype) < 0)
+ {
+ fprintf(stderr,"\nGB parameter(s) missing or negative for atom type '%s'\n",
+ get_atomtype_name(i,atype));
+ nmiss++;
+ }
+ }
+
+ if (nmiss > 0)
+ {
+ 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);
+ }
+
+}
+
+int main (int argc, char *argv[])
+{
+ static const char *desc[] = {
+ "The gromacs preprocessor",
+ "reads a molecular topology file, checks the validity of the",
+ "file, expands the topology from a molecular description to an atomic",
+ "description. The topology file contains information about",
+ "molecule types and the number of molecules, the preprocessor",
+ "copies each molecule as needed. ",
+ "There is no limitation on the number of molecule types. ",
+ "Bonds and bond-angles can be converted into constraints, separately",
+ "for hydrogens and heavy atoms.",
+ "Then a coordinate file is read and velocities can be generated",
+ "from a Maxwellian distribution if requested.",
+ "grompp also reads parameters for the mdrun ",
+ "(eg. number of MD steps, time step, cut-off), and others such as",
+ "NEMD parameters, which are corrected so that the net acceleration",
+ "is zero.",
+ "Eventually a binary file is produced that can serve as the sole input",
+ "file for the MD program.[PAR]",
+
+ "grompp uses the atom names from the topology file. The atom names",
+ "in the coordinate file (option [TT]-c[tt]) are only read to generate",
+ "warnings when they do not match the atom names in the topology.",
+ "Note that the atom names are irrelevant for the simulation as",
+ "only the atom types are used for generating interaction parameters.[PAR]",
+
+ "grompp uses a built-in preprocessor to resolve includes, macros ",
+ "etcetera. The preprocessor supports the following keywords:[BR]",
+ "#ifdef VARIABLE[BR]",
+ "#ifndef VARIABLE[BR]",
+ "#else[BR]",
+ "#endif[BR]",
+ "#define VARIABLE[BR]",
+ "#undef VARIABLE[BR]"
+ "#include \"filename\"[BR]",
+ "#include <filename>[BR]",
+ "The functioning of these statements in your topology may be modulated by",
+ "using the following two flags in your [TT]mdp[tt] file:[BR]",
+ "define = -DVARIABLE1 -DVARIABLE2[BR]",
+ "include = -I/home/john/doe[BR]",
+ "For further information a C-programming textbook may help you out.",
+ "Specifying the [TT]-pp[tt] flag will get the pre-processed",
+ "topology file written out so that you can verify its contents.[PAR]",
+
+ "If your system does not have a c-preprocessor, you can still",
+ "use grompp, but you do not have access to the features ",
+ "from the cpp. Command line options to the c-preprocessor can be given",
+ "in the [TT].mdp[tt] file. See your local manual (man cpp).[PAR]",
+
+ "When using position restraints a file with restraint coordinates",
+ "can be supplied with [TT]-r[tt], otherwise restraining will be done",
+ "with respect to the conformation from the [TT]-c[tt] option.",
+ "For free energy calculation the the coordinates for the B topology",
+ "can be supplied with [TT]-rb[tt], otherwise they will be equal to",
+ "those of the A topology.[PAR]",
+
+ "Starting coordinates can be read from trajectory with [TT]-t[tt].",
+ "The last frame with coordinates and velocities will be read,",
+ "unless the [TT]-time[tt] option is used.",
+ "Note that these velocities will not be used when [TT]gen_vel = yes[tt]",
+ "in your [TT].mdp[tt] file. An energy file can be supplied with",
+ "[TT]-e[tt] to read Nose-Hoover and/or Parrinello-Rahman coupling",
+ "variables. Note that for continuation it is better and easier to supply",
+ "a checkpoint file directly to mdrun, since that always contains",
+ "the complete state of the system and you don't need to generate",
+ "a new run input file. Note that if you only want to change the number",
+ "of run steps tpbconv is more convenient than grompp.[PAR]",
+
+ "By default all bonded interactions which have constant energy due to",
+ "virtual site constructions will be removed. If this constant energy is",
+ "not zero, this will result in a shift in the total energy. All bonded",
+ "interactions can be kept by turning off [TT]-rmvsbds[tt]. Additionally,",
+ "all constraints for distances which will be constant anyway because",
+ "of virtual site constructions will be removed. If any constraints remain",
+ "which involve virtual sites, a fatal error will result.[PAR]"
+
+ "To verify your run input file, please make notice of all warnings",
+ "on the screen, and correct where necessary. Do also look at the contents",
+ "of the [TT]mdout.mdp[tt] file, this contains comment lines, as well as",
+ "the input that [TT]grompp[tt] has read. If in doubt you can start grompp",
+ "with the [TT]-debug[tt] option which will give you more information",
+ "in a file called grompp.log (along with real debug info). Finally, you",
+ "can see the contents of the run input file with the [TT]gmxdump[tt]",
- "Number of allowed warnings during input processing" },
++ "program.[PAR]"
++
++ "The [TT]-maxwarn[tt] option can be used to override warnings printed",
++ "by grompp that otherwise halt output. In some cases, warnings are",
++ "harmless, but usually they are not. The user is advised to carefully",
++ "interpret the output messages before attempting to bypass them with",
++ "this option."
+ };
+ t_gromppopts *opts;
+ gmx_mtop_t *sys;
+ int nmi;
+ t_molinfo *mi;
+ gpp_atomtype_t atype;
+ t_inputrec *ir;
+ int natoms,nvsite,comb,mt;
+ t_params *plist;
+ t_state state;
+ matrix box;
+ real max_spacing,fudgeQQ;
+ double reppow;
+ char fn[STRLEN],fnB[STRLEN];
+ const char *mdparin;
+ int ntype;
+ gmx_bool bNeedVel,bGenVel;
+ gmx_bool have_atomnumber;
+ int n12,n13,n14;
+ t_params *gb_plist = NULL;
+ gmx_genborn_t *born = NULL;
+ output_env_t oenv;
+ gmx_bool bVerbose = FALSE;
+ warninp_t wi;
+ char warn_buf[STRLEN];
+
+ t_filenm fnm[] = {
+ { efMDP, NULL, NULL, ffREAD },
+ { efMDP, "-po", "mdout", ffWRITE },
+ { efSTX, "-c", NULL, ffREAD },
+ { efSTX, "-r", NULL, ffOPTRD },
+ { efSTX, "-rb", NULL, ffOPTRD },
+ { efNDX, NULL, NULL, ffOPTRD },
+ { efTOP, NULL, NULL, ffREAD },
+ { efTOP, "-pp", "processed", ffOPTWR },
+ { efTPX, "-o", NULL, ffWRITE },
+ { efTRN, "-t", NULL, ffOPTRD },
+ { efEDR, "-e", NULL, ffOPTRD },
+ { efTRN, "-ref","rotref", ffOPTRW }
+ };
+#define NFILE asize(fnm)
+
+ /* Command line options */
+ static gmx_bool bRenum=TRUE;
+ static gmx_bool bRmVSBds=TRUE,bZero=FALSE;
+ static int i,maxwarn=0;
+ static real fr_time=-1;
+ t_pargs pa[] = {
+ { "-v", FALSE, etBOOL,{&bVerbose},
+ "Be loud and noisy" },
+ { "-time", FALSE, etREAL, {&fr_time},
+ "Take frame at or first after this time." },
+ { "-rmvsbds",FALSE, etBOOL, {&bRmVSBds},
+ "Remove constant bonded interactions with virtual sites" },
+ { "-maxwarn", FALSE, etINT, {&maxwarn},
++ "Number of allowed warnings during input processing. Not for normal use and may generate unstable systems" },
+ { "-zero", FALSE, etBOOL, {&bZero},
+ "Set parameters for bonded interactions without defaults to zero instead of generating an error" },
+ { "-renum", FALSE, etBOOL, {&bRenum},
+ "Renumber atomtypes and minimize number of atomtypes" }
+ };
+
+ CopyRight(stdout,argv[0]);
+
+ /* Initiate some variables */
+ snew(ir,1);
+ snew(opts,1);
+ init_ir(ir,opts);
+
+ /* Parse the command line */
+ parse_common_args(&argc,argv,0,NFILE,fnm,asize(pa),pa,
+ asize(desc),desc,0,NULL,&oenv);
+
+ wi = init_warning(TRUE,maxwarn);
+
+ /* PARAMETER file processing */
+ mdparin = opt2fn("-f",NFILE,fnm);
+ set_warning_line(wi,mdparin,-1);
+ get_ir(mdparin,opt2fn("-po",NFILE,fnm),ir,opts,wi);
+
+ if (bVerbose)
+ fprintf(stderr,"checking input for internal consistency...\n");
+ check_ir(mdparin,ir,opts,wi);
+
+ if (ir->ld_seed == -1) {
+ ir->ld_seed = make_seed();
+ fprintf(stderr,"Setting the LD random seed to %d\n",ir->ld_seed);
+ }
+
+ bNeedVel = EI_STATE_VELOCITY(ir->eI);
+ bGenVel = (bNeedVel && opts->bGenVel);
+
+ snew(plist,F_NRE);
+ init_plist(plist);
+ snew(sys,1);
+ atype = init_atomtype();
+ if (debug)
+ pr_symtab(debug,0,"Just opened",&sys->symtab);
+
+ strcpy(fn,ftp2fn(efTOP,NFILE,fnm));
+ if (!gmx_fexist(fn))
+ gmx_fatal(FARGS,"%s does not exist",fn);
+ new_status(fn,opt2fn_null("-pp",NFILE,fnm),opt2fn("-c",NFILE,fnm),
+ opts,ir,bZero,bGenVel,bVerbose,&state,
+ atype,sys,&nmi,&mi,plist,&comb,&reppow,&fudgeQQ,
+ opts->bMorse,
+ wi);
+
+ if (debug)
+ pr_symtab(debug,0,"After new_status",&sys->symtab);
+
+ if (count_constraints(sys,mi,wi) && (ir->eConstrAlg == econtSHAKE)) {
+ if (ir->eI == eiCG || ir->eI == eiLBFGS) {
+ sprintf(warn_buf,"Can not do %s with %s, use %s",
+ EI(ir->eI),econstr_names[econtSHAKE],econstr_names[econtLINCS]);
+ warning_error(wi,warn_buf);
+ }
+ if (ir->bPeriodicMols) {
+ sprintf(warn_buf,"Can not do periodic molecules with %s, use %s",
+ econstr_names[econtSHAKE],econstr_names[econtLINCS]);
+ warning_error(wi,warn_buf);
+ }
+ }
+
+ /* If we are doing QM/MM, check that we got the atom numbers */
+ have_atomnumber = TRUE;
+ for (i=0; i<get_atomtype_ntypes(atype); i++) {
+ have_atomnumber = have_atomnumber && (get_atomtype_atomnumber(i,atype) >= 0);
+ }
+ if (!have_atomnumber && ir->bQMMM)
+ {
+ warning_error(wi,
+ "\n"
+ "It appears as if you are trying to run a QM/MM calculation, but the force\n"
+ "field you are using does not contain atom numbers fields. This is an\n"
+ "optional field (introduced in Gromacs 3.3) for general runs, but mandatory\n"
+ "for QM/MM. The good news is that it is easy to add - put the atom number as\n"
+ "an integer just before the mass column in ffXXXnb.itp.\n"
+ "NB: United atoms have the same atom numbers as normal ones.\n\n");
+ }
+
+ /* Check for errors in the input now, since they might cause problems
+ * during processing further down.
+ */
+ check_warning_error(wi,FARGS);
+
+ if (opt2bSet("-r",NFILE,fnm))
+ sprintf(fn,"%s",opt2fn("-r",NFILE,fnm));
+ else
+ sprintf(fn,"%s",opt2fn("-c",NFILE,fnm));
+ if (opt2bSet("-rb",NFILE,fnm))
+ sprintf(fnB,"%s",opt2fn("-rb",NFILE,fnm));
+ else
+ strcpy(fnB,fn);
+
+ if (nint_ftype(sys,mi,F_POSRES) > 0)
+ {
+ if (bVerbose)
+ {
+ fprintf(stderr,"Reading position restraint coords from %s",fn);
+ if (strcmp(fn,fnB) == 0)
+ {
+ fprintf(stderr,"\n");
+ }
+ else
+ {
+ fprintf(stderr," and %s\n",fnB);
+ if (ir->efep != efepNO && ir->n_flambda > 0)
+ {
+ warning_error(wi,"Can not change the position restraint reference coordinates with lambda togther with foreign lambda calculation.");
+ }
+ }
+ }
+ gen_posres(sys,mi,fn,fnB,
+ ir->refcoord_scaling,ir->ePBC,
+ ir->posres_com,ir->posres_comB,
+ wi);
+ }
+
+ nvsite = 0;
+ /* set parameters for virtual site construction (not for vsiten) */
+ for(mt=0; mt<sys->nmoltype; mt++) {
+ nvsite +=
+ set_vsites(bVerbose, &sys->moltype[mt].atoms, atype, mi[mt].plist);
+ }
+ /* now throw away all obsolete bonds, angles and dihedrals: */
+ /* note: constraints are ALWAYS removed */
+ if (nvsite) {
+ for(mt=0; mt<sys->nmoltype; mt++) {
+ clean_vsite_bondeds(mi[mt].plist,sys->moltype[mt].atoms.nr,bRmVSBds);
+ }
+ }
+
+ /* If we are using CMAP, setup the pre-interpolation grid */
+ if(plist->ncmap>0)
+ {
+ init_cmap_grid(&sys->ffparams.cmap_grid, plist->nc, plist->grid_spacing);
+ setup_cmap(plist->grid_spacing, plist->nc, plist->cmap,&sys->ffparams.cmap_grid);
+ }
+
+ set_wall_atomtype(atype,opts,ir);
+ if (bRenum) {
+ renum_atype(plist, sys, ir->wall_atomtype, atype, bVerbose);
+ ntype = get_atomtype_ntypes(atype);
+ }
+
+ if (ir->implicit_solvent != eisNO)
+ {
+ /* Now we have renumbered the atom types, we can check the GBSA params */
+ check_gbsa_params(ir,atype);
+
+ /* Check that all atoms that have charge and/or LJ-parameters also have
+ * sensible GB-parameters
+ */
+ check_gbsa_params_charged(sys,atype);
+ }
+
+ /* PELA: Copy the atomtype data to the topology atomtype list */
+ copy_atomtype_atomtypes(atype,&(sys->atomtypes));
+
+ if (debug)
+ pr_symtab(debug,0,"After renum_atype",&sys->symtab);
+
+ if (bVerbose)
+ fprintf(stderr,"converting bonded parameters...\n");
+
+ ntype = get_atomtype_ntypes(atype);
+ convert_params(ntype, plist, mi, comb, reppow, fudgeQQ, sys);
+
+ if (debug)
+ pr_symtab(debug,0,"After convert_params",&sys->symtab);
+
+ /* set ptype to VSite for virtual sites */
+ for(mt=0; mt<sys->nmoltype; mt++) {
+ set_vsites_ptype(FALSE,&sys->moltype[mt]);
+ }
+ if (debug) {
+ pr_symtab(debug,0,"After virtual sites",&sys->symtab);
+ }
+ /* Check velocity for virtual sites and shells */
+ if (bGenVel) {
+ check_vel(sys,state.v);
+ }
+
+ /* check masses */
+ check_mol(sys,wi);
+
+ for(i=0; i<sys->nmoltype; i++) {
+ check_cg_sizes(ftp2fn(efTOP,NFILE,fnm),&sys->moltype[i].cgs,wi);
+ }
+
+ if (EI_DYNAMICS(ir->eI) && ir->eI != eiBD)
+ {
+ check_bonds_timestep(sys,ir->delta_t,wi);
+ }
+
+ check_warning_error(wi,FARGS);
+
+ if (bVerbose)
+ fprintf(stderr,"initialising group options...\n");
+ do_index(mdparin,ftp2fn_null(efNDX,NFILE,fnm),
+ sys,bVerbose,ir,
+ bGenVel ? state.v : NULL,
+ wi);
+
+ /* Init the temperature coupling state */
+ init_gtc_state(&state,ir->opts.ngtc,0,ir->opts.nhchainlength);
+
+ if (bVerbose)
+ fprintf(stderr,"Checking consistency between energy and charge groups...\n");
+ check_eg_vs_cg(sys);
+
+ if (debug)
+ pr_symtab(debug,0,"After index",&sys->symtab);
+ triple_check(mdparin,ir,sys,wi);
+ close_symtab(&sys->symtab);
+ if (debug)
+ pr_symtab(debug,0,"After close",&sys->symtab);
+
+ /* make exclusions between QM atoms */
+ if (ir->bQMMM) {
+ generate_qmexcl(sys,ir);
+ }
+
+ if (ftp2bSet(efTRN,NFILE,fnm)) {
+ if (bVerbose)
+ fprintf(stderr,"getting data from old trajectory ...\n");
+ cont_status(ftp2fn(efTRN,NFILE,fnm),ftp2fn_null(efEDR,NFILE,fnm),
+ bNeedVel,bGenVel,fr_time,ir,&state,sys,oenv);
+ }
+
+ if (ir->ePBC==epbcXY && ir->nwall!=2)
+ {
+ clear_rvec(state.box[ZZ]);
+ }
+
+ if (ir->rlist > 0)
+ {
+ set_warning_line(wi,mdparin,-1);
+ check_chargegroup_radii(sys,ir,state.x,wi);
+ }
+
+ if (EEL_FULL(ir->coulombtype)) {
+ /* Calculate the optimal grid dimensions */
+ copy_mat(state.box,box);
+ if (ir->ePBC==epbcXY && ir->nwall==2)
+ svmul(ir->wall_ewald_zfac,box[ZZ],box[ZZ]);
+ max_spacing = calc_grid(stdout,box,opts->fourierspacing,
+ &(ir->nkx),&(ir->nky),&(ir->nkz));
+ if ((ir->coulombtype == eelPPPM) && (max_spacing > 0.1)) {
+ set_warning_line(wi,mdparin,-1);
+ warning_note(wi,"Grid spacing larger then 0.1 while using PPPM.");
+ }
+ }
+
+ if (ir->ePull != epullNO)
+ set_pull_init(ir,sys,state.x,state.box,oenv,opts->pull_start);
+
+ if (ir->bRot)
+ {
+ set_reference_positions(ir->rot,sys,state.x,state.box,
+ opt2fn("-ref",NFILE,fnm),opt2bSet("-ref",NFILE,fnm),
+ wi);
+ }
+
+ /* reset_multinr(sys); */
+
+ if (EEL_PME(ir->coulombtype)) {
+ float ratio = pme_load_estimate(sys,ir,state.box);
+ fprintf(stderr,"Estimate for the relative computational load of the PME mesh part: %.2f\n",ratio);
+ /* With free energy we might need to do PME both for the A and B state
+ * charges. This will double the cost, but the optimal performance will
+ * then probably be at a slightly larger cut-off and grid spacing.
+ */
+ if ((ir->efep == efepNO && ratio > 1.0/2.0) ||
+ (ir->efep != efepNO && ratio > 2.0/3.0)) {
+ warning_note(wi,
+ "The optimal PME mesh load for parallel simulations is below 0.5\n"
+ "and for highly parallel simulations between 0.25 and 0.33,\n"
+ "for higher performance, increase the cut-off and the PME grid spacing");
+ }
+ }
+
+ {
+ char warn_buf[STRLEN];
+ double cio = compute_io(ir,sys->natoms,&sys->groups,F_NRE,1);
+ sprintf(warn_buf,"This run will generate roughly %.0f Mb of data",cio);
+ if (cio > 2000) {
+ set_warning_line(wi,mdparin,-1);
+ warning_note(wi,warn_buf);
+ } else {
+ printf("%s\n",warn_buf);
+ }
+ }
+
+ if (bVerbose)
+ fprintf(stderr,"writing run input file...\n");
+
+ done_warning(wi,FARGS);
+
+ state.lambda = ir->init_lambda;
+ write_tpx_state(ftp2fn(efTPX,NFILE,fnm),ir,&state,sys);
+
+ thanx(stderr);
+
+ return 0;
+}
--- /dev/null
+/* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
+ *
+ *
+ * This source code is part of
+ *
+ * G R O M A C S
+ *
+ * GROningen MAchine for Chemical Simulations
+ *
+ * VERSION 3.2.0
+ * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
+ * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
+ * Copyright (c) 2001-2004, The GROMACS development team,
+ * check out http://www.gromacs.org for more information.
+
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version.
+ *
+ * If you want to redistribute modifications, please consider that
+ * scientific software is very special. Version control is crucial -
+ * bugs must be traceable. We will be happy to consider code for
+ * inclusion in the official distribution, but derived work must not
+ * be called official GROMACS. Details are found in the README & COPYING
+ * files - if they are missing, get the official version at www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the papers on the package - you can find them in the top README file.
+ *
+ * For more info, check our website at http://www.gromacs.org
+ *
+ * And Hey:
+ * Gallium Rubidium Oxygen Manganese Argon Carbon Silicon
+ */
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <signal.h>
+#include <stdlib.h>
+
+#if ((defined WIN32 || defined _WIN32 || defined WIN64 || defined _WIN64) && !defined __CYGWIN__ && !defined __CYGWIN32__)
+/* _isnan() */
+#include <float.h>
+#endif
+
+#include "typedefs.h"
+#include "smalloc.h"
+#include "sysstuff.h"
+#include "vec.h"
+#include "statutil.h"
+#include "vcm.h"
+#include "mdebin.h"
+#include "nrnb.h"
+#include "calcmu.h"
+#include "index.h"
+#include "vsite.h"
+#include "update.h"
+#include "ns.h"
+#include "trnio.h"
+#include "xtcio.h"
+#include "mdrun.h"
+#include "confio.h"
+#include "network.h"
+#include "pull.h"
+#include "xvgr.h"
+#include "physics.h"
+#include "names.h"
+#include "xmdrun.h"
+#include "ionize.h"
+#include "disre.h"
+#include "orires.h"
+#include "dihre.h"
+#include "pppm.h"
+#include "pme.h"
+#include "mdatoms.h"
+#include "repl_ex.h"
+#include "qmmm.h"
+#include "mpelogging.h"
+#include "domdec.h"
+#include "partdec.h"
+#include "topsort.h"
+#include "coulomb.h"
+#include "constr.h"
+#include "shellfc.h"
+#include "compute_io.h"
+#include "mvdata.h"
+#include "checkpoint.h"
+#include "mtop_util.h"
+#include "sighandler.h"
+#include "membed.h"
+
+#ifdef GMX_LIB_MPI
+#include <mpi.h>
+#endif
+#ifdef GMX_THREADS
+#include "tmpi.h"
+#endif
+
+#ifdef GMX_FAHCORE
+#include "corewrap.h"
+#endif
+
+
+/* simulation conditions to transmit. Keep in mind that they are
+ transmitted to other nodes through an MPI_Reduce after
+ casting them to a real (so the signals can be sent together with other
+ data). This means that the only meaningful values are positive,
+ negative or zero. */
+enum { eglsNABNSB, eglsCHKPT, eglsSTOPCOND, eglsRESETCOUNTERS, eglsNR };
+/* Is the signal in one simulation independent of other simulations? */
+gmx_bool gs_simlocal[eglsNR] = { TRUE, FALSE, FALSE, TRUE };
+
+typedef struct {
+ int nstms; /* The frequency for intersimulation communication */
+ int sig[eglsNR]; /* The signal set by one process in do_md */
+ int set[eglsNR]; /* The communicated signal, equal for all processes */
+} globsig_t;
+
+
+static int multisim_min(const gmx_multisim_t *ms,int nmin,int n)
+{
+ int *buf;
+ gmx_bool bPos,bEqual;
+ int s,d;
+
+ snew(buf,ms->nsim);
+ buf[ms->sim] = n;
+ gmx_sumi_sim(ms->nsim,buf,ms);
+ bPos = TRUE;
+ bEqual = TRUE;
+ for(s=0; s<ms->nsim; s++)
+ {
+ bPos = bPos && (buf[s] > 0);
+ bEqual = bEqual && (buf[s] == buf[0]);
+ }
+ if (bPos)
+ {
+ if (bEqual)
+ {
+ nmin = min(nmin,buf[0]);
+ }
+ else
+ {
+ /* Find the least common multiple */
+ for(d=2; d<nmin; d++)
+ {
+ s = 0;
+ while (s < ms->nsim && d % buf[s] == 0)
+ {
+ s++;
+ }
+ if (s == ms->nsim)
+ {
+ /* We found the LCM and it is less than nmin */
+ nmin = d;
+ break;
+ }
+ }
+ }
+ }
+ sfree(buf);
+
+ return nmin;
+}
+
+static int multisim_nstsimsync(const t_commrec *cr,
+ const t_inputrec *ir,int repl_ex_nst)
+{
+ int nmin;
+
+ if (MASTER(cr))
+ {
+ nmin = INT_MAX;
+ nmin = multisim_min(cr->ms,nmin,ir->nstlist);
+ nmin = multisim_min(cr->ms,nmin,ir->nstcalcenergy);
+ nmin = multisim_min(cr->ms,nmin,repl_ex_nst);
+ if (nmin == INT_MAX)
+ {
+ gmx_fatal(FARGS,"Can not find an appropriate interval for inter-simulation communication, since nstlist, nstcalcenergy and -replex are all <= 0");
+ }
+ /* Avoid inter-simulation communication at every (second) step */
+ if (nmin <= 2)
+ {
+ nmin = 10;
+ }
+ }
+
+ gmx_bcast(sizeof(int),&nmin,cr);
+
+ return nmin;
+}
+
+static void init_global_signals(globsig_t *gs,const t_commrec *cr,
+ const t_inputrec *ir,int repl_ex_nst)
+{
+ int i;
+
+ if (MULTISIM(cr))
+ {
+ gs->nstms = multisim_nstsimsync(cr,ir,repl_ex_nst);
+ if (debug)
+ {
+ fprintf(debug,"Syncing simulations for checkpointing and termination every %d steps\n",gs->nstms);
+ }
+ }
+ else
+ {
+ gs->nstms = 1;
+ }
+
+ for(i=0; i<eglsNR; i++)
+ {
+ gs->sig[i] = 0;
+ gs->set[i] = 0;
+ }
+}
+
+static void copy_coupling_state(t_state *statea,t_state *stateb,
+ gmx_ekindata_t *ekinda,gmx_ekindata_t *ekindb, t_grpopts* opts)
+{
+
+ /* MRS note -- might be able to get rid of some of the arguments. Look over it when it's all debugged */
+
+ int i,j,nc;
+
+ /* Make sure we have enough space for x and v */
+ if (statea->nalloc > stateb->nalloc)
+ {
+ stateb->nalloc = statea->nalloc;
+ srenew(stateb->x,stateb->nalloc);
+ srenew(stateb->v,stateb->nalloc);
+ }
+
+ stateb->natoms = statea->natoms;
+ stateb->ngtc = statea->ngtc;
+ stateb->nnhpres = statea->nnhpres;
+ stateb->veta = statea->veta;
+ if (ekinda)
+ {
+ copy_mat(ekinda->ekin,ekindb->ekin);
+ for (i=0; i<stateb->ngtc; i++)
+ {
+ ekindb->tcstat[i].T = ekinda->tcstat[i].T;
+ ekindb->tcstat[i].Th = ekinda->tcstat[i].Th;
+ copy_mat(ekinda->tcstat[i].ekinh,ekindb->tcstat[i].ekinh);
+ copy_mat(ekinda->tcstat[i].ekinf,ekindb->tcstat[i].ekinf);
+ ekindb->tcstat[i].ekinscalef_nhc = ekinda->tcstat[i].ekinscalef_nhc;
+ ekindb->tcstat[i].ekinscaleh_nhc = ekinda->tcstat[i].ekinscaleh_nhc;
+ ekindb->tcstat[i].vscale_nhc = ekinda->tcstat[i].vscale_nhc;
+ }
+ }
+ copy_rvecn(statea->x,stateb->x,0,stateb->natoms);
+ copy_rvecn(statea->v,stateb->v,0,stateb->natoms);
+ copy_mat(statea->box,stateb->box);
+ copy_mat(statea->box_rel,stateb->box_rel);
+ copy_mat(statea->boxv,stateb->boxv);
+
+ for (i = 0; i<stateb->ngtc; i++)
+ {
+ nc = i*opts->nhchainlength;
+ for (j=0; j<opts->nhchainlength; j++)
+ {
+ stateb->nosehoover_xi[nc+j] = statea->nosehoover_xi[nc+j];
+ stateb->nosehoover_vxi[nc+j] = statea->nosehoover_vxi[nc+j];
+ }
+ }
+ if (stateb->nhpres_xi != NULL)
+ {
+ for (i = 0; i<stateb->nnhpres; i++)
+ {
+ nc = i*opts->nhchainlength;
+ for (j=0; j<opts->nhchainlength; j++)
+ {
+ stateb->nhpres_xi[nc+j] = statea->nhpres_xi[nc+j];
+ stateb->nhpres_vxi[nc+j] = statea->nhpres_vxi[nc+j];
+ }
+ }
+ }
+}
+
+static real compute_conserved_from_auxiliary(t_inputrec *ir, t_state *state, t_extmass *MassQ)
+{
+ real quantity = 0;
+ switch (ir->etc)
+ {
+ case etcNO:
+ break;
+ case etcBERENDSEN:
+ break;
+ case etcNOSEHOOVER:
+ quantity = NPT_energy(ir,state,MassQ);
+ break;
+ case etcVRESCALE:
+ quantity = vrescale_energy(&(ir->opts),state->therm_integral);
+ break;
+ default:
+ break;
+ }
+ return quantity;
+}
+
+static void compute_globals(FILE *fplog, gmx_global_stat_t gstat, t_commrec *cr, t_inputrec *ir,
+ t_forcerec *fr, gmx_ekindata_t *ekind,
+ t_state *state, t_state *state_global, t_mdatoms *mdatoms,
+ t_nrnb *nrnb, t_vcm *vcm, gmx_wallcycle_t wcycle,
+ gmx_enerdata_t *enerd,tensor force_vir, tensor shake_vir, tensor total_vir,
+ tensor pres, rvec mu_tot, gmx_constr_t constr,
+ globsig_t *gs,gmx_bool bInterSimGS,
+ matrix box, gmx_mtop_t *top_global, real *pcurr,
+ int natoms, gmx_bool *bSumEkinhOld, int flags)
+{
+ int i,gsi;
+ real gs_buf[eglsNR];
+ tensor corr_vir,corr_pres,shakeall_vir;
+ gmx_bool bEner,bPres,bTemp, bVV;
+ gmx_bool bRerunMD, bStopCM, bGStat, bIterate,
+ bFirstIterate,bReadEkin,bEkinAveVel,bScaleEkin, bConstrain;
+ real ekin,temp,prescorr,enercorr,dvdlcorr;
+
+ /* translate CGLO flags to gmx_booleans */
+ bRerunMD = flags & CGLO_RERUNMD;
+ bStopCM = flags & CGLO_STOPCM;
+ bGStat = flags & CGLO_GSTAT;
+
+ bReadEkin = (flags & CGLO_READEKIN);
+ bScaleEkin = (flags & CGLO_SCALEEKIN);
+ bEner = flags & CGLO_ENERGY;
+ bTemp = flags & CGLO_TEMPERATURE;
+ bPres = (flags & CGLO_PRESSURE);
+ bConstrain = (flags & CGLO_CONSTRAINT);
+ bIterate = (flags & CGLO_ITERATE);
+ bFirstIterate = (flags & CGLO_FIRSTITERATE);
+
+ /* we calculate a full state kinetic energy either with full-step velocity verlet
+ or half step where we need the pressure */
+
+ bEkinAveVel = (ir->eI==eiVV || (ir->eI==eiVVAK && bPres) || bReadEkin);
+
+ /* in initalization, it sums the shake virial in vv, and to
+ sums ekinh_old in leapfrog (or if we are calculating ekinh_old) for other reasons */
+
+ /* ########## Kinetic energy ############## */
+
+ if (bTemp)
+ {
+ /* Non-equilibrium MD: this is parallellized, but only does communication
+ * when there really is NEMD.
+ */
+
+ if (PAR(cr) && (ekind->bNEMD))
+ {
+ accumulate_u(cr,&(ir->opts),ekind);
+ }
+ debug_gmx();
+ if (bReadEkin)
+ {
+ restore_ekinstate_from_state(cr,ekind,&state_global->ekinstate);
+ }
+ else
+ {
+
+ calc_ke_part(state,&(ir->opts),mdatoms,ekind,nrnb,bEkinAveVel,bIterate);
+ }
+
+ debug_gmx();
+
+ /* Calculate center of mass velocity if necessary, also parallellized */
+ if (bStopCM && !bRerunMD && bEner)
+ {
+ calc_vcm_grp(fplog,mdatoms->start,mdatoms->homenr,mdatoms,
+ state->x,state->v,vcm);
+ }
+ }
+
+ if (bTemp || bPres || bEner || bConstrain)
+ {
+ if (!bGStat)
+ {
+ /* We will not sum ekinh_old,
+ * so signal that we still have to do it.
+ */
+ *bSumEkinhOld = TRUE;
+
+ }
+ else
+ {
+ if (gs != NULL)
+ {
+ for(i=0; i<eglsNR; i++)
+ {
+ gs_buf[i] = gs->sig[i];
+ }
+ }
+ if (PAR(cr))
+ {
+ wallcycle_start(wcycle,ewcMoveE);
+ GMX_MPE_LOG(ev_global_stat_start);
+ global_stat(fplog,gstat,cr,enerd,force_vir,shake_vir,mu_tot,
+ ir,ekind,constr,vcm,
+ gs != NULL ? eglsNR : 0,gs_buf,
+ top_global,state,
+ *bSumEkinhOld,flags);
+ GMX_MPE_LOG(ev_global_stat_finish);
+ wallcycle_stop(wcycle,ewcMoveE);
+ }
+ if (gs != NULL)
+ {
+ if (MULTISIM(cr) && bInterSimGS)
+ {
+ if (MASTER(cr))
+ {
+ /* Communicate the signals between the simulations */
+ gmx_sum_sim(eglsNR,gs_buf,cr->ms);
+ }
+ /* Communicate the signals form the master to the others */
+ gmx_bcast(eglsNR*sizeof(gs_buf[0]),gs_buf,cr);
+ }
+ for(i=0; i<eglsNR; i++)
+ {
+ if (bInterSimGS || gs_simlocal[i])
+ {
+ /* Set the communicated signal only when it is non-zero,
+ * since signals might not be processed at each MD step.
+ */
+ gsi = (gs_buf[i] >= 0 ?
+ (int)(gs_buf[i] + 0.5) :
+ (int)(gs_buf[i] - 0.5));
+ if (gsi != 0)
+ {
+ gs->set[i] = gsi;
+ }
+ /* Turn off the local signal */
+ gs->sig[i] = 0;
+ }
+ }
+ }
+ *bSumEkinhOld = FALSE;
+ }
+ }
+
+ if (!ekind->bNEMD && debug && bTemp && (vcm->nr > 0))
+ {
+ correct_ekin(debug,
+ mdatoms->start,mdatoms->start+mdatoms->homenr,
+ state->v,vcm->group_p[0],
+ mdatoms->massT,mdatoms->tmass,ekind->ekin);
+ }
+
+ if (bEner) {
+ /* Do center of mass motion removal */
+ if (bStopCM && !bRerunMD) /* is this correct? Does it get called too often with this logic? */
+ {
+ check_cm_grp(fplog,vcm,ir,1);
+ do_stopcm_grp(fplog,mdatoms->start,mdatoms->homenr,mdatoms->cVCM,
+ state->x,state->v,vcm);
+ inc_nrnb(nrnb,eNR_STOPCM,mdatoms->homenr);
+ }
++
++ /* Calculate the amplitude of the cosine velocity profile */
++ ekind->cosacc.vcos = ekind->cosacc.mvcos/mdatoms->tmass;
+ }
+
+ if (bTemp)
+ {
+ /* Sum the kinetic energies of the groups & calc temp */
+ /* compute full step kinetic energies if vv, or if vv-avek and we are computing the pressure with IR_NPT_TROTTER */
+ /* three maincase: VV with AveVel (md-vv), vv with AveEkin (md-vv-avek), leap with AveEkin (md).
+ Leap with AveVel is not supported; it's not clear that it will actually work.
+ bEkinAveVel: If TRUE, we simply multiply ekin by ekinscale to get a full step kinetic energy.
+ If FALSE, we average ekinh_old and ekinh*ekinscale_nhc to get an averaged half step kinetic energy.
+ bSaveEkinOld: If TRUE (in the case of iteration = bIterate is TRUE), we don't reset the ekinscale_nhc.
+ If FALSE, we go ahead and erase over it.
+ */
+ enerd->term[F_TEMP] = sum_ekin(&(ir->opts),ekind,&(enerd->term[F_DKDL]),
+ bEkinAveVel,bIterate,bScaleEkin);
+
+ enerd->term[F_EKIN] = trace(ekind->ekin);
+ }
+
+ /* ########## Long range energy information ###### */
+
+ if (bEner || bPres || bConstrain)
+ {
+ calc_dispcorr(fplog,ir,fr,0,top_global->natoms,box,state->lambda,
+ corr_pres,corr_vir,&prescorr,&enercorr,&dvdlcorr);
+ }
+
+ if (bEner && bFirstIterate)
+ {
+ enerd->term[F_DISPCORR] = enercorr;
+ enerd->term[F_EPOT] += enercorr;
+ enerd->term[F_DVDL] += dvdlcorr;
+ if (fr->efep != efepNO) {
+ enerd->dvdl_lin += dvdlcorr;
+ }
+ }
+
+ /* ########## Now pressure ############## */
+ if (bPres || bConstrain)
+ {
+
+ m_add(force_vir,shake_vir,total_vir);
+
+ /* Calculate pressure and apply LR correction if PPPM is used.
+ * Use the box from last timestep since we already called update().
+ */
+
+ enerd->term[F_PRES] = calc_pres(fr->ePBC,ir->nwall,box,ekind->ekin,total_vir,pres,
+ (fr->eeltype==eelPPPM)?enerd->term[F_COUL_RECIP]:0.0);
+
+ /* Calculate long range corrections to pressure and energy */
+ /* this adds to enerd->term[F_PRES] and enerd->term[F_ETOT],
+ and computes enerd->term[F_DISPCORR]. Also modifies the
+ total_vir and pres tesors */
+
+ m_add(total_vir,corr_vir,total_vir);
+ m_add(pres,corr_pres,pres);
+ enerd->term[F_PDISPCORR] = prescorr;
+ enerd->term[F_PRES] += prescorr;
+ *pcurr = enerd->term[F_PRES];
+ /* calculate temperature using virial */
+ enerd->term[F_VTEMP] = calc_temp(trace(total_vir),ir->opts.nrdf[0]);
+
+ }
+}
+
+
+/* Definitions for convergence of iterated constraints */
+
+/* iterate constraints up to 50 times */
+#define MAXITERCONST 50
+
+/* data type */
+typedef struct
+{
+ real f,fprev,x,xprev;
+ int iter_i;
+ gmx_bool bIterate;
+ real allrelerr[MAXITERCONST+2];
+ int num_close; /* number of "close" violations, caused by limited precision. */
+} gmx_iterate_t;
+
+#ifdef GMX_DOUBLE
+#define CONVERGEITER 0.000000001
+#define CLOSE_ENOUGH 0.000001000
+#else
+#define CONVERGEITER 0.0001
+#define CLOSE_ENOUGH 0.0050
+#endif
+
+/* we want to keep track of the close calls. If there are too many, there might be some other issues.
+ so we make sure that it's either less than some predetermined number, or if more than that number,
+ only some small fraction of the total. */
+#define MAX_NUMBER_CLOSE 50
+#define FRACTION_CLOSE 0.001
+
+/* maximum length of cyclic traps to check, emerging from limited numerical precision */
+#define CYCLEMAX 20
+
+static void gmx_iterate_init(gmx_iterate_t *iterate,gmx_bool bIterate)
+{
+ int i;
+
+ iterate->iter_i = 0;
+ iterate->bIterate = bIterate;
+ iterate->num_close = 0;
+ for (i=0;i<MAXITERCONST+2;i++)
+ {
+ iterate->allrelerr[i] = 0;
+ }
+}
+
+static gmx_bool done_iterating(const t_commrec *cr,FILE *fplog, int nsteps, gmx_iterate_t *iterate, gmx_bool bFirstIterate, real fom, real *newf)
+{
+ /* monitor convergence, and use a secant search to propose new
+ values.
+ x_{i} - x_{i-1}
+ The secant method computes x_{i+1} = x_{i} - f(x_{i}) * ---------------------
+ f(x_{i}) - f(x_{i-1})
+
+ The function we are trying to zero is fom-x, where fom is the
+ "figure of merit" which is the pressure (or the veta value) we
+ would get by putting in an old value of the pressure or veta into
+ the incrementor function for the step or half step. I have
+ verified that this gives the same answer as self consistent
+ iteration, usually in many fewer steps, especially for small tau_p.
+
+ We could possibly eliminate an iteration with proper use
+ of the value from the previous step, but that would take a bit
+ more bookkeeping, especially for veta, since tests indicate the
+ function of veta on the last step is not sufficiently close to
+ guarantee convergence this step. This is
+ good enough for now. On my tests, I could use tau_p down to
+ 0.02, which is smaller that would ever be necessary in
+ practice. Generally, 3-5 iterations will be sufficient */
+
+ real relerr,err,xmin;
+ char buf[256];
+ int i;
+ gmx_bool incycle;
+
+ if (bFirstIterate)
+ {
+ iterate->x = fom;
+ iterate->f = fom-iterate->x;
+ iterate->xprev = 0;
+ iterate->fprev = 0;
+ *newf = fom;
+ }
+ else
+ {
+ iterate->f = fom-iterate->x; /* we want to zero this difference */
+ if ((iterate->iter_i > 1) && (iterate->iter_i < MAXITERCONST))
+ {
+ if (iterate->f==iterate->fprev)
+ {
+ *newf = iterate->f;
+ }
+ else
+ {
+ *newf = iterate->x - (iterate->x-iterate->xprev)*(iterate->f)/(iterate->f-iterate->fprev);
+ }
+ }
+ else
+ {
+ /* just use self-consistent iteration the first step to initialize, or
+ if it's not converging (which happens occasionally -- need to investigate why) */
+ *newf = fom;
+ }
+ }
+ /* Consider a slight shortcut allowing us to exit one sooner -- we check the
+ difference between the closest of x and xprev to the new
+ value. To be 100% certain, we should check the difference between
+ the last result, and the previous result, or
+
+ relerr = (fabs((x-xprev)/fom));
+
+ but this is pretty much never necessary under typical conditions.
+ Checking numerically, it seems to lead to almost exactly the same
+ trajectories, but there are small differences out a few decimal
+ places in the pressure, and eventually in the v_eta, but it could
+ save an interation.
+
+ if (fabs(*newf-x) < fabs(*newf - xprev)) { xmin = x;} else { xmin = xprev;}
+ relerr = (fabs((*newf-xmin) / *newf));
+ */
+
+ err = fabs((iterate->f-iterate->fprev));
+ relerr = fabs(err/fom);
+
+ iterate->allrelerr[iterate->iter_i] = relerr;
+
+ if (iterate->iter_i > 0)
+ {
+ if (debug)
+ {
+ fprintf(debug,"Iterating NPT constraints: %6i %20.12f%14.6g%20.12f\n",
+ iterate->iter_i,fom,relerr,*newf);
+ }
+
+ if ((relerr < CONVERGEITER) || (err < CONVERGEITER) || (fom==0) || ((iterate->x == iterate->xprev) && iterate->iter_i > 1))
+ {
+ iterate->bIterate = FALSE;
+ if (debug)
+ {
+ fprintf(debug,"Iterating NPT constraints: CONVERGED\n");
+ }
+ return TRUE;
+ }
+ if (iterate->iter_i > MAXITERCONST)
+ {
+ if (relerr < CLOSE_ENOUGH)
+ {
+ incycle = FALSE;
+ for (i=1;i<CYCLEMAX;i++) {
+ if ((iterate->allrelerr[iterate->iter_i-(1+i)] == iterate->allrelerr[iterate->iter_i-1]) &&
+ (iterate->allrelerr[iterate->iter_i-(1+i)] == iterate->allrelerr[iterate->iter_i-(1+2*i)])) {
+ incycle = TRUE;
+ if (debug)
+ {
+ fprintf(debug,"Exiting from an NPT iterating cycle of length %d\n",i);
+ }
+ break;
+ }
+ }
+
+ if (incycle) {
+ /* step 1: trapped in a numerical attractor */
+ /* we are trapped in a numerical attractor, and can't converge any more, and are close to the final result.
+ Better to give up convergence here than have the simulation die.
+ */
+ iterate->num_close++;
+ return TRUE;
+ }
+ else
+ {
+ /* Step #2: test if we are reasonably close for other reasons, then monitor the number. If not, die */
+
+ /* how many close calls have we had? If less than a few, we're OK */
+ if (iterate->num_close < MAX_NUMBER_CLOSE)
+ {
+ sprintf(buf,"Slight numerical convergence deviation with NPT at step %d, relative error only %10.5g, likely not a problem, continuing\n",nsteps,relerr);
+ md_print_warning(cr,fplog,buf);
+ iterate->num_close++;
+ return TRUE;
+ /* if more than a few, check the total fraction. If too high, die. */
+ } else if (iterate->num_close/(double)nsteps > FRACTION_CLOSE) {
+ gmx_fatal(FARGS,"Could not converge NPT constraints, too many exceptions (%d%%\n",iterate->num_close/(double)nsteps);
+ }
+ }
+ }
+ else
+ {
+ gmx_fatal(FARGS,"Could not converge NPT constraints\n");
+ }
+ }
+ }
+
+ iterate->xprev = iterate->x;
+ iterate->x = *newf;
+ iterate->fprev = iterate->f;
+ iterate->iter_i++;
+
+ return FALSE;
+}
+
+static void check_nst_param(FILE *fplog,t_commrec *cr,
+ const char *desc_nst,int nst,
+ const char *desc_p,int *p)
+{
+ char buf[STRLEN];
+
+ if (*p > 0 && *p % nst != 0)
+ {
+ /* Round up to the next multiple of nst */
+ *p = ((*p)/nst + 1)*nst;
+ sprintf(buf,"NOTE: %s changes %s to %d\n",desc_nst,desc_p,*p);
+ md_print_warning(cr,fplog,buf);
+ }
+}
+
+static void reset_all_counters(FILE *fplog,t_commrec *cr,
+ gmx_large_int_t step,
+ gmx_large_int_t *step_rel,t_inputrec *ir,
+ gmx_wallcycle_t wcycle,t_nrnb *nrnb,
+ gmx_runtime_t *runtime)
+{
+ char buf[STRLEN],sbuf[STEPSTRSIZE];
+
+ /* Reset all the counters related to performance over the run */
+ sprintf(buf,"Step %s: resetting all time and cycle counters\n",
+ gmx_step_str(step,sbuf));
+ md_print_warning(cr,fplog,buf);
+
+ wallcycle_stop(wcycle,ewcRUN);
+ wallcycle_reset_all(wcycle);
+ if (DOMAINDECOMP(cr))
+ {
+ reset_dd_statistics_counters(cr->dd);
+ }
+ init_nrnb(nrnb);
+ ir->init_step += *step_rel;
+ ir->nsteps -= *step_rel;
+ *step_rel = 0;
+ wallcycle_start(wcycle,ewcRUN);
+ runtime_start(runtime);
+ print_date_and_time(fplog,cr->nodeid,"Restarted time",runtime);
+}
+
+static void min_zero(int *n,int i)
+{
+ if (i > 0 && (*n == 0 || i < *n))
+ {
+ *n = i;
+ }
+}
+
+static int lcd4(int i1,int i2,int i3,int i4)
+{
+ int nst;
+
+ nst = 0;
+ min_zero(&nst,i1);
+ min_zero(&nst,i2);
+ min_zero(&nst,i3);
+ min_zero(&nst,i4);
+ if (nst == 0)
+ {
+ gmx_incons("All 4 inputs for determininig nstglobalcomm are <= 0");
+ }
+
+ while (nst > 1 && ((i1 > 0 && i1 % nst != 0) ||
+ (i2 > 0 && i2 % nst != 0) ||
+ (i3 > 0 && i3 % nst != 0) ||
+ (i4 > 0 && i4 % nst != 0)))
+ {
+ nst--;
+ }
+
+ return nst;
+}
+
+static int check_nstglobalcomm(FILE *fplog,t_commrec *cr,
+ int nstglobalcomm,t_inputrec *ir)
+{
+ char buf[STRLEN];
+
+ if (!EI_DYNAMICS(ir->eI))
+ {
+ nstglobalcomm = 1;
+ }
+
+ if (nstglobalcomm == -1)
+ {
+ if (!(ir->nstcalcenergy > 0 ||
+ ir->nstlist > 0 ||
+ ir->etc != etcNO ||
+ ir->epc != epcNO))
+ {
+ nstglobalcomm = 10;
+ if (ir->nstenergy > 0 && ir->nstenergy < nstglobalcomm)
+ {
+ nstglobalcomm = ir->nstenergy;
+ }
+ }
+ else
+ {
+ /* Ensure that we do timely global communication for
+ * (possibly) each of the four following options.
+ */
+ nstglobalcomm = lcd4(ir->nstcalcenergy,
+ ir->nstlist,
+ ir->etc != etcNO ? ir->nsttcouple : 0,
+ ir->epc != epcNO ? ir->nstpcouple : 0);
+ }
+ }
+ else
+ {
+ if (ir->nstlist > 0 &&
+ nstglobalcomm > ir->nstlist && nstglobalcomm % ir->nstlist != 0)
+ {
+ nstglobalcomm = (nstglobalcomm / ir->nstlist)*ir->nstlist;
+ sprintf(buf,"WARNING: nstglobalcomm is larger than nstlist, but not a multiple, setting it to %d\n",nstglobalcomm);
+ md_print_warning(cr,fplog,buf);
+ }
+ if (ir->nstcalcenergy > 0)
+ {
+ check_nst_param(fplog,cr,"-gcom",nstglobalcomm,
+ "nstcalcenergy",&ir->nstcalcenergy);
+ }
+ if (ir->etc != etcNO && ir->nsttcouple > 0)
+ {
+ check_nst_param(fplog,cr,"-gcom",nstglobalcomm,
+ "nsttcouple",&ir->nsttcouple);
+ }
+ if (ir->epc != epcNO && ir->nstpcouple > 0)
+ {
+ check_nst_param(fplog,cr,"-gcom",nstglobalcomm,
+ "nstpcouple",&ir->nstpcouple);
+ }
+
+ check_nst_param(fplog,cr,"-gcom",nstglobalcomm,
+ "nstenergy",&ir->nstenergy);
+
+ check_nst_param(fplog,cr,"-gcom",nstglobalcomm,
+ "nstlog",&ir->nstlog);
+ }
+
+ if (ir->comm_mode != ecmNO && ir->nstcomm < nstglobalcomm)
+ {
+ sprintf(buf,"WARNING: Changing nstcomm from %d to %d\n",
+ ir->nstcomm,nstglobalcomm);
+ md_print_warning(cr,fplog,buf);
+ ir->nstcomm = nstglobalcomm;
+ }
+
+ return nstglobalcomm;
+}
+
+void check_ir_old_tpx_versions(t_commrec *cr,FILE *fplog,
+ t_inputrec *ir,gmx_mtop_t *mtop)
+{
+ /* Check required for old tpx files */
+ if (IR_TWINRANGE(*ir) && ir->nstlist > 1 &&
+ ir->nstcalcenergy % ir->nstlist != 0)
+ {
+ md_print_warning(cr,fplog,"Old tpr file with twin-range settings: modifying energy calculation and/or T/P-coupling frequencies");
+
+ if (gmx_mtop_ftype_count(mtop,F_CONSTR) +
+ gmx_mtop_ftype_count(mtop,F_CONSTRNC) > 0 &&
+ ir->eConstrAlg == econtSHAKE)
+ {
+ md_print_warning(cr,fplog,"With twin-range cut-off's and SHAKE the virial and pressure are incorrect");
+ if (ir->epc != epcNO)
+ {
+ gmx_fatal(FARGS,"Can not do pressure coupling with twin-range cut-off's and SHAKE");
+ }
+ }
+ check_nst_param(fplog,cr,"nstlist",ir->nstlist,
+ "nstcalcenergy",&ir->nstcalcenergy);
+ if (ir->epc != epcNO)
+ {
+ check_nst_param(fplog,cr,"nstlist",ir->nstlist,
+ "nstpcouple",&ir->nstpcouple);
+ }
+ check_nst_param(fplog,cr,"nstcalcenergy",ir->nstcalcenergy,
+ "nstenergy",&ir->nstenergy);
+ check_nst_param(fplog,cr,"nstcalcenergy",ir->nstcalcenergy,
+ "nstlog",&ir->nstlog);
+ if (ir->efep != efepNO)
+ {
+ check_nst_param(fplog,cr,"nstcalcenergy",ir->nstcalcenergy,
+ "nstdhdl",&ir->nstdhdl);
+ }
+ }
+}
+
+typedef struct {
+ gmx_bool bGStatEveryStep;
+ gmx_large_int_t step_ns;
+ gmx_large_int_t step_nscheck;
+ gmx_large_int_t nns;
+ matrix scale_tot;
+ int nabnsb;
+ double s1;
+ double s2;
+ double ab;
+ double lt_runav;
+ double lt_runav2;
+} gmx_nlheur_t;
+
+static void reset_nlistheuristics(gmx_nlheur_t *nlh,gmx_large_int_t step)
+{
+ nlh->lt_runav = 0;
+ nlh->lt_runav2 = 0;
+ nlh->step_nscheck = step;
+}
+
+static void init_nlistheuristics(gmx_nlheur_t *nlh,
+ gmx_bool bGStatEveryStep,gmx_large_int_t step)
+{
+ nlh->bGStatEveryStep = bGStatEveryStep;
+ nlh->nns = 0;
+ nlh->nabnsb = 0;
+ nlh->s1 = 0;
+ nlh->s2 = 0;
+ nlh->ab = 0;
+
+ reset_nlistheuristics(nlh,step);
+}
+
+static void update_nliststatistics(gmx_nlheur_t *nlh,gmx_large_int_t step)
+{
+ gmx_large_int_t nl_lt;
+ char sbuf[STEPSTRSIZE],sbuf2[STEPSTRSIZE];
+
+ /* Determine the neighbor list life time */
+ nl_lt = step - nlh->step_ns;
+ if (debug)
+ {
+ fprintf(debug,"%d atoms beyond ns buffer, updating neighbor list after %s steps\n",nlh->nabnsb,gmx_step_str(nl_lt,sbuf));
+ }
+ nlh->nns++;
+ nlh->s1 += nl_lt;
+ nlh->s2 += nl_lt*nl_lt;
+ nlh->ab += nlh->nabnsb;
+ if (nlh->lt_runav == 0)
+ {
+ nlh->lt_runav = nl_lt;
+ /* Initialize the fluctuation average
+ * such that at startup we check after 0 steps.
+ */
+ nlh->lt_runav2 = sqr(nl_lt/2.0);
+ }
+ /* Running average with 0.9 gives an exp. history of 9.5 */
+ nlh->lt_runav2 = 0.9*nlh->lt_runav2 + 0.1*sqr(nlh->lt_runav - nl_lt);
+ nlh->lt_runav = 0.9*nlh->lt_runav + 0.1*nl_lt;
+ if (nlh->bGStatEveryStep)
+ {
+ /* Always check the nlist validity */
+ nlh->step_nscheck = step;
+ }
+ else
+ {
+ /* We check after: <life time> - 2*sigma
+ * The factor 2 is quite conservative,
+ * but we assume that with nstlist=-1 the user
+ * prefers exact integration over performance.
+ */
+ nlh->step_nscheck = step
+ + (int)(nlh->lt_runav - 2.0*sqrt(nlh->lt_runav2)) - 1;
+ }
+ if (debug)
+ {
+ fprintf(debug,"nlist life time %s run av. %4.1f sig %3.1f check %s check with -gcom %d\n",
+ gmx_step_str(nl_lt,sbuf),nlh->lt_runav,sqrt(nlh->lt_runav2),
+ gmx_step_str(nlh->step_nscheck-step+1,sbuf2),
+ (int)(nlh->lt_runav - 2.0*sqrt(nlh->lt_runav2)));
+ }
+}
+
+static void set_nlistheuristics(gmx_nlheur_t *nlh,gmx_bool bReset,gmx_large_int_t step)
+{
+ int d;
+
+ if (bReset)
+ {
+ reset_nlistheuristics(nlh,step);
+ }
+ else
+ {
+ update_nliststatistics(nlh,step);
+ }
+
+ nlh->step_ns = step;
+ /* Initialize the cumulative coordinate scaling matrix */
+ clear_mat(nlh->scale_tot);
+ for(d=0; d<DIM; d++)
+ {
+ nlh->scale_tot[d][d] = 1.0;
+ }
+}
+
+static void rerun_parallel_comm(t_commrec *cr,t_trxframe *fr,
+ gmx_bool *bNotLastFrame)
+{
+ gmx_bool bAlloc;
+ rvec *xp,*vp;
+
+ bAlloc = (fr->natoms == 0);
+
+ if (MASTER(cr) && !*bNotLastFrame)
+ {
+ fr->natoms = -1;
+ }
+ xp = fr->x;
+ vp = fr->v;
+ gmx_bcast(sizeof(*fr),fr,cr);
+ fr->x = xp;
+ fr->v = vp;
+
+ *bNotLastFrame = (fr->natoms >= 0);
+
+ if (*bNotLastFrame && PARTDECOMP(cr))
+ {
+ /* x and v are the only variable size quantities stored in trr
+ * that are required for rerun (f is not needed).
+ */
+ if (bAlloc)
+ {
+ snew(fr->x,fr->natoms);
+ snew(fr->v,fr->natoms);
+ }
+ if (fr->bX)
+ {
+ gmx_bcast(fr->natoms*sizeof(fr->x[0]),fr->x[0],cr);
+ }
+ if (fr->bV)
+ {
+ gmx_bcast(fr->natoms*sizeof(fr->v[0]),fr->v[0],cr);
+ }
+ }
+}
+
+double do_md(FILE *fplog,t_commrec *cr,int nfile,const t_filenm fnm[],
+ const output_env_t oenv, gmx_bool bVerbose,gmx_bool bCompact,
+ int nstglobalcomm,
+ gmx_vsite_t *vsite,gmx_constr_t constr,
+ int stepout,t_inputrec *ir,
+ gmx_mtop_t *top_global,
+ t_fcdata *fcd,
+ t_state *state_global,
+ t_mdatoms *mdatoms,
+ t_nrnb *nrnb,gmx_wallcycle_t wcycle,
+ gmx_edsam_t ed,t_forcerec *fr,
+ int repl_ex_nst,int repl_ex_seed,gmx_membed_t *membed,
+ real cpt_period,real max_hours,
+ const char *deviceOptions,
+ unsigned long Flags,
+ gmx_runtime_t *runtime)
+{
+ gmx_mdoutf_t *outf;
+ gmx_large_int_t step,step_rel;
+ double run_time;
+ double t,t0,lam0;
+ gmx_bool bGStatEveryStep,bGStat,bNstEner,bCalcEnerPres;
+ gmx_bool bNS,bNStList,bSimAnn,bStopCM,bRerunMD,bNotLastFrame=FALSE,
+ bFirstStep,bStateFromTPX,bInitStep,bLastStep,
+ bBornRadii,bStartingFromCpt;
+ gmx_bool bDoDHDL=FALSE;
+ gmx_bool do_ene,do_log,do_verbose,bRerunWarnNoV=TRUE,
+ bForceUpdate=FALSE,bCPT;
+ int mdof_flags;
+ gmx_bool bMasterState;
+ int force_flags,cglo_flags;
+ tensor force_vir,shake_vir,total_vir,tmp_vir,pres;
+ int i,m;
+ t_trxstatus *status;
+ rvec mu_tot;
+ t_vcm *vcm;
+ t_state *bufstate=NULL;
+ matrix *scale_tot,pcoupl_mu,M,ebox;
+ gmx_nlheur_t nlh;
+ t_trxframe rerun_fr;
+ gmx_repl_ex_t repl_ex=NULL;
+ int nchkpt=1;
+
+ gmx_localtop_t *top;
+ t_mdebin *mdebin=NULL;
+ t_state *state=NULL;
+ rvec *f_global=NULL;
+ int n_xtc=-1;
+ rvec *x_xtc=NULL;
+ gmx_enerdata_t *enerd;
+ rvec *f=NULL;
+ gmx_global_stat_t gstat;
+ gmx_update_t upd=NULL;
+ t_graph *graph=NULL;
+ globsig_t gs;
+
+ gmx_bool bFFscan;
+ gmx_groups_t *groups;
+ gmx_ekindata_t *ekind, *ekind_save;
+ gmx_shellfc_t shellfc;
+ int count,nconverged=0;
+ real timestep=0;
+ double tcount=0;
+ gmx_bool bIonize=FALSE;
+ gmx_bool bTCR=FALSE,bConverged=TRUE,bOK,bSumEkinhOld,bExchanged;
+ gmx_bool bAppend;
+ gmx_bool bResetCountersHalfMaxH=FALSE;
+ gmx_bool bVV,bIterations,bFirstIterate,bTemp,bPres,bTrotter;
+ real temp0,mu_aver=0,dvdl;
+ int a0,a1,gnx=0,ii;
+ atom_id *grpindex=NULL;
+ char *grpname;
+ t_coupl_rec *tcr=NULL;
+ rvec *xcopy=NULL,*vcopy=NULL,*cbuf=NULL;
+ matrix boxcopy={{0}},lastbox;
+ tensor tmpvir;
+ real fom,oldfom,veta_save,pcurr,scalevir,tracevir;
+ real vetanew = 0;
+ double cycles;
+ real saved_conserved_quantity = 0;
+ real last_ekin = 0;
+ int iter_i;
+ t_extmass MassQ;
+ int **trotter_seq;
+ char sbuf[STEPSTRSIZE],sbuf2[STEPSTRSIZE];
+ int handled_stop_condition=gmx_stop_cond_none; /* compare to get_stop_condition*/
+ gmx_iterate_t iterate;
+#ifdef GMX_FAHCORE
+ /* Temporary addition for FAHCORE checkpointing */
+ int chkpt_ret;
+#endif
+
+ /* Check for special mdrun options */
+ bRerunMD = (Flags & MD_RERUN);
+ bIonize = (Flags & MD_IONIZE);
+ bFFscan = (Flags & MD_FFSCAN);
+ bAppend = (Flags & MD_APPENDFILES);
+ if (Flags & MD_RESETCOUNTERSHALFWAY)
+ {
+ if (ir->nsteps > 0)
+ {
+ /* Signal to reset the counters half the simulation steps. */
+ wcycle_set_reset_counters(wcycle,ir->nsteps/2);
+ }
+ /* Signal to reset the counters halfway the simulation time. */
+ bResetCountersHalfMaxH = (max_hours > 0);
+ }
+
+ /* md-vv uses averaged full step velocities for T-control
+ md-vv-avek uses averaged half step velocities for T-control (but full step ekin for P control)
+ md uses averaged half step kinetic energies to determine temperature unless defined otherwise by GMX_EKIN_AVE_VEL; */
+ bVV = EI_VV(ir->eI);
+ if (bVV) /* to store the initial velocities while computing virial */
+ {
+ snew(cbuf,top_global->natoms);
+ }
+ /* all the iteratative cases - only if there are constraints */
+ bIterations = ((IR_NPT_TROTTER(ir)) && (constr) && (!bRerunMD));
+ bTrotter = (bVV && (IR_NPT_TROTTER(ir) || (IR_NVT_TROTTER(ir))));
+
+ if (bRerunMD)
+ {
+ /* Since we don't know if the frames read are related in any way,
+ * rebuild the neighborlist at every step.
+ */
+ ir->nstlist = 1;
+ ir->nstcalcenergy = 1;
+ nstglobalcomm = 1;
+ }
+
+ check_ir_old_tpx_versions(cr,fplog,ir,top_global);
+
+ nstglobalcomm = check_nstglobalcomm(fplog,cr,nstglobalcomm,ir);
+ bGStatEveryStep = (nstglobalcomm == 1);
+
+ if (!bGStatEveryStep && ir->nstlist == -1 && fplog != NULL)
+ {
+ fprintf(fplog,
+ "To reduce the energy communication with nstlist = -1\n"
+ "the neighbor list validity should not be checked at every step,\n"
+ "this means that exact integration is not guaranteed.\n"
+ "The neighbor list validity is checked after:\n"
+ " <n.list life time> - 2*std.dev.(n.list life time) steps.\n"
+ "In most cases this will result in exact integration.\n"
+ "This reduces the energy communication by a factor of 2 to 3.\n"
+ "If you want less energy communication, set nstlist > 3.\n\n");
+ }
+
+ if (bRerunMD || bFFscan)
+ {
+ ir->nstxtcout = 0;
+ }
+ groups = &top_global->groups;
+
+ /* Initial values */
+ init_md(fplog,cr,ir,oenv,&t,&t0,&state_global->lambda,&lam0,
+ nrnb,top_global,&upd,
+ nfile,fnm,&outf,&mdebin,
+ force_vir,shake_vir,mu_tot,&bSimAnn,&vcm,state_global,Flags);
+
+ clear_mat(total_vir);
+ clear_mat(pres);
+ /* Energy terms and groups */
+ snew(enerd,1);
+ init_enerdata(top_global->groups.grps[egcENER].nr,ir->n_flambda,enerd);
+ if (DOMAINDECOMP(cr))
+ {
+ f = NULL;
+ }
+ else
+ {
+ snew(f,top_global->natoms);
+ }
+
+ /* Kinetic energy data */
+ snew(ekind,1);
+ init_ekindata(fplog,top_global,&(ir->opts),ekind);
+ /* needed for iteration of constraints */
+ snew(ekind_save,1);
+ init_ekindata(fplog,top_global,&(ir->opts),ekind_save);
+ /* Copy the cos acceleration to the groups struct */
+ ekind->cosacc.cos_accel = ir->cos_accel;
+
+ gstat = global_stat_init(ir);
+ debug_gmx();
+
+ /* Check for polarizable models and flexible constraints */
+ shellfc = init_shell_flexcon(fplog,
+ top_global,n_flexible_constraints(constr),
+ (ir->bContinuation ||
+ (DOMAINDECOMP(cr) && !MASTER(cr))) ?
+ NULL : state_global->x);
+
+ if (DEFORM(*ir))
+ {
+#ifdef GMX_THREADS
+ tMPI_Thread_mutex_lock(&deform_init_box_mutex);
+#endif
+ set_deform_reference_box(upd,
+ deform_init_init_step_tpx,
+ deform_init_box_tpx);
+#ifdef GMX_THREADS
+ tMPI_Thread_mutex_unlock(&deform_init_box_mutex);
+#endif
+ }
+
+ {
+ double io = compute_io(ir,top_global->natoms,groups,mdebin->ebin->nener,1);
+ if ((io > 2000) && MASTER(cr))
+ fprintf(stderr,
+ "\nWARNING: This run will generate roughly %.0f Mb of data\n\n",
+ io);
+ }
+
+ if (DOMAINDECOMP(cr)) {
+ top = dd_init_local_top(top_global);
+
+ snew(state,1);
+ dd_init_local_state(cr->dd,state_global,state);
+
+ if (DDMASTER(cr->dd) && ir->nstfout) {
+ snew(f_global,state_global->natoms);
+ }
+ } else {
+ if (PAR(cr)) {
+ /* Initialize the particle decomposition and split the topology */
+ top = split_system(fplog,top_global,ir,cr);
+
+ pd_cg_range(cr,&fr->cg0,&fr->hcg);
+ pd_at_range(cr,&a0,&a1);
+ } else {
+ top = gmx_mtop_generate_local_top(top_global,ir);
+
+ a0 = 0;
+ a1 = top_global->natoms;
+ }
+
+ state = partdec_init_local_state(cr,state_global);
+ f_global = f;
+
+ atoms2md(top_global,ir,0,NULL,a0,a1-a0,mdatoms);
+
+ if (vsite) {
+ set_vsite_top(vsite,top,mdatoms,cr);
+ }
+
+ if (ir->ePBC != epbcNONE && !ir->bPeriodicMols) {
+ graph = mk_graph(fplog,&(top->idef),0,top_global->natoms,FALSE,FALSE);
+ }
+
+ if (shellfc) {
+ make_local_shells(cr,mdatoms,shellfc);
+ }
+
+ if (ir->pull && PAR(cr)) {
+ dd_make_local_pull_groups(NULL,ir->pull,mdatoms);
+ }
+ }
+
+ if (DOMAINDECOMP(cr))
+ {
+ /* Distribute the charge groups over the nodes from the master node */
+ dd_partition_system(fplog,ir->init_step,cr,TRUE,1,
+ state_global,top_global,ir,
+ state,&f,mdatoms,top,fr,
+ vsite,shellfc,constr,
+ nrnb,wcycle,FALSE);
+ }
+
+ update_mdatoms(mdatoms,state->lambda);
+
+ if (MASTER(cr))
+ {
+ if (opt2bSet("-cpi",nfile,fnm))
+ {
+ /* Update mdebin with energy history if appending to output files */
+ if ( Flags & MD_APPENDFILES )
+ {
+ restore_energyhistory_from_state(mdebin,&state_global->enerhist);
+ }
+ else
+ {
+ /* We might have read an energy history from checkpoint,
+ * free the allocated memory and reset the counts.
+ */
+ done_energyhistory(&state_global->enerhist);
+ init_energyhistory(&state_global->enerhist);
+ }
+ }
+ /* Set the initial energy history in state by updating once */
+ update_energyhistory(&state_global->enerhist,mdebin);
+ }
+
+ if ((state->flags & (1<<estLD_RNG)) && (Flags & MD_READ_RNG)) {
+ /* Set the random state if we read a checkpoint file */
+ set_stochd_state(upd,state);
+ }
+
+ /* Initialize constraints */
+ if (constr) {
+ if (!DOMAINDECOMP(cr))
+ set_constraints(constr,top,ir,mdatoms,cr);
+ }
+
+ /* Check whether we have to GCT stuff */
+ bTCR = ftp2bSet(efGCT,nfile,fnm);
+ if (bTCR) {
+ if (MASTER(cr)) {
+ fprintf(stderr,"Will do General Coupling Theory!\n");
+ }
+ gnx = top_global->mols.nr;
+ snew(grpindex,gnx);
+ for(i=0; (i<gnx); i++) {
+ grpindex[i] = i;
+ }
+ }
+
+ if (repl_ex_nst > 0 && MASTER(cr))
+ repl_ex = init_replica_exchange(fplog,cr->ms,state_global,ir,
+ repl_ex_nst,repl_ex_seed);
+
+ if (!ir->bContinuation && !bRerunMD)
+ {
+ if (mdatoms->cFREEZE && (state->flags & (1<<estV)))
+ {
+ /* Set the velocities of frozen particles to zero */
+ for(i=mdatoms->start; i<mdatoms->start+mdatoms->homenr; i++)
+ {
+ for(m=0; m<DIM; m++)
+ {
+ if (ir->opts.nFreeze[mdatoms->cFREEZE[i]][m])
+ {
+ state->v[i][m] = 0;
+ }
+ }
+ }
+ }
+
+ if (constr)
+ {
+ /* Constrain the initial coordinates and velocities */
+ do_constrain_first(fplog,constr,ir,mdatoms,state,f,
+ graph,cr,nrnb,fr,top,shake_vir);
+ }
+ if (vsite)
+ {
+ /* Construct the virtual sites for the initial configuration */
+ construct_vsites(fplog,vsite,state->x,nrnb,ir->delta_t,NULL,
+ top->idef.iparams,top->idef.il,
+ fr->ePBC,fr->bMolPBC,graph,cr,state->box);
+ }
+ }
+
+ debug_gmx();
+
+ /* I'm assuming we need global communication the first time! MRS */
+ cglo_flags = (CGLO_TEMPERATURE | CGLO_GSTAT
+ | (bVV ? CGLO_PRESSURE:0)
+ | (bVV ? CGLO_CONSTRAINT:0)
+ | (bRerunMD ? CGLO_RERUNMD:0)
+ | ((Flags & MD_READ_EKIN) ? CGLO_READEKIN:0));
+
+ bSumEkinhOld = FALSE;
+ compute_globals(fplog,gstat,cr,ir,fr,ekind,state,state_global,mdatoms,nrnb,vcm,
+ wcycle,enerd,force_vir,shake_vir,total_vir,pres,mu_tot,
+ constr,NULL,FALSE,state->box,
+ top_global,&pcurr,top_global->natoms,&bSumEkinhOld,cglo_flags);
+ if (ir->eI == eiVVAK) {
+ /* a second call to get the half step temperature initialized as well */
+ /* we do the same call as above, but turn the pressure off -- internally to
+ compute_globals, this is recognized as a velocity verlet half-step
+ kinetic energy calculation. This minimized excess variables, but
+ perhaps loses some logic?*/
+
+ compute_globals(fplog,gstat,cr,ir,fr,ekind,state,state_global,mdatoms,nrnb,vcm,
+ wcycle,enerd,force_vir,shake_vir,total_vir,pres,mu_tot,
+ constr,NULL,FALSE,state->box,
+ top_global,&pcurr,top_global->natoms,&bSumEkinhOld,
+ cglo_flags &~ CGLO_PRESSURE);
+ }
+
+ /* Calculate the initial half step temperature, and save the ekinh_old */
+ if (!(Flags & MD_STARTFROMCPT))
+ {
+ for(i=0; (i<ir->opts.ngtc); i++)
+ {
+ copy_mat(ekind->tcstat[i].ekinh,ekind->tcstat[i].ekinh_old);
+ }
+ }
+ if (ir->eI != eiVV)
+ {
+ enerd->term[F_TEMP] *= 2; /* result of averages being done over previous and current step,
+ and there is no previous step */
+ }
+ temp0 = enerd->term[F_TEMP];
+
+ /* if using an iterative algorithm, we need to create a working directory for the state. */
+ if (bIterations)
+ {
+ bufstate = init_bufstate(state);
+ }
+ if (bFFscan)
+ {
+ snew(xcopy,state->natoms);
+ snew(vcopy,state->natoms);
+ copy_rvecn(state->x,xcopy,0,state->natoms);
+ copy_rvecn(state->v,vcopy,0,state->natoms);
+ copy_mat(state->box,boxcopy);
+ }
+
+ /* need to make an initiation call to get the Trotter variables set, as well as other constants for non-trotter
+ temperature control */
+ trotter_seq = init_npt_vars(ir,state,&MassQ,bTrotter);
+
+ if (MASTER(cr))
+ {
+ if (constr && !ir->bContinuation && ir->eConstrAlg == econtLINCS)
+ {
+ fprintf(fplog,
+ "RMS relative constraint deviation after constraining: %.2e\n",
+ constr_rmsd(constr,FALSE));
+ }
+ fprintf(fplog,"Initial temperature: %g K\n",enerd->term[F_TEMP]);
+ if (bRerunMD)
+ {
+ fprintf(stderr,"starting md rerun '%s', reading coordinates from"
+ " input trajectory '%s'\n\n",
+ *(top_global->name),opt2fn("-rerun",nfile,fnm));
+ if (bVerbose)
+ {
+ fprintf(stderr,"Calculated time to finish depends on nsteps from "
+ "run input file,\nwhich may not correspond to the time "
+ "needed to process input trajectory.\n\n");
+ }
+ }
+ else
+ {
+ char tbuf[20];
+ fprintf(stderr,"starting mdrun '%s'\n",
+ *(top_global->name));
+ if (ir->nsteps >= 0)
+ {
+ sprintf(tbuf,"%8.1f",(ir->init_step+ir->nsteps)*ir->delta_t);
+ }
+ else
+ {
+ sprintf(tbuf,"%s","infinite");
+ }
+ if (ir->init_step > 0)
+ {
+ fprintf(stderr,"%s steps, %s ps (continuing from step %s, %8.1f ps).\n",
+ gmx_step_str(ir->init_step+ir->nsteps,sbuf),tbuf,
+ gmx_step_str(ir->init_step,sbuf2),
+ ir->init_step*ir->delta_t);
+ }
+ else
+ {
+ fprintf(stderr,"%s steps, %s ps.\n",
+ gmx_step_str(ir->nsteps,sbuf),tbuf);
+ }
+ }
+ fprintf(fplog,"\n");
+ }
+
+ /* Set and write start time */
+ runtime_start(runtime);
+ print_date_and_time(fplog,cr->nodeid,"Started mdrun",runtime);
+ wallcycle_start(wcycle,ewcRUN);
+ if (fplog)
+ fprintf(fplog,"\n");
+
+ /* safest point to do file checkpointing is here. More general point would be immediately before integrator call */
+#ifdef GMX_FAHCORE
+ chkpt_ret=fcCheckPointParallel( cr->nodeid,
+ NULL,0);
+ if ( chkpt_ret == 0 )
+ gmx_fatal( 3,__FILE__,__LINE__, "Checkpoint error on step %d\n", 0 );
+#endif
+
+ debug_gmx();
+ /***********************************************************
+ *
+ * Loop over MD steps
+ *
+ ************************************************************/
+
+ /* if rerunMD then read coordinates and velocities from input trajectory */
+ if (bRerunMD)
+ {
+ if (getenv("GMX_FORCE_UPDATE"))
+ {
+ bForceUpdate = TRUE;
+ }
+
+ rerun_fr.natoms = 0;
+ if (MASTER(cr))
+ {
+ bNotLastFrame = read_first_frame(oenv,&status,
+ opt2fn("-rerun",nfile,fnm),
+ &rerun_fr,TRX_NEED_X | TRX_READ_V);
+ if (rerun_fr.natoms != top_global->natoms)
+ {
+ gmx_fatal(FARGS,
+ "Number of atoms in trajectory (%d) does not match the "
+ "run input file (%d)\n",
+ rerun_fr.natoms,top_global->natoms);
+ }
+ if (ir->ePBC != epbcNONE)
+ {
+ if (!rerun_fr.bBox)
+ {
+ gmx_fatal(FARGS,"Rerun trajectory frame step %d time %f does not contain a box, while pbc is used",rerun_fr.step,rerun_fr.time);
+ }
+ if (max_cutoff2(ir->ePBC,rerun_fr.box) < sqr(fr->rlistlong))
+ {
+ gmx_fatal(FARGS,"Rerun trajectory frame step %d time %f has too small box dimensions",rerun_fr.step,rerun_fr.time);
+ }
+ }
+ }
+
+ if (PAR(cr))
+ {
+ rerun_parallel_comm(cr,&rerun_fr,&bNotLastFrame);
+ }
+
+ if (ir->ePBC != epbcNONE)
+ {
+ /* Set the shift vectors.
+ * Necessary here when have a static box different from the tpr box.
+ */
+ calc_shifts(rerun_fr.box,fr->shift_vec);
+ }
+ }
+
+ /* loop over MD steps or if rerunMD to end of input trajectory */
+ bFirstStep = TRUE;
+ /* Skip the first Nose-Hoover integration when we get the state from tpx */
+ bStateFromTPX = !opt2bSet("-cpi",nfile,fnm);
+ bInitStep = bFirstStep && (bStateFromTPX || bVV);
+ bStartingFromCpt = (Flags & MD_STARTFROMCPT) && bInitStep;
+ bLastStep = FALSE;
+ bSumEkinhOld = FALSE;
+ bExchanged = FALSE;
+
+ init_global_signals(&gs,cr,ir,repl_ex_nst);
+
+ step = ir->init_step;
+ step_rel = 0;
+
+ if (ir->nstlist == -1)
+ {
+ init_nlistheuristics(&nlh,bGStatEveryStep,step);
+ }
+
+ bLastStep = (bRerunMD || (ir->nsteps >= 0 && step_rel > ir->nsteps));
+ while (!bLastStep || (bRerunMD && bNotLastFrame)) {
+
+ wallcycle_start(wcycle,ewcSTEP);
+
+ GMX_MPE_LOG(ev_timestep1);
+
+ if (bRerunMD) {
+ if (rerun_fr.bStep) {
+ step = rerun_fr.step;
+ step_rel = step - ir->init_step;
+ }
+ if (rerun_fr.bTime) {
+ t = rerun_fr.time;
+ }
+ else
+ {
+ t = step;
+ }
+ }
+ else
+ {
+ bLastStep = (step_rel == ir->nsteps);
+ t = t0 + step*ir->delta_t;
+ }
+
+ if (ir->efep != efepNO)
+ {
+ if (bRerunMD && rerun_fr.bLambda && (ir->delta_lambda!=0))
+ {
+ state_global->lambda = rerun_fr.lambda;
+ }
+ else
+ {
+ state_global->lambda = lam0 + step*ir->delta_lambda;
+ }
+ state->lambda = state_global->lambda;
+ bDoDHDL = do_per_step(step,ir->nstdhdl);
+ }
+
+ if (bSimAnn)
+ {
+ update_annealing_target_temp(&(ir->opts),t);
+ }
+
+ if (bRerunMD)
+ {
+ if (!(DOMAINDECOMP(cr) && !MASTER(cr)))
+ {
+ for(i=0; i<state_global->natoms; i++)
+ {
+ copy_rvec(rerun_fr.x[i],state_global->x[i]);
+ }
+ if (rerun_fr.bV)
+ {
+ for(i=0; i<state_global->natoms; i++)
+ {
+ copy_rvec(rerun_fr.v[i],state_global->v[i]);
+ }
+ }
+ else
+ {
+ for(i=0; i<state_global->natoms; i++)
+ {
+ clear_rvec(state_global->v[i]);
+ }
+ if (bRerunWarnNoV)
+ {
+ fprintf(stderr,"\nWARNING: Some frames do not contain velocities.\n"
+ " Ekin, temperature and pressure are incorrect,\n"
+ " the virial will be incorrect when constraints are present.\n"
+ "\n");
+ bRerunWarnNoV = FALSE;
+ }
+ }
+ }
+ copy_mat(rerun_fr.box,state_global->box);
+ copy_mat(state_global->box,state->box);
+
+ if (vsite && (Flags & MD_RERUN_VSITE))
+ {
+ if (DOMAINDECOMP(cr))
+ {
+ gmx_fatal(FARGS,"Vsite recalculation with -rerun is not implemented for domain decomposition, use particle decomposition");
+ }
+ if (graph)
+ {
+ /* Following is necessary because the graph may get out of sync
+ * with the coordinates if we only have every N'th coordinate set
+ */
+ mk_mshift(fplog,graph,fr->ePBC,state->box,state->x);
+ shift_self(graph,state->box,state->x);
+ }
+ construct_vsites(fplog,vsite,state->x,nrnb,ir->delta_t,state->v,
+ top->idef.iparams,top->idef.il,
+ fr->ePBC,fr->bMolPBC,graph,cr,state->box);
+ if (graph)
+ {
+ unshift_self(graph,state->box,state->x);
+ }
+ }
+ }
+
+ /* Stop Center of Mass motion */
+ bStopCM = (ir->comm_mode != ecmNO && do_per_step(step,ir->nstcomm));
+
+ /* Copy back starting coordinates in case we're doing a forcefield scan */
+ if (bFFscan)
+ {
+ for(ii=0; (ii<state->natoms); ii++)
+ {
+ copy_rvec(xcopy[ii],state->x[ii]);
+ copy_rvec(vcopy[ii],state->v[ii]);
+ }
+ copy_mat(boxcopy,state->box);
+ }
+
+ if (bRerunMD)
+ {
+ /* for rerun MD always do Neighbour Searching */
+ bNS = (bFirstStep || ir->nstlist != 0);
+ bNStList = bNS;
+ }
+ else
+ {
+ /* Determine whether or not to do Neighbour Searching and LR */
+ bNStList = (ir->nstlist > 0 && step % ir->nstlist == 0);
+
+ bNS = (bFirstStep || bExchanged || bNStList ||
+ (ir->nstlist == -1 && nlh.nabnsb > 0));
+
+ if (bNS && ir->nstlist == -1)
+ {
+ set_nlistheuristics(&nlh,bFirstStep || bExchanged,step);
+ }
+ }
+
+ /* < 0 means stop at next step, > 0 means stop at next NS step */
+ if ( (gs.set[eglsSTOPCOND] < 0 ) ||
+ ( (gs.set[eglsSTOPCOND] > 0 ) && ( bNS || ir->nstlist==0)) )
+ {
+ bLastStep = TRUE;
+ }
+
+ /* Determine whether or not to update the Born radii if doing GB */
+ bBornRadii=bFirstStep;
+ if (ir->implicit_solvent && (step % ir->nstgbradii==0))
+ {
+ bBornRadii=TRUE;
+ }
+
+ do_log = do_per_step(step,ir->nstlog) || bFirstStep || bLastStep;
+ do_verbose = bVerbose &&
+ (step % stepout == 0 || bFirstStep || bLastStep);
+
+ if (bNS && !(bFirstStep && ir->bContinuation && !bRerunMD))
+ {
+ if (bRerunMD)
+ {
+ bMasterState = TRUE;
+ }
+ else
+ {
+ bMasterState = FALSE;
+ /* Correct the new box if it is too skewed */
+ if (DYNAMIC_BOX(*ir))
+ {
+ if (correct_box(fplog,step,state->box,graph))
+ {
+ bMasterState = TRUE;
+ }
+ }
+ if (DOMAINDECOMP(cr) && bMasterState)
+ {
+ dd_collect_state(cr->dd,state,state_global);
+ }
+ }
+
+ if (DOMAINDECOMP(cr))
+ {
+ /* Repartition the domain decomposition */
+ wallcycle_start(wcycle,ewcDOMDEC);
+ dd_partition_system(fplog,step,cr,
+ bMasterState,nstglobalcomm,
+ state_global,top_global,ir,
+ state,&f,mdatoms,top,fr,
+ vsite,shellfc,constr,
+ nrnb,wcycle,do_verbose);
+ wallcycle_stop(wcycle,ewcDOMDEC);
+ /* If using an iterative integrator, reallocate space to match the decomposition */
+ }
+ }
+
+ if (MASTER(cr) && do_log && !bFFscan)
+ {
+ print_ebin_header(fplog,step,t,state->lambda);
+ }
+
+ if (ir->efep != efepNO)
+ {
+ update_mdatoms(mdatoms,state->lambda);
+ }
+
+ if (bRerunMD && rerun_fr.bV)
+ {
+
+ /* We need the kinetic energy at minus the half step for determining
+ * the full step kinetic energy and possibly for T-coupling.*/
+ /* This may not be quite working correctly yet . . . . */
+ compute_globals(fplog,gstat,cr,ir,fr,ekind,state,state_global,mdatoms,nrnb,vcm,
+ wcycle,enerd,NULL,NULL,NULL,NULL,mu_tot,
+ constr,NULL,FALSE,state->box,
+ top_global,&pcurr,top_global->natoms,&bSumEkinhOld,
+ CGLO_RERUNMD | CGLO_GSTAT | CGLO_TEMPERATURE);
+ }
+ clear_mat(force_vir);
+
+ /* Ionize the atoms if necessary */
+ if (bIonize)
+ {
+ ionize(fplog,oenv,mdatoms,top_global,t,ir,state->x,state->v,
+ mdatoms->start,mdatoms->start+mdatoms->homenr,state->box,cr);
+ }
+
+ /* Update force field in ffscan program */
+ if (bFFscan)
+ {
+ if (update_forcefield(fplog,
+ nfile,fnm,fr,
+ mdatoms->nr,state->x,state->box)) {
+ if (gmx_parallel_env_initialized())
+ {
+ gmx_finalize();
+ }
+ exit(0);
+ }
+ }
+
+ GMX_MPE_LOG(ev_timestep2);
+
+ /* We write a checkpoint at this MD step when:
+ * either at an NS step when we signalled through gs,
+ * or at the last step (but not when we do not want confout),
+ * but never at the first step or with rerun.
+ */
+ bCPT = (((gs.set[eglsCHKPT] && (bNS || ir->nstlist == 0)) ||
+ (bLastStep && (Flags & MD_CONFOUT))) &&
+ step > ir->init_step && !bRerunMD);
+ if (bCPT)
+ {
+ gs.set[eglsCHKPT] = 0;
+ }
+
+ /* Determine the energy and pressure:
+ * at nstcalcenergy steps and at energy output steps (set below).
+ */
+ bNstEner = do_per_step(step,ir->nstcalcenergy);
+ bCalcEnerPres =
+ (bNstEner ||
+ (ir->epc != epcNO && do_per_step(step,ir->nstpcouple)));
+
+ /* Do we need global communication ? */
+ bGStat = (bCalcEnerPres || bStopCM ||
+ do_per_step(step,nstglobalcomm) ||
+ (ir->nstlist == -1 && !bRerunMD && step >= nlh.step_nscheck));
+
+ do_ene = (do_per_step(step,ir->nstenergy) || bLastStep);
+
+ if (do_ene || do_log)
+ {
+ bCalcEnerPres = TRUE;
+ bGStat = TRUE;
+ }
+
+ /* these CGLO_ options remain the same throughout the iteration */
+ cglo_flags = ((bRerunMD ? CGLO_RERUNMD : 0) |
+ (bStopCM ? CGLO_STOPCM : 0) |
+ (bGStat ? CGLO_GSTAT : 0)
+ );
+
+ force_flags = (GMX_FORCE_STATECHANGED |
+ ((DYNAMIC_BOX(*ir) || bRerunMD) ? GMX_FORCE_DYNAMICBOX : 0) |
+ GMX_FORCE_ALLFORCES |
+ (bNStList ? GMX_FORCE_DOLR : 0) |
+ GMX_FORCE_SEPLRF |
+ (bCalcEnerPres ? GMX_FORCE_VIRIAL : 0) |
+ (bDoDHDL ? GMX_FORCE_DHDL : 0)
+ );
+
+ if (shellfc)
+ {
+ /* Now is the time to relax the shells */
+ count=relax_shell_flexcon(fplog,cr,bVerbose,bFFscan ? step+1 : step,
+ ir,bNS,force_flags,
+ bStopCM,top,top_global,
+ constr,enerd,fcd,
+ state,f,force_vir,mdatoms,
+ nrnb,wcycle,graph,groups,
+ shellfc,fr,bBornRadii,t,mu_tot,
+ state->natoms,&bConverged,vsite,
+ outf->fp_field);
+ tcount+=count;
+
+ if (bConverged)
+ {
+ nconverged++;
+ }
+ }
+ else
+ {
+ /* The coordinates (x) are shifted (to get whole molecules)
+ * in do_force.
+ * This is parallellized as well, and does communication too.
+ * Check comments in sim_util.c
+ */
+
+ do_force(fplog,cr,ir,step,nrnb,wcycle,top,top_global,groups,
+ state->box,state->x,&state->hist,
+ f,force_vir,mdatoms,enerd,fcd,
+ state->lambda,graph,
+ fr,vsite,mu_tot,t,outf->fp_field,ed,bBornRadii,
+ (bNS ? GMX_FORCE_NS : 0) | force_flags);
+ }
+
+ GMX_BARRIER(cr->mpi_comm_mygroup);
+
+ if (bTCR)
+ {
+ mu_aver = calc_mu_aver(cr,state->x,mdatoms->chargeA,
+ mu_tot,&top_global->mols,mdatoms,gnx,grpindex);
+ }
+
+ if (bTCR && bFirstStep)
+ {
+ tcr=init_coupling(fplog,nfile,fnm,cr,fr,mdatoms,&(top->idef));
+ fprintf(fplog,"Done init_coupling\n");
+ fflush(fplog);
+ }
+
+ if (bVV && !bStartingFromCpt && !bRerunMD)
+ /* ############### START FIRST UPDATE HALF-STEP FOR VV METHODS############### */
+ {
+ if (ir->eI==eiVV && bInitStep)
+ {
+ /* if using velocity verlet with full time step Ekin,
+ * take the first half step only to compute the
+ * virial for the first step. From there,
+ * revert back to the initial coordinates
+ * so that the input is actually the initial step.
+ */
+ copy_rvecn(state->v,cbuf,0,state->natoms); /* should make this better for parallelizing? */
+ } else {
+ /* this is for NHC in the Ekin(t+dt/2) version of vv */
+ trotter_update(ir,step,ekind,enerd,state,total_vir,mdatoms,&MassQ,trotter_seq,ettTSEQ1);
+ }
+
+ update_coords(fplog,step,ir,mdatoms,state,
+ f,fr->bTwinRange && bNStList,fr->f_twin,fcd,
+ ekind,M,wcycle,upd,bInitStep,etrtVELOCITY1,
+ cr,nrnb,constr,&top->idef);
+
+ if (bIterations)
+ {
+ gmx_iterate_init(&iterate,bIterations && !bInitStep);
+ }
+ /* for iterations, we save these vectors, as we will be self-consistently iterating
+ the calculations */
+
+ /*#### UPDATE EXTENDED VARIABLES IN TROTTER FORMULATION */
+
+ /* save the state */
+ if (bIterations && iterate.bIterate) {
+ copy_coupling_state(state,bufstate,ekind,ekind_save,&(ir->opts));
+ }
+
+ bFirstIterate = TRUE;
+ while (bFirstIterate || (bIterations && iterate.bIterate))
+ {
+ if (bIterations && iterate.bIterate)
+ {
+ copy_coupling_state(bufstate,state,ekind_save,ekind,&(ir->opts));
+ if (bFirstIterate && bTrotter)
+ {
+ /* The first time through, we need a decent first estimate
+ of veta(t+dt) to compute the constraints. Do
+ this by computing the box volume part of the
+ trotter integration at this time. Nothing else
+ should be changed by this routine here. If
+ !(first time), we start with the previous value
+ of veta. */
+
+ veta_save = state->veta;
+ trotter_update(ir,step,ekind,enerd,state,total_vir,mdatoms,&MassQ,trotter_seq,ettTSEQ0);
+ vetanew = state->veta;
+ state->veta = veta_save;
+ }
+ }
+
+ bOK = TRUE;
+ if ( !bRerunMD || rerun_fr.bV || bForceUpdate) { /* Why is rerun_fr.bV here? Unclear. */
+ dvdl = 0;
+
+ update_constraints(fplog,step,&dvdl,ir,ekind,mdatoms,state,graph,f,
+ &top->idef,shake_vir,NULL,
+ cr,nrnb,wcycle,upd,constr,
+ bInitStep,TRUE,bCalcEnerPres,vetanew);
+
+ if (!bOK && !bFFscan)
+ {
+ gmx_fatal(FARGS,"Constraint error: Shake, Lincs or Settle could not solve the constrains");
+ }
+
+ }
+ else if (graph)
+ { /* Need to unshift here if a do_force has been
+ called in the previous step */
+ unshift_self(graph,state->box,state->x);
+ }
+
+
+ /* if VV, compute the pressure and constraints */
+ /* For VV2, we strictly only need this if using pressure
+ * control, but we really would like to have accurate pressures
+ * printed out.
+ * Think about ways around this in the future?
+ * For now, keep this choice in comments.
+ */
+ /*bPres = (ir->eI==eiVV || IR_NPT_TROTTER(ir)); */
+ /*bTemp = ((ir->eI==eiVV &&(!bInitStep)) || (ir->eI==eiVVAK && IR_NPT_TROTTER(ir)));*/
+ bPres = TRUE;
+ bTemp = ((ir->eI==eiVV &&(!bInitStep)) || (ir->eI==eiVVAK));
+ compute_globals(fplog,gstat,cr,ir,fr,ekind,state,state_global,mdatoms,nrnb,vcm,
+ wcycle,enerd,force_vir,shake_vir,total_vir,pres,mu_tot,
+ constr,NULL,FALSE,state->box,
+ top_global,&pcurr,top_global->natoms,&bSumEkinhOld,
+ cglo_flags
+ | CGLO_ENERGY
+ | (bTemp ? CGLO_TEMPERATURE:0)
+ | (bPres ? CGLO_PRESSURE : 0)
+ | (bPres ? CGLO_CONSTRAINT : 0)
+ | ((bIterations && iterate.bIterate) ? CGLO_ITERATE : 0)
+ | (bFirstIterate ? CGLO_FIRSTITERATE : 0)
+ | CGLO_SCALEEKIN
+ );
+ /* explanation of above:
+ a) We compute Ekin at the full time step
+ if 1) we are using the AveVel Ekin, and it's not the
+ initial step, or 2) if we are using AveEkin, but need the full
+ time step kinetic energy for the pressure (always true now, since we want accurate statistics).
+ b) If we are using EkinAveEkin for the kinetic energy for the temperture control, we still feed in
+ EkinAveVel because it's needed for the pressure */
+
+ /* temperature scaling and pressure scaling to produce the extended variables at t+dt */
+ if (!bInitStep)
+ {
+ if (bTrotter)
+ {
+ trotter_update(ir,step,ekind,enerd,state,total_vir,mdatoms,&MassQ,trotter_seq,ettTSEQ2);
+ }
+ else
+ {
+ update_tcouple(fplog,step,ir,state,ekind,wcycle,upd,&MassQ,mdatoms);
+ }
+ }
+
+ if (bIterations &&
+ done_iterating(cr,fplog,step,&iterate,bFirstIterate,
+ state->veta,&vetanew))
+ {
+ break;
+ }
+ bFirstIterate = FALSE;
+ }
+
+ if (bTrotter && !bInitStep) {
+ copy_mat(shake_vir,state->svir_prev);
+ copy_mat(force_vir,state->fvir_prev);
+ if (IR_NVT_TROTTER(ir) && ir->eI==eiVV) {
+ /* update temperature and kinetic energy now that step is over - this is the v(t+dt) point */
+ enerd->term[F_TEMP] = sum_ekin(&(ir->opts),ekind,NULL,(ir->eI==eiVV),FALSE,FALSE);
+ enerd->term[F_EKIN] = trace(ekind->ekin);
+ }
+ }
+ /* if it's the initial step, we performed this first step just to get the constraint virial */
+ if (bInitStep && ir->eI==eiVV) {
+ copy_rvecn(cbuf,state->v,0,state->natoms);
+ }
+
+ if (fr->bSepDVDL && fplog && do_log)
+ {
+ fprintf(fplog,sepdvdlformat,"Constraint",0.0,dvdl);
+ }
+ enerd->term[F_DHDL_CON] += dvdl;
+
+ GMX_MPE_LOG(ev_timestep1);
+ }
+
+ /* MRS -- now done iterating -- compute the conserved quantity */
+ if (bVV) {
+ saved_conserved_quantity = compute_conserved_from_auxiliary(ir,state,&MassQ);
+ if (ir->eI==eiVV)
+ {
+ last_ekin = enerd->term[F_EKIN]; /* does this get preserved through checkpointing? */
+ }
+ if ((ir->eDispCorr != edispcEnerPres) && (ir->eDispCorr != edispcAllEnerPres))
+ {
+ saved_conserved_quantity -= enerd->term[F_DISPCORR];
+ }
+ }
+
+ /* ######## END FIRST UPDATE STEP ############## */
+ /* ######## If doing VV, we now have v(dt) ###### */
+
+ /* ################## START TRAJECTORY OUTPUT ################# */
+
+ /* Now we have the energies and forces corresponding to the
+ * coordinates at time t. We must output all of this before
+ * the update.
+ * for RerunMD t is read from input trajectory
+ */
+ GMX_MPE_LOG(ev_output_start);
+
+ mdof_flags = 0;
+ if (do_per_step(step,ir->nstxout)) { mdof_flags |= MDOF_X; }
+ if (do_per_step(step,ir->nstvout)) { mdof_flags |= MDOF_V; }
+ if (do_per_step(step,ir->nstfout)) { mdof_flags |= MDOF_F; }
+ if (do_per_step(step,ir->nstxtcout)) { mdof_flags |= MDOF_XTC; }
+ if (bCPT) { mdof_flags |= MDOF_CPT; };
+
+#if defined(GMX_FAHCORE) || defined(GMX_WRITELASTSTEP)
+ if (bLastStep)
+ {
+ /* Enforce writing positions and velocities at end of run */
+ mdof_flags |= (MDOF_X | MDOF_V);
+ }
+#endif
+#ifdef GMX_FAHCORE
+ if (MASTER(cr))
+ fcReportProgress( ir->nsteps, step );
+
+ /* sync bCPT and fc record-keeping */
+ if (bCPT && MASTER(cr))
+ fcRequestCheckPoint();
+#endif
+
+ if (mdof_flags != 0)
+ {
+ wallcycle_start(wcycle,ewcTRAJ);
+ if (bCPT)
+ {
+ if (state->flags & (1<<estLD_RNG))
+ {
+ get_stochd_state(upd,state);
+ }
+ if (MASTER(cr))
+ {
+ if (bSumEkinhOld)
+ {
+ state_global->ekinstate.bUpToDate = FALSE;
+ }
+ else
+ {
+ update_ekinstate(&state_global->ekinstate,ekind);
+ state_global->ekinstate.bUpToDate = TRUE;
+ }
+ update_energyhistory(&state_global->enerhist,mdebin);
+ }
+ }
+ write_traj(fplog,cr,outf,mdof_flags,top_global,
+ step,t,state,state_global,f,f_global,&n_xtc,&x_xtc);
+ if (bCPT)
+ {
+ nchkpt++;
+ bCPT = FALSE;
+ }
+ debug_gmx();
+ if (bLastStep && step_rel == ir->nsteps &&
+ (Flags & MD_CONFOUT) && MASTER(cr) &&
+ !bRerunMD && !bFFscan)
+ {
+ /* x and v have been collected in write_traj,
+ * because a checkpoint file will always be written
+ * at the last step.
+ */
+ fprintf(stderr,"\nWriting final coordinates.\n");
+ if (ir->ePBC != epbcNONE && !ir->bPeriodicMols &&
+ DOMAINDECOMP(cr))
+ {
+ /* Make molecules whole only for confout writing */
+ do_pbc_mtop(fplog,ir->ePBC,state->box,top_global,state_global->x);
+ }
+ write_sto_conf_mtop(ftp2fn(efSTO,nfile,fnm),
+ *top_global->name,top_global,
+ state_global->x,state_global->v,
+ ir->ePBC,state->box);
+ debug_gmx();
+ }
+ wallcycle_stop(wcycle,ewcTRAJ);
+ }
+ GMX_MPE_LOG(ev_output_finish);
+
+ /* kludge -- virial is lost with restart for NPT control. Must restart */
+ if (bStartingFromCpt && bVV)
+ {
+ copy_mat(state->svir_prev,shake_vir);
+ copy_mat(state->fvir_prev,force_vir);
+ }
+ /* ################## END TRAJECTORY OUTPUT ################ */
+
+ /* Determine the pressure:
+ * always when we want exact averages in the energy file,
+ * at ns steps when we have pressure coupling,
+ * otherwise only at energy output steps (set below).
+ */
+
+ bNstEner = (bGStatEveryStep || do_per_step(step,ir->nstcalcenergy));
+ bCalcEnerPres = bNstEner;
+
+ /* Do we need global communication ? */
+ bGStat = (bGStatEveryStep || bStopCM || bNS ||
+ (ir->nstlist == -1 && !bRerunMD && step >= nlh.step_nscheck));
+
+ do_ene = (do_per_step(step,ir->nstenergy) || bLastStep);
+
+ if (do_ene || do_log)
+ {
+ bCalcEnerPres = TRUE;
+ bGStat = TRUE;
+ }
+
+ /* Determine the wallclock run time up till now */
+ run_time = gmx_gettime() - (double)runtime->real;
+
+ /* Check whether everything is still allright */
+ if (((int)gmx_get_stop_condition() > handled_stop_condition)
+#ifdef GMX_THREADS
+ && MASTER(cr)
+#endif
+ )
+ {
+ /* this is just make gs.sig compatible with the hack
+ of sending signals around by MPI_Reduce with together with
+ other floats */
+ if ( gmx_get_stop_condition() == gmx_stop_cond_next_ns )
+ gs.sig[eglsSTOPCOND]=1;
+ if ( gmx_get_stop_condition() == gmx_stop_cond_next )
+ gs.sig[eglsSTOPCOND]=-1;
+ /* < 0 means stop at next step, > 0 means stop at next NS step */
+ if (fplog)
+ {
+ fprintf(fplog,
+ "\n\nReceived the %s signal, stopping at the next %sstep\n\n",
+ gmx_get_signal_name(),
+ gs.sig[eglsSTOPCOND]==1 ? "NS " : "");
+ fflush(fplog);
+ }
+ fprintf(stderr,
+ "\n\nReceived the %s signal, stopping at the next %sstep\n\n",
+ gmx_get_signal_name(),
+ gs.sig[eglsSTOPCOND]==1 ? "NS " : "");
+ fflush(stderr);
+ handled_stop_condition=(int)gmx_get_stop_condition();
+ }
+ else if (MASTER(cr) && (bNS || ir->nstlist <= 0) &&
+ (max_hours > 0 && run_time > max_hours*60.0*60.0*0.99) &&
+ gs.sig[eglsSTOPCOND] == 0 && gs.set[eglsSTOPCOND] == 0)
+ {
+ /* Signal to terminate the run */
+ gs.sig[eglsSTOPCOND] = 1;
+ if (fplog)
+ {
+ fprintf(fplog,"\nStep %s: Run time exceeded %.3f hours, will terminate the run\n",gmx_step_str(step,sbuf),max_hours*0.99);
+ }
+ fprintf(stderr, "\nStep %s: Run time exceeded %.3f hours, will terminate the run\n",gmx_step_str(step,sbuf),max_hours*0.99);
+ }
+
+ if (bResetCountersHalfMaxH && MASTER(cr) &&
+ run_time > max_hours*60.0*60.0*0.495)
+ {
+ gs.sig[eglsRESETCOUNTERS] = 1;
+ }
+
+ if (ir->nstlist == -1 && !bRerunMD)
+ {
+ /* When bGStatEveryStep=FALSE, global_stat is only called
+ * when we check the atom displacements, not at NS steps.
+ * This means that also the bonded interaction count check is not
+ * performed immediately after NS. Therefore a few MD steps could
+ * be performed with missing interactions.
+ * But wrong energies are never written to file,
+ * since energies are only written after global_stat
+ * has been called.
+ */
+ if (step >= nlh.step_nscheck)
+ {
+ nlh.nabnsb = natoms_beyond_ns_buffer(ir,fr,&top->cgs,
+ nlh.scale_tot,state->x);
+ }
+ else
+ {
+ /* This is not necessarily true,
+ * but step_nscheck is determined quite conservatively.
+ */
+ nlh.nabnsb = 0;
+ }
+ }
+
+ /* In parallel we only have to check for checkpointing in steps
+ * where we do global communication,
+ * otherwise the other nodes don't know.
+ */
+ if (MASTER(cr) && ((bGStat || !PAR(cr)) &&
+ cpt_period >= 0 &&
+ (cpt_period == 0 ||
+ run_time >= nchkpt*cpt_period*60.0)) &&
+ gs.set[eglsCHKPT] == 0)
+ {
+ gs.sig[eglsCHKPT] = 1;
+ }
+
+ if (bIterations)
+ {
+ gmx_iterate_init(&iterate,bIterations);
+ }
+
+ /* for iterations, we save these vectors, as we will be redoing the calculations */
+ if (bIterations && iterate.bIterate)
+ {
+ copy_coupling_state(state,bufstate,ekind,ekind_save,&(ir->opts));
+ }
+ bFirstIterate = TRUE;
+ while (bFirstIterate || (bIterations && iterate.bIterate))
+ {
+ /* We now restore these vectors to redo the calculation with improved extended variables */
+ if (bIterations)
+ {
+ copy_coupling_state(bufstate,state,ekind_save,ekind,&(ir->opts));
+ }
+
+ /* We make the decision to break or not -after- the calculation of Ekin and Pressure,
+ so scroll down for that logic */
+
+ /* ######### START SECOND UPDATE STEP ################# */
+ GMX_MPE_LOG(ev_update_start);
+ /* Box is changed in update() when we do pressure coupling,
+ * but we should still use the old box for energy corrections and when
+ * writing it to the energy file, so it matches the trajectory files for
+ * the same timestep above. Make a copy in a separate array.
+ */
+ copy_mat(state->box,lastbox);
+
+ bOK = TRUE;
+ if (!(bRerunMD && !rerun_fr.bV && !bForceUpdate))
+ {
+ wallcycle_start(wcycle,ewcUPDATE);
+ dvdl = 0;
+ /* UPDATE PRESSURE VARIABLES IN TROTTER FORMULATION WITH CONSTRAINTS */
+ if (bTrotter)
+ {
+ if (bIterations && iterate.bIterate)
+ {
+ if (bFirstIterate)
+ {
+ scalevir = 1;
+ }
+ else
+ {
+ /* we use a new value of scalevir to converge the iterations faster */
+ scalevir = tracevir/trace(shake_vir);
+ }
+ msmul(shake_vir,scalevir,shake_vir);
+ m_add(force_vir,shake_vir,total_vir);
+ clear_mat(shake_vir);
+ }
+ trotter_update(ir,step,ekind,enerd,state,total_vir,mdatoms,&MassQ,trotter_seq,ettTSEQ3);
+ /* We can only do Berendsen coupling after we have summed
+ * the kinetic energy or virial. Since the happens
+ * in global_state after update, we should only do it at
+ * step % nstlist = 1 with bGStatEveryStep=FALSE.
+ */
+ }
+ else
+ {
+ update_tcouple(fplog,step,ir,state,ekind,wcycle,upd,&MassQ,mdatoms);
+ update_pcouple(fplog,step,ir,state,pcoupl_mu,M,wcycle,
+ upd,bInitStep);
+ }
+
+ if (bVV)
+ {
+ /* velocity half-step update */
+ update_coords(fplog,step,ir,mdatoms,state,f,
+ fr->bTwinRange && bNStList,fr->f_twin,fcd,
+ ekind,M,wcycle,upd,FALSE,etrtVELOCITY2,
+ cr,nrnb,constr,&top->idef);
+ }
+
+ /* Above, initialize just copies ekinh into ekin,
+ * it doesn't copy position (for VV),
+ * and entire integrator for MD.
+ */
+
+ if (ir->eI==eiVVAK)
+ {
+ copy_rvecn(state->x,cbuf,0,state->natoms);
+ }
+
+ update_coords(fplog,step,ir,mdatoms,state,f,fr->bTwinRange && bNStList,fr->f_twin,fcd,
+ ekind,M,wcycle,upd,bInitStep,etrtPOSITION,cr,nrnb,constr,&top->idef);
+ wallcycle_stop(wcycle,ewcUPDATE);
+
+ update_constraints(fplog,step,&dvdl,ir,ekind,mdatoms,state,graph,f,
+ &top->idef,shake_vir,force_vir,
+ cr,nrnb,wcycle,upd,constr,
+ bInitStep,FALSE,bCalcEnerPres,state->veta);
+
+ if (ir->eI==eiVVAK)
+ {
+ /* erase F_EKIN and F_TEMP here? */
+ /* just compute the kinetic energy at the half step to perform a trotter step */
+ compute_globals(fplog,gstat,cr,ir,fr,ekind,state,state_global,mdatoms,nrnb,vcm,
+ wcycle,enerd,force_vir,shake_vir,total_vir,pres,mu_tot,
+ constr,NULL,FALSE,lastbox,
+ top_global,&pcurr,top_global->natoms,&bSumEkinhOld,
+ cglo_flags | CGLO_TEMPERATURE
+ );
+ wallcycle_start(wcycle,ewcUPDATE);
+ trotter_update(ir,step,ekind,enerd,state,total_vir,mdatoms,&MassQ,trotter_seq,ettTSEQ4);
+ /* now we know the scaling, we can compute the positions again again */
+ copy_rvecn(cbuf,state->x,0,state->natoms);
+
+ update_coords(fplog,step,ir,mdatoms,state,f,fr->bTwinRange && bNStList,fr->f_twin,fcd,
+ ekind,M,wcycle,upd,bInitStep,etrtPOSITION,cr,nrnb,constr,&top->idef);
+ wallcycle_stop(wcycle,ewcUPDATE);
+
+ /* do we need an extra constraint here? just need to copy out of state->v to upd->xp? */
+ /* are the small terms in the shake_vir here due
+ * to numerical errors, or are they important
+ * physically? I'm thinking they are just errors, but not completely sure.
+ * For now, will call without actually constraining, constr=NULL*/
+ update_constraints(fplog,step,&dvdl,ir,ekind,mdatoms,state,graph,f,
+ &top->idef,tmp_vir,force_vir,
+ cr,nrnb,wcycle,upd,NULL,
+ bInitStep,FALSE,bCalcEnerPres,
+ state->veta);
+ }
+ if (!bOK && !bFFscan)
+ {
+ gmx_fatal(FARGS,"Constraint error: Shake, Lincs or Settle could not solve the constrains");
+ }
+
+ if (fr->bSepDVDL && fplog && do_log)
+ {
+ fprintf(fplog,sepdvdlformat,"Constraint",0.0,dvdl);
+ }
+ enerd->term[F_DHDL_CON] += dvdl;
+ }
+ else if (graph)
+ {
+ /* Need to unshift here */
+ unshift_self(graph,state->box,state->x);
+ }
+
+ GMX_BARRIER(cr->mpi_comm_mygroup);
+ GMX_MPE_LOG(ev_update_finish);
+
+ if (vsite != NULL)
+ {
+ wallcycle_start(wcycle,ewcVSITECONSTR);
+ if (graph != NULL)
+ {
+ shift_self(graph,state->box,state->x);
+ }
+ construct_vsites(fplog,vsite,state->x,nrnb,ir->delta_t,state->v,
+ top->idef.iparams,top->idef.il,
+ fr->ePBC,fr->bMolPBC,graph,cr,state->box);
+
+ if (graph != NULL)
+ {
+ unshift_self(graph,state->box,state->x);
+ }
+ wallcycle_stop(wcycle,ewcVSITECONSTR);
+ }
+
+ /* ############## IF NOT VV, Calculate globals HERE, also iterate constraints ############ */
+ if (ir->nstlist == -1 && bFirstIterate)
+ {
+ gs.sig[eglsNABNSB] = nlh.nabnsb;
+ }
+ compute_globals(fplog,gstat,cr,ir,fr,ekind,state,state_global,mdatoms,nrnb,vcm,
+ wcycle,enerd,force_vir,shake_vir,total_vir,pres,mu_tot,
+ constr,
+ bFirstIterate ? &gs : NULL,(step % gs.nstms == 0),
+ lastbox,
+ top_global,&pcurr,top_global->natoms,&bSumEkinhOld,
+ cglo_flags
+ | (!EI_VV(ir->eI) ? CGLO_ENERGY : 0)
+ | (!EI_VV(ir->eI) ? CGLO_TEMPERATURE : 0)
+ | (!EI_VV(ir->eI) || bRerunMD ? CGLO_PRESSURE : 0)
+ | (bIterations && iterate.bIterate ? CGLO_ITERATE : 0)
+ | (bFirstIterate ? CGLO_FIRSTITERATE : 0)
+ | CGLO_CONSTRAINT
+ );
+ if (ir->nstlist == -1 && bFirstIterate)
+ {
+ nlh.nabnsb = gs.set[eglsNABNSB];
+ gs.set[eglsNABNSB] = 0;
+ }
+ /* bIterate is set to keep it from eliminating the old ekin kinetic energy terms */
+ /* ############# END CALC EKIN AND PRESSURE ################# */
+
+ /* Note: this is OK, but there are some numerical precision issues with using the convergence of
+ the virial that should probably be addressed eventually. state->veta has better properies,
+ but what we actually need entering the new cycle is the new shake_vir value. Ideally, we could
+ generate the new shake_vir, but test the veta value for convergence. This will take some thought. */
+
+ if (bIterations &&
+ done_iterating(cr,fplog,step,&iterate,bFirstIterate,
+ trace(shake_vir),&tracevir))
+ {
+ break;
+ }
+ bFirstIterate = FALSE;
+ }
+
+ update_box(fplog,step,ir,mdatoms,state,graph,f,
+ ir->nstlist==-1 ? &nlh.scale_tot : NULL,pcoupl_mu,nrnb,wcycle,upd,bInitStep,FALSE);
+
+ /* ################# END UPDATE STEP 2 ################# */
+ /* #### We now have r(t+dt) and v(t+dt/2) ############# */
+
+ /* The coordinates (x) were unshifted in update */
+ if (bFFscan && (shellfc==NULL || bConverged))
+ {
+ if (print_forcefield(fplog,enerd->term,mdatoms->homenr,
+ f,NULL,xcopy,
+ &(top_global->mols),mdatoms->massT,pres))
+ {
+ if (gmx_parallel_env_initialized())
+ {
+ gmx_finalize();
+ }
+ fprintf(stderr,"\n");
+ exit(0);
+ }
+ }
+ if (!bGStat)
+ {
+ /* We will not sum ekinh_old,
+ * so signal that we still have to do it.
+ */
+ bSumEkinhOld = TRUE;
+ }
+
+ if (bTCR)
+ {
+ /* Only do GCT when the relaxation of shells (minimization) has converged,
+ * otherwise we might be coupling to bogus energies.
+ * In parallel we must always do this, because the other sims might
+ * update the FF.
+ */
+
+ /* Since this is called with the new coordinates state->x, I assume
+ * we want the new box state->box too. / EL 20040121
+ */
+ do_coupling(fplog,oenv,nfile,fnm,tcr,t,step,enerd->term,fr,
+ ir,MASTER(cr),
+ mdatoms,&(top->idef),mu_aver,
+ top_global->mols.nr,cr,
+ state->box,total_vir,pres,
+ mu_tot,state->x,f,bConverged);
+ debug_gmx();
+ }
+
+ /* ######### BEGIN PREPARING EDR OUTPUT ########### */
+
+ /* sum up the foreign energy and dhdl terms */
+ sum_dhdl(enerd,state->lambda,ir);
+
+ /* use the directly determined last velocity, not actually the averaged half steps */
+ if (bTrotter && ir->eI==eiVV)
+ {
+ enerd->term[F_EKIN] = last_ekin;
+ }
+ enerd->term[F_ETOT] = enerd->term[F_EPOT] + enerd->term[F_EKIN];
+
+ if (bVV)
+ {
+ enerd->term[F_ECONSERVED] = enerd->term[F_ETOT] + saved_conserved_quantity;
+ }
+ else
+ {
+ enerd->term[F_ECONSERVED] = enerd->term[F_ETOT] + compute_conserved_from_auxiliary(ir,state,&MassQ);
+ }
+ /* Check for excessively large energies */
+ if (bIonize)
+ {
+#ifdef GMX_DOUBLE
+ real etot_max = 1e200;
+#else
+ real etot_max = 1e30;
+#endif
+ if (fabs(enerd->term[F_ETOT]) > etot_max)
+ {
+ fprintf(stderr,"Energy too large (%g), giving up\n",
+ enerd->term[F_ETOT]);
+ }
+ }
+ /* ######### END PREPARING EDR OUTPUT ########### */
+
+ /* Time for performance */
+ if (((step % stepout) == 0) || bLastStep)
+ {
+ runtime_upd_proc(runtime);
+ }
+
+ /* Output stuff */
+ if (MASTER(cr))
+ {
+ gmx_bool do_dr,do_or;
+
+ if (!(bStartingFromCpt && (EI_VV(ir->eI))))
+ {
+ if (bNstEner)
+ {
+ upd_mdebin(mdebin,bDoDHDL, TRUE,
+ t,mdatoms->tmass,enerd,state,lastbox,
+ shake_vir,force_vir,total_vir,pres,
+ ekind,mu_tot,constr);
+ }
+ else
+ {
+ upd_mdebin_step(mdebin);
+ }
+
+ do_dr = do_per_step(step,ir->nstdisreout);
+ do_or = do_per_step(step,ir->nstorireout);
+
+ print_ebin(outf->fp_ene,do_ene,do_dr,do_or,do_log?fplog:NULL,
+ step,t,
+ eprNORMAL,bCompact,mdebin,fcd,groups,&(ir->opts));
+ }
+ if (ir->ePull != epullNO)
+ {
+ pull_print_output(ir->pull,step,t);
+ }
+
+ if (do_per_step(step,ir->nstlog))
+ {
+ if(fflush(fplog) != 0)
+ {
+ gmx_fatal(FARGS,"Cannot flush logfile - maybe you are out of quota?");
+ }
+ }
+ }
+
+
+ /* Remaining runtime */
+ if (MULTIMASTER(cr) && (do_verbose || gmx_got_usr_signal() ))
+ {
+ if (shellfc)
+ {
+ fprintf(stderr,"\n");
+ }
+ print_time(stderr,runtime,step,ir,cr);
+ }
+
+ /* Replica exchange */
+ bExchanged = FALSE;
+ if ((repl_ex_nst > 0) && (step > 0) && !bLastStep &&
+ do_per_step(step,repl_ex_nst))
+ {
+ bExchanged = replica_exchange(fplog,cr,repl_ex,
+ state_global,enerd->term,
+ state,step,t);
+
+ if (bExchanged && DOMAINDECOMP(cr))
+ {
+ dd_partition_system(fplog,step,cr,TRUE,1,
+ state_global,top_global,ir,
+ state,&f,mdatoms,top,fr,
+ vsite,shellfc,constr,
+ nrnb,wcycle,FALSE);
+ }
+ }
+
+ bFirstStep = FALSE;
+ bInitStep = FALSE;
+ bStartingFromCpt = FALSE;
+
+ /* ####### SET VARIABLES FOR NEXT ITERATION IF THEY STILL NEED IT ###### */
+ /* With all integrators, except VV, we need to retain the pressure
+ * at the current step for coupling at the next step.
+ */
+ if ((state->flags & (1<<estPRES_PREV)) &&
+ (bGStatEveryStep ||
+ (ir->nstpcouple > 0 && step % ir->nstpcouple == 0)))
+ {
+ /* Store the pressure in t_state for pressure coupling
+ * at the next MD step.
+ */
+ copy_mat(pres,state->pres_prev);
+ }
+
+ /* ####### END SET VARIABLES FOR NEXT ITERATION ###### */
+
+ if ( (membed!=NULL) && (!bLastStep) )
+ rescale_membed(step_rel,membed,state_global->x);
+
+ if (bRerunMD)
+ {
+ if (MASTER(cr))
+ {
+ /* read next frame from input trajectory */
+ bNotLastFrame = read_next_frame(oenv,status,&rerun_fr);
+ }
+
+ if (PAR(cr))
+ {
+ rerun_parallel_comm(cr,&rerun_fr,&bNotLastFrame);
+ }
+ }
+
+ if (!bRerunMD || !rerun_fr.bStep)
+ {
+ /* increase the MD step number */
+ step++;
+ step_rel++;
+ }
+
+ cycles = wallcycle_stop(wcycle,ewcSTEP);
+ if (DOMAINDECOMP(cr) && wcycle)
+ {
+ dd_cycles_add(cr->dd,cycles,ddCyclStep);
+ }
+
+ if (step_rel == wcycle_get_reset_counters(wcycle) ||
+ gs.set[eglsRESETCOUNTERS] != 0)
+ {
+ /* Reset all the counters related to performance over the run */
+ reset_all_counters(fplog,cr,step,&step_rel,ir,wcycle,nrnb,runtime);
+ wcycle_set_reset_counters(wcycle,-1);
+ /* Correct max_hours for the elapsed time */
+ max_hours -= run_time/(60.0*60.0);
+ bResetCountersHalfMaxH = FALSE;
+ gs.set[eglsRESETCOUNTERS] = 0;
+ }
+ }
+ /* End of main MD loop */
+ debug_gmx();
+
+ /* Stop the time */
+ runtime_end(runtime);
+
+ if (bRerunMD && MASTER(cr))
+ {
+ close_trj(status);
+ }
+
+ if (!(cr->duty & DUTY_PME))
+ {
+ /* Tell the PME only node to finish */
+ gmx_pme_finish(cr);
+ }
+
+ if (MASTER(cr))
+ {
+ if (ir->nstcalcenergy > 0 && !bRerunMD)
+ {
+ print_ebin(outf->fp_ene,FALSE,FALSE,FALSE,fplog,step,t,
+ eprAVER,FALSE,mdebin,fcd,groups,&(ir->opts));
+ }
+ }
+
+ done_mdoutf(outf);
+
+ debug_gmx();
+
+ if (ir->nstlist == -1 && nlh.nns > 0 && fplog)
+ {
+ fprintf(fplog,"Average neighborlist lifetime: %.1f steps, std.dev.: %.1f steps\n",nlh.s1/nlh.nns,sqrt(nlh.s2/nlh.nns - sqr(nlh.s1/nlh.nns)));
+ fprintf(fplog,"Average number of atoms that crossed the half buffer length: %.1f\n\n",nlh.ab/nlh.nns);
+ }
+
+ if (shellfc && fplog)
+ {
+ fprintf(fplog,"Fraction of iterations that converged: %.2f %%\n",
+ (nconverged*100.0)/step_rel);
+ fprintf(fplog,"Average number of force evaluations per MD step: %.2f\n\n",
+ tcount/step_rel);
+ }
+
+ if (repl_ex_nst > 0 && MASTER(cr))
+ {
+ print_replica_exchange_statistics(fplog,repl_ex);
+ }
+
+ runtime->nsteps_done = step_rel;
+
+ return 0;
+}
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff