From: Roland Schulz Date: Fri, 18 May 2012 01:00:16 +0000 (-0400) Subject: Merge remote-tracking branch 'origin/release-4-6' X-Git-Url: http://biod.pnpi.spb.ru/gitweb/?a=commitdiff_plain;h=1698224bd6251c2b0ece0fd375c0c095ce5eddb8;p=alexxy%2Fgromacs.git Merge remote-tracking branch 'origin/release-4-6' Conflicts: include/membed.h src/gromacs/legacyheaders/membed.h src/kernel/membed.h All 3 conflicts old names. Renamed to src/programs/mdrun/membed.h with content of src/kernel/membed.h from 4.6 Change-Id: I64654ce18162fa7887f17a7c565bc5783c3a1238 --- 1698224bd6251c2b0ece0fd375c0c095ce5eddb8 diff --cc src/gromacs/gmxlib/index.c index 6fff42416e,0000000000..c92b580903 mode 100644,000000..100644 --- a/src/gromacs/gmxlib/index.c +++ b/src/gromacs/gmxlib/index.c @@@ -1,1156 -1,0 +1,1160 @@@ +/* + * + * 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 +#endif + +#include +#include +#include +#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; (inr); i++) { + fprintf(out,"[ %s ]\n",gnames[i]); + for(k=0,j=b->index[i]; jindex[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; (ia[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; ia[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 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; (inr); 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; (inres); 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; (knr); 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; (lnr); + naid=0; + for(j=0; (jnr); 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; (katomname[aid[k]]; + for(l=0; (l 1) { + for(l=0; (latomname[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; (inres); 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; (nnr); n++) { + if (0 == gmx_strcasecmp(restype[atoms->atom[n].resind],"Protein")) { + match=FALSE; + for(j=0; (jatomname[n]; + while (isdigit(atnm[0])) { + atnm++; + } + if ( (constructing_data[i].wholename==-1) || (jatom[n].resind < npres) && (nnr));) { + resind = atoms->atom[n].resind; + for(;((atoms->atom[n].resind==resind) && (nnr));n++) { + match=FALSE; + for(j=0;(jatomname[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) && (nnr));) { + resind = atoms->atom[n].resind; + hold = -1; + for(;((atoms->atom[n].resind==resind) && (nnr));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;in && 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;in;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;in;i++) + { + p=rt->restype[i]; + found=0; + for(j=0;jn; +} + +/* 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;in && 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;inr;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;inres;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;knres,p_typename,ntypes); + } + + for(k=0;k0) + { + 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;inr;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];iindex[iwater+1];i++) + { + gb->a[gb->nra++] = gb->a[i]; + } + } + if(nion>0) + { + for(i=gb->index[iion];iindex[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= 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; (inr); 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; (ja[b->index[i]+j]=a; + } + } + } + gmx_fio_fclose(in); + + for(i=0; (inr); i++) { + for(j=b->index[i]; (jindex[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= 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; (inr); i++) + fprintf(stderr,"Group %5d (%15s) has %5d elements\n",i,grpname[i], + grps->index[i+1]-grps->index[i]); + for(i=0; (inr > 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; (ja[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; (iclust->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; (iclust->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; +} + diff --cc src/gromacs/gmxlib/oenv.c index edbc6328fb,0000000000..015b25b94d mode 100644,000000..100644 --- a/src/gromacs/gmxlib/oenv.c +++ b/src/gromacs/gmxlib/oenv.c @@@ -1,261 -1,0 +1,264 @@@ +/* -*- 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 +#endif + + +#include ++#include +#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; icmd_line,cmdlength+argc+1); + if (argv) + { + for (i=0; icmd_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; iview; +} + +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; +} + + diff --cc src/gromacs/gmxlib/string2.c index 9048a2a673,0000000000..a15a961672 mode 100644,000000..100644 --- a/src/gromacs/gmxlib/string2.c +++ b/src/gromacs/gmxlib/string2.c @@@ -1,609 -1,0 +1,611 @@@ +/* + * + * 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 +#endif +#include "gromacs/utility/gmx_header_config.h" + +#ifdef GMX_CRAY_XT3 +#undef HAVE_PWD_H +#endif + +#include +#include +#include +#include +#include +#include + +#ifdef HAVE_SYS_TIME_H +#include +#endif + +#ifdef HAVE_UNISTD_H +#include +#endif + +#ifdef HAVE_PWD_H +#include +#endif +#include ++#include + +#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) + { + 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) ) { + /* 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 +#endif ++#include + +#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; iQPinv[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 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 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 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;dcompress[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;dref_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;depct) { + case epctANISOTROPIC: + for(d=0;depct)); + break; + } + + maxchange=0; + for(d=0;dmaxchange) + 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;dcompress[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; dref_p[d][d] - scalar_pressure) /DIM; + } + break; + case epctSEMIISOTROPIC: + for(d=0; dref_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; dref_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; dref_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; nopts; + + for(i=0; (ingtc); 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; (ingtc); 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;iveta),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; itcstat[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;nstart+md->homenr;n++) + { + if (md->cTC) + { + gc = md->cTC[n]; + } + for (d=0;dv[n][d] *= scalefac[gc]; + } + + if (debug) + { + for (d=0;dv[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;dnrdf[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; (itau_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;dWinvm[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; itau_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;jQinv[i*nh+j] = 1.0/(sqr(opts->tau_t[i]/M_2PI)*ndj*kT); + } + } + else + { + reft=0.0; + for (j=0;jQinv[i*nh+j] = 0.0; + } + } + } + } + + /* first, initialize clear all the trotter calls */ + snew(trotter_seq,ettTSEQMAX); + for (i=0;ieI==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;iQPinv[i*opts->nhchainlength+j] = 1.0/(sqr(opts->tau_t[0]/M_2PI)*qmass*kT); + } + } + } + else + { + for (i=0;iQPinv[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;innhpres;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 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; iopts.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 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; (ingtc); 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; ingtc; 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; nngtc;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]; + } + } +} diff --cc src/gromacs/mdlib/force.c index fc73494ad0,0000000000..0c8751ec86 mode 100644,000000..100644 --- a/src/gromacs/mdlib/force.c +++ b/src/gromacs/mdlib/force.c @@@ -1,727 -1,0 +1,729 @@@ +/* -*- 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 +#endif + +#include +#include ++#include +#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; (ibSepDVDL && 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;inr;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; in_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; igrpp.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; in_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; iterm[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; (igrpp.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; (igrpp.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; (igrpp; + 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; (iterm[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; in_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; (igrpp.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; in_lambda; i++) + { + enerd->enerpart_lambda[i] = 0.0; + } + } +} diff --cc src/programs/mdrun/md.c index 3ffb6fd138,0000000000..d78534049c mode 100644,000000..100644 --- a/src/programs/mdrun/md.c +++ b/src/programs/mdrun/md.c @@@ -1,1866 -1,0 +1,1866 @@@ +/* -*- 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 +#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 +#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, ++ 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" + " - 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<mols.nr; + snew(grpindex,gnx); + for(i=0; (i 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<start; istart+mdatoms->homenr; i++) + { + for(m=0; mopts.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; (iopts.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; inatoms; i++) + { + copy_rvec(rerun_fr.x[i],state_global->x[i]); + } + if (rerun_fr.bV) + { + for(i=0; inatoms; i++) + { + copy_rvec(rerun_fr.v[i],state_global->v[i]); + } + } + else + { + for(i=0; inatoms; 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; (iinatoms); 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<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_relnatoms,&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<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; +} diff --cc src/programs/mdrun/md_openmm.c index 01d2e74524,0000000000..e16f809e27 mode 100644,000000..100644 --- a/src/programs/mdrun/md_openmm.c +++ b/src/programs/mdrun/md_openmm.c @@@ -1,572 -1,0 +1,572 @@@ +/* -*- 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 +#endif + +#include +#include + +#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, ++ 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:\""); + } + 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<start; istart+mdatoms->homenr; i++) + { + for (m=0; mopts.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; +} diff --cc src/programs/mdrun/runner.c index be455a7805,0000000000..1025372f6a mode 100644,000000..100644 --- a/src/programs/mdrun/runner.c +++ b/src/programs/mdrun/runner.c @@@ -1,980 -1,0 +1,982 @@@ +/* -*- 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 +#endif +#ifdef __linux +#define _GNU_SOURCE +#include +#include +#endif +#include +#include +#ifdef HAVE_UNISTD_H +#include +#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 +#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 +#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; - gmx_membed_t *membed=NULL; ++ 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)) + { - fprintf(stderr,"Initializing membed"); - snew(membed,1); - init_membed(fplog,membed,nfile,fnm,mtop,inputrec,state,cr,&cpt_period); ++ 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; +}