--- /dev/null
+/*
+ *
+ * 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:
+ * GROningen Mixture of Alchemy and Childrens' Stories
+ */
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <ctype.h>
+#include <string.h>
+#include <assert.h>
+#include "sysstuff.h"
+#include "strdb.h"
+#include "futil.h"
+#include "macros.h"
+#include "names.h"
+#include "string2.h"
+#include "statutil.h"
+#include "confio.h"
+#include "main.h"
+#include "copyrite.h"
+#include "typedefs.h"
+#include "smalloc.h"
+#include "invblock.h"
+#include "macros.h"
+#include "index.h"
+#include "txtdump.h"
+#include "gmxfio.h"
+
+
+
+const char gmx_residuetype_undefined[]="Other";
+
+struct gmx_residuetype
+{
+ int n;
+ char ** resname;
+ char ** restype;
+
+};
+
+
+static gmx_bool gmx_ask_yesno(gmx_bool bASK)
+{
+ char c;
+
+ if (bASK) {
+ do {
+ c=toupper(fgetc(stdin));
+ } while ((c != 'Y') && (c != 'N'));
+
+ return (c == 'Y');
+ }
+ else
+ return FALSE;
+}
+
+t_blocka *new_blocka(void)
+{
+ t_blocka *block;
+
+ snew(block,1);
+ snew(block->index,1);
+
+ return block;
+}
+
+void write_index(const char *outf, t_blocka *b,char **gnames)
+{
+ FILE *out;
+ int i,j,k;
+
+ out=gmx_fio_fopen(outf,"w");
+ /* fprintf(out,"%5d %5d\n",b->nr,b->nra); */
+ for(i=0; (i<b->nr); i++) {
+ fprintf(out,"[ %s ]\n",gnames[i]);
+ for(k=0,j=b->index[i]; j<b->index[i+1]; j++,k++) {
+ fprintf(out,"%4d ",b->a[j]+1);
+ if ((k % 15) == 14)
+ fprintf(out,"\n");
+ }
+ fprintf(out,"\n");
+ }
+ gmx_fio_fclose(out);
+}
+
+void add_grp(t_blocka *b,char ***gnames,int nra,atom_id a[],const char *name)
+{
+ int i;
+
+ srenew(b->index,b->nr+2);
+ srenew(*gnames,b->nr+1);
+ (*gnames)[b->nr]=strdup(name);
+
+ srenew(b->a,b->nra+nra);
+ for(i=0; (i<nra); i++)
+ b->a[b->nra++]=a[i];
+ b->nr++;
+ b->index[b->nr]=b->nra;
+}
+
+/* compare index in `a' with group in `b' at `index',
+ when `index'<0 it is relative to end of `b' */
+static gmx_bool grp_cmp(t_blocka *b, int nra, atom_id a[], int index)
+{
+ int i;
+
+ if (index < 0)
+ index = b->nr-1+index;
+ if (index >= b->nr)
+ gmx_fatal(FARGS,"no such index group %d in t_blocka (nr=%d)",index,b->nr);
+ /* compare sizes */
+ if ( nra != b->index[index+1] - b->index[index] )
+ return FALSE;
+ for(i=0; i<nra; i++)
+ if ( a[i] != b->a[b->index[index]+i] )
+ return FALSE;
+ return TRUE;
+}
+
+static void
+p_status(const char **restype, int nres, const char **typenames, int ntypes)
+{
+ int i,j;
+ int found;
+
+ int * counter;
+
+ snew(counter,ntypes);
+ for(i=0;i<ntypes;i++)
+ {
+ counter[i]=0;
+ }
+ for(i=0;i<nres;i++)
+ {
+ found=0;
+ for(j=0;j<ntypes;j++)
+ {
+ if(!gmx_strcasecmp(restype[i],typenames[j]))
+ {
+ counter[j]++;
+ }
+ }
+ }
+
+ for(i=0; (i<ntypes); i++)
+ {
+ if (counter[i] > 0)
+ {
+ printf("There are: %5d %10s residues\n",counter[i],typenames[i]);
+ }
+ }
+
+ sfree(counter);
+}
+
+
+atom_id *
+mk_aid(t_atoms *atoms,const char ** restype,const char * typestring,int *nra,gmx_bool bMatch)
+/* Make an array of atom_ids for all atoms with residuetypes matching typestring, or the opposite if bMatch is false */
+{
+ atom_id *a;
+ int i;
+ int res;
+
+ snew(a,atoms->nr);
+ *nra=0;
+ for(i=0; (i<atoms->nr); i++)
+ {
+ res=!gmx_strcasecmp(restype[atoms->atom[i].resind],typestring);
+ if(bMatch==FALSE)
+ {
+ res=!res;
+ }
+ if(res)
+ {
+ a[(*nra)++]=i;
+ }
+ }
+
+ return a;
+}
+
+typedef struct {
+ char *rname;
+ gmx_bool bNeg;
+ char *gname;
+} restp_t;
+
+static void analyse_other(const char ** restype,t_atoms *atoms,
+ t_blocka *gb,char ***gn,gmx_bool bASK,gmx_bool bVerb)
+{
+ restp_t *restp=NULL;
+ char **attp=NULL;
+ char *rname,*aname;
+ atom_id *other_ndx,*aid,*aaid;
+ int i,j,k,l,resind,naid,naaid,natp,nrestp=0;
+
+ for(i=0; (i<atoms->nres); i++)
+ {
+ if (gmx_strcasecmp(restype[i],"Protein") && gmx_strcasecmp(restype[i],"DNA") && gmx_strcasecmp(restype[i],"RNA") && gmx_strcasecmp(restype[i],"Water"))
+ {
+ break;
+ }
+ }
+ if (i < atoms->nres) {
+ /* we have others */
+ if (bVerb)
+ printf("Analysing residues not classified as Protein/DNA/RNA/Water and splitting into groups...\n");
+ snew(other_ndx,atoms->nr);
+ for(k=0; (k<atoms->nr); k++) {
+ resind = atoms->atom[k].resind;
+ rname = *atoms->resinfo[resind].name;
+ if (gmx_strcasecmp(restype[resind],"Protein") && gmx_strcasecmp(restype[resind],"DNA") &&
+ gmx_strcasecmp(restype[resind],"RNA") && gmx_strcasecmp(restype[resind],"Water"))
+ {
+
+ for(l=0; (l<nrestp); l++)
+ if (strcmp(restp[l].rname,rname) == 0)
+ break;
+ if (l==nrestp) {
+ srenew(restp,nrestp+1);
+ restp[nrestp].rname = strdup(rname);
+ restp[nrestp].bNeg = FALSE;
+ restp[nrestp].gname = strdup(rname);
+ nrestp++;
+
+ }
+ }
+ }
+ for(i=0; (i<nrestp); i++) {
+ snew(aid,atoms->nr);
+ naid=0;
+ for(j=0; (j<atoms->nr); j++) {
+ rname = *atoms->resinfo[atoms->atom[j].resind].name;
+ if ((strcmp(restp[i].rname,rname) == 0 && !restp[i].bNeg) ||
+ (strcmp(restp[i].rname,rname) != 0 && restp[i].bNeg)) {
+ aid[naid++] = j;
+ }
+ }
+ add_grp(gb,gn,naid,aid,restp[i].gname);
+ if (bASK) {
+ printf("split %s into atoms (y/n) ? ",restp[i].gname);
+ fflush(stdout);
+ if (gmx_ask_yesno(bASK)) {
+ natp=0;
+ for(k=0; (k<naid); k++) {
+ aname=*atoms->atomname[aid[k]];
+ for(l=0; (l<natp); l++)
+ if (strcmp(aname,attp[l]) == 0)
+ break;
+ if (l == natp) {
+ srenew(attp,++natp);
+ attp[natp-1]=aname;
+ }
+ }
+ if (natp > 1) {
+ for(l=0; (l<natp); l++) {
+ snew(aaid,naid);
+ naaid=0;
+ for(k=0; (k<naid); k++) {
+ aname=*atoms->atomname[aid[k]];
+ if (strcmp(aname,attp[l])==0)
+ aaid[naaid++]=aid[k];
+ }
+ add_grp(gb,gn,naaid,aaid,attp[l]);
+ sfree(aaid);
+ }
+ }
+ sfree(attp);
+ attp=NULL;
+ }
+ sfree(aid);
+ }
+ }
+ sfree(other_ndx);
+ }
+}
+
+/*! /brief Instances of this struct contain the data necessary to
+ * construct a single (protein) index group in
+ * analyse_prot(). */
+typedef struct gmx_help_make_index_group
+{
+ /** The set of atom names that will be used to form this index group */
+ const char **defining_atomnames;
+ /** Size of the defining_atomnames array */
+ const int num_defining_atomnames;
+ /** Name of this index group */
+ const char *group_name;
+ /** Whether the above atom names name the atoms in the group, or
+ those not in the group */
+ gmx_bool bTakeComplement;
+ /** The index in wholename gives the first item in the arrays of
+ atomnames that should be tested with 'gmx_strncasecmp' in stead of
+ gmx_strcasecmp, or -1 if all items should be tested with strcasecmp
+ This is comparable to using a '*' wildcard at the end of specific
+ atom names, but that is more involved to implement...
+ */
+ int wholename;
+ /** Only create this index group if it differs from the one specified in compareto,
+ where -1 means to always create this group. */
+ int compareto;
+} t_gmx_help_make_index_group;
+
+static void analyse_prot(const char ** restype,t_atoms *atoms,
+ t_blocka *gb,char ***gn,gmx_bool bASK,gmx_bool bVerb)
+{
+ /* lists of atomnames to be used in constructing index groups: */
+ static const char *pnoh[] = { "H", "HN" };
+ static const char *pnodum[] = { "MN1", "MN2", "MCB1", "MCB2", "MCG1", "MCG2",
+ "MCD1", "MCD2", "MCE1", "MCE2", "MNZ1", "MNZ2" };
+ static const char *calpha[] = { "CA" };
+ static const char *bb[] = { "N","CA","C" };
+ static const char *mc[] = { "N","CA","C","O","O1","O2","OC1","OC2","OT","OXT" };
+ static const char *mcb[] = { "N","CA","CB","C","O","O1","O2","OC1","OC2","OT","OXT" };
+ static const char *mch[] = { "N","CA","C","O","O1","O2","OC1","OC2","OT","OXT",
+ "H1","H2","H3","H","HN" };
+
+ static const t_gmx_help_make_index_group constructing_data[] =
+ {{ NULL, 0, "Protein", TRUE, -1, -1},
+ { pnoh, asize(pnoh), "Protein-H", TRUE, 0, -1},
+ { calpha, asize(calpha), "C-alpha", FALSE, -1, -1},
+ { bb, asize(bb), "Backbone", FALSE, -1, -1},
+ { mc, asize(mc), "MainChain", FALSE, -1, -1},
+ { mcb, asize(mcb), "MainChain+Cb", FALSE, -1, -1},
+ { mch, asize(mch), "MainChain+H", FALSE, -1, -1},
+ { mch, asize(mch), "SideChain", TRUE, -1, -1},
+ { mch, asize(mch), "SideChain-H", TRUE, 11, -1},
+ { pnodum, asize(pnodum), "Prot-Masses", TRUE, -1, 0},
+ };
+ const int num_index_groups = asize(constructing_data);
+
+ int n,j;
+ atom_id *aid;
+ int nra,nnpres,npres;
+ gmx_bool match;
+ char ndx_name[STRLEN],*atnm;
+ int i;
+
+ if (bVerb)
+ {
+ printf("Analysing Protein...\n");
+ }
+ snew(aid,atoms->nr);
+
+ /* calculate the number of protein residues */
+ npres=0;
+ for(i=0; (i<atoms->nres); i++) {
+ if (0 == gmx_strcasecmp(restype[i],"Protein")) {
+ npres++;
+ }
+ }
+ /* find matching or complement atoms */
+ for(i=0; (i<(int)num_index_groups); i++) {
+ nra=0;
+ for(n=0; (n<atoms->nr); n++) {
+ if (0 == gmx_strcasecmp(restype[atoms->atom[n].resind],"Protein")) {
+ match=FALSE;
+ for(j=0; (j<constructing_data[i].num_defining_atomnames); j++) {
+ /* skip digits at beginning of atomname, e.g. 1H */
+ atnm=*atoms->atomname[n];
+ while (isdigit(atnm[0])) {
+ atnm++;
+ }
+ if ( (constructing_data[i].wholename==-1) || (j<constructing_data[i].wholename) ) {
+ if (0 == gmx_strcasecmp(constructing_data[i].defining_atomnames[j],atnm)) {
+ match=TRUE;
+ }
+ } else {
+ if (0 == gmx_strncasecmp(constructing_data[i].defining_atomnames[j],atnm,strlen(constructing_data[i].defining_atomnames[j]))) {
+ match=TRUE;
+ }
+ }
+ }
+ if (constructing_data[i].bTakeComplement != match) {
+ aid[nra++]=n;
+ }
+ }
+ }
+ /* if we want to add this group always or it differs from previous
+ group, add it: */
+ if ( -1 == constructing_data[i].compareto || !grp_cmp(gb,nra,aid,constructing_data[i].compareto-i) ) {
+ add_grp(gb,gn,nra,aid,constructing_data[i].group_name);
+ }
+ }
+
+ if (bASK) {
+ for(i=0; (i<(int)num_index_groups); i++) {
+ printf("Split %12s into %5d residues (y/n) ? ",constructing_data[i].group_name,npres);
+ if (gmx_ask_yesno(bASK)) {
+ int resind;
+ nra = 0;
+ for(n=0;((atoms->atom[n].resind < npres) && (n<atoms->nr));) {
+ resind = atoms->atom[n].resind;
+ for(;((atoms->atom[n].resind==resind) && (n<atoms->nr));n++) {
+ match=FALSE;
+ for(j=0;(j<constructing_data[i].num_defining_atomnames); j++) {
+ if (0 == gmx_strcasecmp(constructing_data[i].defining_atomnames[j],*atoms->atomname[n])) {
+ match=TRUE;
+ }
+ }
+ if (constructing_data[i].bTakeComplement != match) {
+ aid[nra++]=n;
+ }
+ }
+ /* copy the residuename to the tail of the groupname */
+ if (nra > 0) {
+ t_resinfo *ri;
+ ri = &atoms->resinfo[resind];
+ sprintf(ndx_name,"%s_%s%d%c",
+ constructing_data[i].group_name,*ri->name,ri->nr,ri->ic==' ' ? '\0' : ri->ic);
+ add_grp(gb,gn,nra,aid,ndx_name);
+ nra = 0;
+ }
+ }
+ }
+ }
+ printf("Make group with sidechain and C=O swapped (y/n) ? ");
+ if (gmx_ask_yesno(bASK)) {
+ /* Make swap sidechain C=O index */
+ int resind,hold;
+ nra = 0;
+ for(n=0;((atoms->atom[n].resind < npres) && (n<atoms->nr));) {
+ resind = atoms->atom[n].resind;
+ hold = -1;
+ for(;((atoms->atom[n].resind==resind) && (n<atoms->nr));n++)
+ if (strcmp("CA",*atoms->atomname[n]) == 0) {
+ aid[nra++]=n;
+ hold=nra;
+ nra+=2;
+ } else if (strcmp("C",*atoms->atomname[n]) == 0) {
+ if (hold == -1) {
+ gmx_incons("Atom naming problem");
+ }
+ aid[hold]=n;
+ } else if (strcmp("O",*atoms->atomname[n]) == 0) {
+ if (hold == -1) {
+ gmx_incons("Atom naming problem");
+ }
+ aid[hold+1]=n;
+ } else if (strcmp("O1",*atoms->atomname[n]) == 0) {
+ if (hold == -1) {
+ gmx_incons("Atom naming problem");
+ }
+ aid[hold+1]=n;
+ } else
+ aid[nra++]=n;
+ }
+ /* copy the residuename to the tail of the groupname */
+ if (nra > 0) {
+ add_grp(gb,gn,nra,aid,"SwapSC-CO");
+ nra = 0;
+ }
+ }
+ }
+ sfree(aid);
+}
+
+
+
+
+/* Return 0 if the name was found, otherwise -1.
+ * p_restype is set to a pointer to the type name, or 'Other' if we did not find it.
+ */
+int
+gmx_residuetype_get_type(gmx_residuetype_t rt,const char * resname, const char ** p_restype)
+{
+ int i,rc;
+
+ rc=-1;
+ for(i=0;i<rt->n && rc;i++)
+ {
+ rc=gmx_strcasecmp(rt->resname[i],resname);
+ }
+
+ *p_restype = (rc==0) ? rt->restype[i-1] : gmx_residuetype_undefined;
+
+ return rc;
+}
+
+int
+gmx_residuetype_add(gmx_residuetype_t rt,const char *newresname, const char *newrestype)
+{
+ int i;
+ int found;
+ const char * p_oldtype;
+
+ found = !gmx_residuetype_get_type(rt,newresname,&p_oldtype);
+
+ if(found && gmx_strcasecmp(p_oldtype,newrestype))
+ {
+ fprintf(stderr,"Warning: Residue '%s' already present with type '%s' in database, ignoring new type '%s'.",
+ newresname,p_oldtype,newrestype);
+ }
+
+ if(found==0)
+ {
+ srenew(rt->resname,rt->n+1);
+ srenew(rt->restype,rt->n+1);
+ rt->resname[rt->n]=strdup(newresname);
+ rt->restype[rt->n]=strdup(newrestype);
+ rt->n++;
+ }
+
+ return 0;
+}
+
+
+int
+gmx_residuetype_init(gmx_residuetype_t *prt)
+{
+ FILE * db;
+ char line[STRLEN];
+ char resname[STRLEN],restype[STRLEN],dum[STRLEN];
+ char * p;
+ int i;
+ struct gmx_residuetype *rt;
+
+ snew(rt,1);
+ *prt=rt;
+
+ rt->n = 0;
+ rt->resname = NULL;
+ rt->restype = NULL;
+
+ db=libopen("residuetypes.dat");
+
+ while(get_a_line(db,line,STRLEN))
+ {
+ strip_comment(line);
+ trim(line);
+ if(line[0]!='\0')
+ {
+ if(sscanf(line,"%s %s %s",resname,restype,dum)!=2)
+ {
+ gmx_fatal(FARGS,"Incorrect number of columns (2 expected) for line in residuetypes.dat");
+ }
+ gmx_residuetype_add(rt,resname,restype);
+ }
+ }
+
+ fclose(db);
+
+ return 0;
+}
+
+
+
+int
+gmx_residuetype_destroy(gmx_residuetype_t rt)
+{
+ int i;
+
+ for(i=0;i<rt->n;i++)
+ {
+ free(rt->resname[i]);
+ free(rt->restype[i]);
+ }
+ free(rt);
+
+ return 0;
+}
+
+int
+gmx_residuetype_get_alltypes(gmx_residuetype_t rt,
+ const char *** p_typenames,
+ int * ntypes)
+{
+ int i,j,n;
+ int found;
+ const char ** my_typename;
+ char * p;
+
+ n=0;
+
+ my_typename=NULL;
+ for(i=0;i<rt->n;i++)
+ {
+ p=rt->restype[i];
+ found=0;
+ for(j=0;j<n && !found;j++)
+ {
+ found=!gmx_strcasecmp(p,my_typename[j]);
+ }
+
+ if(!found)
+ {
+ srenew(my_typename,n+1);
+ my_typename[n]=p;
+ n++;
+ }
+ }
+ *ntypes=n;
+ *p_typenames=my_typename;
+
+ return 0;
+}
+
+
+
+gmx_bool
+gmx_residuetype_is_protein(gmx_residuetype_t rt, const char *resnm)
+{
+ gmx_bool rc;
+ const char *p_type;
+
+ if(gmx_residuetype_get_type(rt,resnm,&p_type)==0 &&
+ gmx_strcasecmp(p_type,"Protein")==0)
+ {
+ rc=TRUE;
+ }
+ else
+ {
+ rc=FALSE;
+ }
+ return rc;
+}
+
+gmx_bool
+gmx_residuetype_is_dna(gmx_residuetype_t rt, const char *resnm)
+{
+ gmx_bool rc;
+ const char *p_type;
+
+ if(gmx_residuetype_get_type(rt,resnm,&p_type)==0 &&
+ gmx_strcasecmp(p_type,"DNA")==0)
+ {
+ rc=TRUE;
+ }
+ else
+ {
+ rc=FALSE;
+ }
+ return rc;
+}
+
+gmx_bool
+gmx_residuetype_is_rna(gmx_residuetype_t rt, const char *resnm)
+{
+ gmx_bool rc;
+ const char *p_type;
+
+ if(gmx_residuetype_get_type(rt,resnm,&p_type)==0 &&
+ gmx_strcasecmp(p_type,"RNA")==0)
+ {
+ rc=TRUE;
+ }
+ else
+ {
+ rc=FALSE;
+ }
+ return rc;
+}
+
+/* Return the size of the arrays */
+int
+gmx_residuetype_get_size(gmx_residuetype_t rt)
+{
+ return rt->n;
+}
+
+/* Search for a residuetype with name resnm within the
+ * gmx_residuetype database. Return the index if found,
+ * otherwise -1.
+ */
+int
+gmx_residuetype_get_index(gmx_residuetype_t rt, const char *resnm)
+{
+ int i,rc;
+
+ rc=-1;
+ for(i=0;i<rt->n && rc;i++)
+ {
+ rc=gmx_strcasecmp(rt->resname[i],resnm);
+ }
+
+ return (0 == rc) ? i-1 : -1;
+}
+
+/* Return the name of the residuetype with the given index, or
+ * NULL if not found. */
+const char *
+gmx_residuetype_get_name(gmx_residuetype_t rt, int index)
+{
+ if(index >= 0 && index < rt->n) {
+ return rt->resname[index];
+ } else {
+ return NULL;
+ }
+}
+
+
+
+void analyse(t_atoms *atoms,t_blocka *gb,char ***gn,gmx_bool bASK,gmx_bool bVerb)
+{
+ gmx_residuetype_t rt=NULL;
+ char *resnm;
+ atom_id *aid;
+ const char ** restype;
+ int nra;
+ int i,k;
+ size_t j;
+ int ntypes;
+ char * p;
+ const char ** p_typename;
+ int iwater,iion;
+ int nwater,nion;
+ int found;
+
+ if (bVerb)
+ {
+ printf("Analysing residue names:\n");
+ }
+ /* Create system group, every single atom */
+ snew(aid,atoms->nr);
+ for(i=0;i<atoms->nr;i++)
+ {
+ aid[i]=i;
+ }
+ add_grp(gb,gn,atoms->nr,aid,"System");
+ sfree(aid);
+
+ /* For every residue, get a pointer to the residue type name */
+ gmx_residuetype_init(&rt);
+ assert(rt);
+
+ snew(restype,atoms->nres);
+ ntypes = 0;
+ p_typename = NULL;
+ for(i=0;i<atoms->nres;i++)
+ {
+ resnm = *atoms->resinfo[i].name;
+ gmx_residuetype_get_type(rt,resnm,&(restype[i]));
+
+ /* Note that this does not lead to a N*N loop, but N*K, where
+ * K is the number of residue _types_, which is small and independent of N.
+ */
+ found = 0;
+ for(k=0;k<ntypes && !found;k++)
+ {
+ found = !strcmp(restype[i],p_typename[k]);
+ }
+ if(!found)
+ {
+ srenew(p_typename,ntypes+1);
+ p_typename[ntypes++] = strdup(restype[i]);
+ }
+ }
+
+ if (bVerb)
+ {
+ p_status(restype,atoms->nres,p_typename,ntypes);
+ }
+
+ for(k=0;k<ntypes;k++)
+ {
+ aid=mk_aid(atoms,restype,p_typename[k],&nra,TRUE);
+
+ /* Check for special types to do fancy stuff with */
+
+ if(!gmx_strcasecmp(p_typename[k],"Protein") && nra>0)
+ {
+ sfree(aid);
+ /* PROTEIN */
+ analyse_prot(restype,atoms,gb,gn,bASK,bVerb);
+
+ /* Create a Non-Protein group */
+ aid=mk_aid(atoms,restype,"Protein",&nra,FALSE);
+ if ((nra > 0) && (nra < atoms->nr))
+ {
+ add_grp(gb,gn,nra,aid,"non-Protein");
+ }
+ sfree(aid);
+ }
+ else if(!gmx_strcasecmp(p_typename[k],"Water") && nra>0)
+ {
+ add_grp(gb,gn,nra,aid,p_typename[k]);
+ /* Add this group as 'SOL' too, for backward compatibility with older gromacs versions */
+ add_grp(gb,gn,nra,aid,"SOL");
+
+ sfree(aid);
+
+ /* Solvent, create a negated group too */
+ aid=mk_aid(atoms,restype,"Water",&nra,FALSE);
+ if ((nra > 0) && (nra < atoms->nr))
+ {
+ add_grp(gb,gn,nra,aid,"non-Water");
+ }
+ sfree(aid);
+ }
+ else if(nra>0)
+ {
+ /* Other groups */
+ add_grp(gb,gn,nra,aid,p_typename[k]);
+ sfree(aid);
+ analyse_other(restype,atoms,gb,gn,bASK,bVerb);
+ }
+ }
+
+ sfree(p_typename);
+ sfree(restype);
+ gmx_residuetype_destroy(rt);
+
+ /* Create a merged water_and_ions group */
+ iwater = -1;
+ iion = -1;
+ nwater = 0;
+ nion = 0;
+
+ for(i=0;i<gb->nr;i++)
+ {
+ if(!gmx_strcasecmp((*gn)[i],"Water"))
+ {
+ iwater = i;
+ nwater = gb->index[i+1]-gb->index[i];
+ }
+ else if(!gmx_strcasecmp((*gn)[i],"Ion"))
+ {
+ iion = i;
+ nion = gb->index[i+1]-gb->index[i];
+ }
+ }
+
+ if(nwater>0 && nion>0)
+ {
+ srenew(gb->index,gb->nr+2);
+ srenew(*gn,gb->nr+1);
+ (*gn)[gb->nr] = strdup("Water_and_ions");
+ srenew(gb->a,gb->nra+nwater+nion);
+ if(nwater>0)
+ {
+ for(i=gb->index[iwater];i<gb->index[iwater+1];i++)
+ {
+ gb->a[gb->nra++] = gb->a[i];
+ }
+ }
+ if(nion>0)
+ {
+ for(i=gb->index[iion];i<gb->index[iion+1];i++)
+ {
+ gb->a[gb->nra++] = gb->a[i];
+ }
+ }
+ gb->nr++;
+ gb->index[gb->nr]=gb->nra;
+ }
+}
+
+
+void check_index(char *gname,int n,atom_id index[],char *traj,int natoms)
+{
+ int i;
+
+ for(i=0; i<n; i++)
+ if (index[i] >= natoms)
+ gmx_fatal(FARGS,"%s atom number (index[%d]=%d) is larger than the number of atoms in %s (%d)",
+ gname ? gname : "Index",i+1, index[i]+1,
+ traj ? traj : "the trajectory",natoms);
+ else if (index[i] < 0)
+ gmx_fatal(FARGS,"%s atom number (index[%d]=%d) is less than zero",
+ gname ? gname : "Index",i+1, index[i]+1);
+}
+
+t_blocka *init_index(const char *gfile, char ***grpname)
+{
+ FILE *in;
+ t_blocka *b;
+ int a,maxentries;
+ int i,j,ng,nread;
+ char line[STRLEN],*pt,str[STRLEN];
+
+ in=gmx_fio_fopen(gfile,"r");
+ snew(b,1);
+ get_a_line(in,line,STRLEN);
+ if ( line[0]=='[' ) {
+ /* new format */
+ b->nr=0;
+ b->index=NULL;
+ b->nra=0;
+ b->a=NULL;
+ *grpname=NULL;
+ maxentries=0;
+ do {
+ if (get_header(line,str)) {
+ b->nr++;
+ srenew(b->index,b->nr+1);
+ srenew(*grpname,b->nr);
+ if (b->nr==1)
+ b->index[0]=0;
+ b->index[b->nr]=b->index[b->nr-1];
+ (*grpname)[b->nr-1]=strdup(str);
+ } else {
++ if (b->nr==0)
++ {
++ gmx_fatal(FARGS,"The first header of your indexfile is invalid");
++ }
+ pt=line;
+ while (sscanf(pt,"%s",str) == 1) {
+ i=b->index[b->nr];
+ if (i>=maxentries) {
+ maxentries+=1024;
+ srenew(b->a,maxentries);
+ }
+ b->a[i]=strtol(str, NULL, 10)-1;
+ b->index[b->nr]++;
+ (b->nra)++;
+ pt=strstr(pt,str)+strlen(str);
+ }
+ }
+ } while (get_a_line(in,line,STRLEN));
+ }
+ else {
+ /* old format */
+ sscanf(line,"%d%d",&b->nr,&b->nra);
+ snew(b->index,b->nr+1);
+ snew(*grpname,b->nr);
+ b->index[0]=0;
+ snew(b->a,b->nra);
+ for (i=0; (i<b->nr); i++) {
+ nread=fscanf(in,"%s%d",str,&ng);
+ (*grpname)[i]=strdup(str);
+ b->index[i+1]=b->index[i]+ng;
+ if (b->index[i+1] > b->nra)
+ gmx_fatal(FARGS,"Something wrong in your indexfile at group %s",str);
+ for(j=0; (j<ng); j++) {
+ nread=fscanf(in,"%d",&a);
+ b->a[b->index[i]+j]=a;
+ }
+ }
+ }
+ gmx_fio_fclose(in);
+
+ for(i=0; (i<b->nr); i++) {
+ for(j=b->index[i]; (j<b->index[i+1]); j++) {
+ if (b->a[j] < 0)
+ fprintf(stderr,"\nWARNING: negative index %d in group %s\n\n",
+ b->a[j],(*grpname)[i]);
+ }
+ }
+
+ return b;
+}
+
+static void minstring(char *str)
+{
+ int i;
+
+ for (i=0; (i < (int)strlen(str)); i++)
+ if (str[i]=='-')
+ str[i]='_';
+}
+
+int find_group(char s[], int ngrps, char **grpname)
+{
+ int aa, i, n;
+ char string[STRLEN];
+ gmx_bool bMultiple;
+
+ bMultiple = FALSE;
+ n = strlen(s);
+ aa=NOTSET;
+ /* first look for whole name match */
+ if (aa==NOTSET)
+ for(i=0; i<ngrps; i++)
+ if (gmx_strcasecmp_min(s,grpname[i])==0) {
+ if(aa!=NOTSET)
+ bMultiple = TRUE;
+ aa=i;
+ }
+ /* second look for first string match */
+ if (aa==NOTSET)
+ for(i=0; i<ngrps; i++)
+ if (gmx_strncasecmp_min(s,grpname[i],n)==0) {
+ if(aa!=NOTSET)
+ bMultiple = TRUE;
+ aa=i;
+ }
+ /* last look for arbitrary substring match */
+ if (aa==NOTSET) {
+ upstring(s);
+ minstring(s);
+ for(i=0; i<ngrps; i++) {
+ strcpy(string, grpname[i]);
+ upstring(string);
+ minstring(string);
+ if (strstr(string,s)!=NULL) {
+ if(aa!=NOTSET)
+ bMultiple = TRUE;
+ aa=i;
+ }
+ }
+ }
+ if (bMultiple) {
+ printf("Error: Multiple groups '%s' selected\n", s);
+ aa=NOTSET;
+ }
+ return aa;
+}
+
+static int qgroup(int *a, int ngrps, char **grpname)
+{
+ char s[STRLEN];
+ int aa;
+ gmx_bool bInRange;
+ char *end;
+
+ do {
+ fprintf(stderr,"Select a group: ");
+ do {
+ if ( scanf("%s",s)!=1 )
+ gmx_fatal(FARGS,"Cannot read from input");
+ trim(s); /* remove spaces */
+ } while (strlen(s)==0);
+ aa = strtol(s, &end, 10);
+ if (aa==0 && end[0] != '\0') /* string entered */
+ aa = find_group(s, ngrps, grpname);
+ bInRange = (aa >= 0 && aa < ngrps);
+ if (!bInRange)
+ printf("Error: No such group '%s'\n", s);
+ } while (!bInRange);
+ printf("Selected %d: '%s'\n", aa, grpname[aa]);
+ *a = aa;
+ return aa;
+}
+
+static void rd_groups(t_blocka *grps,char **grpname,char *gnames[],
+ int ngrps,int isize[],atom_id *index[],int grpnr[])
+{
+ int i,j,gnr1;
+
+ if (grps->nr==0)
+ gmx_fatal(FARGS,"Error: no groups in indexfile");
+ for(i=0; (i<grps->nr); i++)
+ fprintf(stderr,"Group %5d (%15s) has %5d elements\n",i,grpname[i],
+ grps->index[i+1]-grps->index[i]);
+ for(i=0; (i<ngrps); i++) {
+ if (grps->nr > 1)
+ do {
+ gnr1=qgroup(&grpnr[i], grps->nr, grpname);
+ if ((gnr1<0) || (gnr1>=grps->nr))
+ fprintf(stderr,"Select between %d and %d.\n",0,grps->nr-1);
+ } while ((gnr1<0) || (gnr1>=grps->nr));
+ else {
+ fprintf(stderr,"There is one group in the index\n");
+ gnr1=0;
+ }
+ gnames[i]=strdup(grpname[gnr1]);
+ isize[i]=grps->index[gnr1+1]-grps->index[gnr1];
+ snew(index[i],isize[i]);
+ for(j=0; (j<isize[i]); j++)
+ index[i][j]=grps->a[grps->index[gnr1]+j];
+ }
+}
+
+void rd_index(const char *statfile,int ngrps,int isize[],
+ atom_id *index[],char *grpnames[])
+{
+ char **gnames;
+ t_blocka *grps;
+ int *grpnr;
+
+ snew(grpnr,ngrps);
+ if (!statfile)
+ gmx_fatal(FARGS,"No index file specified");
+ grps=init_index(statfile,&gnames);
+ rd_groups(grps,gnames,grpnames,ngrps,isize,index,grpnr);
+}
+
+void rd_index_nrs(char *statfile,int ngrps,int isize[],
+ atom_id *index[],char *grpnames[],int grpnr[])
+{
+ char **gnames;
+ t_blocka *grps;
+
+ if (!statfile)
+ gmx_fatal(FARGS,"No index file specified");
+ grps=init_index(statfile,&gnames);
+
+ rd_groups(grps,gnames,grpnames,ngrps,isize,index,grpnr);
+}
+
+void get_index(t_atoms *atoms, const char *fnm, int ngrps,
+ int isize[], atom_id *index[],char *grpnames[])
+{
+ char ***gnames;
+ t_blocka *grps = NULL;
+ int *grpnr;
+
+ snew(grpnr,ngrps);
+ snew(gnames,1);
+ if (fnm != NULL) {
+ grps=init_index(fnm,gnames);
+ }
+ else if (atoms) {
+ snew(grps,1);
+ snew(grps->index,1);
+ analyse(atoms,grps,gnames,FALSE,FALSE);
+ }
+ else
+ gmx_incons("You need to supply a valid atoms structure or a valid index file name");
+
+ rd_groups(grps,*gnames,grpnames,ngrps,isize,index,grpnr);
+}
+
+t_cluster_ndx *cluster_index(FILE *fplog,const char *ndx)
+{
+ t_cluster_ndx *c;
+ int i;
+
+ snew(c,1);
+ c->clust = init_index(ndx,&c->grpname);
+ c->maxframe = -1;
+ for(i=0; (i<c->clust->nra); i++)
+ c->maxframe = max(c->maxframe,c->clust->a[i]);
+ fprintf(fplog ? fplog : stdout,
+ "There are %d clusters containing %d structures, highest framenr is %d\n",
+ c->clust->nr,c->clust->nra,c->maxframe);
+ if (debug) {
+ pr_blocka(debug,0,"clust",c->clust,TRUE);
+ for(i=0; (i<c->clust->nra); i++)
+ if ((c->clust->a[i] < 0) || (c->clust->a[i] > c->maxframe))
+ gmx_fatal(FARGS,"Range check error for c->clust->a[%d] = %d\n"
+ "should be within 0 and %d",i,c->clust->a[i],c->maxframe+1);
+ }
+ c->inv_clust=make_invblocka(c->clust,c->maxframe);
+
+ return c;
+}
+
--- /dev/null
-
+/* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
+ *
+ *
+ * This source code is part of
+ *
+ * G R O M A C S
+ *
+ * GROningen MAchine for Chemical Simulations
+ *
+ * VERSION 3.2.0
+ * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
+ * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
+ * Copyright (c) 2001-2004, The GROMACS development team,
+ * check out http://www.gromacs.org for more information.
+
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version.
+ *
+ * If you want to redistribute modifications, please consider that
+ * scientific software is very special. Version control is crucial -
+ * bugs must be traceable. We will be happy to consider code for
+ * inclusion in the official distribution, but derived work must not
+ * be called official GROMACS. Details are found in the README & COPYING
+ * files - if they are missing, get the official version at www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the papers on the package - you can find them in the top README file.
+ *
+ * For more info, check our website at http://www.gromacs.org
+ *
+ * And Hey:
+ * GROningen Mixture of Alchemy and Childrens' Stories
+ */
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+
+#include <ctype.h>
++#include <assert.h>
+#include "sysstuff.h"
+#include "macros.h"
+#include "string2.h"
+#include "smalloc.h"
+#include "pbc.h"
+#include "statutil.h"
+#include "names.h"
+#include "vec.h"
+#include "futil.h"
+#include "wman.h"
+#include "tpxio.h"
+#include "gmx_fatal.h"
+#include "network.h"
+#include "vec.h"
+#include "mtop_util.h"
+#include "gmxfio.h"
+#include "oenv.h"
+
+#ifdef GMX_THREAD_MPI
+#include "thread_mpi.h"
+#endif
+
+struct output_env
+{
+ time_unit_t time_unit; /* the time unit, enum defined in oenv.h */
+ gmx_bool view; /* view of file requested */
+ xvg_format_t xvg_format; /* xvg output format, enum defined in oenv.h */
+ int verbosity; /* The level of verbosity for this program */
+ int debug_level; /* the debug level */
+
+ char *program_name; /* the program name */
+ char *cmd_line; /* the re-assembled command line */
+};
+
+/* The source code in this file should be thread-safe.
+ Please keep it that way. */
+
+/******************************************************************
+ *
+ * T R A J E C T O R Y S T U F F
+ *
+ ******************************************************************/
+
+/* read only time names */
+/* These must correspond to the time units type time_unit_t in oenv.h */
+static const real timefactors[] = { 0, 1e3, 1, 1e-3, 1e-6, 1e-9, 1e-12, 0 };
+static const real timeinvfactors[] ={ 0, 1e-3, 1, 1e3, 1e6, 1e9, 1e12, 0 };
+static const char *time_units_str[] = { NULL, "fs", "ps", "ns", "us",
+ "\\mus", "ms", "s" };
+static const char *time_units_xvgr[] = { NULL, "fs", "ps", "ns",
+ "ms", "s", NULL };
+
+
+
+/***** OUTPUT_ENV MEMBER FUNCTIONS ******/
+
+void output_env_init(output_env_t *oenvp, int argc, char *argv[],
+ time_unit_t tmu, gmx_bool view, xvg_format_t xvg_format,
+ int verbosity, int debug_level)
+{
+ int i;
+ int cmdlength=0;
+ char *argvzero=NULL;
+ output_env_t oenv;
+
+ snew(oenv, 1);
+ *oenvp = oenv;
+ oenv->time_unit = tmu;
+ oenv->view=view;
+ oenv->xvg_format = xvg_format;
+ oenv->verbosity=verbosity;
+ oenv->debug_level=debug_level;
+ oenv->program_name=NULL;
+
+ if (argv)
++ {
+ argvzero=argv[0];
++ assert(argvzero);
++ }
+ /* set program name */
+ if (argvzero)
+ {
+ oenv->program_name=strdup(argvzero);
+ }
+ if (oenv->program_name == NULL)
+ oenv->program_name = strdup("GROMACS");
+
+ /* copy command line */
+ if (argv)
+ {
+ cmdlength = strlen(argvzero);
+ for (i=1; i<argc; i++)
+ {
+ cmdlength += strlen(argv[i]);
+ }
+ }
+
+ /* Fill the cmdline string */
+ snew(oenv->cmd_line,cmdlength+argc+1);
+ if (argv)
+ {
+ for (i=0; i<argc; i++)
+ {
+ strcat(oenv->cmd_line,argv[i]);
+ strcat(oenv->cmd_line," ");
+ }
+ }
+}
+
+
+void output_env_init_default(output_env_t *oenvp)
+{
+ output_env_init(oenvp, 0, NULL, time_ps, FALSE, exvgNONE, 0, 0);
+}
+
+
+void output_env_done(output_env_t oenv)
+{
+ sfree(oenv->program_name);
+ sfree(oenv->cmd_line);
+ sfree(oenv);
+}
+
+
+
+int output_env_get_verbosity(const output_env_t oenv)
+{
+ return oenv->verbosity;
+}
+
+int output_env_get_debug_level(const output_env_t oenv)
+{
+ return oenv->debug_level;
+}
+
+
+const char *output_env_get_time_unit(const output_env_t oenv)
+{
+ return time_units_str[oenv->time_unit];
+}
+
+const char *output_env_get_time_label(const output_env_t oenv)
+{
+ char *label;
+ snew(label, 20);
+
+ sprintf(label,"Time (%s)",time_units_str[oenv->time_unit] ?
+ time_units_str[oenv->time_unit]: "ps");
+
+ return label;
+}
+
+const char *output_env_get_xvgr_tlabel(const output_env_t oenv)
+{
+ char *label;
+ snew(label, 20);
+
+ sprintf(label,"Time (%s)", time_units_xvgr[oenv->time_unit] ?
+ time_units_xvgr[oenv->time_unit] : "ps");
+
+ return label;
+}
+
+
+real output_env_get_time_factor(const output_env_t oenv)
+{
+ return timefactors[oenv->time_unit];
+}
+
+real output_env_get_time_invfactor(const output_env_t oenv)
+{
+ return timeinvfactors[oenv->time_unit];
+}
+
+real output_env_conv_time(const output_env_t oenv, real time)
+{
+ return time*timefactors[oenv->time_unit];
+}
+
+
+void output_env_conv_times(const output_env_t oenv, int n, real *time)
+{
+ int i;
+ double fact=timefactors[oenv->time_unit];
+
+ if (fact!=1.)
+ for(i=0; i<n; i++)
+ time[i] *= fact;
+}
+
+gmx_bool output_env_get_view(const output_env_t oenv)
+{
+ return oenv->view;
+}
+
+xvg_format_t output_env_get_xvg_format(const output_env_t oenv)
+{
+ return oenv->xvg_format;
+}
+
+const char *output_env_get_program_name(const output_env_t oenv)
+{
+ return oenv->program_name;
+}
+
+const char *output_env_get_short_program_name(const output_env_t oenv)
+{
+ const char *pr,*ret;
+ pr=ret=oenv->program_name;
+ if ((pr=strrchr(ret,DIR_SEPARATOR)) != NULL)
+ ret=pr+1;
+ return ret;
+}
+
+
+
+const char *output_env_get_cmd_line(const output_env_t oenv)
+{
+ return oenv->cmd_line;
+}
+
+
--- /dev/null
+/*
+ *
+ * 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:
+ * GROningen Mixture of Alchemy and Childrens' Stories
+ */
+/* This file is completely threadsafe - keep it that way! */
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+#include "gromacs/utility/gmx_header_config.h"
+
+#ifdef GMX_CRAY_XT3
+#undef HAVE_PWD_H
+#endif
+
+#include <stdio.h>
+#include <ctype.h>
+#include <stdlib.h>
+#include <errno.h>
+#include <sys/types.h>
+#include <time.h>
+
+#ifdef HAVE_SYS_TIME_H
+#include <sys/time.h>
+#endif
+
+#ifdef HAVE_UNISTD_H
+#include <unistd.h>
+#endif
+
+#ifdef HAVE_PWD_H
+#include <pwd.h>
+#endif
+#include <time.h>
++#include <assert.h>
+
+#include "typedefs.h"
+#include "smalloc.h"
+#include "gmx_fatal.h"
+#include "macros.h"
+#include "string2.h"
+#include "futil.h"
+
+int continuing(char *s)
+/* strip trailing spaces and if s ends with a CONTINUE remove that too.
+ * returns TRUE if s ends with a CONTINUE, FALSE otherwise.
+ */
+{
+ int sl;
++ assert(s);
+
+ rtrim(s);
+ sl = strlen(s);
+ if ((sl > 0) && (s[sl-1] == CONTINUE)) {
+ s[sl-1] = 0;
+ return TRUE;
+ }
+ else
+ return FALSE;
+}
+
+
+
+char *fgets2(char *line, int n, FILE *stream)
+/* This routine reads a string from stream of max length n
+ * and zero terminated, without newlines
+ * line should be long enough (>= n)
+ */
+{
+ char *c;
+ if (fgets(line,n,stream) == NULL) {
+ return NULL;
+ }
+ if ((c=strchr(line,'\n')) != NULL) {
+ *c = '\0';
+ } else {
+ /* A line not ending in a newline can only occur at the end of a file,
+ * or because of n being too small.
+ * Since both cases occur very infrequently, we can check for EOF.
+ */
+ if (!gmx_eof(stream)) {
+ gmx_fatal(FARGS,"An input file contains a line longer than %d characters, while the buffer passed to fgets2 has size %d. The line starts with: '%20.20s'",n,n,line);
+ }
+ }
+ if ((c=strchr(line,'\r')) != NULL) {
+ *c = '\0';
+ }
+
+ return line;
+}
+
+void strip_comment (char *line)
+{
+ char *c;
+
+ if (!line)
+ return;
+
+ /* search for a comment mark and replace it by a zero */
+ if ((c = strchr(line,COMMENTSIGN)) != NULL)
+ (*c) = 0;
+}
+
+void upstring (char *str)
+{
+ int i;
+
+ for (i=0; (i < (int)strlen(str)); i++)
+ str[i] = toupper(str[i]);
+}
+
+void ltrim (char *str)
+{
+ char *tr;
+ int i,c;
+
+ if (NULL == str)
+ return;
+
+ c = 0;
+ while (('\0' != str[c]) && isspace(str[c]))
+ c++;
+ if (c > 0)
+ {
+ for(i=c; ('\0' != str[i]); i++)
+ str[i-c] = str[i];
+ str[i-c] = '\0';
+ }
+}
+
+void rtrim (char *str)
+{
+ int nul;
+
+ if (NULL == str)
+ return;
+
+ nul = strlen(str)-1;
+ while ((nul > 0) && ((str[nul] == ' ') || (str[nul] == '\t')) ) {
+ str[nul] = '\0';
+ nul--;
+ }
+}
+
+void trim (char *str)
+{
+ ltrim (str);
+ rtrim (str);
+}
+
+char *
+gmx_ctime_r(const time_t *clock,char *buf, int n)
+{
+ char tmpbuf[STRLEN];
+
+#ifdef GMX_NATIVE_WINDOWS
+ /* Windows */
+ ctime_s( tmpbuf, STRLEN, clock );
+#elif (defined(__sun))
+ /*Solaris*/
+ ctime_r(clock, tmpbuf, n);
+#else
+ ctime_r(clock,tmpbuf);
+#endif
+ strncpy(buf,tmpbuf,n-1);
+ buf[n-1]='\0';
+
+ return buf;
+}
+
+void nice_header (FILE *out,const char *fn)
+{
+ const char *unk = "onbekend";
+ time_t clock;
+ const char *user=unk;
+ int gh;
+ uid_t uid;
+ char buf[256]="";
+ char timebuf[STRLEN];
+#ifdef HAVE_PWD_H
+ struct passwd *pw;
+#endif
+
+ /* Print a nice header above the file */
+ time(&clock);
+ fprintf (out,"%c\n",COMMENTSIGN);
+ fprintf (out,"%c\tFile '%s' was generated\n",COMMENTSIGN,fn ? fn : unk);
+
+#ifdef HAVE_PWD_H
+ uid = getuid();
+ pw = getpwuid(uid);
+ gh = gethostname(buf,255);
+ user= pw->pw_name;
+#else
+ uid = 0;
+ gh = -1;
+#endif
+
+ gmx_ctime_r(&clock,timebuf,STRLEN);
+ fprintf (out,"%c\tBy user: %s (%d)\n",COMMENTSIGN,
+ user ? user : unk,(int) uid);
+ fprintf(out,"%c\tOn host: %s\n",COMMENTSIGN,(gh == 0) ? buf : unk);
+
+ fprintf (out,"%c\tAt date: %s",COMMENTSIGN,timebuf);
+ fprintf (out,"%c\n",COMMENTSIGN);
+}
+
+int gmx_strcasecmp_min(const char *str1, const char *str2)
+{
+ char ch1,ch2;
+
+ do
+ {
+ do
+ ch1=toupper(*(str1++));
+ while ((ch1=='-') || (ch1=='_'));
+ do
+ ch2=toupper(*(str2++));
+ while ((ch2=='-') || (ch2=='_'));
+ if (ch1!=ch2) return (ch1-ch2);
+ }
+ while (ch1);
+ return 0;
+}
+
+int gmx_strncasecmp_min(const char *str1, const char *str2, int n)
+{
+ char ch1,ch2;
+ char *stri1, *stri2;
+
+ stri1=(char *)str1;
+ stri2=(char *)str2;
+ do
+ {
+ do
+ ch1=toupper(*(str1++));
+ while ((ch1=='-') || (ch1=='_'));
+ do
+ ch2=toupper(*(str2++));
+ while ((ch2=='-') || (ch2=='_'));
+ if (ch1!=ch2) return (ch1-ch2);
+ }
+ while (ch1 && (str1-stri1<n) && (str2-stri2<n));
+ return 0;
+}
+
+int gmx_strcasecmp(const char *str1, const char *str2)
+{
+ char ch1,ch2;
+
+ do
+ {
+ ch1=toupper(*(str1++));
+ ch2=toupper(*(str2++));
+ if (ch1!=ch2) return (ch1-ch2);
+ }
+ while (ch1);
+ return 0;
+}
+
+int gmx_strncasecmp(const char *str1, const char *str2, int n)
+{
+ char ch1,ch2;
+
+ if(n==0)
+ return 0;
+
+ do
+ {
+ ch1=toupper(*(str1++));
+ ch2=toupper(*(str2++));
+ if (ch1!=ch2) return (ch1-ch2);
+ n--;
+ }
+ while (ch1 && n);
+ return 0;
+}
+
+char *gmx_strdup(const char *src)
+{
+ char *dest;
+
+ snew(dest,strlen(src)+1);
+ strcpy(dest,src);
+
+ return dest;
+}
+
+char *
+gmx_strndup(const char *src, int n)
+{
+ int len;
+ char *dest;
+
+ len = strlen(src);
+ if (len > n)
+ {
+ len = n;
+ }
+ snew(dest, len+1);
+ strncpy(dest, src, len);
+ dest[len] = 0;
+ return dest;
+}
+
+/*!
+ * \param[in] pattern Pattern to match against.
+ * \param[in] str String to match.
+ * \returns 0 on match, GMX_NO_WCMATCH if there is no match.
+ *
+ * Matches \p str against \p pattern, which may contain * and ? wildcards.
+ * All other characters are matched literally.
+ * Currently, it is not possible to match literal * or ?.
+ */
+int
+gmx_wcmatch(const char *pattern, const char *str)
+{
+ while (*pattern)
+ {
+ if (*pattern == '*')
+ {
+ /* Skip multiple wildcards in a sequence */
+ while (*pattern == '*' || *pattern == '?')
+ {
+ ++pattern;
+ /* For ?, we need to check that there are characters left
+ * in str. */
+ if (*pattern == '?')
+ {
+ if (*str == 0)
+ {
+ return GMX_NO_WCMATCH;
+ }
+ else
+ {
+ ++str;
+ }
+ }
+ }
+ /* If the pattern ends after the star, we have a match */
+ if (*pattern == 0)
+ {
+ return 0;
+ }
+ /* Match the rest against each possible suffix of str */
+ while (*str)
+ {
+ /* Only do the recursive call if the first character
+ * matches. We don't have to worry about wildcards here,
+ * since we have processed them above. */
+ if (*pattern == *str)
+ {
+ int rc;
+ /* Match the suffix, and return if a match or an error */
+ rc = gmx_wcmatch(pattern, str);
+ if (rc != GMX_NO_WCMATCH)
+ {
+ return rc;
+ }
+ }
+ ++str;
+ }
+ /* If no suffix of str matches, we don't have a match */
+ return GMX_NO_WCMATCH;
+ }
+ else if ((*pattern == '?' && *str != 0) || *pattern == *str)
+ {
+ ++str;
+ }
+ else
+ {
+ return GMX_NO_WCMATCH;
+ }
+ ++pattern;
+ }
+ /* When the pattern runs out, we have a match if the string has ended. */
+ return (*str == 0) ? 0 : GMX_NO_WCMATCH;
+}
+
+char *wrap_lines(const char *buf,int line_width, int indent,gmx_bool bIndentFirst)
+{
+ char *b2;
+ int i,i0,i2,j,b2len,lspace=0,l2space=0;
+ gmx_bool bFirst,bFitsOnLine;
+
+ /* characters are copied from buf to b2 with possible spaces changed
+ * into newlines and extra space added for indentation.
+ * i indexes buf (source buffer) and i2 indexes b2 (destination buffer)
+ * i0 points to the beginning of the current line (in buf, source)
+ * lspace and l2space point to the last space on the current line
+ * bFirst is set to prevent indentation of first line
+ * bFitsOnLine says if the first space occurred before line_width, if
+ * that is not the case, we have a word longer than line_width which
+ * will also not fit on the next line, so we might as well keep it on
+ * the current line (where it also won't fit, but looks better)
+ */
+
+ b2=NULL;
+ b2len=strlen(buf)+1+indent;
+ snew(b2,b2len);
+ i0=i2=0;
+ if (bIndentFirst)
+ for(i2=0; (i2<indent); i2++)
+ b2[i2] = ' ';
+ bFirst=TRUE;
+ do {
+ l2space = -1;
+ /* find the last space before end of line */
+ for(i=i0; ((i-i0 < line_width) || (l2space==-1)) && (buf[i]); i++) {
+ b2[i2++] = buf[i];
+ /* remember the position of a space */
+ if (buf[i] == ' ') {
+ lspace = i;
+ l2space = i2-1;
+ }
+ /* if we have a newline before the line is full, reset counters */
+ if (buf[i]=='\n' && buf[i+1]) {
+ i0=i+1;
+ b2len+=indent;
+ srenew(b2, b2len);
+ /* add indentation after the newline */
+ for(j=0; (j<indent); j++)
+ b2[i2++]=' ';
+ }
+ }
+ /* If we are at the last newline, copy it */
+ if (buf[i]=='\n' && !buf[i+1]) {
+ b2[i2++] = buf[i++];
+ }
+ /* if we're not at the end of the string */
+ if (buf[i]) {
+ /* check if one word does not fit on the line */
+ bFitsOnLine = (i-i0 <= line_width);
+ /* reset line counters to just after the space */
+ i0 = lspace+1;
+ i2 = l2space+1;
+ /* if the words fit on the line, and we're beyond the indentation part */
+ if ( (bFitsOnLine) && (l2space >= indent) ) {
+ /* start a new line */
+ b2[l2space] = '\n';
+ /* and add indentation */
+ if (indent) {
+ if (bFirst) {
+ line_width-=indent;
+ bFirst=FALSE;
+ }
+ b2len+=indent;
+ srenew(b2, b2len);
+ for(j=0; (j<indent); j++)
+ b2[i2++]=' ';
+ /* no extra spaces after indent; */
+ while(buf[i0]==' ')
+ i0++;
+ }
+ }
+ }
+ } while (buf[i]);
+ b2[i2] = '\0';
+
+ return b2;
+}
+
+char **split(char sep,const char *str)
+{
+ char **ptr = NULL;
+ int n,nn,nptr = 0;
+
+ if (str == NULL)
+ return NULL;
+ nn = strlen(str);
+ for(n=0; (n<nn); n++)
+ if (str[n] == sep)
+ nptr++;
+ snew(ptr,nptr+2);
+ nptr = 0;
+ while (*str != '\0') {
+ while ((*str != '\0') && (*str == sep))
+ str++;
+ if (*str != '\0') {
+ snew(ptr[nptr],1+strlen(str));
+ n = 0;
+ while ((*str != '\0') && (*str != sep)) {
+ ptr[nptr][n] = *str;
+ str++;
+ n++;
+ }
+ ptr[nptr][n] = '\0';
+ nptr++;
+ }
+ }
+ ptr[nptr] = NULL;
+
+ return ptr;
+}
+
+
+gmx_large_int_t
+str_to_large_int_t(const char *str, char **endptr)
+{
+ int sign = 1;
+ gmx_large_int_t val = 0;
+ char ch;
+ const char *p;
+
+ p = str;
+ if(p==NULL)
+ {
+ *endptr=NULL;
+ return 0;
+ }
+
+ /* Strip off initial white space */
+ while(isspace(*p))
+ {
+ p++;
+ }
+ /* Conform to ISO C99 - return original pointer if string does not contain a number */
+ if(*str=='\0')
+ {
+ *endptr=(char *)str;
+ }
+
+ if(*p=='-')
+ {
+ p++;
+ sign *= -1;
+ }
+
+ while( ((ch=*p) != '\0') && isdigit(ch) )
+ {
+ /* Important to add sign here, so we dont overflow in final multiplication */
+ ch = (ch-'0')*sign;
+ val = val*10 + ch;
+ if(ch != val%10)
+ {
+ /* Some sort of overflow has occured, set endptr to original string */
+ *endptr=(char *)str;
+ errno = ERANGE;
+ return(0);
+ }
+ p++;
+ }
+
+ *endptr=(char *)p;
+
+ return val;
+}
+
+char *gmx_strsep(char **stringp, const char *delim)
+{
+ char *ret;
+ int len=strlen(delim);
+ int i,j=0;
+ int found=0;
+
+ if (! *stringp)
+ return NULL;
+ ret=*stringp;
+ do
+ {
+ if ( (*stringp)[j] == '\0')
+ {
+ found=1;
+ *stringp=NULL;
+ break;
+ }
+ for (i=0;i<len;i++)
+ {
+ if ( (*stringp)[j]==delim[i])
+ {
+ (*stringp)[j]='\0';
+ *stringp=*stringp+j+1;
+ found=1;
+ break;
+ }
+ }
+ j++;
+ } while (!found);
+
+ return ret;
+}
+
--- /dev/null
+/* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
+ *
+ * This source code is part of
+ *
+ * G R O M A C S
+ *
+ * GROningen MAchine for Chemical Simulations
+ *
+ * VERSION 3.2.0
+ * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
+ * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
+ * Copyright (c) 2001-2004, The GROMACS development team,
+ * check out http://www.gromacs.org for more information.
+
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version.
+ *
+ * If you want to redistribute modifications, please consider that
+ * scientific software is very special. Version control is crucial -
+ * bugs must be traceable. We will be happy to consider code for
+ * inclusion in the official distribution, but derived work must not
+ * be called official GROMACS. Details are found in the README & COPYING
+ * files - if they are missing, get the official version at www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the papers on the package - you can find them in the top README file.
+ *
+ * For more info, check our website at http://www.gromacs.org
+ *
+ * And Hey:
+ * GROwing Monsters And Cloning Shrimps
+ */
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
++#include <assert.h>
+
+#include "typedefs.h"
+#include "smalloc.h"
+#include "update.h"
+#include "vec.h"
+#include "macros.h"
+#include "physics.h"
+#include "names.h"
+#include "gmx_fatal.h"
+#include "txtdump.h"
+#include "nrnb.h"
+#include "gmx_random.h"
+#include "update.h"
+#include "mdrun.h"
+
+#define NTROTTERPARTS 3
+
+/* Suzuki-Yoshida Constants, for n=3 and n=5, for symplectic integration */
+/* for n=1, w0 = 1 */
+/* for n=3, w0 = w2 = 1/(2-2^-(1/3)), w1 = 1-2*w0 */
+/* for n=5, w0 = w1 = w3 = w4 = 1/(4-4^-(1/3)), w1 = 1-4*w0 */
+
+#define MAX_SUZUKI_YOSHIDA_NUM 5
+#define SUZUKI_YOSHIDA_NUM 5
+
+static const double sy_const_1[] = { 1. };
+static const double sy_const_3[] = { 0.828981543588751,-0.657963087177502,0.828981543588751 };
+static const double sy_const_5[] = { 0.2967324292201065,0.2967324292201065,-0.186929716880426,0.2967324292201065,0.2967324292201065 };
+
+static const double* sy_const[] = {
+ NULL,
+ sy_const_1,
+ NULL,
+ sy_const_3,
+ NULL,
+ sy_const_5
+};
+
+/*
+static const double sy_const[MAX_SUZUKI_YOSHIDA_NUM+1][MAX_SUZUKI_YOSHIDA_NUM+1] = {
+ {},
+ {1},
+ {},
+ {0.828981543588751,-0.657963087177502,0.828981543588751},
+ {},
+ {0.2967324292201065,0.2967324292201065,-0.186929716880426,0.2967324292201065,0.2967324292201065}
+};*/
+
+/* these integration routines are only referenced inside this file */
+static void NHC_trotter(t_grpopts *opts,int nvar, gmx_ekindata_t *ekind,real dtfull,
+ double xi[],double vxi[], double scalefac[], real *veta, t_extmass *MassQ, gmx_bool bEkinAveVel)
+
+{
+ /* general routine for both barostat and thermostat nose hoover chains */
+
+ int i,j,mi,mj,jmax;
+ double Ekin,Efac,reft,kT,nd;
+ double dt;
+ t_grp_tcstat *tcstat;
+ double *ivxi,*ixi;
+ double *iQinv;
+ double *GQ;
+ gmx_bool bBarostat;
+ int mstepsi, mstepsj;
+ int ns = SUZUKI_YOSHIDA_NUM; /* set the degree of integration in the types/state.h file */
+ int nh = opts->nhchainlength;
+
+ snew(GQ,nh);
+ mstepsi = mstepsj = ns;
+
+/* if scalefac is NULL, we are doing the NHC of the barostat */
+
+ bBarostat = FALSE;
+ if (scalefac == NULL) {
+ bBarostat = TRUE;
+ }
+
+ for (i=0; i<nvar; i++)
+ {
+
+ /* make it easier to iterate by selecting
+ out the sub-array that corresponds to this T group */
+
+ ivxi = &vxi[i*nh];
+ ixi = &xi[i*nh];
+ if (bBarostat) {
+ iQinv = &(MassQ->QPinv[i*nh]);
+ nd = 1; /* THIS WILL CHANGE IF NOT ISOTROPIC */
+ reft = max(0.0,opts->ref_t[0]);
+ Ekin = sqr(*veta)/MassQ->Winv;
+ } else {
+ iQinv = &(MassQ->Qinv[i*nh]);
+ tcstat = &ekind->tcstat[i];
+ nd = opts->nrdf[i];
+ reft = max(0.0,opts->ref_t[i]);
+ if (bEkinAveVel)
+ {
+ Ekin = 2*trace(tcstat->ekinf)*tcstat->ekinscalef_nhc;
+ } else {
+ Ekin = 2*trace(tcstat->ekinh)*tcstat->ekinscaleh_nhc;
+ }
+ }
+ kT = BOLTZ*reft;
+
+ for(mi=0;mi<mstepsi;mi++)
+ {
+ for(mj=0;mj<mstepsj;mj++)
+ {
+ /* weighting for this step using Suzuki-Yoshida integration - fixed at 5 */
+ dt = sy_const[ns][mj] * dtfull / mstepsi;
+
+ /* compute the thermal forces */
+ GQ[0] = iQinv[0]*(Ekin - nd*kT);
+
+ for (j=0;j<nh-1;j++)
+ {
+ if (iQinv[j+1] > 0) {
+ /* we actually don't need to update here if we save the
+ state of the GQ, but it's easier to just recompute*/
+ GQ[j+1] = iQinv[j+1]*((sqr(ivxi[j])/iQinv[j])-kT);
+ } else {
+ GQ[j+1] = 0;
+ }
+ }
+
+ ivxi[nh-1] += 0.25*dt*GQ[nh-1];
+ for (j=nh-1;j>0;j--)
+ {
+ Efac = exp(-0.125*dt*ivxi[j]);
+ ivxi[j-1] = Efac*(ivxi[j-1]*Efac + 0.25*dt*GQ[j-1]);
+ }
+
+ Efac = exp(-0.5*dt*ivxi[0]);
+ if (bBarostat) {
+ *veta *= Efac;
+ } else {
+ scalefac[i] *= Efac;
+ }
+ Ekin *= (Efac*Efac);
+
+ /* Issue - if the KE is an average of the last and the current temperatures, then we might not be
+ able to scale the kinetic energy directly with this factor. Might take more bookkeeping -- have to
+ think about this a bit more . . . */
+
+ GQ[0] = iQinv[0]*(Ekin - nd*kT);
+
+ /* update thermostat positions */
+ for (j=0;j<nh;j++)
+ {
+ ixi[j] += 0.5*dt*ivxi[j];
+ }
+
+ for (j=0;j<nh-1;j++)
+ {
+ Efac = exp(-0.125*dt*ivxi[j+1]);
+ ivxi[j] = Efac*(ivxi[j]*Efac + 0.25*dt*GQ[j]);
+ if (iQinv[j+1] > 0) {
+ GQ[j+1] = iQinv[j+1]*((sqr(ivxi[j])/iQinv[j])-kT);
+ } else {
+ GQ[j+1] = 0;
+ }
+ }
+ ivxi[nh-1] += 0.25*dt*GQ[nh-1];
+ }
+ }
+ }
+ sfree(GQ);
+}
+
+static void boxv_trotter(t_inputrec *ir, real *veta, real dt, tensor box,
+ gmx_ekindata_t *ekind, tensor vir, real pcorr, real ecorr, t_extmass *MassQ)
+{
+
+ real pscal;
+ double alpha;
+ int i,j,d,n,nwall;
+ real T,GW,vol;
+ tensor Winvm,ekinmod,localpres;
+
+ /* The heat bath is coupled to a separate barostat, the last temperature group. In the
+ 2006 Tuckerman et al paper., the order is iL_{T_baro} iL {T_part}
+ */
+
+ if (ir->epct==epctSEMIISOTROPIC)
+ {
+ nwall = 2;
+ }
+ else
+ {
+ nwall = 3;
+ }
+
+ /* eta is in pure units. veta is in units of ps^-1. GW is in
+ units of ps^-2. However, eta has a reference of 1 nm^3, so care must be
+ taken to use only RATIOS of eta in updating the volume. */
+
+ /* we take the partial pressure tensors, modify the
+ kinetic energy tensor, and recovert to pressure */
+
+ if (ir->opts.nrdf[0]==0)
+ {
+ gmx_fatal(FARGS,"Barostat is coupled to a T-group with no degrees of freedom\n");
+ }
+ /* alpha factor for phase space volume, then multiply by the ekin scaling factor. */
+ alpha = 1.0 + DIM/((double)ir->opts.nrdf[0]);
+ alpha *= ekind->tcstat[0].ekinscalef_nhc;
+ msmul(ekind->ekin,alpha,ekinmod);
+
+ pscal = calc_pres(ir->ePBC,nwall,box,ekinmod,vir,localpres) + pcorr;
+
+ vol = det(box);
+ GW = (vol*(MassQ->Winv/PRESFAC))*(DIM*pscal - trace(ir->ref_p)); /* W is in ps^2 * bar * nm^3 */
+
+ *veta += 0.5*dt*GW;
+}
+
+/*
+ * This file implements temperature and pressure coupling algorithms:
+ * For now only the Weak coupling and the modified weak coupling.
+ *
+ * Furthermore computation of pressure and temperature is done here
+ *
+ */
+
+real calc_pres(int ePBC,int nwall,matrix box,tensor ekin,tensor vir,
+ tensor pres)
+{
+ int n,m;
+ real fac;
+
+ if (ePBC==epbcNONE || (ePBC==epbcXY && nwall!=2))
+ clear_mat(pres);
+ else {
+ /* Uitzoeken welke ekin hier van toepassing is, zie Evans & Morris - E.
+ * Wrs. moet de druktensor gecorrigeerd worden voor de netto stroom in
+ * het systeem...
+ */
+
+ fac=PRESFAC*2.0/det(box);
+ for(n=0; (n<DIM); n++)
+ for(m=0; (m<DIM); m++)
+ pres[n][m] = (ekin[n][m] - vir[n][m])*fac;
+
+ if (debug) {
+ pr_rvecs(debug,0,"PC: pres",pres,DIM);
+ pr_rvecs(debug,0,"PC: ekin",ekin,DIM);
+ pr_rvecs(debug,0,"PC: vir ",vir, DIM);
+ pr_rvecs(debug,0,"PC: box ",box, DIM);
+ }
+ }
+ return trace(pres)/DIM;
+}
+
+real calc_temp(real ekin,real nrdf)
+{
+ if (nrdf > 0)
+ return (2.0*ekin)/(nrdf*BOLTZ);
+ else
+ return 0;
+}
+
+void parrinellorahman_pcoupl(FILE *fplog,gmx_large_int_t step,
+ t_inputrec *ir,real dt,tensor pres,
+ tensor box,tensor box_rel,tensor boxv,
+ tensor M,matrix mu,gmx_bool bFirstStep)
+{
+ /* This doesn't do any coordinate updating. It just
+ * integrates the box vector equations from the calculated
+ * acceleration due to pressure difference. We also compute
+ * the tensor M which is used in update to couple the particle
+ * coordinates to the box vectors.
+ *
+ * In Nose and Klein (Mol.Phys 50 (1983) no 5., p 1055) this is
+ * given as
+ * -1 . . -1
+ * M_nk = (h') * (h' * h + h' h) * h
+ *
+ * with the dots denoting time derivatives and h is the transformation from
+ * the scaled frame to the real frame, i.e. the TRANSPOSE of the box.
+ * This also goes for the pressure and M tensors - they are transposed relative
+ * to ours. Our equation thus becomes:
+ *
+ * -1 . . -1
+ * M_gmx = M_nk' = b * (b * b' + b * b') * b'
+ *
+ * where b is the gromacs box matrix.
+ * Our box accelerations are given by
+ * .. ..
+ * b = vol/W inv(box') * (P-ref_P) (=h')
+ */
+
+ int d,n;
+ tensor winv;
+ real vol=box[XX][XX]*box[YY][YY]*box[ZZ][ZZ];
+ real atot,arel,change,maxchange,xy_pressure;
+ tensor invbox,pdiff,t1,t2;
+
+ real maxl;
+
+ m_inv_ur0(box,invbox);
+
+ if (!bFirstStep) {
+ /* Note that PRESFAC does not occur here.
+ * The pressure and compressibility always occur as a product,
+ * therefore the pressure unit drops out.
+ */
+ maxl=max(box[XX][XX],box[YY][YY]);
+ maxl=max(maxl,box[ZZ][ZZ]);
+ for(d=0;d<DIM;d++)
+ for(n=0;n<DIM;n++)
+ winv[d][n]=
+ (4*M_PI*M_PI*ir->compress[d][n])/(3*ir->tau_p*ir->tau_p*maxl);
+
+ m_sub(pres,ir->ref_p,pdiff);
+
+ if(ir->epct==epctSURFACETENSION) {
+ /* Unlike Berendsen coupling it might not be trivial to include a z
+ * pressure correction here? On the other hand we don't scale the
+ * box momentarily, but change accelerations, so it might not be crucial.
+ */
+ xy_pressure=0.5*(pres[XX][XX]+pres[YY][YY]);
+ for(d=0;d<ZZ;d++)
+ pdiff[d][d]=(xy_pressure-(pres[ZZ][ZZ]-ir->ref_p[d][d]/box[d][d]));
+ }
+
+ tmmul(invbox,pdiff,t1);
+ /* Move the off-diagonal elements of the 'force' to one side to ensure
+ * that we obey the box constraints.
+ */
+ for(d=0;d<DIM;d++) {
+ for(n=0;n<d;n++) {
+ t1[d][n] += t1[n][d];
+ t1[n][d] = 0;
+ }
+ }
+
+ switch (ir->epct) {
+ case epctANISOTROPIC:
+ for(d=0;d<DIM;d++)
+ for(n=0;n<=d;n++)
+ t1[d][n] *= winv[d][n]*vol;
+ break;
+ case epctISOTROPIC:
+ /* calculate total volume acceleration */
+ atot=box[XX][XX]*box[YY][YY]*t1[ZZ][ZZ]+
+ box[XX][XX]*t1[YY][YY]*box[ZZ][ZZ]+
+ t1[XX][XX]*box[YY][YY]*box[ZZ][ZZ];
+ arel=atot/(3*vol);
+ /* set all RELATIVE box accelerations equal, and maintain total V
+ * change speed */
+ for(d=0;d<DIM;d++)
+ for(n=0;n<=d;n++)
+ t1[d][n] = winv[0][0]*vol*arel*box[d][n];
+ break;
+ case epctSEMIISOTROPIC:
+ case epctSURFACETENSION:
+ /* Note the correction to pdiff above for surftens. coupling */
+
+ /* calculate total XY volume acceleration */
+ atot=box[XX][XX]*t1[YY][YY]+t1[XX][XX]*box[YY][YY];
+ arel=atot/(2*box[XX][XX]*box[YY][YY]);
+ /* set RELATIVE XY box accelerations equal, and maintain total V
+ * change speed. Dont change the third box vector accelerations */
+ for(d=0;d<ZZ;d++)
+ for(n=0;n<=d;n++)
+ t1[d][n] = winv[d][n]*vol*arel*box[d][n];
+ for(n=0;n<DIM;n++)
+ t1[ZZ][n] *= winv[d][n]*vol;
+ break;
+ default:
+ gmx_fatal(FARGS,"Parrinello-Rahman pressure coupling type %s "
+ "not supported yet\n",EPCOUPLTYPETYPE(ir->epct));
+ break;
+ }
+
+ maxchange=0;
+ for(d=0;d<DIM;d++)
+ for(n=0;n<=d;n++) {
+ boxv[d][n] += dt*t1[d][n];
+
+ /* We do NOT update the box vectors themselves here, since
+ * we need them for shifting later. It is instead done last
+ * in the update() routine.
+ */
+
+ /* Calculate the change relative to diagonal elements-
+ since it's perfectly ok for the off-diagonal ones to
+ be zero it doesn't make sense to check the change relative
+ to its current size.
+ */
+
+ change=fabs(dt*boxv[d][n]/box[d][d]);
+
+ if (change>maxchange)
+ maxchange=change;
+ }
+
+ if (maxchange > 0.01 && fplog) {
+ char buf[22];
+ fprintf(fplog,
+ "\nStep %s Warning: Pressure scaling more than 1%%. "
+ "This may mean your system\n is not yet equilibrated. "
+ "Use of Parrinello-Rahman pressure coupling during\n"
+ "equilibration can lead to simulation instability, "
+ "and is discouraged.\n",
+ gmx_step_str(step,buf));
+ }
+ }
+
+ preserve_box_shape(ir,box_rel,boxv);
+
+ mtmul(boxv,box,t1); /* t1=boxv * b' */
+ mmul(invbox,t1,t2);
+ mtmul(t2,invbox,M);
+
+ /* Determine the scaling matrix mu for the coordinates */
+ for(d=0;d<DIM;d++)
+ for(n=0;n<=d;n++)
+ t1[d][n] = box[d][n] + dt*boxv[d][n];
+ preserve_box_shape(ir,box_rel,t1);
+ /* t1 is the box at t+dt, determine mu as the relative change */
+ mmul_ur0(invbox,t1,mu);
+}
+
+void berendsen_pcoupl(FILE *fplog,gmx_large_int_t step,
+ t_inputrec *ir,real dt, tensor pres,matrix box,
+ matrix mu)
+{
+ int d,n;
+ real scalar_pressure, xy_pressure, p_corr_z;
+ char *ptr,buf[STRLEN];
+
+ /*
+ * Calculate the scaling matrix mu
+ */
+ scalar_pressure=0;
+ xy_pressure=0;
+ for(d=0; d<DIM; d++) {
+ scalar_pressure += pres[d][d]/DIM;
+ if (d != ZZ)
+ xy_pressure += pres[d][d]/(DIM-1);
+ }
+ /* Pressure is now in bar, everywhere. */
+#define factor(d,m) (ir->compress[d][m]*dt/ir->tau_p)
+
+ /* mu has been changed from pow(1+...,1/3) to 1+.../3, since this is
+ * necessary for triclinic scaling
+ */
+ clear_mat(mu);
+ switch (ir->epct) {
+ case epctISOTROPIC:
+ for(d=0; d<DIM; d++)
+ {
+ mu[d][d] = 1.0 - factor(d,d)*(ir->ref_p[d][d] - scalar_pressure) /DIM;
+ }
+ break;
+ case epctSEMIISOTROPIC:
+ for(d=0; d<ZZ; d++)
+ mu[d][d] = 1.0 - factor(d,d)*(ir->ref_p[d][d]-xy_pressure)/DIM;
+ mu[ZZ][ZZ] =
+ 1.0 - factor(ZZ,ZZ)*(ir->ref_p[ZZ][ZZ] - pres[ZZ][ZZ])/DIM;
+ break;
+ case epctANISOTROPIC:
+ for(d=0; d<DIM; d++)
+ for(n=0; n<DIM; n++)
+ mu[d][n] = (d==n ? 1.0 : 0.0)
+ -factor(d,n)*(ir->ref_p[d][n] - pres[d][n])/DIM;
+ break;
+ case epctSURFACETENSION:
+ /* ir->ref_p[0/1] is the reference surface-tension times *
+ * the number of surfaces */
+ if (ir->compress[ZZ][ZZ])
+ p_corr_z = dt/ir->tau_p*(ir->ref_p[ZZ][ZZ] - pres[ZZ][ZZ]);
+ else
+ /* when the compressibity is zero, set the pressure correction *
+ * in the z-direction to zero to get the correct surface tension */
+ p_corr_z = 0;
+ mu[ZZ][ZZ] = 1.0 - ir->compress[ZZ][ZZ]*p_corr_z;
+ for(d=0; d<DIM-1; d++)
+ mu[d][d] = 1.0 + factor(d,d)*(ir->ref_p[d][d]/(mu[ZZ][ZZ]*box[ZZ][ZZ])
+ - (pres[ZZ][ZZ]+p_corr_z - xy_pressure))/(DIM-1);
+ break;
+ default:
+ gmx_fatal(FARGS,"Berendsen pressure coupling type %s not supported yet\n",
+ EPCOUPLTYPETYPE(ir->epct));
+ break;
+ }
+ /* To fullfill the orientation restrictions on triclinic boxes
+ * we will set mu_yx, mu_zx and mu_zy to 0 and correct
+ * the other elements of mu to first order.
+ */
+ mu[YY][XX] += mu[XX][YY];
+ mu[ZZ][XX] += mu[XX][ZZ];
+ mu[ZZ][YY] += mu[YY][ZZ];
+ mu[XX][YY] = 0;
+ mu[XX][ZZ] = 0;
+ mu[YY][ZZ] = 0;
+
+ if (debug) {
+ pr_rvecs(debug,0,"PC: pres ",pres,3);
+ pr_rvecs(debug,0,"PC: mu ",mu,3);
+ }
+
+ if (mu[XX][XX]<0.99 || mu[XX][XX]>1.01 ||
+ mu[YY][YY]<0.99 || mu[YY][YY]>1.01 ||
+ mu[ZZ][ZZ]<0.99 || mu[ZZ][ZZ]>1.01) {
+ char buf2[22];
+ sprintf(buf,"\nStep %s Warning: pressure scaling more than 1%%, "
+ "mu: %g %g %g\n",
+ gmx_step_str(step,buf2),mu[XX][XX],mu[YY][YY],mu[ZZ][ZZ]);
+ if (fplog)
+ fprintf(fplog,"%s",buf);
+ fprintf(stderr,"%s",buf);
+ }
+}
+
+void berendsen_pscale(t_inputrec *ir,matrix mu,
+ matrix box,matrix box_rel,
+ int start,int nr_atoms,
+ rvec x[],unsigned short cFREEZE[],
+ t_nrnb *nrnb)
+{
+ ivec *nFreeze=ir->opts.nFreeze;
+ int n,d,g=0;
+
+ /* Scale the positions */
+ for (n=start; n<start+nr_atoms; n++) {
+ if (cFREEZE)
+ g = cFREEZE[n];
+
+ if (!nFreeze[g][XX])
+ x[n][XX] = mu[XX][XX]*x[n][XX]+mu[YY][XX]*x[n][YY]+mu[ZZ][XX]*x[n][ZZ];
+ if (!nFreeze[g][YY])
+ x[n][YY] = mu[YY][YY]*x[n][YY]+mu[ZZ][YY]*x[n][ZZ];
+ if (!nFreeze[g][ZZ])
+ x[n][ZZ] = mu[ZZ][ZZ]*x[n][ZZ];
+ }
+ /* compute final boxlengths */
+ for (d=0; d<DIM; d++) {
+ box[d][XX] = mu[XX][XX]*box[d][XX]+mu[YY][XX]*box[d][YY]+mu[ZZ][XX]*box[d][ZZ];
+ box[d][YY] = mu[YY][YY]*box[d][YY]+mu[ZZ][YY]*box[d][ZZ];
+ box[d][ZZ] = mu[ZZ][ZZ]*box[d][ZZ];
+ }
+
+ preserve_box_shape(ir,box_rel,box);
+
+ /* (un)shifting should NOT be done after this,
+ * since the box vectors might have changed
+ */
+ inc_nrnb(nrnb,eNR_PCOUPL,nr_atoms);
+}
+
+void berendsen_tcoupl(t_inputrec *ir,gmx_ekindata_t *ekind,real dt)
+{
+ t_grpopts *opts;
+ int i;
+ real T,reft=0,lll;
+
+ opts = &ir->opts;
+
+ for(i=0; (i<opts->ngtc); i++)
+ {
+ if (ir->eI == eiVV)
+ {
+ T = ekind->tcstat[i].T;
+ }
+ else
+ {
+ T = ekind->tcstat[i].Th;
+ }
+
+ if ((opts->tau_t[i] > 0) && (T > 0.0)) {
+
+ reft = max(0.0,opts->ref_t[i]);
+ lll = sqrt(1.0 + (dt/opts->tau_t[i])*(reft/T-1.0));
+ ekind->tcstat[i].lambda = max(min(lll,1.25),0.8);
+ }
+ else {
+ ekind->tcstat[i].lambda = 1.0;
+ }
+
+ if (debug)
+ fprintf(debug,"TC: group %d: T: %g, Lambda: %g\n",
+ i,T,ekind->tcstat[i].lambda);
+ }
+}
+
+void nosehoover_tcoupl(t_grpopts *opts,gmx_ekindata_t *ekind,real dt,
+ double xi[],double vxi[], t_extmass *MassQ)
+{
+ int i;
+ real reft,oldvxi;
+
+ /* note that this routine does not include Nose-hoover chains yet. Should be easy to add. */
+
+ for(i=0; (i<opts->ngtc); i++) {
+ reft = max(0.0,opts->ref_t[i]);
+ oldvxi = vxi[i];
+ vxi[i] += dt*MassQ->Qinv[i]*(ekind->tcstat[i].Th - reft);
+ xi[i] += dt*(oldvxi + vxi[i])*0.5;
+ }
+}
+
+t_state *init_bufstate(const t_state *template_state)
+{
+ t_state *state;
+ int nc = template_state->nhchainlength;
+ snew(state,1);
+ snew(state->nosehoover_xi,nc*template_state->ngtc);
+ snew(state->nosehoover_vxi,nc*template_state->ngtc);
+ snew(state->therm_integral,template_state->ngtc);
+ snew(state->nhpres_xi,nc*template_state->nnhpres);
+ snew(state->nhpres_vxi,nc*template_state->nnhpres);
+
+ return state;
+}
+
+void destroy_bufstate(t_state *state)
+{
+ sfree(state->x);
+ sfree(state->v);
+ sfree(state->nosehoover_xi);
+ sfree(state->nosehoover_vxi);
+ sfree(state->therm_integral);
+ sfree(state->nhpres_xi);
+ sfree(state->nhpres_vxi);
+ sfree(state);
+}
+
+void trotter_update(t_inputrec *ir,gmx_large_int_t step, gmx_ekindata_t *ekind,
+ gmx_enerdata_t *enerd, t_state *state,
+ tensor vir, t_mdatoms *md,
+ t_extmass *MassQ, int **trotter_seqlist, int trotter_seqno)
+{
+
+ int n,i,j,d,ntgrp,ngtc,gc=0;
+ t_grp_tcstat *tcstat;
+ t_grpopts *opts;
+ gmx_large_int_t step_eff;
+ real ecorr,pcorr,dvdlcorr;
+ real bmass,qmass,reft,kT,dt,nd;
+ tensor dumpres,dumvir;
+ double *scalefac,dtc;
+ int *trotter_seq;
+ rvec sumv={0,0,0},consk;
+ gmx_bool bCouple;
+
+ if (trotter_seqno <= ettTSEQ2)
+ {
+ step_eff = step-1; /* the velocity verlet calls are actually out of order -- the first half step
+ is actually the last half step from the previous step. Thus the first half step
+ actually corresponds to the n-1 step*/
+
+ } else {
+ step_eff = step;
+ }
+
+ bCouple = (ir->nsttcouple == 1 ||
+ do_per_step(step_eff+ir->nsttcouple,ir->nsttcouple));
+
+ trotter_seq = trotter_seqlist[trotter_seqno];
+
+ /* signal we are returning if nothing is going to be done in this routine */
+ if ((trotter_seq[0] == etrtSKIPALL) || !(bCouple))
+ {
+ return;
+ }
+
+ dtc = ir->nsttcouple*ir->delta_t;
+ opts = &(ir->opts); /* just for ease of referencing */
+ ngtc = opts->ngtc;
++ assert(ngtc>0);
+ snew(scalefac,opts->ngtc);
+ for (i=0;i<ngtc;i++)
+ {
+ scalefac[i] = 1;
+ }
+ /* execute the series of trotter updates specified in the trotterpart array */
+
+ for (i=0;i<NTROTTERPARTS;i++){
+ /* allow for doubled intgrators by doubling dt instead of making 2 calls */
+ if ((trotter_seq[i] == etrtBAROV2) || (trotter_seq[i] == etrtBARONHC2) || (trotter_seq[i] == etrtNHC2))
+ {
+ dt = 2 * dtc;
+ }
+ else
+ {
+ dt = dtc;
+ }
+
+ switch (trotter_seq[i])
+ {
+ case etrtBAROV:
+ case etrtBAROV2:
+ boxv_trotter(ir,&(state->veta),dt,state->box,ekind,vir,
+ enerd->term[F_PDISPCORR],enerd->term[F_DISPCORR],MassQ);
+ break;
+ case etrtBARONHC:
+ case etrtBARONHC2:
+ NHC_trotter(opts,state->nnhpres,ekind,dt,state->nhpres_xi,
+ state->nhpres_vxi,NULL,&(state->veta),MassQ,FALSE);
+ break;
+ case etrtNHC:
+ case etrtNHC2:
+ NHC_trotter(opts,opts->ngtc,ekind,dt,state->nosehoover_xi,
+ state->nosehoover_vxi,scalefac,NULL,MassQ,(ir->eI==eiVV));
+ /* need to rescale the kinetic energies and velocities here. Could
+ scale the velocities later, but we need them scaled in order to
+ produce the correct outputs, so we'll scale them here. */
+
+ for (i=0; i<ngtc;i++)
+ {
+ tcstat = &ekind->tcstat[i];
+ tcstat->vscale_nhc = scalefac[i];
+ tcstat->ekinscaleh_nhc *= (scalefac[i]*scalefac[i]);
+ tcstat->ekinscalef_nhc *= (scalefac[i]*scalefac[i]);
+ }
+ /* now that we've scaled the groupwise velocities, we can add them up to get the total */
+ /* but do we actually need the total? */
+
+ /* modify the velocities as well */
+ for (n=md->start;n<md->start+md->homenr;n++)
+ {
+ if (md->cTC)
+ {
+ gc = md->cTC[n];
+ }
+ for (d=0;d<DIM;d++)
+ {
+ state->v[n][d] *= scalefac[gc];
+ }
+
+ if (debug)
+ {
+ for (d=0;d<DIM;d++)
+ {
+ sumv[d] += (state->v[n][d])/md->invmass[n];
+ }
+ }
+ }
+ break;
+ default:
+ break;
+ }
+ }
+ /* check for conserved momentum -- worth looking at this again eventually, but not working right now.*/
+#if 0
+ if (debug)
+ {
+ if (bFirstHalf)
+ {
+ for (d=0;d<DIM;d++)
+ {
+ consk[d] = sumv[d]*exp((1 + 1.0/opts->nrdf[0])*((1.0/DIM)*log(det(state->box)/state->vol0)) + state->nosehoover_xi[0]);
+ }
+ fprintf(debug,"Conserved kappa: %15.8f %15.8f %15.8f\n",consk[0],consk[1],consk[2]);
+ }
+ }
+#endif
+ sfree(scalefac);
+}
+
+int **init_npt_vars(t_inputrec *ir, t_state *state, t_extmass *MassQ, gmx_bool bTrotter)
+{
+ int n,i,j,d,ntgrp,ngtc,nnhpres,nh,gc=0;
+ t_grp_tcstat *tcstat;
+ t_grpopts *opts;
+ real ecorr,pcorr,dvdlcorr;
+ real bmass,qmass,reft,kT,dt,ndj,nd;
+ tensor dumpres,dumvir;
+ int **trotter_seq;
+
+ opts = &(ir->opts); /* just for ease of referencing */
+ ngtc = state->ngtc;
+ nnhpres = state->nnhpres;
+ nh = state->nhchainlength;
+
+ if (ir->eI == eiMD) {
+ snew(MassQ->Qinv,ngtc);
+ for(i=0; (i<ngtc); i++)
+ {
+ if ((opts->tau_t[i] > 0) && (opts->ref_t[i] > 0))
+ {
+ MassQ->Qinv[i]=1.0/(sqr(opts->tau_t[i]/M_2PI)*opts->ref_t[i]);
+ }
+ else
+ {
+ MassQ->Qinv[i]=0.0;
+ }
+ }
+ }
+ else if (EI_VV(ir->eI))
+ {
+ /* Set pressure variables */
+
+ if (state->vol0 == 0)
+ {
+ state->vol0 = det(state->box); /* because we start by defining a fixed compressibility,
+ we need the volume at this compressibility to solve the problem */
+ }
+
+ /* units are nm^3 * ns^2 / (nm^3 * bar / kJ/mol) = kJ/mol */
+ /* Investigate this more -- is this the right mass to make it? */
+ MassQ->Winv = (PRESFAC*trace(ir->compress)*BOLTZ*opts->ref_t[0])/(DIM*state->vol0*sqr(ir->tau_p/M_2PI));
+ /* An alternate mass definition, from Tuckerman et al. */
+ /* MassQ->Winv = 1.0/(sqr(ir->tau_p/M_2PI)*(opts->nrdf[0]+DIM)*BOLTZ*opts->ref_t[0]); */
+ for (d=0;d<DIM;d++)
+ {
+ for (n=0;n<DIM;n++)
+ {
+ MassQ->Winvm[d][n]= PRESFAC*ir->compress[d][n]/(state->vol0*sqr(ir->tau_p/M_2PI));
+ /* not clear this is correct yet for the anisotropic case*/
+ }
+ }
+ /* Allocate space for thermostat variables */
+ snew(MassQ->Qinv,ngtc*nh);
+
+ /* now, set temperature variables */
+ for(i=0; i<ngtc; i++)
+ {
+ if ((opts->tau_t[i] > 0) && (opts->ref_t[i] > 0))
+ {
+ reft = max(0.0,opts->ref_t[i]);
+ nd = opts->nrdf[i];
+ kT = BOLTZ*reft;
+ for (j=0;j<nh;j++)
+ {
+ if (j==0)
+ {
+ ndj = nd;
+ }
+ else
+ {
+ ndj = 1;
+ }
+ MassQ->Qinv[i*nh+j] = 1.0/(sqr(opts->tau_t[i]/M_2PI)*ndj*kT);
+ }
+ }
+ else
+ {
+ reft=0.0;
+ for (j=0;j<nh;j++)
+ {
+ MassQ->Qinv[i*nh+j] = 0.0;
+ }
+ }
+ }
+ }
+
+ /* first, initialize clear all the trotter calls */
+ snew(trotter_seq,ettTSEQMAX);
+ for (i=0;i<ettTSEQMAX;i++)
+ {
+ snew(trotter_seq[i],NTROTTERPARTS);
+ for (j=0;j<NTROTTERPARTS;j++) {
+ trotter_seq[i][j] = etrtNONE;
+ }
+ trotter_seq[i][0] = etrtSKIPALL;
+ }
+
+ if (!bTrotter)
+ {
+ /* no trotter calls, so we never use the values in the array.
+ * We access them (so we need to define them, but ignore
+ * then.*/
+
+ return trotter_seq;
+ }
+
+ /* compute the kinetic energy by using the half step velocities or
+ * the kinetic energies, depending on the order of the trotter calls */
+
+ if (ir->eI==eiVV)
+ {
+ if (IR_NPT_TROTTER(ir))
+ {
+ /* This is the complicated version - there are 4 possible calls, depending on ordering.
+ We start with the initial one. */
+ /* first, a round that estimates veta. */
+ trotter_seq[0][0] = etrtBAROV;
+
+ /* trotter_seq[1] is etrtNHC for 1/2 step velocities - leave zero */
+
+ /* The first half trotter update */
+ trotter_seq[2][0] = etrtBAROV;
+ trotter_seq[2][1] = etrtNHC;
+ trotter_seq[2][2] = etrtBARONHC;
+
+ /* The second half trotter update */
+ trotter_seq[3][0] = etrtBARONHC;
+ trotter_seq[3][1] = etrtNHC;
+ trotter_seq[3][2] = etrtBAROV;
+
+ /* trotter_seq[4] is etrtNHC for second 1/2 step velocities - leave zero */
+
+ }
+ else
+ {
+ if (IR_NVT_TROTTER(ir))
+ {
+ /* This is the easy version - there are only two calls, both the same.
+ Otherwise, even easier -- no calls */
+ trotter_seq[2][0] = etrtNHC;
+ trotter_seq[3][0] = etrtNHC;
+ }
+ }
+ } else if (ir->eI==eiVVAK) {
+ if (IR_NPT_TROTTER(ir))
+ {
+ /* This is the complicated version - there are 4 possible calls, depending on ordering.
+ We start with the initial one. */
+ /* first, a round that estimates veta. */
+ trotter_seq[0][0] = etrtBAROV;
+
+ /* The first half trotter update, part 1 -- double update, because it commutes */
+ trotter_seq[1][0] = etrtNHC;
+
+ /* The first half trotter update, part 2 */
+ trotter_seq[2][0] = etrtBAROV;
+ trotter_seq[2][1] = etrtBARONHC;
+
+ /* The second half trotter update, part 1 */
+ trotter_seq[3][0] = etrtBARONHC;
+ trotter_seq[3][1] = etrtBAROV;
+
+ /* The second half trotter update -- blank for now */
+ trotter_seq[4][0] = etrtNHC;
+ }
+ else
+ {
+ if (IR_NVT_TROTTER(ir))
+ {
+ /* This is the easy version - there is only one call, both the same.
+ Otherwise, even easier -- no calls */
+ trotter_seq[1][0] = etrtNHC;
+ trotter_seq[4][0] = etrtNHC;
+ }
+ }
+ }
+
+ switch (ir->epct)
+ {
+ case epctISOTROPIC:
+ default:
+ bmass = DIM*DIM; /* recommended mass parameters for isotropic barostat */
+ }
+
+ snew(MassQ->QPinv,nnhpres*opts->nhchainlength);
+
+ /* barostat temperature */
+ if ((ir->tau_p > 0) && (opts->ref_t[0] > 0))
+ {
+ reft = max(0.0,opts->ref_t[0]);
+ kT = BOLTZ*reft;
+ for (i=0;i<nnhpres;i++) {
+ for (j=0;j<nh;j++)
+ {
+ if (j==0) {
+ qmass = bmass;
+ }
+ else
+ {
+ qmass = 1;
+ }
+ MassQ->QPinv[i*opts->nhchainlength+j] = 1.0/(sqr(opts->tau_t[0]/M_2PI)*qmass*kT);
+ }
+ }
+ }
+ else
+ {
+ for (i=0;i<nnhpres;i++) {
+ for (j=0;j<nh;j++)
+ {
+ MassQ->QPinv[i*nh+j]=0.0;
+ }
+ }
+ }
+ return trotter_seq;
+}
+
+real NPT_energy(t_inputrec *ir, t_state *state, t_extmass *MassQ)
+{
+ int i,j,nd,ndj,bmass,qmass,ngtcall;
+ real ener_npt,reft,eta,kT,tau;
+ double *ivxi, *ixi;
+ double *iQinv;
+ real vol,dbaro,W,Q;
+ int nh = state->nhchainlength;
+
+ ener_npt = 0;
+
+ /* now we compute the contribution of the pressure to the conserved quantity*/
+
+ if (ir->epc==epcMTTK)
+ {
+ /* find the volume, and the kinetic energy of the volume */
+
+ switch (ir->epct) {
+
+ case epctISOTROPIC:
+ /* contribution from the pressure momenenta */
+ ener_npt += 0.5*sqr(state->veta)/MassQ->Winv;
+
+ /* contribution from the PV term */
+ vol = det(state->box);
+ ener_npt += vol*trace(ir->ref_p)/(DIM*PRESFAC);
+
+ break;
+ case epctANISOTROPIC:
+
+ break;
+
+ case epctSURFACETENSION:
+
+ break;
+ case epctSEMIISOTROPIC:
+
+ break;
+ default:
+ break;
+ }
+ }
+
+ if (IR_NPT_TROTTER(ir))
+ {
+ /* add the energy from the barostat thermostat chain */
+ for (i=0;i<state->nnhpres;i++) {
+
+ /* note -- assumes only one degree of freedom that is thermostatted in barostat */
+ ivxi = &state->nhpres_vxi[i*nh];
+ ixi = &state->nhpres_xi[i*nh];
+ iQinv = &(MassQ->QPinv[i*nh]);
+ reft = max(ir->opts.ref_t[0],0); /* using 'System' temperature */
+ kT = BOLTZ * reft;
+
+ for (j=0;j<nh;j++)
+ {
+ if (iQinv[j] > 0)
+ {
+ ener_npt += 0.5*sqr(ivxi[j])/iQinv[j];
+ /* contribution from the thermal variable of the NH chain */
+ ener_npt += ixi[j]*kT;
+ }
+ if (debug)
+ {
+ fprintf(debug,"P-T-group: %10d Chain %4d ThermV: %15.8f ThermX: %15.8f",i,j,ivxi[j],ixi[j]);
+ }
+ }
+ }
+ }
+
+ if (ir->etc)
+ {
+ for(i=0; i<ir->opts.ngtc; i++)
+ {
+ ixi = &state->nosehoover_xi[i*nh];
+ ivxi = &state->nosehoover_vxi[i*nh];
+ iQinv = &(MassQ->Qinv[i*nh]);
+
+ nd = ir->opts.nrdf[i];
+ reft = max(ir->opts.ref_t[i],0);
+ kT = BOLTZ * reft;
+
+ if (nd > 0)
+ {
+ if (IR_NVT_TROTTER(ir))
+ {
+ /* contribution from the thermal momenta of the NH chain */
+ for (j=0;j<nh;j++)
+ {
+ if (iQinv[j] > 0) {
+ ener_npt += 0.5*sqr(ivxi[j])/iQinv[j];
+ /* contribution from the thermal variable of the NH chain */
+ if (j==0) {
+ ndj = nd;
+ }
+ else
+ {
+ ndj = 1;
+ }
+ ener_npt += ndj*ixi[j]*kT;
+ }
+ }
+ }
+ else /* Other non Trotter temperature NH control -- no chains yet. */
+ {
+ ener_npt += 0.5*BOLTZ*nd*sqr(ivxi[0])/iQinv[0];
+ ener_npt += nd*ixi[0]*kT;
+ }
+ }
+ }
+ }
+ return ener_npt;
+}
+
+static real vrescale_gamdev(int ia, gmx_rng_t rng)
+/* Gamma distribution, adapted from numerical recipes */
+{
+ int j;
+ real am,e,s,v1,v2,x,y;
+
+ if (ia < 6)
+ {
+ do
+ {
+ x = 1.0;
+ for(j=1; j<=ia; j++)
+ {
+ x *= gmx_rng_uniform_real(rng);
+ }
+ }
+ while (x == 0);
+ x = -log(x);
+ }
+ else
+ {
+ do
+ {
+ do
+ {
+ do
+ {
+ v1 = gmx_rng_uniform_real(rng);
+ v2 = 2.0*gmx_rng_uniform_real(rng)-1.0;
+ }
+ while (v1*v1 + v2*v2 > 1.0 ||
+ v1*v1*GMX_REAL_MAX < 3.0*ia);
+ /* The last check above ensures that both x (3.0 > 2.0 in s)
+ * and the pre-factor for e do not go out of range.
+ */
+ y = v2/v1;
+ am = ia - 1;
+ s = sqrt(2.0*am + 1.0);
+ x = s*y + am;
+ }
+ while (x <= 0.0);
+ e = (1.0 + y*y)*exp(am*log(x/am) - s*y);
+ }
+ while (gmx_rng_uniform_real(rng) > e);
+ }
+
+ return x;
+}
+
+static real vrescale_sumnoises(int nn,gmx_rng_t rng)
+{
+/*
+ * Returns the sum of n independent gaussian noises squared
+ * (i.e. equivalent to summing the square of the return values
+ * of nn calls to gmx_rng_gaussian_real).xs
+ */
+ real rr;
+
+ if (nn == 0) {
+ return 0.0;
+ } else if (nn == 1) {
+ rr = gmx_rng_gaussian_real(rng);
+ return rr*rr;
+ } else if (nn % 2 == 0) {
+ return 2.0*vrescale_gamdev(nn/2,rng);
+ } else {
+ rr = gmx_rng_gaussian_real(rng);
+ return 2.0*vrescale_gamdev((nn-1)/2,rng) + rr*rr;
+ }
+}
+
+static real vrescale_resamplekin(real kk,real sigma, int ndeg, real taut,
+ gmx_rng_t rng)
+{
+/*
+ * Generates a new value for the kinetic energy,
+ * according to Bussi et al JCP (2007), Eq. (A7)
+ * kk: present value of the kinetic energy of the atoms to be thermalized (in arbitrary units)
+ * sigma: target average value of the kinetic energy (ndeg k_b T/2) (in the same units as kk)
+ * ndeg: number of degrees of freedom of the atoms to be thermalized
+ * taut: relaxation time of the thermostat, in units of 'how often this routine is called'
+ */
+ real factor,rr;
+
+ if (taut > 0.1) {
+ factor = exp(-1.0/taut);
+ } else {
+ factor = 0.0;
+ }
+ rr = gmx_rng_gaussian_real(rng);
+ return
+ kk +
+ (1.0 - factor)*(sigma*(vrescale_sumnoises(ndeg-1,rng) + rr*rr)/ndeg - kk) +
+ 2.0*rr*sqrt(kk*sigma/ndeg*(1.0 - factor)*factor);
+}
+
+void vrescale_tcoupl(t_inputrec *ir,gmx_ekindata_t *ekind,real dt,
+ double therm_integral[],gmx_rng_t rng)
+{
+ t_grpopts *opts;
+ int i;
+ real Ek,Ek_ref1,Ek_ref,Ek_new;
+
+ opts = &ir->opts;
+
+ for(i=0; (i<opts->ngtc); i++)
+ {
+ if (ir->eI == eiVV)
+ {
+ Ek = trace(ekind->tcstat[i].ekinf);
+ }
+ else
+ {
+ Ek = trace(ekind->tcstat[i].ekinh);
+ }
+
+ if (opts->tau_t[i] >= 0 && opts->nrdf[i] > 0 && Ek > 0)
+ {
+ Ek_ref1 = 0.5*opts->ref_t[i]*BOLTZ;
+ Ek_ref = Ek_ref1*opts->nrdf[i];
+
+ Ek_new = vrescale_resamplekin(Ek,Ek_ref,opts->nrdf[i],
+ opts->tau_t[i]/dt,rng);
+
+ /* Analytically Ek_new>=0, but we check for rounding errors */
+ if (Ek_new <= 0)
+ {
+ ekind->tcstat[i].lambda = 0.0;
+ }
+ else
+ {
+ ekind->tcstat[i].lambda = sqrt(Ek_new/Ek);
+ }
+
+ therm_integral[i] -= Ek_new - Ek;
+
+ if (debug)
+ {
+ fprintf(debug,"TC: group %d: Ekr %g, Ek %g, Ek_new %g, Lambda: %g\n",
+ i,Ek_ref,Ek,Ek_new,ekind->tcstat[i].lambda);
+ }
+ }
+ else
+ {
+ ekind->tcstat[i].lambda = 1.0;
+ }
+ }
+}
+
+real vrescale_energy(t_grpopts *opts,double therm_integral[])
+{
+ int i;
+ real ener;
+
+ ener = 0;
+ for(i=0; i<opts->ngtc; i++) {
+ ener += therm_integral[i];
+ }
+
+ return ener;
+}
+
+void rescale_velocities(gmx_ekindata_t *ekind,t_mdatoms *mdatoms,
+ int start,int end,rvec v[])
+{
+ t_grp_acc *gstat;
+ t_grp_tcstat *tcstat;
+ unsigned short *cACC,*cTC;
+ int ga,gt,n,d;
+ real lg;
+ rvec vrel;
+
+ tcstat = ekind->tcstat;
+ cTC = mdatoms->cTC;
+
+ if (ekind->bNEMD)
+ {
+ gstat = ekind->grpstat;
+ cACC = mdatoms->cACC;
+
+ ga = 0;
+ gt = 0;
+ for(n=start; n<end; n++)
+ {
+ if (cACC)
+ {
+ ga = cACC[n];
+ }
+ if (cTC)
+ {
+ gt = cTC[n];
+ }
+ /* Only scale the velocity component relative to the COM velocity */
+ rvec_sub(v[n],gstat[ga].u,vrel);
+ lg = tcstat[gt].lambda;
+ for(d=0; d<DIM; d++)
+ {
+ v[n][d] = gstat[ga].u[d] + lg*vrel[d];
+ }
+ }
+ }
+ else
+ {
+ gt = 0;
+ for(n=start; n<end; n++)
+ {
+ if (cTC)
+ {
+ gt = cTC[n];
+ }
+ lg = tcstat[gt].lambda;
+ for(d=0; d<DIM; d++)
+ {
+ v[n][d] *= lg;
+ }
+ }
+ }
+}
+
+
+/* set target temperatures if we are annealing */
+void update_annealing_target_temp(t_grpopts *opts,real t)
+{
+ int i,j,n,npoints;
+ real pert,thist=0,x;
+
+ for(i=0;i<opts->ngtc;i++) {
+ npoints = opts->anneal_npoints[i];
+ switch (opts->annealing[i]) {
+ case eannNO:
+ continue;
+ case eannPERIODIC:
+ /* calculate time modulo the period */
+ pert = opts->anneal_time[i][npoints-1];
+ n = t / pert;
+ thist = t - n*pert; /* modulo time */
+ /* Make sure rounding didn't get us outside the interval */
+ if (fabs(thist-pert) < GMX_REAL_EPS*100)
+ thist=0;
+ break;
+ case eannSINGLE:
+ thist = t;
+ break;
+ default:
+ gmx_fatal(FARGS,"Death horror in update_annealing_target_temp (i=%d/%d npoints=%d)",i,opts->ngtc,npoints);
+ }
+ /* We are doing annealing for this group if we got here,
+ * and we have the (relative) time as thist.
+ * calculate target temp */
+ j=0;
+ while ((j < npoints-1) && (thist>(opts->anneal_time[i][j+1])))
+ j++;
+ if (j < npoints-1) {
+ /* Found our position between points j and j+1.
+ * Interpolate: x is the amount from j+1, (1-x) from point j
+ * First treat possible jumps in temperature as a special case.
+ */
+ if ((opts->anneal_time[i][j+1]-opts->anneal_time[i][j]) < GMX_REAL_EPS*100)
+ opts->ref_t[i]=opts->anneal_temp[i][j+1];
+ else {
+ x = ((thist-opts->anneal_time[i][j])/
+ (opts->anneal_time[i][j+1]-opts->anneal_time[i][j]));
+ opts->ref_t[i] = x*opts->anneal_temp[i][j+1]+(1-x)*opts->anneal_temp[i][j];
+ }
+ }
+ else {
+ opts->ref_t[i] = opts->anneal_temp[i][npoints-1];
+ }
+ }
+}
--- /dev/null
+/* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
+ *
+ *
+ * This source code is part of
+ *
+ * G R O M A C S
+ *
+ * GROningen MAchine for Chemical Simulations
+ *
+ * VERSION 3.2.0
+ * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
+ * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
+ * Copyright (c) 2001-2004, The GROMACS development team,
+ * check out http://www.gromacs.org for more information.
+
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version.
+ *
+ * If you want to redistribute modifications, please consider that
+ * scientific software is very special. Version control is crucial -
+ * bugs must be traceable. We will be happy to consider code for
+ * inclusion in the official distribution, but derived work must not
+ * be called official GROMACS. Details are found in the README & COPYING
+ * files - if they are missing, get the official version at www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the papers on the package - you can find them in the top README file.
+ *
+ * For more info, check our website at http://www.gromacs.org
+ *
+ * And Hey:
+ * GROwing Monsters And Cloning Shrimps
+ */
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <math.h>
+#include <string.h>
++#include <assert.h>
+#include "sysstuff.h"
+#include "typedefs.h"
+#include "macros.h"
+#include "smalloc.h"
+#include "macros.h"
+#include "physics.h"
+#include "force.h"
+#include "nonbonded.h"
+#include "names.h"
+#include "network.h"
+#include "pbc.h"
+#include "ns.h"
+#include "nrnb.h"
+#include "bondf.h"
+#include "mshift.h"
+#include "txtdump.h"
+#include "coulomb.h"
+#include "pme.h"
+#include "mdrun.h"
+#include "domdec.h"
+#include "partdec.h"
+#include "qmmm.h"
+
+
+void ns(FILE *fp,
+ t_forcerec *fr,
+ rvec x[],
+ matrix box,
+ gmx_groups_t *groups,
+ t_grpopts *opts,
+ gmx_localtop_t *top,
+ t_mdatoms *md,
+ t_commrec *cr,
+ t_nrnb *nrnb,
+ real lambda,
+ real *dvdlambda,
+ gmx_grppairener_t *grppener,
+ gmx_bool bFillGrid,
+ gmx_bool bDoLongRange,
+ gmx_bool bDoForces,
+ rvec *f)
+{
+ char *ptr;
+ int nsearch;
+
+
+ if (!fr->ns.nblist_initialized)
+ {
+ init_neighbor_list(fp, fr, md->homenr);
+ }
+
+ if (fr->bTwinRange)
+ fr->nlr=0;
+
+ nsearch = search_neighbours(fp,fr,x,box,top,groups,cr,nrnb,md,
+ lambda,dvdlambda,grppener,
+ bFillGrid,bDoLongRange,
+ bDoForces,f);
+ if (debug)
+ fprintf(debug,"nsearch = %d\n",nsearch);
+
+ /* Check whether we have to do dynamic load balancing */
+ /*if ((nsb->nstDlb > 0) && (mod(step,nsb->nstDlb) == 0))
+ count_nb(cr,nsb,&(top->blocks[ebCGS]),nns,fr->nlr,
+ &(top->idef),opts->ngener);
+ */
+ if (fr->ns.dump_nl > 0)
+ dump_nblist(fp,cr,fr,fr->ns.dump_nl);
+}
+
+void do_force_lowlevel(FILE *fplog, gmx_large_int_t step,
+ t_forcerec *fr, t_inputrec *ir,
+ t_idef *idef, t_commrec *cr,
+ t_nrnb *nrnb, gmx_wallcycle_t wcycle,
+ t_mdatoms *md,
+ t_grpopts *opts,
+ rvec x[], history_t *hist,
+ rvec f[],
+ gmx_enerdata_t *enerd,
+ t_fcdata *fcd,
+ gmx_mtop_t *mtop,
+ gmx_localtop_t *top,
+ gmx_genborn_t *born,
+ t_atomtypes *atype,
+ gmx_bool bBornRadii,
+ matrix box,
+ real lambda,
+ t_graph *graph,
+ t_blocka *excl,
+ rvec mu_tot[],
+ int flags,
+ float *cycles_pme)
+{
+ int i,status;
+ int donb_flags;
+ gmx_bool bDoEpot,bSepDVDL,bSB;
+ int pme_flags;
+ matrix boxs;
+ rvec box_size;
+ real dvdlambda,Vsr,Vlr,Vcorr=0,vdip,vcharge;
+ t_pbc pbc;
+ real dvdgb;
+ char buf[22];
+ gmx_enerdata_t ed_lam;
+ double lam_i;
+ real dvdl_dum;
+
+#ifdef GMX_MPI
+ double t0=0.0,t1,t2,t3; /* time measurement for coarse load balancing */
+#endif
+
+#define PRINT_SEPDVDL(s,v,dvdl) if (bSepDVDL) fprintf(fplog,sepdvdlformat,s,v,dvdl);
+
+
+ set_pbc(&pbc,fr->ePBC,box);
+
+ /* Reset box */
+ for(i=0; (i<DIM); i++)
+ {
+ box_size[i]=box[i][i];
+ }
+
+ bSepDVDL=(fr->bSepDVDL && do_per_step(step,ir->nstlog));
+ debug_gmx();
+
+ /* do QMMM first if requested */
+ if(fr->bQMMM)
+ {
+ enerd->term[F_EQM] = calculate_QMMM(cr,x,f,fr,md);
+ }
+
+ if (bSepDVDL)
+ {
+ fprintf(fplog,"Step %s: non-bonded V and dVdl for node %d:\n",
+ gmx_step_str(step,buf),cr->nodeid);
+ }
+
+ /* Call the short range functions all in one go. */
+
+ dvdlambda = 0;
+
+#ifdef GMX_MPI
+ /*#define TAKETIME ((cr->npmenodes) && (fr->timesteps < 12))*/
+#define TAKETIME FALSE
+ if (TAKETIME)
+ {
+ MPI_Barrier(cr->mpi_comm_mygroup);
+ t0=MPI_Wtime();
+ }
+#endif
+
+ if (ir->nwall)
+ {
+ dvdlambda = do_walls(ir,fr,box,md,x,f,lambda,
+ enerd->grpp.ener[egLJSR],nrnb);
+ PRINT_SEPDVDL("Walls",0.0,dvdlambda);
+ enerd->dvdl_lin += dvdlambda;
+ }
+
+ /* If doing GB, reset dvda and calculate the Born radii */
+ if (ir->implicit_solvent)
+ {
+ /* wallcycle_start(wcycle,ewcGB); */
+
+ for(i=0;i<born->nr;i++)
+ {
+ fr->dvda[i]=0;
+ }
+
+ if(bBornRadii)
+ {
+ calc_gb_rad(cr,fr,ir,top,atype,x,&(fr->gblist),born,md,nrnb);
+ }
+
+ /* wallcycle_stop(wcycle, ewcGB); */
+ }
+
+ where();
+ donb_flags = 0;
+ if (flags & GMX_FORCE_FORCES)
+ {
+ donb_flags |= GMX_DONB_FORCES;
+ }
+ do_nonbonded(cr,fr,x,f,md,excl,
+ fr->bBHAM ?
+ enerd->grpp.ener[egBHAMSR] :
+ enerd->grpp.ener[egLJSR],
+ enerd->grpp.ener[egCOULSR],
+ enerd->grpp.ener[egGB],box_size,nrnb,
+ lambda,&dvdlambda,-1,-1,donb_flags);
+ /* If we do foreign lambda and we have soft-core interactions
+ * we have to recalculate the (non-linear) energies contributions.
+ */
+ if (ir->n_flambda > 0 && (flags & GMX_FORCE_DHDL) && ir->sc_alpha != 0)
+ {
+ init_enerdata(mtop->groups.grps[egcENER].nr,ir->n_flambda,&ed_lam);
+
+ for(i=0; i<enerd->n_lambda; i++)
+ {
+ lam_i = (i==0 ? lambda : ir->flambda[i-1]);
+ dvdl_dum = 0;
+ reset_enerdata(&ir->opts,fr,TRUE,&ed_lam,FALSE);
+ do_nonbonded(cr,fr,x,f,md,excl,
+ fr->bBHAM ?
+ ed_lam.grpp.ener[egBHAMSR] :
+ ed_lam.grpp.ener[egLJSR],
+ ed_lam.grpp.ener[egCOULSR],
+ enerd->grpp.ener[egGB], box_size,nrnb,
+ lam_i,&dvdl_dum,-1,-1,
+ GMX_DONB_FOREIGNLAMBDA);
+ sum_epot(&ir->opts,&ed_lam);
+ enerd->enerpart_lambda[i] += ed_lam.term[F_EPOT];
+ }
+ destroy_enerdata(&ed_lam);
+ }
+ where();
+
+ /* If we are doing GB, calculate bonded forces and apply corrections
+ * to the solvation forces */
+ if (ir->implicit_solvent) {
+ calc_gb_forces(cr,md,born,top,atype,x,f,fr,idef,
+ ir->gb_algorithm,ir->sa_algorithm,nrnb,bBornRadii,&pbc,graph,enerd);
+ }
+
+#ifdef GMX_MPI
+ if (TAKETIME)
+ {
+ t1=MPI_Wtime();
+ fr->t_fnbf += t1-t0;
+ }
+#endif
+
+ if (ir->sc_alpha != 0)
+ {
+ enerd->dvdl_nonlin += dvdlambda;
+ }
+ else
+ {
+ enerd->dvdl_lin += dvdlambda;
+ }
+ Vsr = 0;
+ if (bSepDVDL)
+ {
+ for(i=0; i<enerd->grpp.nener; i++)
+ {
+ Vsr +=
+ (fr->bBHAM ?
+ enerd->grpp.ener[egBHAMSR][i] :
+ enerd->grpp.ener[egLJSR][i])
+ + enerd->grpp.ener[egCOULSR][i] + enerd->grpp.ener[egGB][i];
+ }
+ }
+ PRINT_SEPDVDL("VdW and Coulomb SR particle-p.",Vsr,dvdlambda);
+ debug_gmx();
+
+ if (debug)
+ {
+ pr_rvecs(debug,0,"fshift after SR",fr->fshift,SHIFTS);
+ }
+
+ /* Shift the coordinates. Must be done before bonded forces and PPPM,
+ * but is also necessary for SHAKE and update, therefore it can NOT
+ * go when no bonded forces have to be evaluated.
+ */
+
+ /* Here sometimes we would not need to shift with NBFonly,
+ * but we do so anyhow for consistency of the returned coordinates.
+ */
+ if (graph)
+ {
+ shift_self(graph,box,x);
+ if (TRICLINIC(box))
+ {
+ inc_nrnb(nrnb,eNR_SHIFTX,2*graph->nnodes);
+ }
+ else
+ {
+ inc_nrnb(nrnb,eNR_SHIFTX,graph->nnodes);
+ }
+ }
+ /* Check whether we need to do bondeds or correct for exclusions */
+ if (fr->bMolPBC &&
+ ((flags & GMX_FORCE_BONDED)
+ || EEL_RF(fr->eeltype) || EEL_FULL(fr->eeltype)))
+ {
+ /* Since all atoms are in the rectangular or triclinic unit-cell,
+ * only single box vector shifts (2 in x) are required.
+ */
+ set_pbc_dd(&pbc,fr->ePBC,cr->dd,TRUE,box);
+ }
+ debug_gmx();
+
+ if (flags & GMX_FORCE_BONDED)
+ {
+ calc_bonds(fplog,cr->ms,
+ idef,x,hist,f,fr,&pbc,graph,enerd,nrnb,lambda,md,fcd,
+ DOMAINDECOMP(cr) ? cr->dd->gatindex : NULL, atype, born,
+ fr->bSepDVDL && do_per_step(step,ir->nstlog),step);
+
+ /* Check if we have to determine energy differences
+ * at foreign lambda's.
+ */
+ if (ir->n_flambda > 0 && (flags & GMX_FORCE_DHDL) &&
+ idef->ilsort != ilsortNO_FE)
+ {
+ if (idef->ilsort != ilsortFE_SORTED)
+ {
+ gmx_incons("The bonded interactions are not sorted for free energy");
+ }
+ init_enerdata(mtop->groups.grps[egcENER].nr,ir->n_flambda,&ed_lam);
+
+ for(i=0; i<enerd->n_lambda; i++)
+ {
+ lam_i = (i==0 ? lambda : ir->flambda[i-1]);
+ dvdl_dum = 0;
+ reset_enerdata(&ir->opts,fr,TRUE,&ed_lam,FALSE);
+ calc_bonds_lambda(fplog,
+ idef,x,fr,&pbc,graph,&ed_lam,nrnb,lam_i,md,
+ fcd,
+ DOMAINDECOMP(cr) ? cr->dd->gatindex : NULL);
+ sum_epot(&ir->opts,&ed_lam);
+ enerd->enerpart_lambda[i] += ed_lam.term[F_EPOT];
+ }
+ destroy_enerdata(&ed_lam);
+ }
+ debug_gmx();
+ }
+
+ where();
+
+ *cycles_pme = 0;
+ if (EEL_FULL(fr->eeltype))
+ {
+ bSB = (ir->nwall == 2);
+ if (bSB)
+ {
+ copy_mat(box,boxs);
+ svmul(ir->wall_ewald_zfac,boxs[ZZ],boxs[ZZ]);
+ box_size[ZZ] *= ir->wall_ewald_zfac;
+ }
+
+ clear_mat(fr->vir_el_recip);
+
+ if (fr->bEwald)
+ {
+ if (fr->n_tpi == 0)
+ {
+ dvdlambda = 0;
+ Vcorr = ewald_LRcorrection(fplog,md->start,md->start+md->homenr,
+ cr,fr,
+ md->chargeA,
+ md->nChargePerturbed ? md->chargeB : NULL,
+ excl,x,bSB ? boxs : box,mu_tot,
+ ir->ewald_geometry,
+ ir->epsilon_surface,
+ lambda,&dvdlambda,&vdip,&vcharge);
+ PRINT_SEPDVDL("Ewald excl./charge/dip. corr.",Vcorr,dvdlambda);
+ enerd->dvdl_lin += dvdlambda;
+ }
+ else
+ {
+ if (ir->ewald_geometry != eewg3D || ir->epsilon_surface != 0)
+ {
+ gmx_fatal(FARGS,"TPI with PME currently only works in a 3D geometry with tin-foil boundary conditions");
+ }
+ /* The TPI molecule does not have exclusions with the rest
+ * of the system and no intra-molecular PME grid contributions
+ * will be calculated in gmx_pme_calc_energy.
+ */
+ Vcorr = 0;
+ }
+ }
+
+ dvdlambda = 0;
+ status = 0;
+ switch (fr->eeltype)
+ {
+ case eelPME:
+ case eelPMESWITCH:
+ case eelPMEUSER:
+ case eelPMEUSERSWITCH:
+ case eelP3M_AD:
+ if (cr->duty & DUTY_PME)
+ {
++ assert(fr->n_tpi >= 0);
+ if (fr->n_tpi == 0 || (flags & GMX_FORCE_STATECHANGED))
+ {
+ pme_flags = GMX_PME_SPREAD_Q | GMX_PME_SOLVE;
+ if (flags & GMX_FORCE_FORCES)
+ {
+ pme_flags |= GMX_PME_CALC_F;
+ }
+ if (flags & GMX_FORCE_VIRIAL)
+ {
+ pme_flags |= GMX_PME_CALC_ENER_VIR;
+ }
+ if (fr->n_tpi > 0)
+ {
+ /* We don't calculate f, but we do want the potential */
+ pme_flags |= GMX_PME_CALC_POT;
+ }
+ wallcycle_start(wcycle,ewcPMEMESH);
+ status = gmx_pme_do(fr->pmedata,
+ md->start,md->homenr - fr->n_tpi,
+ x,fr->f_novirsum,
+ md->chargeA,md->chargeB,
+ bSB ? boxs : box,cr,
+ DOMAINDECOMP(cr) ? dd_pme_maxshift_x(cr->dd) : 0,
+ DOMAINDECOMP(cr) ? dd_pme_maxshift_y(cr->dd) : 0,
+ nrnb,wcycle,
+ fr->vir_el_recip,fr->ewaldcoeff,
+ &Vlr,lambda,&dvdlambda,
+ pme_flags);
+ *cycles_pme = wallcycle_stop(wcycle,ewcPMEMESH);
+
+ /* We should try to do as little computation after
+ * this as possible, because parallel PME synchronizes
+ * the nodes, so we want all load imbalance of the rest
+ * of the force calculation to be before the PME call.
+ * DD load balancing is done on the whole time of
+ * the force call (without PME).
+ */
+ }
+ if (fr->n_tpi > 0)
+ {
+ /* Determine the PME grid energy of the test molecule
+ * with the PME grid potential of the other charges.
+ */
+ gmx_pme_calc_energy(fr->pmedata,fr->n_tpi,
+ x + md->homenr - fr->n_tpi,
+ md->chargeA + md->homenr - fr->n_tpi,
+ &Vlr);
+ }
+ PRINT_SEPDVDL("PME mesh",Vlr,dvdlambda);
+ }
+ else
+ {
+ /* Energies and virial are obtained later from the PME nodes */
+ /* but values have to be zeroed out here */
+ Vlr=0.0;
+ }
+ break;
+ case eelEWALD:
+ Vlr = do_ewald(fplog,FALSE,ir,x,fr->f_novirsum,
+ md->chargeA,md->chargeB,
+ box_size,cr,md->homenr,
+ fr->vir_el_recip,fr->ewaldcoeff,
+ lambda,&dvdlambda,fr->ewald_table);
+ PRINT_SEPDVDL("Ewald long-range",Vlr,dvdlambda);
+ break;
+ default:
+ Vlr = 0;
+ gmx_fatal(FARGS,"No such electrostatics method implemented %s",
+ eel_names[fr->eeltype]);
+ }
+ if (status != 0)
+ {
+ gmx_fatal(FARGS,"Error %d in long range electrostatics routine %s",
+ status,EELTYPE(fr->eeltype));
+ }
+ enerd->dvdl_lin += dvdlambda;
+ enerd->term[F_COUL_RECIP] = Vlr + Vcorr;
+ if (debug)
+ {
+ fprintf(debug,"Vlr = %g, Vcorr = %g, Vlr_corr = %g\n",
+ Vlr,Vcorr,enerd->term[F_COUL_RECIP]);
+ pr_rvecs(debug,0,"vir_el_recip after corr",fr->vir_el_recip,DIM);
+ pr_rvecs(debug,0,"fshift after LR Corrections",fr->fshift,SHIFTS);
+ }
+ }
+ else
+ {
+ if (EEL_RF(fr->eeltype))
+ {
+ dvdlambda = 0;
+
+ if (fr->eeltype != eelRF_NEC)
+ {
+ enerd->term[F_RF_EXCL] =
+ RF_excl_correction(fplog,fr,graph,md,excl,x,f,
+ fr->fshift,&pbc,lambda,&dvdlambda);
+ }
+
+ enerd->dvdl_lin += dvdlambda;
+ PRINT_SEPDVDL("RF exclusion correction",
+ enerd->term[F_RF_EXCL],dvdlambda);
+ }
+ }
+ where();
+ debug_gmx();
+
+ if (debug)
+ {
+ print_nrnb(debug,nrnb);
+ }
+ debug_gmx();
+
+#ifdef GMX_MPI
+ if (TAKETIME)
+ {
+ t2=MPI_Wtime();
+ MPI_Barrier(cr->mpi_comm_mygroup);
+ t3=MPI_Wtime();
+ fr->t_wait += t3-t2;
+ if (fr->timesteps == 11)
+ {
+ fprintf(stderr,"* PP load balancing info: node %d, step %s, rel wait time=%3.0f%% , load string value: %7.2f\n",
+ cr->nodeid, gmx_step_str(fr->timesteps,buf),
+ 100*fr->t_wait/(fr->t_wait+fr->t_fnbf),
+ (fr->t_fnbf+fr->t_wait)/fr->t_fnbf);
+ }
+ fr->timesteps++;
+ }
+#endif
+
+ if (debug)
+ {
+ pr_rvecs(debug,0,"fshift after bondeds",fr->fshift,SHIFTS);
+ }
+
+}
+
+void init_enerdata(int ngener,int n_flambda,gmx_enerdata_t *enerd)
+{
+ int i,n2;
+
+ for(i=0; i<F_NRE; i++)
+ {
+ enerd->term[i] = 0;
+ }
+
+ n2=ngener*ngener;
+ if (debug)
+ {
+ fprintf(debug,"Creating %d sized group matrix for energies\n",n2);
+ }
+ enerd->grpp.nener = n2;
+ for(i=0; (i<egNR); i++)
+ {
+ snew(enerd->grpp.ener[i],n2);
+ }
+
+ if (n_flambda)
+ {
+ enerd->n_lambda = 1 + n_flambda;
+ snew(enerd->enerpart_lambda,enerd->n_lambda);
+ }
+ else
+ {
+ enerd->n_lambda = 0;
+ }
+}
+
+void destroy_enerdata(gmx_enerdata_t *enerd)
+{
+ int i;
+
+ for(i=0; (i<egNR); i++)
+ {
+ sfree(enerd->grpp.ener[i]);
+ }
+
+ if (enerd->n_lambda)
+ {
+ sfree(enerd->enerpart_lambda);
+ }
+}
+
+static real sum_v(int n,real v[])
+{
+ real t;
+ int i;
+
+ t = 0.0;
+ for(i=0; (i<n); i++)
+ t = t + v[i];
+
+ return t;
+}
+
+void sum_epot(t_grpopts *opts,gmx_enerdata_t *enerd)
+{
+ gmx_grppairener_t *grpp;
+ real *epot;
+ int i;
+
+ grpp = &enerd->grpp;
+ epot = enerd->term;
+
+ /* Accumulate energies */
+ epot[F_COUL_SR] = sum_v(grpp->nener,grpp->ener[egCOULSR]);
+ epot[F_LJ] = sum_v(grpp->nener,grpp->ener[egLJSR]);
+ epot[F_LJ14] = sum_v(grpp->nener,grpp->ener[egLJ14]);
+ epot[F_COUL14] = sum_v(grpp->nener,grpp->ener[egCOUL14]);
+ epot[F_COUL_LR] = sum_v(grpp->nener,grpp->ener[egCOULLR]);
+ epot[F_LJ_LR] = sum_v(grpp->nener,grpp->ener[egLJLR]);
+ /* We have already added 1-2,1-3, and 1-4 terms to F_GBPOL */
+ epot[F_GBPOL] += sum_v(grpp->nener,grpp->ener[egGB]);
+
+/* lattice part of LR doesnt belong to any group
+ * and has been added earlier
+ */
+ epot[F_BHAM] = sum_v(grpp->nener,grpp->ener[egBHAMSR]);
+ epot[F_BHAM_LR] = sum_v(grpp->nener,grpp->ener[egBHAMLR]);
+
+ epot[F_EPOT] = 0;
+ for(i=0; (i<F_EPOT); i++)
+ if (i != F_DISRESVIOL && i != F_ORIRESDEV && i != F_DIHRESVIOL)
+ epot[F_EPOT] += epot[i];
+}
+
+void sum_dhdl(gmx_enerdata_t *enerd,double lambda,t_inputrec *ir)
+{
+ int i;
+ double dlam,dhdl_lin;
+
+ enerd->term[F_DVDL] = enerd->dvdl_lin + enerd->dvdl_nonlin;
+
+ if (debug)
+ {
+ fprintf(debug,"dvdl: %f, non-linear %f + linear %f\n",
+ enerd->term[F_DVDL],enerd->dvdl_nonlin,enerd->dvdl_lin);
+ }
+
+ /* Notes on the foreign lambda free energy difference evaluation:
+ * Adding the potential and ekin terms that depend linearly on lambda
+ * as delta lam * dvdl to the energy differences is exact.
+ * For the constraint dvdl this is not exact, but we have no other option.
+ * For the non-bonded LR term we assume that the soft-core (if present)
+ * no longer affects the energy beyond the short-range cut-off,
+ * which is a very good approximation (except for exotic settings).
+ */
+ for(i=1; i<enerd->n_lambda; i++)
+ {
+ dlam = (ir->flambda[i-1] - lambda);
+ dhdl_lin =
+ enerd->dvdl_lin + enerd->term[F_DKDL] + enerd->term[F_DHDL_CON];
+ if (debug)
+ {
+ fprintf(debug,"enerdiff lam %g: non-linear %f linear %f*%f\n",
+ ir->flambda[i-1],
+ enerd->enerpart_lambda[i] - enerd->enerpart_lambda[0],
+ dlam,dhdl_lin);
+ }
+ enerd->enerpart_lambda[i] += dlam*dhdl_lin;
+
+ }
+}
+
+void reset_enerdata(t_grpopts *opts,
+ t_forcerec *fr,gmx_bool bNS,
+ gmx_enerdata_t *enerd,
+ gmx_bool bMaster)
+{
+ gmx_bool bKeepLR;
+ int i,j;
+
+ /* First reset all energy components, except for the long range terms
+ * on the master at non neighbor search steps, since the long range
+ * terms have already been summed at the last neighbor search step.
+ */
+ bKeepLR = (fr->bTwinRange && !bNS);
+ for(i=0; (i<egNR); i++) {
+ if (!(bKeepLR && bMaster && (i == egCOULLR || i == egLJLR))) {
+ for(j=0; (j<enerd->grpp.nener); j++)
+ enerd->grpp.ener[i][j] = 0.0;
+ }
+ }
+ enerd->dvdl_lin = 0.0;
+ enerd->dvdl_nonlin = 0.0;
+
+ /* Normal potential energy components */
+ for(i=0; (i<=F_EPOT); i++) {
+ enerd->term[i] = 0.0;
+ }
+ /* Initialize the dVdlambda term with the long range contribution */
+ enerd->term[F_DVDL] = 0.0;
+ enerd->term[F_DKDL] = 0.0;
+ enerd->term[F_DHDL_CON] = 0.0;
+ if (enerd->n_lambda > 0)
+ {
+ for(i=0; i<enerd->n_lambda; i++)
+ {
+ enerd->enerpart_lambda[i] = 0.0;
+ }
+ }
+}
--- /dev/null
- int repl_ex_nst,int repl_ex_seed,gmx_membed_t *membed,
+/* -*- 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 "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 "pme.h"
+#include "mdatoms.h"
+#include "repl_ex.h"
+#include "qmmm.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 "string2.h"
+#include "membed.h"
+
+#ifdef GMX_LIB_MPI
+#include <mpi.h>
+#endif
+#ifdef GMX_THREAD_MPI
+#include "tmpi.h"
+#endif
+
+#ifdef GMX_FAHCORE
+#include "corewrap.h"
+#endif
+
+
+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 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;
+ gmx_large_int_t multisim_nsteps=-1; /* number of steps to do before first multisim
+ simulation stops. If equal to zero, don't
+ communicate any more between multisims.*/
+#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_THREAD_MPI
+ 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_THREAD_MPI
+ 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
+ | ((ir->comm_mode != ecmNO) ? CGLO_STOPCM:0)
+ | (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,
+ NULL,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,
+ NULL,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_STOPCM | 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 */
+ }
+
+ /* 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));
+ }
+ if (EI_STATE_VELOCITY(ir->eI))
+ {
+ 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);
+ }
+
+ if (MULTISIM(cr) && (repl_ex_nst <=0 ))
+ {
+ /* check how many steps are left in other sims */
+ multisim_nsteps=get_multisim_nsteps(cr, ir->nsteps);
+ }
+
+
+ /* and stop now if we should */
+ bLastStep = (bRerunMD || (ir->nsteps >= 0 && step_rel > ir->nsteps) ||
+ ((multisim_nsteps >= 0) && (step_rel >= multisim_nsteps )));
+ while (!bLastStep || (bRerunMD && bNotLastFrame)) {
+
+ wallcycle_start(wcycle,ewcSTEP);
+
+ 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);
+ }
+ }
+
+ /* check whether we should stop because another simulation has
+ stopped. */
+ if (MULTISIM(cr))
+ {
+ if ( (multisim_nsteps >= 0) && (step_rel >= multisim_nsteps) &&
+ (multisim_nsteps != ir->nsteps) )
+ {
+ if (bNS)
+ {
+ if (MASTER(cr))
+ {
+ fprintf(stderr,
+ "Stopping simulation %d because another one has finished\n",
+ cr->ms->sim);
+ }
+ bLastStep=TRUE;
+ gs.sig[eglsCHKPT] = 1;
+ }
+ }
+ }
+
+ /* < 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);
+ }
+ }
+
+ /* 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) |
+ (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);
+ }
+
+ 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
+ | (bStopCM ? CGLO_STOPCM : 0)
+ | (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;
+ }
+
+ /* 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
+ */
+ 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);
+ }
+
+ /* 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_THREAD_MPI
+ && 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 ################# */
+ /* 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);
+ }
+
+ 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_rel % gs.nstms == 0) &&
+ (multisim_nsteps<0 || (step_rel<multisim_nsteps)),
+ lastbox,
+ top_global,&pcurr,top_global->natoms,&bSumEkinhOld,
+ cglo_flags
+ | (!EI_VV(ir->eI) ? CGLO_ENERGY : 0)
+ | (!EI_VV(ir->eI) && bStopCM ? CGLO_STOPCM : 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 disk space?");
+ }
+ }
+ }
+
+
+ /* 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;
+}
--- /dev/null
- gmx_membed_t *membed,
+/* -*- 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
+ *
+ * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
+ * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
+ * Copyright (c) 2001-2010, 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>
+
+#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 "pme.h"
+#include "mdatoms.h"
+#include "qmmm.h"
+#include "domdec.h"
+#include "partdec.h"
+#include "topsort.h"
+#include "coulomb.h"
+#include "constr.h"
+#include "compute_io.h"
+#include "mvdata.h"
+#include "checkpoint.h"
+#include "mtop_util.h"
+#include "sighandler.h"
+#include "genborn.h"
+#include "string2.h"
+#include "copyrite.h"
+#include "membed.h"
+
+#ifdef GMX_THREAD_MPI
+#include "tmpi.h"
+#endif
+
+/* include even when OpenMM not used to force compilation of do_md_openmm */
+#include "openmm_wrapper.h"
+
+double do_md_openmm(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 bSimAnn,
+ bFirstStep,bStateFromTPX,bLastStep,bStartingFromCpt;
+ gmx_bool bInitStep=TRUE;
+ gmx_bool do_ene,do_log, do_verbose,
+ bX,bV,bF,bCPT;
+ tensor force_vir,shake_vir,total_vir,pres;
+ int i,m;
+ int mdof_flags;
+ rvec mu_tot;
+ t_vcm *vcm;
+ int nchkpt=1;
+ gmx_localtop_t *top;
+ t_mdebin *mdebin;
+ 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_groups_t *groups;
+ gmx_ekindata_t *ekind, *ekind_save;
+ gmx_bool bAppend;
+ int a0,a1;
+ matrix lastbox;
+ real reset_counters=0,reset_counters_now=0;
+ char sbuf[STEPSTRSIZE],sbuf2[STEPSTRSIZE];
+ int handled_stop_condition=gmx_stop_cond_none;
+
+ const char *ommOptions = NULL;
+ void *openmmData;
+
+#ifdef GMX_DOUBLE
+ /* Checks in cmake should prevent the compilation in double precision
+ * with OpenMM, but just to be sure we check here.
+ */
+ gmx_fatal(FARGS,"Compilation was performed in double precision, but OpenMM only supports single precision. If you want to use to OpenMM, compile in single precision.");
+#endif
+
+ bAppend = (Flags & MD_APPENDFILES);
+ check_ir_old_tpx_versions(cr,fplog,ir,top_global);
+
+ 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);
+ 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();
+
+ {
+ 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);
+ }
+
+ 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);
+ }
+
+ update_mdatoms(mdatoms,state->lambda);
+
+ if (deviceOptions[0]=='\0')
+ {
+ /* empty options, which should default to OpenMM in this build */
+ ommOptions=deviceOptions;
+ }
+ else
+ {
+ if (gmx_strncasecmp(deviceOptions,"OpenMM",6)!=0)
+ {
+ gmx_fatal(FARGS, "This Gromacs version currently only works with OpenMM. Use -device \"OpenMM:<options>\"");
+ }
+ else
+ {
+ ommOptions=strchr(deviceOptions,':');
+ if (NULL!=ommOptions)
+ {
+ /* Increase the pointer to skip the colon */
+ ommOptions++;
+ }
+ }
+ }
+
+ openmmData = openmm_init(fplog, ommOptions, ir, top_global, top, mdatoms, fr, state);
+ please_cite(fplog,"Friedrichs2009");
+
+ if (MASTER(cr))
+ {
+ /* Update mdebin with energy history if appending to output files */
+ if ( Flags & MD_APPENDFILES )
+ {
+ restore_energyhistory_from_state(mdebin,&state_global->enerhist);
+ }
+ /* Set the initial energy history in state to zero by updating once */
+ update_energyhistory(&state_global->enerhist,mdebin);
+ }
+
+ if (constr)
+ {
+ set_constraints(constr,top,ir,mdatoms,cr);
+ }
+
+ if (!ir->bContinuation)
+ {
+ 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();
+
+ if (MASTER(cr))
+ {
+ 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 */
+
+ debug_gmx();
+ /***********************************************************
+ *
+ * Loop over MD steps
+ *
+ ************************************************************/
+
+ /* 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;
+ bStartingFromCpt = (Flags & MD_STARTFROMCPT) && bInitStep;
+ bLastStep = FALSE;
+
+ init_global_signals(&gs,cr,ir,repl_ex_nst);
+
+ step = ir->init_step;
+ step_rel = 0;
+
+ while (!bLastStep)
+ {
+ wallcycle_start(wcycle,ewcSTEP);
+
+ bLastStep = (step_rel == ir->nsteps);
+ t = t0 + step*ir->delta_t;
+
+ if (gs.set[eglsSTOPCOND] != 0)
+ {
+ bLastStep = TRUE;
+ }
+
+ do_log = do_per_step(step,ir->nstlog) || bFirstStep || bLastStep;
+ do_verbose = bVerbose &&
+ (step % stepout == 0 || bFirstStep || bLastStep);
+
+ if (MASTER(cr) && do_log)
+ {
+ print_ebin_header(fplog,step,t,state->lambda);
+ }
+
+ clear_mat(force_vir);
+
+ /* We write a checkpoint at this MD step when:
+ * either when we signalled through gs (in OpenMM NS works different),
+ * or at the last step (but not when we do not want confout),
+ * but never at the first step.
+ */
+ bCPT = ((gs.set[eglsCHKPT] ||
+ (bLastStep && (Flags & MD_CONFOUT))) &&
+ step > ir->init_step );
+ if (bCPT)
+ {
+ gs.set[eglsCHKPT] = 0;
+ }
+
+ /* 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
+ */
+ 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;
+ };
+ do_ene = (do_per_step(step,ir->nstenergy) || bLastStep);
+
+ if (mdof_flags != 0 || do_ene || do_log)
+ {
+ wallcycle_start(wcycle,ewcTRAJ);
+ bF = (mdof_flags & MDOF_F);
+ bX = (mdof_flags & (MDOF_X | MDOF_XTC | MDOF_CPT));
+ bV = (mdof_flags & (MDOF_V | MDOF_CPT));
+
+ openmm_copy_state(openmmData, state, &t, f, enerd, bX, bV, bF, do_ene);
+
+ upd_mdebin(mdebin,FALSE,TRUE,
+ t,mdatoms->tmass,enerd,state,lastbox,
+ shake_vir,force_vir,total_vir,pres,
+ ekind,mu_tot,constr);
+ print_ebin(outf->fp_ene,do_ene,FALSE,FALSE,do_log?fplog:NULL,
+ step,t,
+ eprNORMAL,bCompact,mdebin,fcd,groups,&(ir->opts));
+ 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))
+ {
+ /* 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)
+ {
+ /* 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);
+ }
+
+ /* 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_THREAD_MPI
+ && 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 */
+ /* NOTE: this only works for serial code. For code that allows
+ MPI nodes to propagate their condition, see kernel/md.c*/
+ if ( gmx_get_stop_condition() == gmx_stop_cond_next_ns )
+ gs.set[eglsSTOPCOND]=1;
+ if ( gmx_get_stop_condition() == gmx_stop_cond_next )
+ gs.set[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) &&
+ (max_hours > 0 && run_time > max_hours*60.0*60.0*0.99) &&
+ gs.set[eglsSTOPCOND] == 0)
+ {
+ /* Signal to terminate the run */
+ gs.set[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);
+ }
+
+ /* checkpoints */
+ if (MASTER(cr) && (cpt_period >= 0 &&
+ (cpt_period == 0 ||
+ run_time >= nchkpt*cpt_period*60.0)) &&
+ gs.set[eglsCHKPT] == 0)
+ {
+ gs.set[eglsCHKPT] = 1;
+ }
+
+ /* Time for performance */
+ if (((step % stepout) == 0) || bLastStep)
+ {
+ runtime_upd_proc(runtime);
+ }
+
+ if (do_per_step(step,ir->nstlog))
+ {
+ if (fflush(fplog) != 0)
+ {
+ gmx_fatal(FARGS,"Cannot flush logfile - maybe you are out of disk space?");
+ }
+ }
+
+ /* Remaining runtime */
+ if (MULTIMASTER(cr) && (do_verbose || gmx_got_usr_signal() ))
+ {
+ print_time(stderr,runtime,step,ir,cr);
+ }
+
+ bFirstStep = FALSE;
+ bInitStep = FALSE;
+ bStartingFromCpt = FALSE;
+ step++;
+ step_rel++;
+
+ openmm_take_one_step(openmmData);
+ }
+ /* End of main MD loop */
+ debug_gmx();
+
+ /* Stop the time */
+ runtime_end(runtime);
+
+ if (MASTER(cr))
+ {
+ if (ir->nstcalcenergy > 0)
+ {
+ print_ebin(outf->fp_ene,FALSE,FALSE,FALSE,fplog,step,t,
+ eprAVER,FALSE,mdebin,fcd,groups,&(ir->opts));
+ }
+ }
+
+ openmm_cleanup(fplog, openmmData);
+
+ done_mdoutf(outf);
+
+ debug_gmx();
+
+ runtime->nsteps_done = step_rel;
+
+ return 0;
+}
--- /dev/null
- gmx_membed_t *membed=NULL;
+/* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
+ *
+ *
+ * This source code is part of
+ *
+ * G R O M A C S
+ *
+ * GROningen MAchine for Chemical Simulations
+ *
+ * VERSION 3.2.0
+ * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
+ * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
+ * Copyright (c) 2001-2004, The GROMACS development team,
+ * check out http://www.gromacs.org for more information.
+
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version.
+ *
+ * If you want to redistribute modifications, please consider that
+ * scientific software is very special. Version control is crucial -
+ * bugs must be traceable. We will be happy to consider code for
+ * inclusion in the official distribution, but derived work must not
+ * be called official GROMACS. Details are found in the README & COPYING
+ * files - if they are missing, get the official version at www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the papers on the package - you can find them in the top README file.
+ *
+ * For more info, check our website at http://www.gromacs.org
+ *
+ * And Hey:
+ * Gallium Rubidium Oxygen Manganese Argon Carbon Silicon
+ */
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+#ifdef __linux
+#define _GNU_SOURCE
+#include <sched.h>
+#include <sys/syscall.h>
+#endif
+#include <signal.h>
+#include <stdlib.h>
+#ifdef HAVE_UNISTD_H
+#include <unistd.h>
+#endif
+
+#include "typedefs.h"
+#include "smalloc.h"
+#include "sysstuff.h"
+#include "statutil.h"
+#include "mdrun.h"
+#include "network.h"
+#include "pull.h"
+#include "pull_rotation.h"
+#include "names.h"
+#include "disre.h"
+#include "orires.h"
+#include "dihre.h"
+#include "pme.h"
+#include "mdatoms.h"
+#include "repl_ex.h"
+#include "qmmm.h"
+#include "domdec.h"
+#include "partdec.h"
+#include "coulomb.h"
+#include "constr.h"
+#include "mvdata.h"
+#include "checkpoint.h"
+#include "mtop_util.h"
+#include "sighandler.h"
+#include "tpxio.h"
+#include "txtdump.h"
+#include "pull_rotation.h"
+#include "membed.h"
+#include "macros.h"
+
+#ifdef GMX_LIB_MPI
+#include <mpi.h>
+#endif
+#ifdef GMX_THREAD_MPI
+#include "tmpi.h"
+#endif
+
+#ifdef GMX_FAHCORE
+#include "corewrap.h"
+#endif
+
+#ifdef GMX_OPENMM
+#include "md_openmm.h"
+#endif
+
+#ifdef GMX_OPENMP
+#include <omp.h>
+#endif
+
+
+typedef struct {
+ gmx_integrator_t *func;
+} gmx_intp_t;
+
+/* The array should match the eI array in include/types/enums.h */
+#ifdef GMX_OPENMM /* FIXME do_md_openmm needs fixing */
+const gmx_intp_t integrator[eiNR] = { {do_md_openmm}, {do_md_openmm}, {do_md_openmm}, {do_md_openmm}, {do_md_openmm}, {do_md_openmm}, {do_md_openmm}, {do_md_openmm}, {do_md_openmm}, {do_md_openmm}, {do_md_openmm},{do_md_openmm}};
+#else
+const gmx_intp_t integrator[eiNR] = { {do_md}, {do_steep}, {do_cg}, {do_md}, {do_md}, {do_nm}, {do_lbfgs}, {do_tpi}, {do_tpi}, {do_md}, {do_md},{do_md}};
+#endif
+
+gmx_large_int_t deform_init_init_step_tpx;
+matrix deform_init_box_tpx;
+#ifdef GMX_THREAD_MPI
+tMPI_Thread_mutex_t deform_init_box_mutex=TMPI_THREAD_MUTEX_INITIALIZER;
+#endif
+
+
+#ifdef GMX_THREAD_MPI
+struct mdrunner_arglist
+{
+ FILE *fplog;
+ t_commrec *cr;
+ int nfile;
+ const t_filenm *fnm;
+ output_env_t oenv;
+ gmx_bool bVerbose;
+ gmx_bool bCompact;
+ int nstglobalcomm;
+ ivec ddxyz;
+ int dd_node_order;
+ real rdd;
+ real rconstr;
+ const char *dddlb_opt;
+ real dlb_scale;
+ const char *ddcsx;
+ const char *ddcsy;
+ const char *ddcsz;
+ int nstepout;
+ int resetstep;
+ int nmultisim;
+ int repl_ex_nst;
+ int repl_ex_seed;
+ real pforce;
+ real cpt_period;
+ real max_hours;
+ const char *deviceOptions;
+ unsigned long Flags;
+ int ret; /* return value */
+};
+
+
+/* The function used for spawning threads. Extracts the mdrunner()
+ arguments from its one argument and calls mdrunner(), after making
+ a commrec. */
+static void mdrunner_start_fn(void *arg)
+{
+ struct mdrunner_arglist *mda=(struct mdrunner_arglist*)arg;
+ struct mdrunner_arglist mc=*mda; /* copy the arg list to make sure
+ that it's thread-local. This doesn't
+ copy pointed-to items, of course,
+ but those are all const. */
+ t_commrec *cr; /* we need a local version of this */
+ FILE *fplog=NULL;
+ t_filenm *fnm;
+
+ fnm = dup_tfn(mc.nfile, mc.fnm);
+
+ cr = init_par_threads(mc.cr);
+
+ if (MASTER(cr))
+ {
+ fplog=mc.fplog;
+ }
+
+ mda->ret=mdrunner(cr->nnodes, fplog, cr, mc.nfile, fnm, mc.oenv,
+ mc.bVerbose, mc.bCompact, mc.nstglobalcomm,
+ mc.ddxyz, mc.dd_node_order, mc.rdd,
+ mc.rconstr, mc.dddlb_opt, mc.dlb_scale,
+ mc.ddcsx, mc.ddcsy, mc.ddcsz, mc.nstepout, mc.resetstep,
+ mc.nmultisim, mc.repl_ex_nst, mc.repl_ex_seed, mc.pforce,
+ mc.cpt_period, mc.max_hours, mc.deviceOptions, mc.Flags);
+}
+
+/* called by mdrunner() to start a specific number of threads (including
+ the main thread) for thread-parallel runs. This in turn calls mdrunner()
+ for each thread.
+ All options besides nthreads are the same as for mdrunner(). */
+static t_commrec *mdrunner_start_threads(int nthreads,
+ FILE *fplog,t_commrec *cr,int nfile,
+ const t_filenm fnm[], const output_env_t oenv, gmx_bool bVerbose,
+ gmx_bool bCompact, int nstglobalcomm,
+ ivec ddxyz,int dd_node_order,real rdd,real rconstr,
+ const char *dddlb_opt,real dlb_scale,
+ const char *ddcsx,const char *ddcsy,const char *ddcsz,
+ int nstepout,int resetstep,int nmultisim,int repl_ex_nst,
+ int repl_ex_seed, real pforce,real cpt_period, real max_hours,
+ const char *deviceOptions, unsigned long Flags)
+{
+ int ret;
+ struct mdrunner_arglist *mda;
+ t_commrec *crn; /* the new commrec */
+ t_filenm *fnmn;
+
+ /* first check whether we even need to start tMPI */
+ if (nthreads<2)
+ return cr;
+
+ /* a few small, one-time, almost unavoidable memory leaks: */
+ snew(mda,1);
+ fnmn=dup_tfn(nfile, fnm);
+
+ /* fill the data structure to pass as void pointer to thread start fn */
+ mda->fplog=fplog;
+ mda->cr=cr;
+ mda->nfile=nfile;
+ mda->fnm=fnmn;
+ mda->oenv=oenv;
+ mda->bVerbose=bVerbose;
+ mda->bCompact=bCompact;
+ mda->nstglobalcomm=nstglobalcomm;
+ mda->ddxyz[XX]=ddxyz[XX];
+ mda->ddxyz[YY]=ddxyz[YY];
+ mda->ddxyz[ZZ]=ddxyz[ZZ];
+ mda->dd_node_order=dd_node_order;
+ mda->rdd=rdd;
+ mda->rconstr=rconstr;
+ mda->dddlb_opt=dddlb_opt;
+ mda->dlb_scale=dlb_scale;
+ mda->ddcsx=ddcsx;
+ mda->ddcsy=ddcsy;
+ mda->ddcsz=ddcsz;
+ mda->nstepout=nstepout;
+ mda->resetstep=resetstep;
+ mda->nmultisim=nmultisim;
+ mda->repl_ex_nst=repl_ex_nst;
+ mda->repl_ex_seed=repl_ex_seed;
+ mda->pforce=pforce;
+ mda->cpt_period=cpt_period;
+ mda->max_hours=max_hours;
+ mda->deviceOptions=deviceOptions;
+ mda->Flags=Flags;
+
+ fprintf(stderr, "Starting %d threads\n",nthreads);
+ fflush(stderr);
+ /* now spawn new threads that start mdrunner_start_fn(), while
+ the main thread returns */
+ ret=tMPI_Init_fn(TRUE, nthreads, mdrunner_start_fn, (void*)(mda) );
+ if (ret!=TMPI_SUCCESS)
+ return NULL;
+
+ /* make a new comm_rec to reflect the new situation */
+ crn=init_par_threads(cr);
+ return crn;
+}
+
+
+/* Get the number of threads to use for thread-MPI based on how many
+ * were requested, which algorithms we're using,
+ * and how many particles there are.
+ */
+static int get_nthreads_mpi(int nthreads_requested, t_inputrec *inputrec,
+ gmx_mtop_t *mtop)
+{
+ int nthreads,nthreads_new;
+ int min_atoms_per_thread;
+ char *env;
+
+ nthreads = nthreads_requested;
+
+ /* determine # of hardware threads. */
+ if (nthreads_requested < 1)
+ {
+ if ((env = getenv("GMX_MAX_THREADS")) != NULL)
+ {
+ nthreads = 0;
+ sscanf(env,"%d",&nthreads);
+ if (nthreads < 1)
+ {
+ gmx_fatal(FARGS,"GMX_MAX_THREADS (%d) should be larger than 0",
+ nthreads);
+ }
+ }
+ else
+ {
+ nthreads = tMPI_Thread_get_hw_number();
+ }
+ }
+
+ if (inputrec->eI == eiNM || EI_TPI(inputrec->eI))
+ {
+ /* Steps are divided over the nodes iso splitting the atoms */
+ min_atoms_per_thread = 0;
+ }
+ else
+ {
+ min_atoms_per_thread = MIN_ATOMS_PER_THREAD;
+ }
+
+ /* Check if an algorithm does not support parallel simulation. */
+ if (nthreads != 1 &&
+ ( inputrec->eI == eiLBFGS ||
+ inputrec->coulombtype == eelEWALD ) )
+ {
+ fprintf(stderr,"\nThe integration or electrostatics algorithm doesn't support parallel runs. Not starting any threads.\n");
+ nthreads = 1;
+ }
+ else if (nthreads_requested < 1 &&
+ mtop->natoms/nthreads < min_atoms_per_thread)
+ {
+ /* the thread number was chosen automatically, but there are too many
+ threads (too few atoms per thread) */
+ nthreads_new = max(1,mtop->natoms/min_atoms_per_thread);
+
+ if (nthreads_new > 8 || (nthreads == 8 && nthreads_new > 4))
+ {
+ /* Use only multiples of 4 above 8 threads
+ * or with an 8-core processor
+ * (to avoid 6 threads on 8 core processors with 4 real cores).
+ */
+ nthreads_new = (nthreads_new/4)*4;
+ }
+ else if (nthreads_new > 4)
+ {
+ /* Avoid 5 or 7 threads */
+ nthreads_new = (nthreads_new/2)*2;
+ }
+
+ nthreads = nthreads_new;
+
+ fprintf(stderr,"\n");
+ fprintf(stderr,"NOTE: Parallelization is limited by the small number of atoms,\n");
+ fprintf(stderr," only starting %d threads.\n",nthreads);
+ fprintf(stderr," You can use the -nt option to optimize the number of threads.\n\n");
+ }
+ return nthreads;
+}
+#endif
+
+
+int mdrunner(int nthreads_requested, FILE *fplog,t_commrec *cr,int nfile,
+ const t_filenm fnm[], const output_env_t oenv, gmx_bool bVerbose,
+ gmx_bool bCompact, int nstglobalcomm,
+ ivec ddxyz,int dd_node_order,real rdd,real rconstr,
+ const char *dddlb_opt,real dlb_scale,
+ const char *ddcsx,const char *ddcsy,const char *ddcsz,
+ int nstepout,int resetstep,int nmultisim,int repl_ex_nst,
+ int repl_ex_seed, real pforce,real cpt_period,real max_hours,
+ const char *deviceOptions, unsigned long Flags)
+{
+ double nodetime=0,realtime;
+ t_inputrec *inputrec;
+ t_state *state=NULL;
+ matrix box;
+ gmx_ddbox_t ddbox={0};
+ int npme_major,npme_minor;
+ real tmpr1,tmpr2;
+ t_nrnb *nrnb;
+ gmx_mtop_t *mtop=NULL;
+ t_mdatoms *mdatoms=NULL;
+ t_forcerec *fr=NULL;
+ t_fcdata *fcd=NULL;
+ real ewaldcoeff=0;
+ gmx_pme_t *pmedata=NULL;
+ gmx_vsite_t *vsite=NULL;
+ gmx_constr_t constr;
+ int i,m,nChargePerturbed=-1,status,nalloc;
+ char *gro;
+ gmx_wallcycle_t wcycle;
+ gmx_bool bReadRNG,bReadEkin;
+ int list;
+ gmx_runtime_t runtime;
+ int rc;
+ gmx_large_int_t reset_counters;
+ gmx_edsam_t ed=NULL;
+ t_commrec *cr_old=cr;
+ int nthreads_mpi=1;
+ int nthreads_pme=1;
- fprintf(stderr,"Initializing membed");
- snew(membed,1);
- init_membed(fplog,membed,nfile,fnm,mtop,inputrec,state,cr,&cpt_period);
++ gmx_membed_t membed=NULL;
+
+ /* CAUTION: threads may be started later on in this function, so
+ cr doesn't reflect the final parallel state right now */
+ snew(inputrec,1);
+ snew(mtop,1);
+
+ if (bVerbose && SIMMASTER(cr))
+ {
+ fprintf(stderr,"Getting Loaded...\n");
+ }
+
+ if (Flags & MD_APPENDFILES)
+ {
+ fplog = NULL;
+ }
+
+ snew(state,1);
+ if (MASTER(cr))
+ {
+ /* Read (nearly) all data required for the simulation */
+ read_tpx_state(ftp2fn(efTPX,nfile,fnm),inputrec,state,NULL,mtop);
+
+ /* NOW the threads will be started: */
+#ifdef GMX_THREAD_MPI
+ nthreads_mpi = get_nthreads_mpi(nthreads_requested, inputrec, mtop);
+
+ if (nthreads_mpi > 1)
+ {
+ /* now start the threads. */
+ cr=mdrunner_start_threads(nthreads_mpi, fplog, cr_old, nfile, fnm,
+ oenv, bVerbose, bCompact, nstglobalcomm,
+ ddxyz, dd_node_order, rdd, rconstr,
+ dddlb_opt, dlb_scale, ddcsx, ddcsy, ddcsz,
+ nstepout, resetstep, nmultisim,
+ repl_ex_nst, repl_ex_seed, pforce,
+ cpt_period, max_hours, deviceOptions,
+ Flags);
+ /* the main thread continues here with a new cr. We don't deallocate
+ the old cr because other threads may still be reading it. */
+ if (cr == NULL)
+ {
+ gmx_comm("Failed to spawn threads");
+ }
+ }
+#endif
+ }
+ /* END OF CAUTION: cr is now reliable */
+
+ /* g_membed initialisation *
+ * Because we change the mtop, init_membed is called before the init_parallel *
+ * (in case we ever want to make it run in parallel) */
+ if (opt2bSet("-membed",nfile,fnm))
+ {
++ if (MASTER(cr))
++ {
++ fprintf(stderr,"Initializing membed");
++ }
++ membed = init_membed(fplog,nfile,fnm,mtop,inputrec,state,cr,&cpt_period);
+ }
+
+ if (PAR(cr))
+ {
+ /* now broadcast everything to the non-master nodes/threads: */
+ init_parallel(fplog, cr, inputrec, mtop);
+ }
+ if (fplog != NULL)
+ {
+ pr_inputrec(fplog,0,"Input Parameters",inputrec,FALSE);
+ }
+
+ /* now make sure the state is initialized and propagated */
+ set_state_entries(state,inputrec,cr->nnodes);
+
+ /* A parallel command line option consistency check that we can
+ only do after any threads have started. */
+ if (!PAR(cr) &&
+ (ddxyz[XX] > 1 || ddxyz[YY] > 1 || ddxyz[ZZ] > 1 || cr->npmenodes > 0))
+ {
+ gmx_fatal(FARGS,
+ "The -dd or -npme option request a parallel simulation, "
+#ifndef GMX_MPI
+ "but mdrun was compiled without threads or MPI enabled"
+#else
+#ifdef GMX_THREAD_MPI
+ "but the number of threads (option -nt) is 1"
+#else
+ "but mdrun was not started through mpirun/mpiexec or only one process was requested through mpirun/mpiexec"
+#endif
+#endif
+ );
+ }
+
+ if ((Flags & MD_RERUN) &&
+ (EI_ENERGY_MINIMIZATION(inputrec->eI) || eiNM == inputrec->eI))
+ {
+ gmx_fatal(FARGS, "The .mdp file specified an energy mininization or normal mode algorithm, and these are not compatible with mdrun -rerun");
+ }
+
+ if (can_use_allvsall(inputrec,mtop,TRUE,cr,fplog))
+ {
+ /* All-vs-all loops do not work with domain decomposition */
+ Flags |= MD_PARTDEC;
+ }
+
+ if (!EEL_PME(inputrec->coulombtype) || (Flags & MD_PARTDEC))
+ {
+ if (cr->npmenodes > 0)
+ {
+ if (!EEL_PME(inputrec->coulombtype))
+ {
+ gmx_fatal_collective(FARGS,cr,NULL,
+ "PME nodes are requested, but the system does not use PME electrostatics");
+ }
+ if (Flags & MD_PARTDEC)
+ {
+ gmx_fatal_collective(FARGS,cr,NULL,
+ "PME nodes are requested, but particle decomposition does not support separate PME nodes");
+ }
+ }
+
+ cr->npmenodes = 0;
+ }
+
+#ifdef GMX_FAHCORE
+ fcRegisterSteps(inputrec->nsteps,inputrec->init_step);
+#endif
+
+ /* NMR restraints must be initialized before load_checkpoint,
+ * since with time averaging the history is added to t_state.
+ * For proper consistency check we therefore need to extend
+ * t_state here.
+ * So the PME-only nodes (if present) will also initialize
+ * the distance restraints.
+ */
+ snew(fcd,1);
+
+ /* This needs to be called before read_checkpoint to extend the state */
+ init_disres(fplog,mtop,inputrec,cr,Flags & MD_PARTDEC,fcd,state);
+
+ if (gmx_mtop_ftype_count(mtop,F_ORIRES) > 0)
+ {
+ if (PAR(cr) && !(Flags & MD_PARTDEC))
+ {
+ gmx_fatal(FARGS,"Orientation restraints do not work (yet) with domain decomposition, use particle decomposition (mdrun option -pd)");
+ }
+ /* Orientation restraints */
+ if (MASTER(cr))
+ {
+ init_orires(fplog,mtop,state->x,inputrec,cr->ms,&(fcd->orires),
+ state);
+ }
+ }
+
+ if (DEFORM(*inputrec))
+ {
+ /* Store the deform reference box before reading the checkpoint */
+ if (SIMMASTER(cr))
+ {
+ copy_mat(state->box,box);
+ }
+ if (PAR(cr))
+ {
+ gmx_bcast(sizeof(box),box,cr);
+ }
+ /* Because we do not have the update struct available yet
+ * in which the reference values should be stored,
+ * we store them temporarily in static variables.
+ * This should be thread safe, since they are only written once
+ * and with identical values.
+ */
+#ifdef GMX_THREAD_MPI
+ tMPI_Thread_mutex_lock(&deform_init_box_mutex);
+#endif
+ deform_init_init_step_tpx = inputrec->init_step;
+ copy_mat(box,deform_init_box_tpx);
+#ifdef GMX_THREAD_MPI
+ tMPI_Thread_mutex_unlock(&deform_init_box_mutex);
+#endif
+ }
+
+ if (opt2bSet("-cpi",nfile,fnm))
+ {
+ /* Check if checkpoint file exists before doing continuation.
+ * This way we can use identical input options for the first and subsequent runs...
+ */
+ if( gmx_fexist_master(opt2fn_master("-cpi",nfile,fnm,cr),cr) )
+ {
+ load_checkpoint(opt2fn_master("-cpi",nfile,fnm,cr),&fplog,
+ cr,Flags & MD_PARTDEC,ddxyz,
+ inputrec,state,&bReadRNG,&bReadEkin,
+ (Flags & MD_APPENDFILES),
+ (Flags & MD_APPENDFILESSET));
+
+ if (bReadRNG)
+ {
+ Flags |= MD_READ_RNG;
+ }
+ if (bReadEkin)
+ {
+ Flags |= MD_READ_EKIN;
+ }
+ }
+ }
+
+ if (((MASTER(cr) || (Flags & MD_SEPPOT)) && (Flags & MD_APPENDFILES))
+#ifdef GMX_THREAD_MPI
+ /* With thread MPI only the master node/thread exists in mdrun.c,
+ * therefore non-master nodes need to open the "seppot" log file here.
+ */
+ || (!MASTER(cr) && (Flags & MD_SEPPOT))
+#endif
+ )
+ {
+ gmx_log_open(ftp2fn(efLOG,nfile,fnm),cr,!(Flags & MD_SEPPOT),
+ Flags,&fplog);
+ }
+
+ if (SIMMASTER(cr))
+ {
+ copy_mat(state->box,box);
+ }
+
+ if (PAR(cr))
+ {
+ gmx_bcast(sizeof(box),box,cr);
+ }
+
+ /* Essential dynamics */
+ if (opt2bSet("-ei",nfile,fnm))
+ {
+ /* Open input and output files, allocate space for ED data structure */
+ ed = ed_open(nfile,fnm,Flags,cr);
+ }
+
+ if (bVerbose && SIMMASTER(cr))
+ {
+ fprintf(stderr,"Loaded with Money\n\n");
+ }
+
+ if (PAR(cr) && !((Flags & MD_PARTDEC) ||
+ EI_TPI(inputrec->eI) ||
+ inputrec->eI == eiNM))
+ {
+ cr->dd = init_domain_decomposition(fplog,cr,Flags,ddxyz,rdd,rconstr,
+ dddlb_opt,dlb_scale,
+ ddcsx,ddcsy,ddcsz,
+ mtop,inputrec,
+ box,state->x,
+ &ddbox,&npme_major,&npme_minor);
+
+ make_dd_communicators(fplog,cr,dd_node_order);
+
+ /* Set overallocation to avoid frequent reallocation of arrays */
+ set_over_alloc_dd(TRUE);
+ }
+ else
+ {
+ /* PME, if used, is done on all nodes with 1D decomposition */
+ cr->npmenodes = 0;
+ cr->duty = (DUTY_PP | DUTY_PME);
+ npme_major = 1;
+ npme_minor = 1;
+ if (!EI_TPI(inputrec->eI))
+ {
+ npme_major = cr->nnodes;
+ }
+
+ if (inputrec->ePBC == epbcSCREW)
+ {
+ gmx_fatal(FARGS,
+ "pbc=%s is only implemented with domain decomposition",
+ epbc_names[inputrec->ePBC]);
+ }
+ }
+
+ if (PAR(cr))
+ {
+ /* After possible communicator splitting in make_dd_communicators.
+ * we can set up the intra/inter node communication.
+ */
+ gmx_setup_nodecomm(fplog,cr);
+ }
+
+ /* get number of OpenMP/PME threads
+ * env variable should be read only on one node to make sure it is identical everywhere */
+#ifdef GMX_OPENMP
+ if (EEL_PME(inputrec->coulombtype))
+ {
+ if (MASTER(cr))
+ {
+ char *ptr;
+ if ((ptr=getenv("GMX_PME_NTHREADS")) != NULL)
+ {
+ sscanf(ptr,"%d",&nthreads_pme);
+ }
+ if (fplog != NULL && nthreads_pme > 1)
+ {
+ fprintf(fplog,"Using %d threads for PME\n",nthreads_pme);
+ }
+ }
+ if (PAR(cr))
+ {
+ gmx_bcast_sim(sizeof(nthreads_pme),&nthreads_pme,cr);
+ }
+ }
+#endif
+
+ wcycle = wallcycle_init(fplog,resetstep,cr,nthreads_pme);
+ if (PAR(cr))
+ {
+ /* Master synchronizes its value of reset_counters with all nodes
+ * including PME only nodes */
+ reset_counters = wcycle_get_reset_counters(wcycle);
+ gmx_bcast_sim(sizeof(reset_counters),&reset_counters,cr);
+ wcycle_set_reset_counters(wcycle, reset_counters);
+ }
+
+
+ snew(nrnb,1);
+ if (cr->duty & DUTY_PP)
+ {
+ /* For domain decomposition we allocate dynamically
+ * in dd_partition_system.
+ */
+ if (DOMAINDECOMP(cr))
+ {
+ bcast_state_setup(cr,state);
+ }
+ else
+ {
+ if (PAR(cr))
+ {
+ bcast_state(cr,state,TRUE);
+ }
+ }
+
+ /* Dihedral Restraints */
+ if (gmx_mtop_ftype_count(mtop,F_DIHRES) > 0)
+ {
+ init_dihres(fplog,mtop,inputrec,fcd);
+ }
+
+ /* Initiate forcerecord */
+ fr = mk_forcerec();
+ init_forcerec(fplog,oenv,fr,fcd,inputrec,mtop,cr,box,FALSE,
+ opt2fn("-table",nfile,fnm),
+ opt2fn("-tabletf",nfile,fnm),
+ opt2fn("-tablep",nfile,fnm),
+ opt2fn("-tableb",nfile,fnm),FALSE,pforce);
+
+ /* version for PCA_NOT_READ_NODE (see md.c) */
+ /*init_forcerec(fplog,fr,fcd,inputrec,mtop,cr,box,FALSE,
+ "nofile","nofile","nofile","nofile",FALSE,pforce);
+ */
+ fr->bSepDVDL = ((Flags & MD_SEPPOT) == MD_SEPPOT);
+
+ /* Initialize QM-MM */
+ if(fr->bQMMM)
+ {
+ init_QMMMrec(cr,box,mtop,inputrec,fr);
+ }
+
+ /* Initialize the mdatoms structure.
+ * mdatoms is not filled with atom data,
+ * as this can not be done now with domain decomposition.
+ */
+ mdatoms = init_mdatoms(fplog,mtop,inputrec->efep!=efepNO);
+
+ /* Initialize the virtual site communication */
+ vsite = init_vsite(mtop,cr);
+
+ calc_shifts(box,fr->shift_vec);
+
+ /* With periodic molecules the charge groups should be whole at start up
+ * and the virtual sites should not be far from their proper positions.
+ */
+ if (!inputrec->bContinuation && MASTER(cr) &&
+ !(inputrec->ePBC != epbcNONE && inputrec->bPeriodicMols))
+ {
+ /* Make molecules whole at start of run */
+ if (fr->ePBC != epbcNONE)
+ {
+ do_pbc_first_mtop(fplog,inputrec->ePBC,box,mtop,state->x);
+ }
+ if (vsite)
+ {
+ /* Correct initial vsite positions are required
+ * for the initial distribution in the domain decomposition
+ * and for the initial shell prediction.
+ */
+ construct_vsites_mtop(fplog,vsite,mtop,state->x);
+ }
+ }
+
+ if (EEL_PME(fr->eeltype))
+ {
+ ewaldcoeff = fr->ewaldcoeff;
+ pmedata = &fr->pmedata;
+ }
+ else
+ {
+ pmedata = NULL;
+ }
+ }
+ else
+ {
+ /* This is a PME only node */
+
+ /* We don't need the state */
+ done_state(state);
+
+ ewaldcoeff = calc_ewaldcoeff(inputrec->rcoulomb, inputrec->ewald_rtol);
+ snew(pmedata,1);
+ }
+
+ /* Initiate PME if necessary,
+ * either on all nodes or on dedicated PME nodes only. */
+ if (EEL_PME(inputrec->coulombtype))
+ {
+ if (mdatoms)
+ {
+ nChargePerturbed = mdatoms->nChargePerturbed;
+ }
+ if (cr->npmenodes > 0)
+ {
+ /* The PME only nodes need to know nChargePerturbed */
+ gmx_bcast_sim(sizeof(nChargePerturbed),&nChargePerturbed,cr);
+ }
+
+
+ /* Set CPU affinity. Can be important for performance.
+ On some systems (e.g. Cray) CPU Affinity is set by default.
+ But default assigning doesn't work (well) with only some ranks
+ having threads. This causes very low performance.
+ External tools have cumbersome syntax for setting affinity
+ in the case that only some ranks have threads.
+ Thus it is important that GROMACS sets the affinity internally at
+ if only PME is using threads.
+ */
+
+#ifdef GMX_OPENMP
+#ifdef __linux
+#ifdef GMX_LIB_MPI
+ {
+ int core;
+ MPI_Comm comm_intra; /* intra communicator (but different to nc.comm_intra includes PME nodes) */
+ MPI_Comm_split(MPI_COMM_WORLD,gmx_hostname_num(),gmx_node_rank(),&comm_intra);
+ int local_omp_nthreads = (cr->duty & DUTY_PME) ? nthreads_pme : 1; /* threads on this node */
+ MPI_Scan(&local_omp_nthreads,&core, 1, MPI_INT, MPI_SUM, comm_intra);
+ core-=local_omp_nthreads; /* make exclusive scan */
+#pragma omp parallel firstprivate(core) num_threads(local_omp_nthreads)
+ {
+ cpu_set_t mask;
+ CPU_ZERO(&mask);
+ core+=omp_get_thread_num();
+ CPU_SET(core,&mask);
+ sched_setaffinity((pid_t) syscall (SYS_gettid),sizeof(cpu_set_t),&mask);
+ }
+ }
+#endif /*GMX_MPI*/
+#endif /*__linux*/
+#endif /*GMX_OPENMP*/
+
+ if (cr->duty & DUTY_PME)
+ {
+ status = gmx_pme_init(pmedata,cr,npme_major,npme_minor,inputrec,
+ mtop ? mtop->natoms : 0,nChargePerturbed,
+ (Flags & MD_REPRODUCIBLE),nthreads_pme);
+ if (status != 0)
+ {
+ gmx_fatal(FARGS,"Error %d initializing PME",status);
+ }
+ }
+ }
+
+
+ if (integrator[inputrec->eI].func == do_md
+#ifdef GMX_OPENMM
+ ||
+ integrator[inputrec->eI].func == do_md_openmm
+#endif
+ )
+ {
+ /* Turn on signal handling on all nodes */
+ /*
+ * (A user signal from the PME nodes (if any)
+ * is communicated to the PP nodes.
+ */
+ signal_handler_install();
+ }
+
+ if (cr->duty & DUTY_PP)
+ {
+ if (inputrec->ePull != epullNO)
+ {
+ /* Initialize pull code */
+ init_pull(fplog,inputrec,nfile,fnm,mtop,cr,oenv,
+ EI_DYNAMICS(inputrec->eI) && MASTER(cr),Flags);
+ }
+
+ if (inputrec->bRot)
+ {
+ /* Initialize enforced rotation code */
+ init_rot(fplog,inputrec,nfile,fnm,cr,state->x,state->box,mtop,oenv,
+ bVerbose,Flags);
+ }
+
+ constr = init_constraints(fplog,mtop,inputrec,ed,state,cr);
+
+ if (DOMAINDECOMP(cr))
+ {
+ dd_init_bondeds(fplog,cr->dd,mtop,vsite,constr,inputrec,
+ Flags & MD_DDBONDCHECK,fr->cginfo_mb);
+
+ set_dd_parameters(fplog,cr->dd,dlb_scale,inputrec,fr,&ddbox);
+
+ setup_dd_grid(fplog,cr->dd);
+ }
+
+ /* Now do whatever the user wants us to do (how flexible...) */
+ integrator[inputrec->eI].func(fplog,cr,nfile,fnm,
+ oenv,bVerbose,bCompact,
+ nstglobalcomm,
+ vsite,constr,
+ nstepout,inputrec,mtop,
+ fcd,state,
+ mdatoms,nrnb,wcycle,ed,fr,
+ repl_ex_nst,repl_ex_seed,
+ membed,
+ cpt_period,max_hours,
+ deviceOptions,
+ Flags,
+ &runtime);
+
+ if (inputrec->ePull != epullNO)
+ {
+ finish_pull(fplog,inputrec->pull);
+ }
+
+ if (inputrec->bRot)
+ {
+ finish_rot(fplog,inputrec->rot);
+ }
+
+ }
+ else
+ {
+ /* do PME only */
+ gmx_pmeonly(*pmedata,cr,nrnb,wcycle,ewaldcoeff,FALSE,inputrec);
+ }
+
+ if (EI_DYNAMICS(inputrec->eI) || EI_TPI(inputrec->eI))
+ {
+ /* Some timing stats */
+ if (SIMMASTER(cr))
+ {
+ if (runtime.proc == 0)
+ {
+ runtime.proc = runtime.real;
+ }
+ }
+ else
+ {
+ runtime.real = 0;
+ }
+ }
+
+ wallcycle_stop(wcycle,ewcRUN);
+
+ /* Finish up, write some stuff
+ * if rerunMD, don't write last frame again
+ */
+ finish_run(fplog,cr,ftp2fn(efSTO,nfile,fnm),
+ inputrec,nrnb,wcycle,&runtime,
+ EI_DYNAMICS(inputrec->eI) && !MULTISIM(cr));
+
+ if (opt2bSet("-membed",nfile,fnm))
+ {
+ sfree(membed);
+ }
+
+ /* Does what it says */
+ print_date_and_time(fplog,cr->nodeid,"Finished mdrun",&runtime);
+
+ /* Close logfile already here if we were appending to it */
+ if (MASTER(cr) && (Flags & MD_APPENDFILES))
+ {
+ gmx_log_close(fplog);
+ }
+
+ rc=(int)gmx_get_stop_condition();
+
+#ifdef GMX_THREAD_MPI
+ /* we need to join all threads. The sub-threads join when they
+ exit this function, but the master thread needs to be told to
+ wait for that. */
+ if (PAR(cr) && MASTER(cr))
+ {
+ tMPI_Finalize();
+ }
+#endif
+
+ return rc;
+}
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff