--- /dev/null
- *vol*sqr(box[ZZ][ZZ]*NANO/(2*M_PI)));
+/* -*- 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:
+ * GROwing Monsters And Cloning Shrimps
+ */
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <string.h>
+#include <float.h>
+#include "typedefs.h"
+#include "string2.h"
+#include "mdebin.h"
+#include "smalloc.h"
+#include "physics.h"
+#include "enxio.h"
+#include "vec.h"
+#include "disre.h"
+#include "main.h"
+#include "network.h"
+#include "names.h"
+#include "orires.h"
+#include "constr.h"
+#include "mtop_util.h"
+#include "xvgr.h"
+#include "gmxfio.h"
+
+#include "mdebin_bar.h"
+
+
+static const char *conrmsd_nm[] = { "Constr. rmsd", "Constr.2 rmsd" };
+
+static const char *boxs_nm[] = { "Box-X", "Box-Y", "Box-Z" };
+
+static const char *tricl_boxs_nm[] = {
+ "Box-XX", "Box-YY", "Box-ZZ",
+ "Box-YX", "Box-ZX", "Box-ZY"
+};
+
+static const char *vol_nm[] = { "Volume" };
+
+static const char *dens_nm[] = {"Density" };
+
+static const char *pv_nm[] = {"pV" };
+
+static const char *enthalpy_nm[] = {"Enthalpy" };
+
+static const char *boxvel_nm[] = {
+ "Box-Vel-XX", "Box-Vel-YY", "Box-Vel-ZZ",
+ "Box-Vel-YX", "Box-Vel-ZX", "Box-Vel-ZY"
+};
+
+#define NBOXS asize(boxs_nm)
+#define NTRICLBOXS asize(tricl_boxs_nm)
+
+static gmx_bool bTricl,bDynBox;
+static int f_nre=0,epc,etc,nCrmsd;
+
+
+
+
+
+t_mdebin *init_mdebin(ener_file_t fp_ene,
+ const gmx_mtop_t *mtop,
+ const t_inputrec *ir,
+ FILE *fp_dhdl)
+{
+ const char *ener_nm[F_NRE];
+ static const char *vir_nm[] = {
+ "Vir-XX", "Vir-XY", "Vir-XZ",
+ "Vir-YX", "Vir-YY", "Vir-YZ",
+ "Vir-ZX", "Vir-ZY", "Vir-ZZ"
+ };
+ static const char *sv_nm[] = {
+ "ShakeVir-XX", "ShakeVir-XY", "ShakeVir-XZ",
+ "ShakeVir-YX", "ShakeVir-YY", "ShakeVir-YZ",
+ "ShakeVir-ZX", "ShakeVir-ZY", "ShakeVir-ZZ"
+ };
+ static const char *fv_nm[] = {
+ "ForceVir-XX", "ForceVir-XY", "ForceVir-XZ",
+ "ForceVir-YX", "ForceVir-YY", "ForceVir-YZ",
+ "ForceVir-ZX", "ForceVir-ZY", "ForceVir-ZZ"
+ };
+ static const char *pres_nm[] = {
+ "Pres-XX","Pres-XY","Pres-XZ",
+ "Pres-YX","Pres-YY","Pres-YZ",
+ "Pres-ZX","Pres-ZY","Pres-ZZ"
+ };
+ static const char *surft_nm[] = {
+ "#Surf*SurfTen"
+ };
+ static const char *mu_nm[] = {
+ "Mu-X", "Mu-Y", "Mu-Z"
+ };
+ static const char *vcos_nm[] = {
+ "2CosZ*Vel-X"
+ };
+ static const char *visc_nm[] = {
+ "1/Viscosity"
+ };
+ static const char *baro_nm[] = {
+ "Barostat"
+ };
+
+ char **grpnms;
+ const gmx_groups_t *groups;
+ char **gnm;
+ char buf[256];
+ const char *bufi;
+ t_mdebin *md;
+ int i,j,ni,nj,n,nh,k,kk,ncon,nset;
+ gmx_bool bBHAM,bNoseHoover,b14;
+
+ snew(md,1);
+
+ if (EI_DYNAMICS(ir->eI))
+ {
+ md->delta_t = ir->delta_t;
+ }
+ else
+ {
+ md->delta_t = 0;
+ }
+
+ groups = &mtop->groups;
+
+ bBHAM = (mtop->ffparams.functype[0] == F_BHAM);
+ b14 = (gmx_mtop_ftype_count(mtop,F_LJ14) > 0 ||
+ gmx_mtop_ftype_count(mtop,F_LJC14_Q) > 0);
+
+ ncon = gmx_mtop_ftype_count(mtop,F_CONSTR);
+ nset = gmx_mtop_ftype_count(mtop,F_SETTLE);
+ md->bConstr = (ncon > 0 || nset > 0);
+ md->bConstrVir = FALSE;
+ if (md->bConstr) {
+ if (ncon > 0 && ir->eConstrAlg == econtLINCS) {
+ if (ir->eI == eiSD2)
+ md->nCrmsd = 2;
+ else
+ md->nCrmsd = 1;
+ }
+ md->bConstrVir = (getenv("GMX_CONSTRAINTVIR") != NULL);
+ } else {
+ md->nCrmsd = 0;
+ }
+
+ /* Energy monitoring */
+ for(i=0;i<egNR;i++)
+ {
+ md->bEInd[i]=FALSE;
+ }
+
+#ifndef GMX_OPENMM
+ for(i=0; i<F_NRE; i++)
+ {
+ md->bEner[i] = FALSE;
+ if (i == F_LJ)
+ md->bEner[i] = !bBHAM;
+ else if (i == F_BHAM)
+ md->bEner[i] = bBHAM;
+ else if (i == F_EQM)
+ md->bEner[i] = ir->bQMMM;
+ else if (i == F_COUL_LR)
+ md->bEner[i] = (ir->rcoulomb > ir->rlist);
+ else if (i == F_LJ_LR)
+ md->bEner[i] = (!bBHAM && ir->rvdw > ir->rlist);
+ else if (i == F_BHAM_LR)
+ md->bEner[i] = (bBHAM && ir->rvdw > ir->rlist);
+ else if (i == F_RF_EXCL)
+ md->bEner[i] = (EEL_RF(ir->coulombtype) && ir->coulombtype != eelRF_NEC);
+ else if (i == F_COUL_RECIP)
+ md->bEner[i] = EEL_FULL(ir->coulombtype);
+ else if (i == F_LJ14)
+ md->bEner[i] = b14;
+ else if (i == F_COUL14)
+ md->bEner[i] = b14;
+ else if (i == F_LJC14_Q || i == F_LJC_PAIRS_NB)
+ md->bEner[i] = FALSE;
+ else if ((i == F_DVDL) || (i == F_DKDL))
+ md->bEner[i] = (ir->efep != efepNO);
+ else if (i == F_DHDL_CON)
+ md->bEner[i] = (ir->efep != efepNO && md->bConstr);
+ else if ((interaction_function[i].flags & IF_VSITE) ||
+ (i == F_CONSTR) || (i == F_CONSTRNC) || (i == F_SETTLE))
+ md->bEner[i] = FALSE;
+ else if ((i == F_COUL_SR) || (i == F_EPOT) || (i == F_PRES) || (i==F_EQM))
+ md->bEner[i] = TRUE;
+ else if ((i == F_GBPOL) && ir->implicit_solvent==eisGBSA)
+ md->bEner[i] = TRUE;
+ else if ((i == F_NPSOLVATION) && ir->implicit_solvent==eisGBSA && (ir->sa_algorithm != esaNO))
+ md->bEner[i] = TRUE;
+ else if ((i == F_GB12) || (i == F_GB13) || (i == F_GB14))
+ md->bEner[i] = FALSE;
+ else if ((i == F_ETOT) || (i == F_EKIN) || (i == F_TEMP))
+ md->bEner[i] = EI_DYNAMICS(ir->eI);
+ else if (i==F_VTEMP)
+ md->bEner[i] = (EI_DYNAMICS(ir->eI) && getenv("GMX_VIRIAL_TEMPERATURE"));
+ else if (i == F_DISPCORR || i == F_PDISPCORR)
+ md->bEner[i] = (ir->eDispCorr != edispcNO);
+ else if (i == F_DISRESVIOL)
+ md->bEner[i] = (gmx_mtop_ftype_count(mtop,F_DISRES) > 0);
+ else if (i == F_ORIRESDEV)
+ md->bEner[i] = (gmx_mtop_ftype_count(mtop,F_ORIRES) > 0);
+ else if (i == F_CONNBONDS)
+ md->bEner[i] = FALSE;
+ else if (i == F_COM_PULL)
+ md->bEner[i] = (ir->ePull == epullUMBRELLA || ir->ePull == epullCONST_F || ir->bRot);
+ else if (i == F_ECONSERVED)
+ md->bEner[i] = ((ir->etc == etcNOSEHOOVER || ir->etc == etcVRESCALE) &&
+ (ir->epc == epcNO || ir->epc==epcMTTK));
+ else
+ md->bEner[i] = (gmx_mtop_ftype_count(mtop,i) > 0);
+ }
+#else
+ /* OpenMM always produces only the following 4 energy terms */
+ md->bEner[F_EPOT] = TRUE;
+ md->bEner[F_EKIN] = TRUE;
+ md->bEner[F_ETOT] = TRUE;
+ md->bEner[F_TEMP] = TRUE;
+#endif
+
+ md->f_nre=0;
+ for(i=0; i<F_NRE; i++)
+ {
+ if (md->bEner[i])
+ {
+ /* FIXME: The constness should not be cast away */
+ /*ener_nm[f_nre]=(char *)interaction_function[i].longname;*/
+ ener_nm[md->f_nre]=interaction_function[i].longname;
+ md->f_nre++;
+ }
+ }
+
+ md->epc = ir->epc;
+ for (i=0;i<DIM;i++)
+ {
+ for (j=0;j<DIM;j++)
+ {
+ md->ref_p[i][j] = ir->ref_p[i][j];
+ }
+ }
+ md->bTricl = TRICLINIC(ir->compress) || TRICLINIC(ir->deform);
+ md->bDynBox = DYNAMIC_BOX(*ir);
+ md->etc = ir->etc;
+ md->bNHC_trotter = IR_NVT_TROTTER(ir);
+ md->bMTTK = IR_NPT_TROTTER(ir);
+
+ md->ebin = mk_ebin();
+ /* Pass NULL for unit to let get_ebin_space determine the units
+ * for interaction_function[i].longname
+ */
+ md->ie = get_ebin_space(md->ebin,md->f_nre,ener_nm,NULL);
+ if (md->nCrmsd)
+ {
+ /* This should be called directly after the call for md->ie,
+ * such that md->iconrmsd follows directly in the list.
+ */
+ md->iconrmsd = get_ebin_space(md->ebin,md->nCrmsd,conrmsd_nm,"");
+ }
+ if (md->bDynBox)
+ {
+ md->ib = get_ebin_space(md->ebin,
+ md->bTricl ? NTRICLBOXS : NBOXS,
+ md->bTricl ? tricl_boxs_nm : boxs_nm,
+ unit_length);
+ md->ivol = get_ebin_space(md->ebin, 1, vol_nm, unit_volume);
+ md->idens = get_ebin_space(md->ebin, 1, dens_nm, unit_density_SI);
+ md->ipv = get_ebin_space(md->ebin, 1, pv_nm, unit_energy);
+ md->ienthalpy = get_ebin_space(md->ebin, 1, enthalpy_nm, unit_energy);
+ }
+ if (md->bConstrVir)
+ {
+ md->isvir = get_ebin_space(md->ebin,asize(sv_nm),sv_nm,unit_energy);
+ md->ifvir = get_ebin_space(md->ebin,asize(fv_nm),fv_nm,unit_energy);
+ }
+ md->ivir = get_ebin_space(md->ebin,asize(vir_nm),vir_nm,unit_energy);
+ md->ipres = get_ebin_space(md->ebin,asize(pres_nm),pres_nm,unit_pres_bar);
+ md->isurft = get_ebin_space(md->ebin,asize(surft_nm),surft_nm,
+ unit_surft_bar);
+ if (md->epc == epcPARRINELLORAHMAN || md->epc == epcMTTK)
+ {
+ md->ipc = get_ebin_space(md->ebin,md->bTricl ? 6 : 3,
+ boxvel_nm,unit_vel);
+ }
+ md->imu = get_ebin_space(md->ebin,asize(mu_nm),mu_nm,unit_dipole_D);
+ if (ir->cos_accel != 0)
+ {
+ md->ivcos = get_ebin_space(md->ebin,asize(vcos_nm),vcos_nm,unit_vel);
+ md->ivisc = get_ebin_space(md->ebin,asize(visc_nm),visc_nm,
+ unit_invvisc_SI);
+ }
+
+ /* Energy monitoring */
+ for(i=0;i<egNR;i++)
+ {
+ md->bEInd[i] = FALSE;
+ }
+ md->bEInd[egCOULSR] = TRUE;
+ md->bEInd[egLJSR ] = TRUE;
+
+ if (ir->rcoulomb > ir->rlist)
+ {
+ md->bEInd[egCOULLR] = TRUE;
+ }
+ if (!bBHAM)
+ {
+ if (ir->rvdw > ir->rlist)
+ {
+ md->bEInd[egLJLR] = TRUE;
+ }
+ }
+ else
+ {
+ md->bEInd[egLJSR] = FALSE;
+ md->bEInd[egBHAMSR] = TRUE;
+ if (ir->rvdw > ir->rlist)
+ {
+ md->bEInd[egBHAMLR] = TRUE;
+ }
+ }
+ if (b14)
+ {
+ md->bEInd[egLJ14] = TRUE;
+ md->bEInd[egCOUL14] = TRUE;
+ }
+ md->nEc=0;
+ for(i=0; (i<egNR); i++)
+ {
+ if (md->bEInd[i])
+ {
+ md->nEc++;
+ }
+ }
+
+ n=groups->grps[egcENER].nr;
+ md->nEg=n;
+ md->nE=(n*(n+1))/2;
+ snew(md->igrp,md->nE);
+ if (md->nE > 1)
+ {
+ n=0;
+ snew(gnm,md->nEc);
+ for(k=0; (k<md->nEc); k++)
+ {
+ snew(gnm[k],STRLEN);
+ }
+ for(i=0; (i<groups->grps[egcENER].nr); i++)
+ {
+ ni=groups->grps[egcENER].nm_ind[i];
+ for(j=i; (j<groups->grps[egcENER].nr); j++)
+ {
+ nj=groups->grps[egcENER].nm_ind[j];
+ for(k=kk=0; (k<egNR); k++)
+ {
+ if (md->bEInd[k])
+ {
+ sprintf(gnm[kk],"%s:%s-%s",egrp_nm[k],
+ *(groups->grpname[ni]),*(groups->grpname[nj]));
+ kk++;
+ }
+ }
+ md->igrp[n]=get_ebin_space(md->ebin,md->nEc,
+ (const char **)gnm,unit_energy);
+ n++;
+ }
+ }
+ for(k=0; (k<md->nEc); k++)
+ {
+ sfree(gnm[k]);
+ }
+ sfree(gnm);
+
+ if (n != md->nE)
+ {
+ gmx_incons("Number of energy terms wrong");
+ }
+ }
+
+ md->nTC=groups->grps[egcTC].nr;
+ md->nNHC = ir->opts.nhchainlength; /* shorthand for number of NH chains */
+ if (md->bMTTK)
+ {
+ md->nTCP = 1; /* assume only one possible coupling system for barostat
+ for now */
+ }
+ else
+ {
+ md->nTCP = 0;
+ }
+
+ if (md->etc == etcNOSEHOOVER) {
+ if (md->bNHC_trotter) {
+ md->mde_n = 2*md->nNHC*md->nTC;
+ }
+ else
+ {
+ md->mde_n = 2*md->nTC;
+ }
+ if (md->epc == epcMTTK)
+ {
+ md->mdeb_n = 2*md->nNHC*md->nTCP;
+ }
+ } else {
+ md->mde_n = md->nTC;
+ md->mdeb_n = 0;
+ }
+
+ snew(md->tmp_r,md->mde_n);
+ snew(md->tmp_v,md->mde_n);
+ snew(md->grpnms,md->mde_n);
+ grpnms = md->grpnms;
+
+ for(i=0; (i<md->nTC); i++)
+ {
+ ni=groups->grps[egcTC].nm_ind[i];
+ sprintf(buf,"T-%s",*(groups->grpname[ni]));
+ grpnms[i]=strdup(buf);
+ }
+ md->itemp=get_ebin_space(md->ebin,md->nTC,(const char **)grpnms,
+ unit_temp_K);
+
+ bNoseHoover = (getenv("GMX_NOSEHOOVER_CHAINS") != NULL); /* whether to print Nose-Hoover chains */
+
+ if (md->etc == etcNOSEHOOVER)
+ {
+ if (bNoseHoover)
+ {
+ if (md->bNHC_trotter)
+ {
+ for(i=0; (i<md->nTC); i++)
+ {
+ ni=groups->grps[egcTC].nm_ind[i];
+ bufi = *(groups->grpname[ni]);
+ for(j=0; (j<md->nNHC); j++)
+ {
+ sprintf(buf,"Xi-%d-%s",j,bufi);
+ grpnms[2*(i*md->nNHC+j)]=strdup(buf);
+ sprintf(buf,"vXi-%d-%s",j,bufi);
+ grpnms[2*(i*md->nNHC+j)+1]=strdup(buf);
+ }
+ }
+ md->itc=get_ebin_space(md->ebin,md->mde_n,
+ (const char **)grpnms,unit_invtime);
+ if (md->bMTTK)
+ {
+ for(i=0; (i<md->nTCP); i++)
+ {
+ bufi = baro_nm[0]; /* All barostat DOF's together for now. */
+ for(j=0; (j<md->nNHC); j++)
+ {
+ sprintf(buf,"Xi-%d-%s",j,bufi);
+ grpnms[2*(i*md->nNHC+j)]=strdup(buf);
+ sprintf(buf,"vXi-%d-%s",j,bufi);
+ grpnms[2*(i*md->nNHC+j)+1]=strdup(buf);
+ }
+ }
+ md->itcb=get_ebin_space(md->ebin,md->mdeb_n,
+ (const char **)grpnms,unit_invtime);
+ }
+ }
+ else
+ {
+ for(i=0; (i<md->nTC); i++)
+ {
+ ni=groups->grps[egcTC].nm_ind[i];
+ bufi = *(groups->grpname[ni]);
+ sprintf(buf,"Xi-%s",bufi);
+ grpnms[2*i]=strdup(buf);
+ sprintf(buf,"vXi-%s",bufi);
+ grpnms[2*i+1]=strdup(buf);
+ }
+ md->itc=get_ebin_space(md->ebin,md->mde_n,
+ (const char **)grpnms,unit_invtime);
+ }
+ }
+ }
+ else if (md->etc == etcBERENDSEN || md->etc == etcYES ||
+ md->etc == etcVRESCALE)
+ {
+ for(i=0; (i<md->nTC); i++)
+ {
+ ni=groups->grps[egcTC].nm_ind[i];
+ sprintf(buf,"Lamb-%s",*(groups->grpname[ni]));
+ grpnms[i]=strdup(buf);
+ }
+ md->itc=get_ebin_space(md->ebin,md->mde_n,(const char **)grpnms,"");
+ }
+
+ sfree(grpnms);
+
+
+ md->nU=groups->grps[egcACC].nr;
+ if (md->nU > 1)
+ {
+ snew(grpnms,3*md->nU);
+ for(i=0; (i<md->nU); i++)
+ {
+ ni=groups->grps[egcACC].nm_ind[i];
+ sprintf(buf,"Ux-%s",*(groups->grpname[ni]));
+ grpnms[3*i+XX]=strdup(buf);
+ sprintf(buf,"Uy-%s",*(groups->grpname[ni]));
+ grpnms[3*i+YY]=strdup(buf);
+ sprintf(buf,"Uz-%s",*(groups->grpname[ni]));
+ grpnms[3*i+ZZ]=strdup(buf);
+ }
+ md->iu=get_ebin_space(md->ebin,3*md->nU,(const char **)grpnms,unit_vel);
+ sfree(grpnms);
+ }
+
+ if ( fp_ene )
+ {
+ do_enxnms(fp_ene,&md->ebin->nener,&md->ebin->enm);
+ }
+
+ md->print_grpnms=NULL;
+
+ /* check whether we're going to write dh histograms */
+ md->dhc=NULL;
+ if (ir->separate_dhdl_file == sepdhdlfileNO )
+ {
+ int i;
+ snew(md->dhc, 1);
+
+ mde_delta_h_coll_init(md->dhc, ir);
+ md->fp_dhdl = NULL;
+ }
+ else
+ {
+ md->fp_dhdl = fp_dhdl;
+ }
+ md->dhdl_derivatives = (ir->dhdl_derivatives==dhdlderivativesYES);
+ return md;
+}
+
+FILE *open_dhdl(const char *filename,const t_inputrec *ir,
+ const output_env_t oenv)
+{
+ FILE *fp;
+ const char *dhdl="dH/d\\lambda",*deltag="\\DeltaH",*lambda="\\lambda";
+ char title[STRLEN],label_x[STRLEN],label_y[STRLEN];
+ char **setname;
+ char buf[STRLEN];
+
+ sprintf(label_x,"%s (%s)","Time",unit_time);
+ if (ir->n_flambda == 0)
+ {
+ sprintf(title,"%s",dhdl);
+ sprintf(label_y,"%s (%s %s)",
+ dhdl,unit_energy,"[\\lambda]\\S-1\\N");
+ }
+ else
+ {
+ sprintf(title,"%s, %s",dhdl,deltag);
+ sprintf(label_y,"(%s)",unit_energy);
+ }
+ fp = gmx_fio_fopen(filename,"w+");
+ xvgr_header(fp,title,label_x,label_y,exvggtXNY,oenv);
+
+ if (ir->delta_lambda == 0)
+ {
+ sprintf(buf,"T = %g (K), %s = %g",
+ ir->opts.ref_t[0],lambda,ir->init_lambda);
+ }
+ else
+ {
+ sprintf(buf,"T = %g (K)",
+ ir->opts.ref_t[0]);
+ }
+ xvgr_subtitle(fp,buf,oenv);
+
+ if (ir->n_flambda > 0)
+ {
+ int nsets,s,nsi=0;
+ /* g_bar has to determine the lambda values used in this simulation
+ * from this xvg legend. */
+ nsets = ( (ir->dhdl_derivatives==dhdlderivativesYES) ? 1 : 0) +
+ ir->n_flambda;
+ snew(setname,nsets);
+ if (ir->dhdl_derivatives == dhdlderivativesYES)
+ {
+ sprintf(buf,"%s %s %g",dhdl,lambda,ir->init_lambda);
+ setname[nsi++] = strdup(buf);
+ }
+ for(s=0; s<ir->n_flambda; s++)
+ {
+ sprintf(buf,"%s %s %g",deltag,lambda,ir->flambda[s]);
+ setname[nsi++] = strdup(buf);
+ }
+ xvgr_legend(fp,nsets,(const char**)setname,oenv);
+
+ for(s=0; s<nsets; s++)
+ {
+ sfree(setname[s]);
+ }
+ sfree(setname);
+ }
+
+ return fp;
+}
+
+static void copy_energy(t_mdebin *md, real e[],real ecpy[])
+{
+ int i,j;
+
+ for(i=j=0; (i<F_NRE); i++)
+ if (md->bEner[i])
+ ecpy[j++] = e[i];
+ if (j != md->f_nre)
+ gmx_incons("Number of energy terms wrong");
+}
+
+void upd_mdebin(t_mdebin *md, gmx_bool write_dhdl,
+ gmx_bool bSum,
+ double time,
+ real tmass,
+ gmx_enerdata_t *enerd,
+ t_state *state,
+ matrix box,
+ tensor svir,
+ tensor fvir,
+ tensor vir,
+ tensor pres,
+ gmx_ekindata_t *ekind,
+ rvec mu_tot,
+ gmx_constr_t constr)
+{
+ int i,j,k,kk,m,n,gid;
+ real crmsd[2],tmp6[6];
+ real bs[NTRICLBOXS],vol,dens,pv,enthalpy;
+ real eee[egNR];
+ real ecopy[F_NRE];
+ real tmp;
+ gmx_bool bNoseHoover;
+
+ /* Do NOT use the box in the state variable, but the separate box provided
+ * as an argument. This is because we sometimes need to write the box from
+ * the last timestep to match the trajectory frames.
+ */
+ copy_energy(md, enerd->term,ecopy);
+ add_ebin(md->ebin,md->ie,md->f_nre,ecopy,bSum);
+ if (md->nCrmsd)
+ {
+ crmsd[0] = constr_rmsd(constr,FALSE);
+ if (md->nCrmsd > 1)
+ {
+ crmsd[1] = constr_rmsd(constr,TRUE);
+ }
+ add_ebin(md->ebin,md->iconrmsd,md->nCrmsd,crmsd,FALSE);
+ }
+ if (md->bDynBox)
+ {
+ int nboxs;
+ if(md->bTricl)
+ {
+ bs[0] = box[XX][XX];
+ bs[1] = box[YY][YY];
+ bs[2] = box[ZZ][ZZ];
+ bs[3] = box[YY][XX];
+ bs[4] = box[ZZ][XX];
+ bs[5] = box[ZZ][YY];
+ nboxs=NTRICLBOXS;
+ }
+ else
+ {
+ bs[0] = box[XX][XX];
+ bs[1] = box[YY][YY];
+ bs[2] = box[ZZ][ZZ];
+ nboxs=NBOXS;
+ }
+ vol = box[XX][XX]*box[YY][YY]*box[ZZ][ZZ];
+ dens = (tmass*AMU)/(vol*NANO*NANO*NANO);
+
+ /* This is pV (in kJ/mol). The pressure is the reference pressure,
+ not the instantaneous pressure */
+ pv = 0;
+ for (i=0;i<DIM;i++)
+ {
+ for (j=0;j<DIM;j++)
+ {
+ if (i>j)
+ {
+ pv += box[i][j]*md->ref_p[i][j]/PRESFAC;
+ }
+ else
+ {
+ pv += box[j][i]*md->ref_p[j][i]/PRESFAC;
+ }
+ }
+ }
+
+ add_ebin(md->ebin,md->ib ,nboxs,bs ,bSum);
+ add_ebin(md->ebin,md->ivol ,1 ,&vol ,bSum);
+ add_ebin(md->ebin,md->idens,1 ,&dens,bSum);
+ add_ebin(md->ebin,md->ipv ,1 ,&pv ,bSum);
+ enthalpy = pv + enerd->term[F_ETOT];
+ add_ebin(md->ebin,md->ienthalpy ,1 ,&enthalpy ,bSum);
+ }
+ if (md->bConstrVir)
+ {
+ add_ebin(md->ebin,md->isvir,9,svir[0],bSum);
+ add_ebin(md->ebin,md->ifvir,9,fvir[0],bSum);
+ }
+ add_ebin(md->ebin,md->ivir,9,vir[0],bSum);
+ add_ebin(md->ebin,md->ipres,9,pres[0],bSum);
+ tmp = (pres[ZZ][ZZ]-(pres[XX][XX]+pres[YY][YY])*0.5)*box[ZZ][ZZ];
+ add_ebin(md->ebin,md->isurft,1,&tmp,bSum);
+ if (md->epc == epcPARRINELLORAHMAN || md->epc == epcMTTK)
+ {
+ tmp6[0] = state->boxv[XX][XX];
+ tmp6[1] = state->boxv[YY][YY];
+ tmp6[2] = state->boxv[ZZ][ZZ];
+ tmp6[3] = state->boxv[YY][XX];
+ tmp6[4] = state->boxv[ZZ][XX];
+ tmp6[5] = state->boxv[ZZ][YY];
+ add_ebin(md->ebin,md->ipc,md->bTricl ? 6 : 3,tmp6,bSum);
+ }
+ add_ebin(md->ebin,md->imu,3,mu_tot,bSum);
+ if (ekind && ekind->cosacc.cos_accel != 0)
+ {
+ vol = box[XX][XX]*box[YY][YY]*box[ZZ][ZZ];
+ dens = (tmass*AMU)/(vol*NANO*NANO*NANO);
+ add_ebin(md->ebin,md->ivcos,1,&(ekind->cosacc.vcos),bSum);
+ /* 1/viscosity, unit 1/(kg m^-1 s^-1) */
+ tmp = 1/(ekind->cosacc.cos_accel/(ekind->cosacc.vcos*PICO)
++ *dens*vol*sqr(box[ZZ][ZZ]*NANO/(2*M_PI)));
+ add_ebin(md->ebin,md->ivisc,1,&tmp,bSum);
+ }
+ if (md->nE > 1)
+ {
+ n=0;
+ for(i=0; (i<md->nEg); i++)
+ {
+ for(j=i; (j<md->nEg); j++)
+ {
+ gid=GID(i,j,md->nEg);
+ for(k=kk=0; (k<egNR); k++)
+ {
+ if (md->bEInd[k])
+ {
+ eee[kk++] = enerd->grpp.ener[k][gid];
+ }
+ }
+ add_ebin(md->ebin,md->igrp[n],md->nEc,eee,bSum);
+ n++;
+ }
+ }
+ }
+
+ if (ekind)
+ {
+ for(i=0; (i<md->nTC); i++)
+ {
+ md->tmp_r[i] = ekind->tcstat[i].T;
+ }
+ add_ebin(md->ebin,md->itemp,md->nTC,md->tmp_r,bSum);
+
+ /* whether to print Nose-Hoover chains: */
+ bNoseHoover = (getenv("GMX_NOSEHOOVER_CHAINS") != NULL);
+
+ if (md->etc == etcNOSEHOOVER)
+ {
+ if (bNoseHoover)
+ {
+ if (md->bNHC_trotter)
+ {
+ for(i=0; (i<md->nTC); i++)
+ {
+ for (j=0;j<md->nNHC;j++)
+ {
+ k = i*md->nNHC+j;
+ md->tmp_r[2*k] = state->nosehoover_xi[k];
+ md->tmp_r[2*k+1] = state->nosehoover_vxi[k];
+ }
+ }
+ add_ebin(md->ebin,md->itc,md->mde_n,md->tmp_r,bSum);
+
+ if (md->bMTTK) {
+ for(i=0; (i<md->nTCP); i++)
+ {
+ for (j=0;j<md->nNHC;j++)
+ {
+ k = i*md->nNHC+j;
+ md->tmp_r[2*k] = state->nhpres_xi[k];
+ md->tmp_r[2*k+1] = state->nhpres_vxi[k];
+ }
+ }
+ add_ebin(md->ebin,md->itcb,md->mdeb_n,md->tmp_r,bSum);
+ }
+
+ }
+ else
+ {
+ for(i=0; (i<md->nTC); i++)
+ {
+ md->tmp_r[2*i] = state->nosehoover_xi[i];
+ md->tmp_r[2*i+1] = state->nosehoover_vxi[i];
+ }
+ add_ebin(md->ebin,md->itc,md->mde_n,md->tmp_r,bSum);
+ }
+ }
+ }
+ else if (md->etc == etcBERENDSEN || md->etc == etcYES ||
+ md->etc == etcVRESCALE)
+ {
+ for(i=0; (i<md->nTC); i++)
+ {
+ md->tmp_r[i] = ekind->tcstat[i].lambda;
+ }
+ add_ebin(md->ebin,md->itc,md->nTC,md->tmp_r,bSum);
+ }
+ }
+
+ if (ekind && md->nU > 1)
+ {
+ for(i=0; (i<md->nU); i++)
+ {
+ copy_rvec(ekind->grpstat[i].u,md->tmp_v[i]);
+ }
+ add_ebin(md->ebin,md->iu,3*md->nU,md->tmp_v[0],bSum);
+ }
+
+ ebin_increase_count(md->ebin,bSum);
+
+ /* BAR + thermodynamic integration values */
+ if (write_dhdl)
+ {
+ if (md->fp_dhdl)
+ {
+ fprintf(md->fp_dhdl,"%.4f", time);
+
+ if (md->dhdl_derivatives)
+ {
+ fprintf(md->fp_dhdl," %g", enerd->term[F_DVDL]+
+ enerd->term[F_DKDL]+
+ enerd->term[F_DHDL_CON]);
+ }
+ for(i=1; i<enerd->n_lambda; i++)
+ {
+ fprintf(md->fp_dhdl," %g",
+ enerd->enerpart_lambda[i]-enerd->enerpart_lambda[0]);
+ }
+ fprintf(md->fp_dhdl,"\n");
+ }
+ /* and the binary BAR output */
+ if (md->dhc)
+ {
+ mde_delta_h_coll_add_dh(md->dhc,
+ enerd->term[F_DVDL]+ enerd->term[F_DKDL]+
+ enerd->term[F_DHDL_CON],
+ enerd->enerpart_lambda, time,
+ state->lambda);
+ }
+ }
+}
+
+void upd_mdebin_step(t_mdebin *md)
+{
+ ebin_increase_count(md->ebin,FALSE);
+}
+
+static void npr(FILE *log,int n,char c)
+{
+ for(; (n>0); n--) fprintf(log,"%c",c);
+}
+
+static void pprint(FILE *log,const char *s,t_mdebin *md)
+{
+ char CHAR='#';
+ int slen;
+ char buf1[22],buf2[22];
+
+ slen = strlen(s);
+ fprintf(log,"\t<====== ");
+ npr(log,slen,CHAR);
+ fprintf(log," ==>\n");
+ fprintf(log,"\t<==== %s ====>\n",s);
+ fprintf(log,"\t<== ");
+ npr(log,slen,CHAR);
+ fprintf(log," ======>\n\n");
+
+ fprintf(log,"\tStatistics over %s steps using %s frames\n",
+ gmx_step_str(md->ebin->nsteps_sim,buf1),
+ gmx_step_str(md->ebin->nsum_sim,buf2));
+ fprintf(log,"\n");
+}
+
+void print_ebin_header(FILE *log,gmx_large_int_t steps,double time,real lamb)
+{
+ char buf[22];
+
+ fprintf(log," %12s %12s %12s\n"
+ " %12s %12.5f %12.5f\n\n",
+ "Step","Time","Lambda",gmx_step_str(steps,buf),time,lamb);
+}
+
+void print_ebin(ener_file_t fp_ene,gmx_bool bEne,gmx_bool bDR,gmx_bool bOR,
+ FILE *log,
+ gmx_large_int_t step,double time,
+ int mode,gmx_bool bCompact,
+ t_mdebin *md,t_fcdata *fcd,
+ gmx_groups_t *groups,t_grpopts *opts)
+{
+ /*static char **grpnms=NULL;*/
+ char buf[246];
+ int i,j,n,ni,nj,ndr,nor,b;
+ int ndisre=0;
+ real *disre_rm3tav, *disre_rt;
+
+ /* these are for the old-style blocks (1 subblock, only reals), because
+ there can be only one per ID for these */
+ int nr[enxNR];
+ int id[enxNR];
+ real *block[enxNR];
+
+ /* temporary arrays for the lambda values to write out */
+ double enxlambda_data[2];
+
+ t_enxframe fr;
+
+ switch (mode)
+ {
+ case eprNORMAL:
+ init_enxframe(&fr);
+ fr.t = time;
+ fr.step = step;
+ fr.nsteps = md->ebin->nsteps;
+ fr.dt = md->delta_t;
+ fr.nsum = md->ebin->nsum;
+ fr.nre = (bEne) ? md->ebin->nener : 0;
+ fr.ener = md->ebin->e;
+ ndisre = bDR ? fcd->disres.npair : 0;
+ disre_rm3tav = fcd->disres.rm3tav;
+ disre_rt = fcd->disres.rt;
+ /* Optional additional old-style (real-only) blocks. */
+ for(i=0; i<enxNR; i++)
+ {
+ nr[i] = 0;
+ }
+ if (fcd->orires.nr > 0 && bOR)
+ {
+ diagonalize_orires_tensors(&(fcd->orires));
+ nr[enxOR] = fcd->orires.nr;
+ block[enxOR] = fcd->orires.otav;
+ id[enxOR] = enxOR;
+ nr[enxORI] = (fcd->orires.oinsl != fcd->orires.otav) ?
+ fcd->orires.nr : 0;
+ block[enxORI] = fcd->orires.oinsl;
+ id[enxORI] = enxORI;
+ nr[enxORT] = fcd->orires.nex*12;
+ block[enxORT] = fcd->orires.eig;
+ id[enxORT] = enxORT;
+ }
+
+ /* whether we are going to wrte anything out: */
+ if (fr.nre || ndisre || nr[enxOR] || nr[enxORI])
+ {
+
+ /* the old-style blocks go first */
+ fr.nblock = 0;
+ for(i=0; i<enxNR; i++)
+ {
+ if (nr[i] > 0)
+ {
+ fr.nblock = i + 1;
+ }
+ }
+ add_blocks_enxframe(&fr, fr.nblock);
+ for(b=0;b<fr.nblock;b++)
+ {
+ add_subblocks_enxblock(&(fr.block[b]), 1);
+ fr.block[b].id=id[b];
+ fr.block[b].sub[0].nr = nr[b];
+#ifndef GMX_DOUBLE
+ fr.block[b].sub[0].type = xdr_datatype_float;
+ fr.block[b].sub[0].fval = block[b];
+#else
+ fr.block[b].sub[0].type = xdr_datatype_double;
+ fr.block[b].sub[0].dval = block[b];
+#endif
+ }
+
+ /* check for disre block & fill it. */
+ if (ndisre>0)
+ {
+ int db = fr.nblock;
+ fr.nblock+=1;
+ add_blocks_enxframe(&fr, fr.nblock);
+
+ add_subblocks_enxblock(&(fr.block[db]), 2);
+ fr.block[db].id=enxDISRE;
+ fr.block[db].sub[0].nr=ndisre;
+ fr.block[db].sub[1].nr=ndisre;
+#ifndef GMX_DOUBLE
+ fr.block[db].sub[0].type=xdr_datatype_float;
+ fr.block[db].sub[1].type=xdr_datatype_float;
+ fr.block[db].sub[0].fval=disre_rt;
+ fr.block[db].sub[1].fval=disre_rm3tav;
+#else
+ fr.block[db].sub[0].type=xdr_datatype_double;
+ fr.block[db].sub[1].type=xdr_datatype_double;
+ fr.block[db].sub[0].dval=disre_rt;
+ fr.block[db].sub[1].dval=disre_rm3tav;
+#endif
+ }
+ /* here we can put new-style blocks */
+
+ /* Free energy perturbation blocks */
+ if (md->dhc)
+ {
+ mde_delta_h_coll_handle_block(md->dhc, &fr, fr.nblock);
+ }
+
+ /* do the actual I/O */
+ do_enx(fp_ene,&fr);
+ gmx_fio_check_file_position(enx_file_pointer(fp_ene));
+ if (fr.nre)
+ {
+ /* We have stored the sums, so reset the sum history */
+ reset_ebin_sums(md->ebin);
+ }
+
+ /* we can now free & reset the data in the blocks */
+ if (md->dhc)
+ mde_delta_h_coll_reset(md->dhc);
+ }
+ free_enxframe(&fr);
+ break;
+ case eprAVER:
+ if (log)
+ {
+ pprint(log,"A V E R A G E S",md);
+ }
+ break;
+ case eprRMS:
+ if (log)
+ {
+ pprint(log,"R M S - F L U C T U A T I O N S",md);
+ }
+ break;
+ default:
+ gmx_fatal(FARGS,"Invalid print mode (%d)",mode);
+ }
+
+ if (log)
+ {
+ for(i=0;i<opts->ngtc;i++)
+ {
+ if(opts->annealing[i]!=eannNO)
+ {
+ fprintf(log,"Current ref_t for group %s: %8.1f\n",
+ *(groups->grpname[groups->grps[egcTC].nm_ind[i]]),
+ opts->ref_t[i]);
+ }
+ }
+ if (mode==eprNORMAL && fcd->orires.nr>0)
+ {
+ print_orires_log(log,&(fcd->orires));
+ }
+ fprintf(log," Energies (%s)\n",unit_energy);
+ pr_ebin(log,md->ebin,md->ie,md->f_nre+md->nCrmsd,5,mode,TRUE);
+ fprintf(log,"\n");
+
+ if (!bCompact)
+ {
+ if (md->bDynBox)
+ {
+ pr_ebin(log,md->ebin,md->ib, md->bTricl ? NTRICLBOXS : NBOXS,5,
+ mode,TRUE);
+ fprintf(log,"\n");
+ }
+ if (md->bConstrVir)
+ {
+ fprintf(log," Constraint Virial (%s)\n",unit_energy);
+ pr_ebin(log,md->ebin,md->isvir,9,3,mode,FALSE);
+ fprintf(log,"\n");
+ fprintf(log," Force Virial (%s)\n",unit_energy);
+ pr_ebin(log,md->ebin,md->ifvir,9,3,mode,FALSE);
+ fprintf(log,"\n");
+ }
+ fprintf(log," Total Virial (%s)\n",unit_energy);
+ pr_ebin(log,md->ebin,md->ivir,9,3,mode,FALSE);
+ fprintf(log,"\n");
+ fprintf(log," Pressure (%s)\n",unit_pres_bar);
+ pr_ebin(log,md->ebin,md->ipres,9,3,mode,FALSE);
+ fprintf(log,"\n");
+ fprintf(log," Total Dipole (%s)\n",unit_dipole_D);
+ pr_ebin(log,md->ebin,md->imu,3,3,mode,FALSE);
+ fprintf(log,"\n");
+
+ if (md->nE > 1)
+ {
+ if (md->print_grpnms==NULL)
+ {
+ snew(md->print_grpnms,md->nE);
+ n=0;
+ for(i=0; (i<md->nEg); i++)
+ {
+ ni=groups->grps[egcENER].nm_ind[i];
+ for(j=i; (j<md->nEg); j++)
+ {
+ nj=groups->grps[egcENER].nm_ind[j];
+ sprintf(buf,"%s-%s",*(groups->grpname[ni]),
+ *(groups->grpname[nj]));
+ md->print_grpnms[n++]=strdup(buf);
+ }
+ }
+ }
+ sprintf(buf,"Epot (%s)",unit_energy);
+ fprintf(log,"%15s ",buf);
+ for(i=0; (i<egNR); i++)
+ {
+ if (md->bEInd[i])
+ {
+ fprintf(log,"%12s ",egrp_nm[i]);
+ }
+ }
+ fprintf(log,"\n");
+ for(i=0; (i<md->nE); i++)
+ {
+ fprintf(log,"%15s",md->print_grpnms[i]);
+ pr_ebin(log,md->ebin,md->igrp[i],md->nEc,md->nEc,mode,
+ FALSE);
+ }
+ fprintf(log,"\n");
+ }
+ if (md->nTC > 1)
+ {
+ pr_ebin(log,md->ebin,md->itemp,md->nTC,4,mode,TRUE);
+ fprintf(log,"\n");
+ }
+ if (md->nU > 1)
+ {
+ fprintf(log,"%15s %12s %12s %12s\n",
+ "Group","Ux","Uy","Uz");
+ for(i=0; (i<md->nU); i++)
+ {
+ ni=groups->grps[egcACC].nm_ind[i];
+ fprintf(log,"%15s",*groups->grpname[ni]);
+ pr_ebin(log,md->ebin,md->iu+3*i,3,3,mode,FALSE);
+ }
+ fprintf(log,"\n");
+ }
+ }
+ }
+
+}
+
+void update_energyhistory(energyhistory_t * enerhist,t_mdebin * mdebin)
+{
+ int i;
+
+ enerhist->nsteps = mdebin->ebin->nsteps;
+ enerhist->nsum = mdebin->ebin->nsum;
+ enerhist->nsteps_sim = mdebin->ebin->nsteps_sim;
+ enerhist->nsum_sim = mdebin->ebin->nsum_sim;
+ enerhist->nener = mdebin->ebin->nener;
+
+ if (mdebin->ebin->nsum > 0)
+ {
+ /* Check if we need to allocate first */
+ if(enerhist->ener_ave == NULL)
+ {
+ snew(enerhist->ener_ave,enerhist->nener);
+ snew(enerhist->ener_sum,enerhist->nener);
+ }
+
+ for(i=0;i<enerhist->nener;i++)
+ {
+ enerhist->ener_ave[i] = mdebin->ebin->e[i].eav;
+ enerhist->ener_sum[i] = mdebin->ebin->e[i].esum;
+ }
+ }
+
+ if (mdebin->ebin->nsum_sim > 0)
+ {
+ /* Check if we need to allocate first */
+ if(enerhist->ener_sum_sim == NULL)
+ {
+ snew(enerhist->ener_sum_sim,enerhist->nener);
+ }
+
+ for(i=0;i<enerhist->nener;i++)
+ {
+ enerhist->ener_sum_sim[i] = mdebin->ebin->e_sim[i].esum;
+ }
+ }
+ if (mdebin->dhc)
+ {
+ mde_delta_h_coll_update_energyhistory(mdebin->dhc, enerhist);
+ }
+}
+
+void restore_energyhistory_from_state(t_mdebin * mdebin,
+ energyhistory_t * enerhist)
+{
+ int i;
+
+ if ((enerhist->nsum > 0 || enerhist->nsum_sim > 0) &&
+ mdebin->ebin->nener != enerhist->nener)
+ {
+ gmx_fatal(FARGS,"Mismatch between number of energies in run input (%d) and checkpoint file (%d).",
+ mdebin->ebin->nener,enerhist->nener);
+ }
+
+ mdebin->ebin->nsteps = enerhist->nsteps;
+ mdebin->ebin->nsum = enerhist->nsum;
+ mdebin->ebin->nsteps_sim = enerhist->nsteps_sim;
+ mdebin->ebin->nsum_sim = enerhist->nsum_sim;
+
+ for(i=0; i<mdebin->ebin->nener; i++)
+ {
+ mdebin->ebin->e[i].eav =
+ (enerhist->nsum > 0 ? enerhist->ener_ave[i] : 0);
+ mdebin->ebin->e[i].esum =
+ (enerhist->nsum > 0 ? enerhist->ener_sum[i] : 0);
+ mdebin->ebin->e_sim[i].esum =
+ (enerhist->nsum_sim > 0 ? enerhist->ener_sum_sim[i] : 0);
+ }
+ if (mdebin->dhc)
+ {
+ mde_delta_h_coll_restore_energyhistory(mdebin->dhc, enerhist);
+ }
+}
--- /dev/null
- struct timezone tz = { 0,0 };
+/* -*- 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:
+ * GROwing Monsters And Cloning Shrimps
+ */
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#ifdef GMX_CRAY_XT3
+#include<catamount/dclock.h>
+#endif
+
+
+#include <stdio.h>
+#include <time.h>
+#ifdef HAVE_SYS_TIME_H
+#include <sys/time.h>
+#endif
+#include <math.h>
+#include "typedefs.h"
+#include "string2.h"
+#include "gmxfio.h"
+#include "smalloc.h"
+#include "names.h"
+#include "confio.h"
+#include "mvdata.h"
+#include "txtdump.h"
+#include "pbc.h"
+#include "chargegroup.h"
+#include "vec.h"
+#include "time.h"
+#include "nrnb.h"
+#include "mshift.h"
+#include "mdrun.h"
+#include "update.h"
+#include "physics.h"
+#include "main.h"
+#include "mdatoms.h"
+#include "force.h"
+#include "bondf.h"
+#include "pme.h"
+#include "pppm.h"
+#include "disre.h"
+#include "orires.h"
+#include "network.h"
+#include "calcmu.h"
+#include "constr.h"
+#include "xvgr.h"
+#include "trnio.h"
+#include "xtcio.h"
+#include "copyrite.h"
+#include "pull_rotation.h"
+#include "mpelogging.h"
+#include "domdec.h"
+#include "partdec.h"
+#include "gmx_wallcycle.h"
+#include "genborn.h"
+
+#ifdef GMX_LIB_MPI
+#include <mpi.h>
+#endif
+#ifdef GMX_THREADS
+#include "tmpi.h"
+#endif
+
+#include "qmmm.h"
+
+#if 0
+typedef struct gmx_timeprint {
+
+} t_gmx_timeprint;
+#endif
+
+/* Portable version of ctime_r implemented in src/gmxlib/string2.c, but we do not want it declared in public installed headers */
+char *
+gmx_ctime_r(const time_t *clock,char *buf, int n);
+
+
+double
+gmx_gettime()
+{
+#ifdef HAVE_GETTIMEOFDAY
+ struct timeval t;
- gettimeofday(&t,&tz);
+ double seconds;
+
++ gettimeofday(&t,NULL);
+
+ seconds = (double) t.tv_sec + 1e-6*(double)t.tv_usec;
+
+ return seconds;
+#else
+ double seconds;
+
+ seconds = time(NULL);
+
+ return seconds;
+#endif
+}
+
+
+#define difftime(end,start) ((double)(end)-(double)(start))
+
+void print_time(FILE *out,gmx_runtime_t *runtime,gmx_large_int_t step,
+ t_inputrec *ir, t_commrec *cr)
+{
+ time_t finish;
+ char timebuf[STRLEN];
+ double dt;
+ char buf[48];
+
+#ifndef GMX_THREADS
+ if (!PAR(cr))
+#endif
+ {
+ fprintf(out,"\r");
+ }
+ fprintf(out,"step %s",gmx_step_str(step,buf));
+ if ((step >= ir->nstlist))
+ {
+ if ((ir->nstlist == 0) || ((step % ir->nstlist) == 0))
+ {
+ /* We have done a full cycle let's update time_per_step */
+ runtime->last = gmx_gettime();
+ dt = difftime(runtime->last,runtime->real);
+ runtime->time_per_step = dt/(step - ir->init_step + 1);
+ }
+ dt = (ir->nsteps + ir->init_step - step)*runtime->time_per_step;
+
+ if (ir->nsteps >= 0)
+ {
+ if (dt >= 300)
+ {
+ finish = (time_t) (runtime->last + dt);
+ gmx_ctime_r(&finish,timebuf,STRLEN);
+ sprintf(buf,"%s",timebuf);
+ buf[strlen(buf)-1]='\0';
+ fprintf(out,", will finish %s",buf);
+ }
+ else
+ fprintf(out,", remaining runtime: %5d s ",(int)dt);
+ }
+ else
+ {
+ fprintf(out," performance: %.1f ns/day ",
+ ir->delta_t/1000*24*60*60/runtime->time_per_step);
+ }
+ }
+#ifndef GMX_THREADS
+ if (PAR(cr))
+ {
+ fprintf(out,"\n");
+ }
+#endif
+
+ fflush(out);
+}
+
+#ifdef NO_CLOCK
+#define clock() -1
+#endif
+
+static double set_proctime(gmx_runtime_t *runtime)
+{
+ double diff;
+#ifdef GMX_CRAY_XT3
+ double prev;
+
+ prev = runtime->proc;
+ runtime->proc = dclock();
+
+ diff = runtime->proc - prev;
+#else
+ clock_t prev;
+
+ prev = runtime->proc;
+ runtime->proc = clock();
+
+ diff = (double)(runtime->proc - prev)/(double)CLOCKS_PER_SEC;
+#endif
+ if (diff < 0)
+ {
+ /* The counter has probably looped, ignore this data */
+ diff = 0;
+ }
+
+ return diff;
+}
+
+void runtime_start(gmx_runtime_t *runtime)
+{
+ runtime->real = gmx_gettime();
+ runtime->proc = 0;
+ set_proctime(runtime);
+ runtime->realtime = 0;
+ runtime->proctime = 0;
+ runtime->last = 0;
+ runtime->time_per_step = 0;
+}
+
+void runtime_end(gmx_runtime_t *runtime)
+{
+ double now;
+
+ now = gmx_gettime();
+
+ runtime->proctime += set_proctime(runtime);
+ runtime->realtime = now - runtime->real;
+ runtime->real = now;
+}
+
+void runtime_upd_proc(gmx_runtime_t *runtime)
+{
+ runtime->proctime += set_proctime(runtime);
+}
+
+void print_date_and_time(FILE *fplog,int nodeid,const char *title,
+ const gmx_runtime_t *runtime)
+{
+ int i;
+ char timebuf[STRLEN];
+ char time_string[STRLEN];
+ time_t tmptime;
+
+ if (fplog)
+ {
+ if (runtime != NULL)
+ {
+ tmptime = (time_t) runtime->real;
+ gmx_ctime_r(&tmptime,timebuf,STRLEN);
+ }
+ else
+ {
+ tmptime = (time_t) gmx_gettime();
+ gmx_ctime_r(&tmptime,timebuf,STRLEN);
+ }
+ for(i=0; timebuf[i]>=' '; i++)
+ {
+ time_string[i]=timebuf[i];
+ }
+ time_string[i]='\0';
+
+ fprintf(fplog,"%s on node %d %s\n",title,nodeid,time_string);
+ }
+}
+
+static void sum_forces(int start,int end,rvec f[],rvec flr[])
+{
+ int i;
+
+ if (gmx_debug_at) {
+ pr_rvecs(debug,0,"fsr",f+start,end-start);
+ pr_rvecs(debug,0,"flr",flr+start,end-start);
+ }
+ for(i=start; (i<end); i++)
+ rvec_inc(f[i],flr[i]);
+}
+
+/*
+ * calc_f_el calculates forces due to an electric field.
+ *
+ * force is kJ mol^-1 nm^-1 = e * kJ mol^-1 nm^-1 / e
+ *
+ * Et[] contains the parameters for the time dependent
+ * part of the field (not yet used).
+ * Ex[] contains the parameters for
+ * the spatial dependent part of the field. You can have cool periodic
+ * fields in principle, but only a constant field is supported
+ * now.
+ * The function should return the energy due to the electric field
+ * (if any) but for now returns 0.
+ *
+ * WARNING:
+ * There can be problems with the virial.
+ * Since the field is not self-consistent this is unavoidable.
+ * For neutral molecules the virial is correct within this approximation.
+ * For neutral systems with many charged molecules the error is small.
+ * But for systems with a net charge or a few charged molecules
+ * the error can be significant when the field is high.
+ * Solution: implement a self-consitent electric field into PME.
+ */
+static void calc_f_el(FILE *fp,int start,int homenr,
+ real charge[],rvec x[],rvec f[],
+ t_cosines Ex[],t_cosines Et[],double t)
+{
+ rvec Ext;
+ real t0;
+ int i,m;
+
+ for(m=0; (m<DIM); m++)
+ {
+ if (Et[m].n > 0)
+ {
+ if (Et[m].n == 3)
+ {
+ t0 = Et[m].a[1];
+ Ext[m] = cos(Et[m].a[0]*(t-t0))*exp(-sqr(t-t0)/(2.0*sqr(Et[m].a[2])));
+ }
+ else
+ {
+ Ext[m] = cos(Et[m].a[0]*t);
+ }
+ }
+ else
+ {
+ Ext[m] = 1.0;
+ }
+ if (Ex[m].n > 0)
+ {
+ /* Convert the field strength from V/nm to MD-units */
+ Ext[m] *= Ex[m].a[0]*FIELDFAC;
+ for(i=start; (i<start+homenr); i++)
+ f[i][m] += charge[i]*Ext[m];
+ }
+ else
+ {
+ Ext[m] = 0;
+ }
+ }
+ if (fp != NULL)
+ {
+ fprintf(fp,"%10g %10g %10g %10g #FIELD\n",t,
+ Ext[XX]/FIELDFAC,Ext[YY]/FIELDFAC,Ext[ZZ]/FIELDFAC);
+ }
+}
+
+static void calc_virial(FILE *fplog,int start,int homenr,rvec x[],rvec f[],
+ tensor vir_part,t_graph *graph,matrix box,
+ t_nrnb *nrnb,const t_forcerec *fr,int ePBC)
+{
+ int i,j;
+ tensor virtest;
+
+ /* The short-range virial from surrounding boxes */
+ clear_mat(vir_part);
+ calc_vir(fplog,SHIFTS,fr->shift_vec,fr->fshift,vir_part,ePBC==epbcSCREW,box);
+ inc_nrnb(nrnb,eNR_VIRIAL,SHIFTS);
+
+ /* Calculate partial virial, for local atoms only, based on short range.
+ * Total virial is computed in global_stat, called from do_md
+ */
+ f_calc_vir(fplog,start,start+homenr,x,f,vir_part,graph,box);
+ inc_nrnb(nrnb,eNR_VIRIAL,homenr);
+
+ /* Add position restraint contribution */
+ for(i=0; i<DIM; i++) {
+ vir_part[i][i] += fr->vir_diag_posres[i];
+ }
+
+ /* Add wall contribution */
+ for(i=0; i<DIM; i++) {
+ vir_part[i][ZZ] += fr->vir_wall_z[i];
+ }
+
+ if (debug)
+ pr_rvecs(debug,0,"vir_part",vir_part,DIM);
+}
+
+static void print_large_forces(FILE *fp,t_mdatoms *md,t_commrec *cr,
+ gmx_large_int_t step,real pforce,rvec *x,rvec *f)
+{
+ int i;
+ real pf2,fn2;
+ char buf[STEPSTRSIZE];
+
+ pf2 = sqr(pforce);
+ for(i=md->start; i<md->start+md->homenr; i++) {
+ fn2 = norm2(f[i]);
+ /* We also catch NAN, if the compiler does not optimize this away. */
+ if (fn2 >= pf2 || fn2 != fn2) {
+ fprintf(fp,"step %s atom %6d x %8.3f %8.3f %8.3f force %12.5e\n",
+ gmx_step_str(step,buf),
+ ddglatnr(cr->dd,i),x[i][XX],x[i][YY],x[i][ZZ],sqrt(fn2));
+ }
+ }
+}
+
+void do_force(FILE *fplog,t_commrec *cr,
+ t_inputrec *inputrec,
+ gmx_large_int_t step,t_nrnb *nrnb,gmx_wallcycle_t wcycle,
+ gmx_localtop_t *top,
+ gmx_mtop_t *mtop,
+ gmx_groups_t *groups,
+ matrix box,rvec x[],history_t *hist,
+ rvec f[],
+ tensor vir_force,
+ t_mdatoms *mdatoms,
+ gmx_enerdata_t *enerd,t_fcdata *fcd,
+ real lambda,t_graph *graph,
+ t_forcerec *fr,gmx_vsite_t *vsite,rvec mu_tot,
+ double t,FILE *field,gmx_edsam_t ed,
+ gmx_bool bBornRadii,
+ int flags)
+{
+ int cg0,cg1,i,j;
+ int start,homenr;
+ double mu[2*DIM];
+ gmx_bool bSepDVDL,bStateChanged,bNS,bFillGrid,bCalcCGCM,bBS;
+ gmx_bool bDoLongRange,bDoForces,bSepLRF;
+ matrix boxs;
+ real e,v,dvdl;
+ t_pbc pbc;
+ float cycles_ppdpme,cycles_pme,cycles_seppme,cycles_force;
+
+ start = mdatoms->start;
+ homenr = mdatoms->homenr;
+
+ bSepDVDL = (fr->bSepDVDL && do_per_step(step,inputrec->nstlog));
+
+ clear_mat(vir_force);
+
+ if (PARTDECOMP(cr))
+ {
+ pd_cg_range(cr,&cg0,&cg1);
+ }
+ else
+ {
+ cg0 = 0;
+ if (DOMAINDECOMP(cr))
+ {
+ cg1 = cr->dd->ncg_tot;
+ }
+ else
+ {
+ cg1 = top->cgs.nr;
+ }
+ if (fr->n_tpi > 0)
+ {
+ cg1--;
+ }
+ }
+
+ bStateChanged = (flags & GMX_FORCE_STATECHANGED);
+ bNS = (flags & GMX_FORCE_NS) && (fr->bAllvsAll==FALSE);
+ bFillGrid = (bNS && bStateChanged);
+ bCalcCGCM = (bFillGrid && !DOMAINDECOMP(cr));
+ bDoLongRange = (fr->bTwinRange && bNS && (flags & GMX_FORCE_DOLR));
+ bDoForces = (flags & GMX_FORCE_FORCES);
+ bSepLRF = (bDoLongRange && bDoForces && (flags & GMX_FORCE_SEPLRF));
+
+ if (bStateChanged)
+ {
+ update_forcerec(fplog,fr,box);
+
+ /* Calculate total (local) dipole moment in a temporary common array.
+ * This makes it possible to sum them over nodes faster.
+ */
+ calc_mu(start,homenr,
+ x,mdatoms->chargeA,mdatoms->chargeB,mdatoms->nChargePerturbed,
+ mu,mu+DIM);
+ }
+
+ if (fr->ePBC != epbcNONE) {
+ /* Compute shift vectors every step,
+ * because of pressure coupling or box deformation!
+ */
+ if ((flags & GMX_FORCE_DYNAMICBOX) && bStateChanged)
+ calc_shifts(box,fr->shift_vec);
+
+ if (bCalcCGCM) {
+ put_charge_groups_in_box(fplog,cg0,cg1,fr->ePBC,box,
+ &(top->cgs),x,fr->cg_cm);
+ inc_nrnb(nrnb,eNR_CGCM,homenr);
+ inc_nrnb(nrnb,eNR_RESETX,cg1-cg0);
+ }
+ else if (EI_ENERGY_MINIMIZATION(inputrec->eI) && graph) {
+ unshift_self(graph,box,x);
+ }
+ }
+ else if (bCalcCGCM) {
+ calc_cgcm(fplog,cg0,cg1,&(top->cgs),x,fr->cg_cm);
+ inc_nrnb(nrnb,eNR_CGCM,homenr);
+ }
+
+ if (bCalcCGCM) {
+ if (PAR(cr)) {
+ move_cgcm(fplog,cr,fr->cg_cm);
+ }
+ if (gmx_debug_at)
+ pr_rvecs(debug,0,"cgcm",fr->cg_cm,top->cgs.nr);
+ }
+
+#ifdef GMX_MPI
+ if (!(cr->duty & DUTY_PME)) {
+ /* Send particle coordinates to the pme nodes.
+ * Since this is only implemented for domain decomposition
+ * and domain decomposition does not use the graph,
+ * we do not need to worry about shifting.
+ */
+
+ wallcycle_start(wcycle,ewcPP_PMESENDX);
+ GMX_MPE_LOG(ev_send_coordinates_start);
+
+ bBS = (inputrec->nwall == 2);
+ if (bBS) {
+ copy_mat(box,boxs);
+ svmul(inputrec->wall_ewald_zfac,boxs[ZZ],boxs[ZZ]);
+ }
+
+ gmx_pme_send_x(cr,bBS ? boxs : box,x,
+ mdatoms->nChargePerturbed,lambda,
+ ( flags & GMX_FORCE_VIRIAL),step);
+
+ GMX_MPE_LOG(ev_send_coordinates_finish);
+ wallcycle_stop(wcycle,ewcPP_PMESENDX);
+ }
+#endif /* GMX_MPI */
+
+ /* Communicate coordinates and sum dipole if necessary */
+ if (PAR(cr))
+ {
+ wallcycle_start(wcycle,ewcMOVEX);
+ if (DOMAINDECOMP(cr))
+ {
+ dd_move_x(cr->dd,box,x);
+ }
+ else
+ {
+ move_x(fplog,cr,GMX_LEFT,GMX_RIGHT,x,nrnb);
+ }
+ /* When we don't need the total dipole we sum it in global_stat */
+ if (bStateChanged && NEED_MUTOT(*inputrec))
+ {
+ gmx_sumd(2*DIM,mu,cr);
+ }
+ wallcycle_stop(wcycle,ewcMOVEX);
+ }
+ if (bStateChanged)
+ {
+ for(i=0; i<2; i++)
+ {
+ for(j=0;j<DIM;j++)
+ {
+ fr->mu_tot[i][j] = mu[i*DIM + j];
+ }
+ }
+ }
+ if (fr->efep == efepNO)
+ {
+ copy_rvec(fr->mu_tot[0],mu_tot);
+ }
+ else
+ {
+ for(j=0; j<DIM; j++)
+ {
+ mu_tot[j] =
+ (1.0 - lambda)*fr->mu_tot[0][j] + lambda*fr->mu_tot[1][j];
+ }
+ }
+
+ /* Reset energies */
+ reset_enerdata(&(inputrec->opts),fr,bNS,enerd,MASTER(cr));
+ clear_rvecs(SHIFTS,fr->fshift);
+
+ if (bNS)
+ {
+ wallcycle_start(wcycle,ewcNS);
+
+ if (graph && bStateChanged)
+ {
+ /* Calculate intramolecular shift vectors to make molecules whole */
+ mk_mshift(fplog,graph,fr->ePBC,box,x);
+ }
+
+ /* Reset long range forces if necessary */
+ if (fr->bTwinRange)
+ {
+ /* Reset the (long-range) forces if necessary */
+ clear_rvecs(fr->natoms_force_constr,bSepLRF ? fr->f_twin : f);
+ }
+
+ /* Do the actual neighbour searching and if twin range electrostatics
+ * also do the calculation of long range forces and energies.
+ */
+ dvdl = 0;
+ ns(fplog,fr,x,box,
+ groups,&(inputrec->opts),top,mdatoms,
+ cr,nrnb,lambda,&dvdl,&enerd->grpp,bFillGrid,
+ bDoLongRange,bDoForces,bSepLRF ? fr->f_twin : f);
+ if (bSepDVDL)
+ {
+ fprintf(fplog,sepdvdlformat,"LR non-bonded",0.0,dvdl);
+ }
+ enerd->dvdl_lin += dvdl;
+
+ wallcycle_stop(wcycle,ewcNS);
+ }
+
+ if (inputrec->implicit_solvent && bNS)
+ {
+ make_gb_nblist(cr,inputrec->gb_algorithm,inputrec->rlist,
+ x,box,fr,&top->idef,graph,fr->born);
+ }
+
+ if (DOMAINDECOMP(cr))
+ {
+ if (!(cr->duty & DUTY_PME))
+ {
+ wallcycle_start(wcycle,ewcPPDURINGPME);
+ dd_force_flop_start(cr->dd,nrnb);
+ }
+ }
+
+ if (inputrec->bRot)
+ {
+ /* Enforced rotation has its own cycle counter that starts after the collective
+ * coordinates have been communicated. It is added to ddCyclF */
+ do_rotation(cr,inputrec,box,x,t,step,wcycle,bNS);
+ }
+
+ /* Start the force cycle counter.
+ * This counter is stopped in do_forcelow_level.
+ * No parallel communication should occur while this counter is running,
+ * since that will interfere with the dynamic load balancing.
+ */
+ wallcycle_start(wcycle,ewcFORCE);
+ GMX_MPE_LOG(ev_forcecycles_start);
+
+ if (bDoForces)
+ {
+ /* Reset forces for which the virial is calculated separately:
+ * PME/Ewald forces if necessary */
+ if (fr->bF_NoVirSum)
+ {
+ if (flags & GMX_FORCE_VIRIAL)
+ {
+ fr->f_novirsum = fr->f_novirsum_alloc;
+ GMX_BARRIER(cr->mpi_comm_mygroup);
+ if (fr->bDomDec)
+ {
+ clear_rvecs(fr->f_novirsum_n,fr->f_novirsum);
+ }
+ else
+ {
+ clear_rvecs(homenr,fr->f_novirsum+start);
+ }
+ GMX_BARRIER(cr->mpi_comm_mygroup);
+ }
+ else
+ {
+ /* We are not calculating the pressure so we do not need
+ * a separate array for forces that do not contribute
+ * to the pressure.
+ */
+ fr->f_novirsum = f;
+ }
+ }
+
+ if (bSepLRF)
+ {
+ /* Add the long range forces to the short range forces */
+ for(i=0; i<fr->natoms_force_constr; i++)
+ {
+ copy_rvec(fr->f_twin[i],f[i]);
+ }
+ }
+ else if (!(fr->bTwinRange && bNS))
+ {
+ /* Clear the short-range forces */
+ clear_rvecs(fr->natoms_force_constr,f);
+ }
+
+ clear_rvec(fr->vir_diag_posres);
+
+ GMX_BARRIER(cr->mpi_comm_mygroup);
+ }
+ if (inputrec->ePull == epullCONSTRAINT)
+ {
+ clear_pull_forces(inputrec->pull);
+ }
+
+ /* update QMMMrec, if necessary */
+ if(fr->bQMMM)
+ {
+ update_QMMMrec(cr,fr,x,mdatoms,box,top);
+ }
+
+ if ((flags & GMX_FORCE_BONDED) && top->idef.il[F_POSRES].nr > 0)
+ {
+ /* Position restraints always require full pbc */
+ set_pbc(&pbc,inputrec->ePBC,box);
+ v = posres(top->idef.il[F_POSRES].nr,top->idef.il[F_POSRES].iatoms,
+ top->idef.iparams_posres,
+ (const rvec*)x,fr->f_novirsum,fr->vir_diag_posres,
+ inputrec->ePBC==epbcNONE ? NULL : &pbc,lambda,&dvdl,
+ fr->rc_scaling,fr->ePBC,fr->posres_com,fr->posres_comB);
+ if (bSepDVDL)
+ {
+ fprintf(fplog,sepdvdlformat,
+ interaction_function[F_POSRES].longname,v,dvdl);
+ }
+ enerd->term[F_POSRES] += v;
+ /* This linear lambda dependence assumption is only correct
+ * when only k depends on lambda,
+ * not when the reference position depends on lambda.
+ * grompp checks for this.
+ */
+ enerd->dvdl_lin += dvdl;
+ inc_nrnb(nrnb,eNR_POSRES,top->idef.il[F_POSRES].nr/2);
+ }
+
+ /* Compute the bonded and non-bonded energies and optionally forces */
+ do_force_lowlevel(fplog,step,fr,inputrec,&(top->idef),
+ cr,nrnb,wcycle,mdatoms,&(inputrec->opts),
+ x,hist,f,enerd,fcd,mtop,top,fr->born,
+ &(top->atomtypes),bBornRadii,box,
+ lambda,graph,&(top->excls),fr->mu_tot,
+ flags,&cycles_pme);
+
+ cycles_force = wallcycle_stop(wcycle,ewcFORCE);
+ GMX_BARRIER(cr->mpi_comm_mygroup);
+
+ if (ed)
+ {
+ do_flood(fplog,cr,x,f,ed,box,step);
+ }
+
+ if (DOMAINDECOMP(cr))
+ {
+ dd_force_flop_stop(cr->dd,nrnb);
+ if (wcycle)
+ {
+ dd_cycles_add(cr->dd,cycles_force-cycles_pme,ddCyclF);
+ }
+ }
+
+ if (bDoForces)
+ {
+ if (IR_ELEC_FIELD(*inputrec))
+ {
+ /* Compute forces due to electric field */
+ calc_f_el(MASTER(cr) ? field : NULL,
+ start,homenr,mdatoms->chargeA,x,fr->f_novirsum,
+ inputrec->ex,inputrec->et,t);
+ }
+
+ /* Communicate the forces */
+ if (PAR(cr))
+ {
+ wallcycle_start(wcycle,ewcMOVEF);
+ if (DOMAINDECOMP(cr))
+ {
+ dd_move_f(cr->dd,f,fr->fshift);
+ /* Do we need to communicate the separate force array
+ * for terms that do not contribute to the single sum virial?
+ * Position restraints and electric fields do not introduce
+ * inter-cg forces, only full electrostatics methods do.
+ * When we do not calculate the virial, fr->f_novirsum = f,
+ * so we have already communicated these forces.
+ */
+ if (EEL_FULL(fr->eeltype) && cr->dd->n_intercg_excl &&
+ (flags & GMX_FORCE_VIRIAL))
+ {
+ dd_move_f(cr->dd,fr->f_novirsum,NULL);
+ }
+ if (bSepLRF)
+ {
+ /* We should not update the shift forces here,
+ * since f_twin is already included in f.
+ */
+ dd_move_f(cr->dd,fr->f_twin,NULL);
+ }
+ }
+ else
+ {
+ pd_move_f(cr,f,nrnb);
+ if (bSepLRF)
+ {
+ pd_move_f(cr,fr->f_twin,nrnb);
+ }
+ }
+ wallcycle_stop(wcycle,ewcMOVEF);
+ }
+
+ /* If we have NoVirSum forces, but we do not calculate the virial,
+ * we sum fr->f_novirum=f later.
+ */
+ if (vsite && !(fr->bF_NoVirSum && !(flags & GMX_FORCE_VIRIAL)))
+ {
+ wallcycle_start(wcycle,ewcVSITESPREAD);
+ spread_vsite_f(fplog,vsite,x,f,fr->fshift,nrnb,
+ &top->idef,fr->ePBC,fr->bMolPBC,graph,box,cr);
+ wallcycle_stop(wcycle,ewcVSITESPREAD);
+
+ if (bSepLRF)
+ {
+ wallcycle_start(wcycle,ewcVSITESPREAD);
+ spread_vsite_f(fplog,vsite,x,fr->f_twin,NULL,
+ nrnb,
+ &top->idef,fr->ePBC,fr->bMolPBC,graph,box,cr);
+ wallcycle_stop(wcycle,ewcVSITESPREAD);
+ }
+ }
+
+ if (flags & GMX_FORCE_VIRIAL)
+ {
+ /* Calculation of the virial must be done after vsites! */
+ calc_virial(fplog,mdatoms->start,mdatoms->homenr,x,f,
+ vir_force,graph,box,nrnb,fr,inputrec->ePBC);
+ }
+ }
+
+ if (inputrec->ePull == epullUMBRELLA || inputrec->ePull == epullCONST_F)
+ {
+ /* Calculate the center of mass forces, this requires communication,
+ * which is why pull_potential is called close to other communication.
+ * The virial contribution is calculated directly,
+ * which is why we call pull_potential after calc_virial.
+ */
+ set_pbc(&pbc,inputrec->ePBC,box);
+ dvdl = 0;
+ enerd->term[F_COM_PULL] =
+ pull_potential(inputrec->ePull,inputrec->pull,mdatoms,&pbc,
+ cr,t,lambda,x,f,vir_force,&dvdl);
+ if (bSepDVDL)
+ {
+ fprintf(fplog,sepdvdlformat,"Com pull",enerd->term[F_COM_PULL],dvdl);
+ }
+ enerd->dvdl_lin += dvdl;
+ }
+ else
+ enerd->term[F_COM_PULL] = 0.0;
+
+ /* Add the forces from enforced rotation potentials (if any) */
+ if (inputrec->bRot)
+ enerd->term[F_COM_PULL] += add_rot_forces(inputrec->rot, f, cr, step, t);
+
+
+ if (PAR(cr) && !(cr->duty & DUTY_PME))
+ {
+ cycles_ppdpme = wallcycle_stop(wcycle,ewcPPDURINGPME);
+ dd_cycles_add(cr->dd,cycles_ppdpme,ddCyclPPduringPME);
+
+ /* In case of node-splitting, the PP nodes receive the long-range
+ * forces, virial and energy from the PME nodes here.
+ */
+ wallcycle_start(wcycle,ewcPP_PMEWAITRECVF);
+ dvdl = 0;
+ gmx_pme_receive_f(cr,fr->f_novirsum,fr->vir_el_recip,&e,&dvdl,
+ &cycles_seppme);
+ if (bSepDVDL)
+ {
+ fprintf(fplog,sepdvdlformat,"PME mesh",e,dvdl);
+ }
+ enerd->term[F_COUL_RECIP] += e;
+ enerd->dvdl_lin += dvdl;
+ if (wcycle)
+ {
+ dd_cycles_add(cr->dd,cycles_seppme,ddCyclPME);
+ }
+ wallcycle_stop(wcycle,ewcPP_PMEWAITRECVF);
+ }
+
+ if (bDoForces && fr->bF_NoVirSum)
+ {
+ if (vsite)
+ {
+ /* Spread the mesh force on virtual sites to the other particles...
+ * This is parallellized. MPI communication is performed
+ * if the constructing atoms aren't local.
+ */
+ wallcycle_start(wcycle,ewcVSITESPREAD);
+ spread_vsite_f(fplog,vsite,x,fr->f_novirsum,NULL,nrnb,
+ &top->idef,fr->ePBC,fr->bMolPBC,graph,box,cr);
+ wallcycle_stop(wcycle,ewcVSITESPREAD);
+ }
+ if (flags & GMX_FORCE_VIRIAL)
+ {
+ /* Now add the forces, this is local */
+ if (fr->bDomDec)
+ {
+ sum_forces(0,fr->f_novirsum_n,f,fr->f_novirsum);
+ }
+ else
+ {
+ sum_forces(start,start+homenr,f,fr->f_novirsum);
+ }
+ if (EEL_FULL(fr->eeltype))
+ {
+ /* Add the mesh contribution to the virial */
+ m_add(vir_force,fr->vir_el_recip,vir_force);
+ }
+ if (debug)
+ {
+ pr_rvecs(debug,0,"vir_force",vir_force,DIM);
+ }
+ }
+ }
+
+ /* Sum the potential energy terms from group contributions */
+ sum_epot(&(inputrec->opts),enerd);
+
+ if (fr->print_force >= 0 && bDoForces)
+ {
+ print_large_forces(stderr,mdatoms,cr,step,fr->print_force,x,f);
+ }
+}
+
+void do_constrain_first(FILE *fplog,gmx_constr_t constr,
+ t_inputrec *ir,t_mdatoms *md,
+ t_state *state,rvec *f,
+ t_graph *graph,t_commrec *cr,t_nrnb *nrnb,
+ t_forcerec *fr, gmx_localtop_t *top, tensor shake_vir)
+{
+ int i,m,start,end;
+ gmx_large_int_t step;
+ double mass,tmass,vcm[4];
+ real dt=ir->delta_t;
+ real dvdlambda;
+ rvec *savex;
+
+ snew(savex,state->natoms);
+
+ start = md->start;
+ end = md->homenr + start;
+
+ if (debug)
+ fprintf(debug,"vcm: start=%d, homenr=%d, end=%d\n",
+ start,md->homenr,end);
+ /* Do a first constrain to reset particles... */
+ step = ir->init_step;
+ if (fplog)
+ {
+ char buf[STEPSTRSIZE];
+ fprintf(fplog,"\nConstraining the starting coordinates (step %s)\n",
+ gmx_step_str(step,buf));
+ }
+ dvdlambda = 0;
+
+ /* constrain the current position */
+ constrain(NULL,TRUE,FALSE,constr,&(top->idef),
+ ir,NULL,cr,step,0,md,
+ state->x,state->x,NULL,
+ state->box,state->lambda,&dvdlambda,
+ NULL,NULL,nrnb,econqCoord,ir->epc==epcMTTK,state->veta,state->veta);
+ if (EI_VV(ir->eI))
+ {
+ /* constrain the inital velocity, and save it */
+ /* also may be useful if we need the ekin from the halfstep for velocity verlet */
+ /* might not yet treat veta correctly */
+ constrain(NULL,TRUE,FALSE,constr,&(top->idef),
+ ir,NULL,cr,step,0,md,
+ state->x,state->v,state->v,
+ state->box,state->lambda,&dvdlambda,
+ NULL,NULL,nrnb,econqVeloc,ir->epc==epcMTTK,state->veta,state->veta);
+ }
+ /* constrain the inital velocities at t-dt/2 */
+ if (EI_STATE_VELOCITY(ir->eI) && ir->eI!=eiVV)
+ {
+ for(i=start; (i<end); i++)
+ {
+ for(m=0; (m<DIM); m++)
+ {
+ /* Reverse the velocity */
+ state->v[i][m] = -state->v[i][m];
+ /* Store the position at t-dt in buf */
+ savex[i][m] = state->x[i][m] + dt*state->v[i][m];
+ }
+ }
+ /* Shake the positions at t=-dt with the positions at t=0
+ * as reference coordinates.
+ */
+ if (fplog)
+ {
+ char buf[STEPSTRSIZE];
+ fprintf(fplog,"\nConstraining the coordinates at t0-dt (step %s)\n",
+ gmx_step_str(step,buf));
+ }
+ dvdlambda = 0;
+ constrain(NULL,TRUE,FALSE,constr,&(top->idef),
+ ir,NULL,cr,step,-1,md,
+ state->x,savex,NULL,
+ state->box,state->lambda,&dvdlambda,
+ state->v,NULL,nrnb,econqCoord,ir->epc==epcMTTK,state->veta,state->veta);
+
+ for(i=start; i<end; i++) {
+ for(m=0; m<DIM; m++) {
+ /* Re-reverse the velocities */
+ state->v[i][m] = -state->v[i][m];
+ }
+ }
+ }
+
+ for(m=0; (m<4); m++)
+ vcm[m] = 0;
+ for(i=start; i<end; i++) {
+ mass = md->massT[i];
+ for(m=0; m<DIM; m++) {
+ vcm[m] += state->v[i][m]*mass;
+ }
+ vcm[3] += mass;
+ }
+
+ if (ir->nstcomm != 0 || debug) {
+ /* Compute the global sum of vcm */
+ if (debug)
+ fprintf(debug,"vcm: %8.3f %8.3f %8.3f,"
+ " total mass = %12.5e\n",vcm[XX],vcm[YY],vcm[ZZ],vcm[3]);
+ if (PAR(cr))
+ gmx_sumd(4,vcm,cr);
+ tmass = vcm[3];
+ for(m=0; (m<DIM); m++)
+ vcm[m] /= tmass;
+ if (debug)
+ fprintf(debug,"vcm: %8.3f %8.3f %8.3f,"
+ " total mass = %12.5e\n",vcm[XX],vcm[YY],vcm[ZZ],tmass);
+ if (ir->nstcomm != 0) {
+ /* Now we have the velocity of center of mass, let's remove it */
+ for(i=start; (i<end); i++) {
+ for(m=0; (m<DIM); m++)
+ state->v[i][m] -= vcm[m];
+ }
+
+ }
+ }
+ sfree(savex);
+}
+
+void calc_enervirdiff(FILE *fplog,int eDispCorr,t_forcerec *fr)
+{
+ double eners[2],virs[2],enersum,virsum,y0,f,g,h;
+ double r0,r1,r,rc3,rc9,ea,eb,ec,pa,pb,pc,pd;
+ double invscale,invscale2,invscale3;
+ int ri0,ri1,ri,i,offstart,offset;
+ real scale,*vdwtab;
+
+ fr->enershiftsix = 0;
+ fr->enershifttwelve = 0;
+ fr->enerdiffsix = 0;
+ fr->enerdifftwelve = 0;
+ fr->virdiffsix = 0;
+ fr->virdifftwelve = 0;
+
+ if (eDispCorr != edispcNO) {
+ for(i=0; i<2; i++) {
+ eners[i] = 0;
+ virs[i] = 0;
+ }
+ if ((fr->vdwtype == evdwSWITCH) || (fr->vdwtype == evdwSHIFT)) {
+ if (fr->rvdw_switch == 0)
+ gmx_fatal(FARGS,
+ "With dispersion correction rvdw-switch can not be zero "
+ "for vdw-type = %s",evdw_names[fr->vdwtype]);
+
+ scale = fr->nblists[0].tab.scale;
+ vdwtab = fr->nblists[0].vdwtab;
+
+ /* Round the cut-offs to exact table values for precision */
+ ri0 = floor(fr->rvdw_switch*scale);
+ ri1 = ceil(fr->rvdw*scale);
+ r0 = ri0/scale;
+ r1 = ri1/scale;
+ rc3 = r0*r0*r0;
+ rc9 = rc3*rc3*rc3;
+
+ if (fr->vdwtype == evdwSHIFT) {
+ /* Determine the constant energy shift below rvdw_switch */
+ fr->enershiftsix = (real)(-1.0/(rc3*rc3)) - vdwtab[8*ri0];
+ fr->enershifttwelve = (real)( 1.0/(rc9*rc3)) - vdwtab[8*ri0 + 4];
+ }
+ /* Add the constant part from 0 to rvdw_switch.
+ * This integration from 0 to rvdw_switch overcounts the number
+ * of interactions by 1, as it also counts the self interaction.
+ * We will correct for this later.
+ */
+ eners[0] += 4.0*M_PI*fr->enershiftsix*rc3/3.0;
+ eners[1] += 4.0*M_PI*fr->enershifttwelve*rc3/3.0;
+
+ invscale = 1.0/(scale);
+ invscale2 = invscale*invscale;
+ invscale3 = invscale*invscale2;
+
+ /* following summation derived from cubic spline definition,
+ Numerical Recipies in C, second edition, p. 113-116. Exact
+ for the cubic spline. We first calculate the negative of
+ the energy from rvdw to rvdw_switch, assuming that g(r)=1,
+ and then add the more standard, abrupt cutoff correction to
+ that result, yielding the long-range correction for a
+ switched function. We perform both the pressure and energy
+ loops at the same time for simplicity, as the computational
+ cost is low. */
+
+ for (i=0;i<2;i++) {
+ enersum = 0.0; virsum = 0.0;
+ if (i==0)
+ offstart = 0;
+ else
+ offstart = 4;
+ for (ri=ri0; ri<ri1; ri++) {
+ r = ri*invscale;
+ ea = invscale3;
+ eb = 2.0*invscale2*r;
+ ec = invscale*r*r;
+
+ pa = invscale3;
+ pb = 3.0*invscale2*r;
+ pc = 3.0*invscale*r*r;
+ pd = r*r*r;
+
+ /* this "8" is from the packing in the vdwtab array - perhaps
+ should be #define'ed? */
+ offset = 8*ri + offstart;
+ y0 = vdwtab[offset];
+ f = vdwtab[offset+1];
+ g = vdwtab[offset+2];
+ h = vdwtab[offset+3];
+
+ enersum += y0*(ea/3 + eb/2 + ec) + f*(ea/4 + eb/3 + ec/2)+
+ g*(ea/5 + eb/4 + ec/3) + h*(ea/6 + eb/5 + ec/4);
+ virsum += f*(pa/4 + pb/3 + pc/2 + pd) +
+ 2*g*(pa/5 + pb/4 + pc/3 + pd/2) + 3*h*(pa/6 + pb/5 + pc/4 + pd/3);
+
+ }
+ enersum *= 4.0*M_PI;
+ virsum *= 4.0*M_PI;
+ eners[i] -= enersum;
+ virs[i] -= virsum;
+ }
+
+ /* now add the correction for rvdw_switch to infinity */
+ eners[0] += -4.0*M_PI/(3.0*rc3);
+ eners[1] += 4.0*M_PI/(9.0*rc9);
+ virs[0] += 8.0*M_PI/rc3;
+ virs[1] += -16.0*M_PI/(3.0*rc9);
+ }
+ else if ((fr->vdwtype == evdwCUT) || (fr->vdwtype == evdwUSER)) {
+ if (fr->vdwtype == evdwUSER && fplog)
+ fprintf(fplog,
+ "WARNING: using dispersion correction with user tables\n");
+ rc3 = fr->rvdw*fr->rvdw*fr->rvdw;
+ rc9 = rc3*rc3*rc3;
+ eners[0] += -4.0*M_PI/(3.0*rc3);
+ eners[1] += 4.0*M_PI/(9.0*rc9);
+ virs[0] += 8.0*M_PI/rc3;
+ virs[1] += -16.0*M_PI/(3.0*rc9);
+ } else {
+ gmx_fatal(FARGS,
+ "Dispersion correction is not implemented for vdw-type = %s",
+ evdw_names[fr->vdwtype]);
+ }
+ fr->enerdiffsix = eners[0];
+ fr->enerdifftwelve = eners[1];
+ /* The 0.5 is due to the Gromacs definition of the virial */
+ fr->virdiffsix = 0.5*virs[0];
+ fr->virdifftwelve = 0.5*virs[1];
+ }
+}
+
+void calc_dispcorr(FILE *fplog,t_inputrec *ir,t_forcerec *fr,
+ gmx_large_int_t step,int natoms,
+ matrix box,real lambda,tensor pres,tensor virial,
+ real *prescorr, real *enercorr, real *dvdlcorr)
+{
+ gmx_bool bCorrAll,bCorrPres;
+ real dvdlambda,invvol,dens,ninter,avcsix,avctwelve,enerdiff,svir=0,spres=0;
+ int m;
+
+ *prescorr = 0;
+ *enercorr = 0;
+ *dvdlcorr = 0;
+
+ clear_mat(virial);
+ clear_mat(pres);
+
+ if (ir->eDispCorr != edispcNO) {
+ bCorrAll = (ir->eDispCorr == edispcAllEner ||
+ ir->eDispCorr == edispcAllEnerPres);
+ bCorrPres = (ir->eDispCorr == edispcEnerPres ||
+ ir->eDispCorr == edispcAllEnerPres);
+
+ invvol = 1/det(box);
+ if (fr->n_tpi)
+ {
+ /* Only correct for the interactions with the inserted molecule */
+ dens = (natoms - fr->n_tpi)*invvol;
+ ninter = fr->n_tpi;
+ }
+ else
+ {
+ dens = natoms*invvol;
+ ninter = 0.5*natoms;
+ }
+
+ if (ir->efep == efepNO)
+ {
+ avcsix = fr->avcsix[0];
+ avctwelve = fr->avctwelve[0];
+ }
+ else
+ {
+ avcsix = (1 - lambda)*fr->avcsix[0] + lambda*fr->avcsix[1];
+ avctwelve = (1 - lambda)*fr->avctwelve[0] + lambda*fr->avctwelve[1];
+ }
+
+ enerdiff = ninter*(dens*fr->enerdiffsix - fr->enershiftsix);
+ *enercorr += avcsix*enerdiff;
+ dvdlambda = 0.0;
+ if (ir->efep != efepNO)
+ {
+ dvdlambda += (fr->avcsix[1] - fr->avcsix[0])*enerdiff;
+ }
+ if (bCorrAll)
+ {
+ enerdiff = ninter*(dens*fr->enerdifftwelve - fr->enershifttwelve);
+ *enercorr += avctwelve*enerdiff;
+ if (fr->efep != efepNO)
+ {
+ dvdlambda += (fr->avctwelve[1] - fr->avctwelve[0])*enerdiff;
+ }
+ }
+
+ if (bCorrPres)
+ {
+ svir = ninter*dens*avcsix*fr->virdiffsix/3.0;
+ if (ir->eDispCorr == edispcAllEnerPres)
+ {
+ svir += ninter*dens*avctwelve*fr->virdifftwelve/3.0;
+ }
+ /* The factor 2 is because of the Gromacs virial definition */
+ spres = -2.0*invvol*svir*PRESFAC;
+
+ for(m=0; m<DIM; m++) {
+ virial[m][m] += svir;
+ pres[m][m] += spres;
+ }
+ *prescorr += spres;
+ }
+
+ /* Can't currently control when it prints, for now, just print when degugging */
+ if (debug)
+ {
+ if (bCorrAll) {
+ fprintf(debug,"Long Range LJ corr.: <C6> %10.4e, <C12> %10.4e\n",
+ avcsix,avctwelve);
+ }
+ if (bCorrPres)
+ {
+ fprintf(debug,
+ "Long Range LJ corr.: Epot %10g, Pres: %10g, Vir: %10g\n",
+ *enercorr,spres,svir);
+ }
+ else
+ {
+ fprintf(debug,"Long Range LJ corr.: Epot %10g\n",*enercorr);
+ }
+ }
+
+ if (fr->bSepDVDL && do_per_step(step,ir->nstlog))
+ {
+ fprintf(fplog,sepdvdlformat,"Dispersion correction",
+ *enercorr,dvdlambda);
+ }
+ if (fr->efep != efepNO)
+ {
+ *dvdlcorr += dvdlambda;
+ }
+ }
+}
+
+void do_pbc_first(FILE *fplog,matrix box,t_forcerec *fr,
+ t_graph *graph,rvec x[])
+{
+ if (fplog)
+ fprintf(fplog,"Removing pbc first time\n");
+ calc_shifts(box,fr->shift_vec);
+ if (graph) {
+ mk_mshift(fplog,graph,fr->ePBC,box,x);
+ if (gmx_debug_at)
+ p_graph(debug,"do_pbc_first 1",graph);
+ shift_self(graph,box,x);
+ /* By doing an extra mk_mshift the molecules that are broken
+ * because they were e.g. imported from another software
+ * will be made whole again. Such are the healing powers
+ * of GROMACS.
+ */
+ mk_mshift(fplog,graph,fr->ePBC,box,x);
+ if (gmx_debug_at)
+ p_graph(debug,"do_pbc_first 2",graph);
+ }
+ if (fplog)
+ fprintf(fplog,"Done rmpbc\n");
+}
+
+static void low_do_pbc_mtop(FILE *fplog,int ePBC,matrix box,
+ gmx_mtop_t *mtop,rvec x[],
+ gmx_bool bFirst)
+{
+ t_graph *graph;
+ int mb,as,mol;
+ gmx_molblock_t *molb;
+
+ if (bFirst && fplog)
+ fprintf(fplog,"Removing pbc first time\n");
+
+ snew(graph,1);
+ as = 0;
+ for(mb=0; mb<mtop->nmolblock; mb++) {
+ molb = &mtop->molblock[mb];
+ if (molb->natoms_mol == 1 ||
+ (!bFirst && mtop->moltype[molb->type].cgs.nr == 1)) {
+ /* Just one atom or charge group in the molecule, no PBC required */
+ as += molb->nmol*molb->natoms_mol;
+ } else {
+ /* Pass NULL iso fplog to avoid graph prints for each molecule type */
+ mk_graph_ilist(NULL,mtop->moltype[molb->type].ilist,
+ 0,molb->natoms_mol,FALSE,FALSE,graph);
+
+ for(mol=0; mol<molb->nmol; mol++) {
+ mk_mshift(fplog,graph,ePBC,box,x+as);
+
+ shift_self(graph,box,x+as);
+ /* The molecule is whole now.
+ * We don't need the second mk_mshift call as in do_pbc_first,
+ * since we no longer need this graph.
+ */
+
+ as += molb->natoms_mol;
+ }
+ done_graph(graph);
+ }
+ }
+ sfree(graph);
+}
+
+void do_pbc_first_mtop(FILE *fplog,int ePBC,matrix box,
+ gmx_mtop_t *mtop,rvec x[])
+{
+ low_do_pbc_mtop(fplog,ePBC,box,mtop,x,TRUE);
+}
+
+void do_pbc_mtop(FILE *fplog,int ePBC,matrix box,
+ gmx_mtop_t *mtop,rvec x[])
+{
+ low_do_pbc_mtop(fplog,ePBC,box,mtop,x,FALSE);
+}
+
+void finish_run(FILE *fplog,t_commrec *cr,const char *confout,
+ t_inputrec *inputrec,
+ t_nrnb nrnb[],gmx_wallcycle_t wcycle,
+ gmx_runtime_t *runtime,
+ gmx_bool bWriteStat)
+{
+ int i,j;
+ t_nrnb *nrnb_tot=NULL;
+ real delta_t;
+ double nbfs,mflop;
+ double cycles[ewcNR];
+
+ wallcycle_sum(cr,wcycle,cycles);
+
+ if (cr->nnodes > 1) {
+ if (SIMMASTER(cr))
+ snew(nrnb_tot,1);
+#ifdef GMX_MPI
+ MPI_Reduce(nrnb->n,nrnb_tot->n,eNRNB,MPI_DOUBLE,MPI_SUM,
+ MASTERRANK(cr),cr->mpi_comm_mysim);
+#endif
+ } else {
+ nrnb_tot = nrnb;
+ }
+
+ if (SIMMASTER(cr)) {
+ print_flop(fplog,nrnb_tot,&nbfs,&mflop);
+ if (cr->nnodes > 1) {
+ sfree(nrnb_tot);
+ }
+ }
+
+ if ((cr->duty & DUTY_PP) && DOMAINDECOMP(cr)) {
+ print_dd_statistics(cr,inputrec,fplog);
+ }
+
+#ifdef GMX_MPI
+ if (PARTDECOMP(cr))
+ {
+ if (MASTER(cr))
+ {
+ t_nrnb *nrnb_all;
+ int s;
+ MPI_Status stat;
+
+ snew(nrnb_all,cr->nnodes);
+ nrnb_all[0] = *nrnb;
+ for(s=1; s<cr->nnodes; s++)
+ {
+ MPI_Recv(nrnb_all[s].n,eNRNB,MPI_DOUBLE,s,0,
+ cr->mpi_comm_mysim,&stat);
+ }
+ pr_load(fplog,cr,nrnb_all);
+ sfree(nrnb_all);
+ }
+ else
+ {
+ MPI_Send(nrnb->n,eNRNB,MPI_DOUBLE,MASTERRANK(cr),0,
+ cr->mpi_comm_mysim);
+ }
+ }
+#endif
+
+ if (SIMMASTER(cr)) {
+ wallcycle_print(fplog,cr->nnodes,cr->npmenodes,runtime->realtime,
+ wcycle,cycles);
+
+ if (EI_DYNAMICS(inputrec->eI)) {
+ delta_t = inputrec->delta_t;
+ } else {
+ delta_t = 0;
+ }
+
+ if (fplog) {
+ print_perf(fplog,runtime->proctime,runtime->realtime,
+ cr->nnodes-cr->npmenodes,
+ runtime->nsteps_done,delta_t,nbfs,mflop);
+ }
+ if (bWriteStat) {
+ print_perf(stderr,runtime->proctime,runtime->realtime,
+ cr->nnodes-cr->npmenodes,
+ runtime->nsteps_done,delta_t,nbfs,mflop);
+ }
+
+ /*
+ runtime=inputrec->nsteps*inputrec->delta_t;
+ if (bWriteStat) {
+ if (cr->nnodes == 1)
+ fprintf(stderr,"\n\n");
+ print_perf(stderr,nodetime,realtime,runtime,&ntot,
+ cr->nnodes-cr->npmenodes,FALSE);
+ }
+ wallcycle_print(fplog,cr->nnodes,cr->npmenodes,realtime,wcycle,cycles);
+ print_perf(fplog,nodetime,realtime,runtime,&ntot,cr->nnodes-cr->npmenodes,
+ TRUE);
+ if (PARTDECOMP(cr))
+ pr_load(fplog,cr,nrnb_all);
+ if (cr->nnodes > 1)
+ sfree(nrnb_all);
+ */
+ }
+}
+
+void init_md(FILE *fplog,
+ t_commrec *cr,t_inputrec *ir,const output_env_t oenv,
+ double *t,double *t0,
+ real *lambda,double *lam0,
+ t_nrnb *nrnb,gmx_mtop_t *mtop,
+ gmx_update_t *upd,
+ int nfile,const t_filenm fnm[],
+ gmx_mdoutf_t **outf,t_mdebin **mdebin,
+ tensor force_vir,tensor shake_vir,rvec mu_tot,
+ gmx_bool *bSimAnn,t_vcm **vcm, t_state *state, unsigned long Flags)
+{
+ int i,j,n;
+ real tmpt,mod;
+
+ /* Initial values */
+ *t = *t0 = ir->init_t;
+ if (ir->efep != efepNO)
+ {
+ *lam0 = ir->init_lambda;
+ *lambda = *lam0 + ir->init_step*ir->delta_lambda;
+ }
+ else
+ {
+ *lambda = *lam0 = 0.0;
+ }
+
+ *bSimAnn=FALSE;
+ for(i=0;i<ir->opts.ngtc;i++)
+ {
+ /* set bSimAnn if any group is being annealed */
+ if(ir->opts.annealing[i]!=eannNO)
+ {
+ *bSimAnn = TRUE;
+ }
+ }
+ if (*bSimAnn)
+ {
+ update_annealing_target_temp(&(ir->opts),ir->init_t);
+ }
+
+ if (upd)
+ {
+ *upd = init_update(fplog,ir);
+ }
+
+ if (vcm != NULL)
+ {
+ *vcm = init_vcm(fplog,&mtop->groups,ir);
+ }
+
+ if (EI_DYNAMICS(ir->eI) && !(Flags & MD_APPENDFILES))
+ {
+ if (ir->etc == etcBERENDSEN)
+ {
+ please_cite(fplog,"Berendsen84a");
+ }
+ if (ir->etc == etcVRESCALE)
+ {
+ please_cite(fplog,"Bussi2007a");
+ }
+ }
+
+ init_nrnb(nrnb);
+
+ if (nfile != -1)
+ {
+ *outf = init_mdoutf(nfile,fnm,Flags,cr,ir,oenv);
+
+ *mdebin = init_mdebin((Flags & MD_APPENDFILES) ? NULL : (*outf)->fp_ene,
+ mtop,ir, (*outf)->fp_dhdl);
+ }
+
+ /* Initiate variables */
+ clear_mat(force_vir);
+ clear_mat(shake_vir);
+ clear_rvec(mu_tot);
+
+ debug_gmx();
+}
+
--- /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.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 [TT].tpx[tt] 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.",
+ "[TT]grompp[tt] also reads parameters for the [TT]mdrun[tt] ",
+ "(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]",
+
+ "[TT]grompp[tt] 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]",
+
+ "[TT]grompp[tt] 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 [TT]grompp[tt], 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. Only if this information",
+ "is absent will the coordinates in the [TT]-c[tt] file be 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.[PAR]",
+
+ "[TT]grompp[tt] can be used to restart simulations preserving",
+ "continuity by supplying just a checkpoint file with [TT]-t[tt].",
+ "However, for simply changing the number of run steps to extend",
+ "a run, using [TT]tpbconv[tt] is more convenient than [TT]grompp[tt].",
+ "You then supply the old checkpoint file directly to [TT]mdrun[tt]",
+ "with [TT]-cpi[tt]. If you wish to change the ensemble or things",
+ "like output frequency, then supplying the checkpoint file to",
+ "[TT]grompp[tt] with [TT]-t[tt] along with a new [TT].mdp[tt] file",
+ "with [TT]-f[tt] is the recommended procedure.[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 [TT]grompp.log[tt] (along with real debug info). You",
+ "can see the contents of the run input file with the [TT]gmxdump[tt]",
+ "program. [TT]gmxcheck[tt] can be used to compare the contents of two",
+ "run input files.[PAR]"
+
+ "The [TT]-maxwarn[tt] option can be used to override warnings printed",
+ "by [TT]grompp[tt] 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);
+ }
+
++ if (EI_ENERGY_MINIMIZATION(ir->eI) && 0 == ir->nsteps)
++ {
++ warning_note(wi,"Zero-step energy minimization will alter the coordinates before calculating the energy. If you just want the energy of a single point, try zero-step MD (with unconstrained_start = yes). To do multiple single-point energy evaluations of different configurations of the same topology, use mdrun -rerun.");
++ }
++
+ 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
- /* 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;
- }
-
+/* -*- 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)
+ {
+ /* We need to be sure replica exchange can only occur
+ * when the energies are current */
+ check_nst_param(fplog,cr,"nstcalcenergy",ir->nstcalcenergy,
+ "repl_ex_nst",&repl_ex_nst);
+ /* This check needs to happen before inter-simulation
+ * signals are initialized, too */
+ }
+ 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 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