#include "pull.h"
#include "update.h"
-#ifdef GMX_THREADS
+
++#ifdef GMX_THREAD_MPI
+ #include "thread_mpi/threads.h"
+ #endif
+
#ifdef __cplusplus
extern "C" {
#endif
*/
extern gmx_large_int_t deform_init_init_step_tpx;
extern matrix deform_init_box_tpx;
-#ifdef GMX_THREADS
+#ifdef GMX_THREAD_MPI
extern tMPI_Thread_mutex_t deform_init_box_mutex;
/* The minimum number of atoms per thread. With fewer atoms than this,
gmx_edsam_t ed,
t_forcerec *fr,
int repl_ex_nst,int repl_ex_seed,
+ gmx_membed_t membed,
real cpt_period,real max_hours,
const char *deviceOptions,
unsigned long Flags,
#ifdef GMX_LIB_MPI
#include <mpi.h>
#else
-#ifdef GMX_THREADS
+#ifdef GMX_THREAD_MPI
- #include "../tmpi.h"
+ #include "../thread_mpi/tmpi.h"
+ #include "../thread_mpi/mpi_bindings.h"
#else
typedef void* MPI_Comm;
typedef void* MPI_Request;
mpi_in_place_buf_t *mpb;
} t_commrec;
-#define MASTERNODE(cr) ((cr)->nodeid == 0)
+#define MASTERNODE(cr) (((cr)->nodeid == 0) || !PAR(cr))
/* #define MASTERTHREAD(cr) ((cr)->threadid == 0) */
/* #define MASTER(cr) (MASTERNODE(cr) && MASTERTHREAD(cr)) */
#define MASTER(cr) MASTERNODE(cr)
-#define SIMMASTER(cr) (MASTER(cr) && ((cr)->duty & DUTY_PP))
+#define SIMMASTER(cr) ((MASTER(cr) && ((cr)->duty & DUTY_PP)) || !PAR(cr))
#define NODEPAR(cr) ((cr)->nnodes > 1)
/* #define THREADPAR(cr) ((cr)->nthreads > 1) */
/* #define PAR(cr) (NODEPAR(cr) || THREADPAR(cr)) */
#define RANK(cr,nodeid) (nodeid)
#define MASTERRANK(cr) (0)
-#define DOMAINDECOMP(cr) ((cr)->dd != NULL)
+#define DOMAINDECOMP(cr) (((cr)->dd != NULL) && PAR(cr))
#define DDMASTER(dd) ((dd)->rank == (dd)->masterrank)
#define PARTDECOMP(cr) ((cr)->pd != NULL)
#include "writeps.h"
#include "macros.h"
#include "xvgr.h"
-#include "pppm.h"
#include "gmxfio.h"
--#ifdef GMX_THREADS
- #include "thread_mpi.h"
++#ifdef GMX_THREAD_MPI
+ #include "thread_mpi/threads.h"
#endif
#define p2(x) ((x)*(x))
#define p4(x) ((x)*(x)*(x)*(x))
static real A,A_3,B,B_4,C,c1,c2,c3,c4,c5,c6,One_4pi,FourPi_V,Vol,N0;
-#ifdef GMX_THREADS
+#ifdef GMX_THREAD_MPI
static tMPI_Thread_mutex_t shift_mutex=TMPI_THREAD_MUTEX_INITIALIZER;
#endif
void set_shift_consts(FILE *log,real r1,real rc,rvec box,t_forcerec *fr)
{
-#ifdef GMX_THREADS
+#ifdef GMX_THREAD_MPI
/* at the very least we shouldn't allow multiple threads to set these
simulataneously */
tMPI_Thread_mutex_lock(&shift_mutex);
}
One_4pi = 1.0/(4.0*M_PI);
-#ifdef GMX_THREADS
+#ifdef GMX_THREAD_MPI
tMPI_Thread_mutex_unlock(&shift_mutex);
#endif
}
-
-real gk(real k,real rc,real r1)
-/* Spread function in Fourier space. Eqs. 56-64 */
-{
- real N,gg;
-
- N = N0*p4(k);
-
- /* c1 thru c6 consts are global variables! */
- gg = (FourPi_V / (N*k)) *
- ( c1*k*cos(k*rc) + (c2+c3*k*k)*sin(k*rc) +
- c4*k*cos(k*r1) + (c5+c6*k*k)*sin(k*r1) );
-
- return gg;
-}
-
-real gknew(real k,real rc,real r1)
-{
- real rck,rck2;
-
- rck = rc*k;
- rck2= rck*rck;
-
- return -15.0*((rck2-3.0)*sin(rck) + 3*rck*cos(rck))/(Vol*rck2*rck2*rck);
-}
-
-real calc_dx2dx(rvec xi,rvec xj,rvec box,rvec dx)
-{
- int m;
- real ddx,dx2;
-
- dx2=0;
- for(m=0; (m<DIM); m++) {
- ddx=xj[m]-xi[m];
- if (ddx < -0.5*box[m])
- ddx+=box[m];
- else if (ddx >= 0.5*box[m])
- ddx-=box[m];
- dx[m]=ddx;
- dx2 += ddx*ddx;
- }
- return dx2;
-}
-
-real calc_dx2(rvec xi,rvec xj,rvec box)
-{
- rvec dx;
-
- return calc_dx2dx(xi,xj,box,dx);
-}
-
-real shiftfunction(real r1,real rc,real R)
-{
- real dr;
-
- if (R <= r1)
- return 0.0;
- else if (R >= rc)
- return -1.0/(R*R);
-
- dr=R-r1;
-
- return A*dr*dr+B*dr*dr*dr;
-}
-
-real new_f(real r,real rc)
-{
- real rrc,rrc2,rrc3;
-
- rrc = r/rc;
- rrc2 = rrc*rrc;
- rrc3 = rrc2*rrc;
- return 1.5*rrc2*rrc3 - 2.5*rrc3 + 1.0;
-}
-
-real new_phi(real r,real rc)
-{
- real rrr;
-
- rrr = sqr(r/rc);
-
- return 1/r-(0.125/rc)*(15 + 3*rrr*rrr - 10*rrr);
-}
-
-real old_f(real r,real rc,real r1)
-{
- real dr,r2;
-
- if (r <= r1)
- return 1.0;
- else if (r >= rc)
- return 0;
-
- dr = r-r1;
- r2 = r*r;
- return 1+A*r2*dr*dr+B*r2*dr*dr*dr;
-}
-
-real old_phi(real r,real rc,real r1)
-{
- real dr;
-
- if (r <= r1)
- return 1/r-C;
- else if (r >= rc)
- return 0.0;
-
- dr = r-r1;
-
- return 1/r-A_3*dr*dr*dr-B_4*dr*dr*dr*dr - C;
-}
-
-real spreadfunction(real r1,real rc,real R)
-{
- real dr;
-
- if (R <= r1)
- return 0.0;
- else if (R >= rc)
- return 0.0;
-
- dr=R-r1;
-
- return -One_4pi*(dr/R)*(2*A*(2*R-r1)+B*dr*(5*R-2*r1));
-}
-
-real potential(real r1,real rc,real R)
-{
- if (R < r1)
- return (1.0/R-C);
- else if (R <= rc)
- return (1.0/R - A_3*p3(R-r1) - B_4*p4(R-r1) - C);
- else
- return 0.0;
-}
-
-
-
-real shift_LRcorrection(FILE *fp,int start,int natoms,
- t_commrec *cr,t_forcerec *fr,
- real charge[],t_blocka *excl,rvec x[],
- gmx_bool bOld,matrix box,matrix lr_vir)
-{
- static gmx_bool bFirst=TRUE;
- static real Vself;
- int i,i1,i2,j,k,m,iv,jv;
- int *AA;
- double qq; /* Necessary for precision */
- double isp=0.564189583547756;
- real qi,dr,dr2,dr_1,dr_3,fscal,Vexcl,qtot=0;
- rvec df,dx;
- real r1=fr->rcoulomb_switch;
- real rc=fr->rcoulomb;
- ivec shift;
-
- if (bFirst) {
- qq =0;
- for(i=start; (i<start+natoms); i++)
- qq += charge[i]*charge[i];
-
- /* Obsolete, until we implement dipole and charge corrections.
- qtot=0;
- for(i=0;i<nsb->natoms;i++)
- qtot+=charge[i];
- */
-
- Vself = 0.5*C*ONE_4PI_EPS0*qq;
- if(debug) {
- fprintf(fp,"Long Range corrections for shifted interactions:\n");
- fprintf(fp,"r1 = %g, rc=%g\n",r1,rc);
- fprintf(fp,"start=%d,natoms=%d\n",start,natoms);
- fprintf(fp,"qq = %g, Vself=%g\n",qq,Vself);
- }
-
- }
- AA = excl->a;
- Vexcl = 0;
-
- for(i=start; (i<start+natoms); i++) {
- /* Initiate local variables (for this i-particle) to 0 */
- i1 = excl->index[i];
- i2 = excl->index[i+1];
- qi = charge[i]*ONE_4PI_EPS0;
-
- /* Loop over excluded neighbours */
- for(j=i1; (j<i2); j++) {
- k = AA[j];
- /*
- * First we must test whether k <> i, and then, because the
- * exclusions are all listed twice i->k and k->i we must select
- * just one of the two.
- * As a minor optimization we only compute forces when the charges
- * are non-zero.
- */
- if (k > i) {
- qq = qi*charge[k];
- if (qq != 0.0) {
- dr2 = 0;
- rvec_sub(x[i],x[k],dx);
- for(m=DIM-1; m>=0; m--) {
- if (dx[m] > 0.5*box[m][m])
- rvec_dec(dx,box[m]);
- else if (dx[m] < -0.5*box[m][m])
- rvec_inc(dx,box[m]);
-
- dr2 += dx[m]*dx[m];
- }
- dr_1 = gmx_invsqrt(dr2);
- dr = 1.0/dr_1;
- dr_3 = dr_1*dr_1*dr_1;
- /* Compute exclusion energy and scalar force */
-
- Vexcl += qq*(dr_1-potential(r1,rc,dr));
- fscal = qq*(-shiftfunction(r1,rc,dr))*dr_3;
-
- if ((fscal != 0) && debug)
- fprintf(debug,"i: %d, k: %d, dr: %.3f fscal: %.3f\n",i,k,dr,fscal);
-
- /* The force vector is obtained by multiplication with the
- * distance vector
- */
- svmul(fscal,dx,df);
- rvec_inc(fr->f_novirsum[k],df);
- rvec_dec(fr->f_novirsum[i],df);
- for(iv=0;iv<DIM;iv++)
- for(jv=0;jv<DIM;jv++)
- lr_vir[iv][jv]+=0.5*dx[iv]*df[jv];
- }
- }
- }
- }
- if (bFirst && debug)
- fprintf(fp,"Long Range correction: Vexcl=%g\n",Vexcl);
-
- bFirst = FALSE;
- /* Return the correction to the energy */
- return (-(Vself+Vexcl));
-}
-
-real phi_aver(int natoms,real phi[])
-{
- real phitot;
- int i;
-
- phitot=0;
- for(i=0; (i<natoms); i++)
- phitot=phitot+phi[i];
-
- return (phitot/natoms);
-}
-
-real symmetrize_phi(FILE *log,int natoms,real phi[],gmx_bool bVerbose)
-{
- real phitot;
- int i;
-
- phitot=phi_aver(natoms,phi);
- if (bVerbose)
- fprintf(log,"phi_aver = %10.3e\n",phitot);
-
- for(i=0; (i<natoms); i++)
- phi[i]-=phitot;
-
- return phitot;
-}
-
-static real rgbset(real col)
-{
- int icol32;
-
- icol32=32.0*col;
- return icol32/32.0;
-}
-
-
-
-real analyse_diff(FILE *log,char *label,const output_env_t oenv,
- int natom,rvec ffour[],rvec fpppm[],
- real phi_f[],real phi_p[],real phi_sr[],
- char *fcorr,char *pcorr,char *ftotcorr,char *ptotcorr)
-/* Analyse difference between forces from fourier (_f) and other (_p)
- * LR solvers (and potential also).
- * If the filenames are given, xvgr files are written.
- * returns the root mean square error in the force.
- */
-{
- int i,m;
- FILE *fp=NULL,*gp=NULL;
- real f2sum=0,p2sum=0;
- real df,fmax,dp,pmax,rmsf;
-
- fmax = fabs(ffour[0][0]-fpppm[0][0]);
- pmax = fabs(phi_f[0] - phi_p[0]);
-
- for(i=0; (i<natom); i++) {
- dp = fabs(phi_f[i] - phi_p[i]);
- pmax = max(dp,pmax);
- p2sum += dp*dp;
- for(m=0; (m<DIM); m++) {
- df = fabs(ffour[i][m] - fpppm[i][m]);
- fmax = max(df,fmax);
- f2sum += df*df;
- }
- }
-
- rmsf = sqrt(f2sum/(3.0*natom));
- fprintf(log,"\n********************************\nERROR ANALYSIS for %s\n",
- label);
- fprintf(log,"%-10s%12s%12s\n","Error:","Max Abs","RMS");
- fprintf(log,"%-10s %10.3f %10.3f\n","Force",fmax,rmsf);
- fprintf(log,"%-10s %10.3f %10.3f\n","Potential",pmax,sqrt(p2sum/(natom)));
-
- if (fcorr) {
- fp = xvgropen(fcorr,"LR Force Correlation","Four-Force","PPPM-Force",oenv);
- for(i=0; (i<natom); i++) {
- for(m=0; (m<DIM); m++) {
- fprintf(fp,"%10.3f %10.3f\n",ffour[i][m],fpppm[i][m]);
- }
- }
- gmx_fio_fclose(fp);
- do_view(oenv,fcorr,NULL);
- }
- if (pcorr)
- fp = xvgropen(pcorr,"LR Potential Correlation","Four-Pot","PPPM-Pot",oenv);
- if (ptotcorr)
- gp = xvgropen(ptotcorr,"Total Potential Correlation","Four-Pot","PPPM-Pot",
- oenv);
- for(i=0; (i<natom); i++) {
- if (pcorr)
- fprintf(fp,"%10.3f %10.3f\n",phi_f[i],phi_p[i]);
- if (ptotcorr)
- fprintf(gp,"%10.3f %10.3f\n",phi_f[i]+phi_sr[i],phi_p[i]+phi_sr[i]);
- }
- if (pcorr) {
- gmx_fio_fclose(fp);
- do_view(oenv,pcorr,NULL);
- }
- if (ptotcorr) {
- gmx_fio_fclose(gp);
- do_view(oenv,ptotcorr,NULL);
- }
-
- return rmsf;
-}
-
#endif
/* This file is completely threadsafe - keep it that way! */
-#ifdef GMX_THREADS
+#ifdef GMX_THREAD_MPI
#include <thread_mpi.h>
#endif
#include "vec.h"
#include "mtop_util.h"
+#define TPX_TAG_RELEASE "release"
+
+/* This is the tag string which is stored in the tpx file.
+ * Change this if you want to change the tpx format in a feature branch.
+ * This ensures that there will not be different tpx formats around which
+ * can not be distinguished.
+ */
+static const char *tpx_tag = TPX_TAG_RELEASE;
+
/* This number should be increased whenever the file format changes! */
-static const int tpx_version = 73;
+static const int tpx_version = 78;
/* This number should only be increased when you edit the TOPOLOGY section
* of the tpx format. This way we can maintain forward compatibility too
* to the end of the tpx file, so we can just skip it if we only
* want the topology.
*/
-static const int tpx_generation = 23;
+static const int tpx_generation = 24;
/* This number should be the most recent backwards incompatible version
* I.e., if this number is 9, we cannot read tpx version 9 with this code.
{ 43, F_TABBONDS },
{ 43, F_TABBONDSNC },
{ 70, F_RESTRBONDS },
+ { 76, F_LINEAR_ANGLES },
{ 30, F_CROSS_BOND_BONDS },
{ 30, F_CROSS_BOND_ANGLES },
{ 30, F_UREY_BRADLEY },
{ 69, F_VTEMP },
{ 66, F_PDISPCORR },
{ 54, F_DHDL_CON },
+ { 76, F_ANHARM_POL }
};
#define NFTUPD asize(ftupd)
do_pullgrp(fio,&pull->grp[g],bRead,file_version);
}
+
+static void do_rotgrp(t_fileio *fio, t_rotgrp *rotg,gmx_bool bRead, int file_version)
+{
+ gmx_bool bDum=TRUE;
+ int i;
+
+ gmx_fio_do_int(fio,rotg->eType);
+ gmx_fio_do_int(fio,rotg->bMassW);
+ gmx_fio_do_int(fio,rotg->nat);
+ if (bRead)
+ snew(rotg->ind,rotg->nat);
+ gmx_fio_ndo_int(fio,rotg->ind,rotg->nat);
+ if (bRead)
+ snew(rotg->x_ref,rotg->nat);
+ gmx_fio_ndo_rvec(fio,rotg->x_ref,rotg->nat);
+ gmx_fio_do_rvec(fio,rotg->vec);
+ gmx_fio_do_rvec(fio,rotg->pivot);
+ gmx_fio_do_real(fio,rotg->rate);
+ gmx_fio_do_real(fio,rotg->k);
+ gmx_fio_do_real(fio,rotg->slab_dist);
+ gmx_fio_do_real(fio,rotg->min_gaussian);
+ gmx_fio_do_real(fio,rotg->eps);
+ gmx_fio_do_int(fio,rotg->eFittype);
+ gmx_fio_do_int(fio,rotg->PotAngle_nstep);
+ gmx_fio_do_real(fio,rotg->PotAngle_step);
+}
+
+static void do_rot(t_fileio *fio, t_rot *rot,gmx_bool bRead, int file_version)
+{
+ int g;
+
+ gmx_fio_do_int(fio,rot->ngrp);
+ gmx_fio_do_int(fio,rot->nstrout);
+ gmx_fio_do_int(fio,rot->nstsout);
+ if (bRead)
+ snew(rot->grp,rot->ngrp);
+ for(g=0; g<rot->ngrp; g++)
+ do_rotgrp(fio, &rot->grp[g],bRead,file_version);
+}
+
+
static void do_inputrec(t_fileio *fio, t_inputrec *ir,gmx_bool bRead,
int file_version, real *fudgeQQ)
{
if (file_version != tpx_version)
{
/* Give a warning about features that are not accessible */
- fprintf(stderr,"Note: tpx file_version %d, software version %d\n",
+ fprintf(stderr,"Note: file tpx version %d, software tpx version %d\n",
file_version,tpx_version);
}
gmx_fio_do_real(fio,ir->userreal3);
gmx_fio_do_real(fio,ir->userreal4);
+ /* AdResS stuff */
+ if (file_version >= 75) {
+ gmx_fio_do_gmx_bool(fio,ir->bAdress);
+ if(ir->bAdress){
+ if (bRead) snew(ir->adress, 1);
+ gmx_fio_do_int(fio,ir->adress->type);
+ gmx_fio_do_real(fio,ir->adress->const_wf);
+ gmx_fio_do_real(fio,ir->adress->ex_width);
+ gmx_fio_do_real(fio,ir->adress->hy_width);
+ gmx_fio_do_int(fio,ir->adress->icor);
+ gmx_fio_do_int(fio,ir->adress->site);
+ gmx_fio_do_rvec(fio,ir->adress->refs);
+ gmx_fio_do_int(fio,ir->adress->n_tf_grps);
+ gmx_fio_do_real(fio, ir->adress->ex_forcecap);
+ gmx_fio_do_int(fio, ir->adress->n_energy_grps);
+ gmx_fio_do_int(fio,ir->adress->do_hybridpairs);
+
+ if (bRead)snew(ir->adress->tf_table_index,ir->adress->n_tf_grps);
+ if (ir->adress->n_tf_grps > 0) {
+ bDum=gmx_fio_ndo_int(fio,ir->adress->tf_table_index,ir->adress->n_tf_grps);
+ }
+ if (bRead)snew(ir->adress->group_explicit,ir->adress->n_energy_grps);
+ if (ir->adress->n_energy_grps > 0) {
+ bDum=gmx_fio_ndo_int(fio, ir->adress->group_explicit,ir->adress->n_energy_grps);
+ }
+ }
+ } else {
+ ir->bAdress = FALSE;
+ }
+
/* pull stuff */
if (file_version >= 48) {
gmx_fio_do_int(fio,ir->ePull);
ir->ePull = epullNO;
}
+ /* Enforced rotation */
+ if (file_version >= 74) {
+ gmx_fio_do_int(fio,ir->bRot);
+ if (ir->bRot == TRUE) {
+ if (bRead)
+ snew(ir->rot,1);
+ do_rot(fio, ir->rot,bRead,file_version);
+ }
+ } else {
+ ir->bRot = FALSE;
+ }
+
/* grpopts stuff */
gmx_fio_do_int(fio,ir->opts.ngtc);
if (file_version >= 69) {
iparams->pdihs.cpB = iparams->pdihs.cpA;
}
break;
+ case F_LINEAR_ANGLES:
+ gmx_fio_do_real(fio,iparams->linangle.klinA);
+ gmx_fio_do_real(fio,iparams->linangle.aA);
+ gmx_fio_do_real(fio,iparams->linangle.klinB);
+ gmx_fio_do_real(fio,iparams->linangle.aB);
+ break;
case F_FENEBONDS:
gmx_fio_do_real(fio,iparams->fene.bm);
gmx_fio_do_real(fio,iparams->fene.kb);
case F_POLARIZATION:
gmx_fio_do_real(fio,iparams->polarize.alpha);
break;
+ case F_ANHARM_POL:
+ gmx_fio_do_real(fio,iparams->anharm_polarize.alpha);
+ gmx_fio_do_real(fio,iparams->anharm_polarize.drcut);
+ gmx_fio_do_real(fio,iparams->anharm_polarize.khyp);
+ break;
case F_WATER_POL:
if (file_version < 31)
gmx_fatal(FARGS,"Old tpr files with water_polarization not supported. Make a new.");
}
}
+static void add_settle_atoms(t_ilist *ilist)
+{
+ int i;
+
+ /* Settle used to only store the first atom: add the other two */
+ srenew(ilist->iatoms,2*ilist->nr);
+ for(i=ilist->nr/2-1; i>=0; i--)
+ {
+ ilist->iatoms[4*i+0] = ilist->iatoms[2*i+0];
+ ilist->iatoms[4*i+1] = ilist->iatoms[2*i+1];
+ ilist->iatoms[4*i+2] = ilist->iatoms[2*i+1] + 1;
+ ilist->iatoms[4*i+3] = ilist->iatoms[2*i+1] + 2;
+ }
+ ilist->nr = 2*ilist->nr;
+}
+
static void do_ilists(t_fileio *fio, t_ilist *ilist,gmx_bool bRead,
int file_version)
{
ilist[j].iatoms = NULL;
} else {
do_ilist(fio, &ilist[j],bRead,file_version,j);
+ if (file_version < 78 && j == F_SETTLE && ilist[j].nr > 0)
+ {
+ add_settle_atoms(&ilist[j]);
+ }
}
/*
if (bRead && gmx_debug_at)
else
{
mtop->ffparams.cmap_grid.ngrid = 0;
- mtop->ffparams.cmap_grid.grid_spacing = 0.1;
+ mtop->ffparams.cmap_grid.grid_spacing = 0;
mtop->ffparams.cmap_grid.cmapdata = NULL;
}
gmx_bool TopOnlyOK, int *file_version,
int *file_generation)
{
- char buf[STRLEN];
+ char buf[STRLEN];
+ char file_tag[STRLEN];
gmx_bool bDouble;
int precision;
int fver,fgen;
gmx_fio_setprecision(fio,bDouble);
gmx_fio_do_int(fio,precision);
fver = tpx_version;
+ sprintf(file_tag,"%s",tpx_tag);
fgen = tpx_generation;
}
- /* Check versions! */
- gmx_fio_do_int(fio,fver);
+ /* Check versions! */
+ gmx_fio_do_int(fio,fver);
- if(fver>=26)
- gmx_fio_do_int(fio,fgen);
- else
- fgen=0;
+ if (fver >= 77)
+ {
+ gmx_fio_do_string(fio,file_tag);
+ }
+ if (bRead)
+ {
+ if (fver < 77)
+ {
+ /* Versions before 77 don't have the tag, set it to release */
+ sprintf(file_tag,"%s",TPX_TAG_RELEASE);
+ }
+
+ if (strcmp(file_tag,tpx_tag) != 0)
+ {
+ fprintf(stderr,"Note: file tpx tag '%s', software tpx tag '%s'\n",
+ file_tag,tpx_tag);
+
+ /* We only support reading tpx files with the same tag as the code
+ * or tpx files with the release tag and with lower version number.
+ */
+ if (!(strcmp(file_tag,TPX_TAG_RELEASE) == 0 && fver < tpx_version))
+ {
+ gmx_fatal(FARGS,"tpx tag/version mismatch: reading tpx file (%s) version %d, tag '%s' with program for tpx version %d, tag '%s'",
+ gmx_fio_getname(fio),fver,file_tag,
+ tpx_version,tpx_tag);
+ }
+ }
+ }
+
+ if (fver >= 26)
+ {
+ gmx_fio_do_int(fio,fgen);
+ }
+ else
+ {
+ fgen=0;
+ }
- if(file_version!=NULL)
- *file_version = fver;
- if(file_generation!=NULL)
- *file_generation = fgen;
+ if (file_version != NULL)
+ {
+ *file_version = fver;
+ }
+ if (file_generation != NULL)
+ {
+ *file_generation = fgen;
+ }
if ((fver <= tpx_incompatible_version) ||
{
t_tpxheader header;
int natoms,i,version,generation;
- gmx_bool bTop,bXNULL;
+ gmx_bool bTop,bXNULL=FALSE;
gmx_mtop_t *mtop;
t_topology *topconv;
gmx_atomprop_t aps;
else {
get_stx_coordnum(infile,&natoms);
init_t_atoms(&top->atoms,natoms,(fn2ftp(infile) == efPDB));
- bXNULL = (x == NULL);
+ if (x == NULL)
+ {
+ snew(x,1);
+ bXNULL = TRUE;
+ }
snew(*x,natoms);
if (v)
snew(*v,natoms);
read_stx_conf(infile,title,&top->atoms,*x,(v==NULL) ? NULL : *v,ePBC,box);
- if (bXNULL) {
+ if (bXNULL)
+ {
sfree(*x);
- x = NULL;
+ sfree(x);
}
if (bMass) {
aps = gmx_atomprop_init();
#include <config.h>
#endif
-#if ((defined WIN32 || defined _WIN32 || defined WIN64 || defined _WIN64) && !defined __CYGWIN__ && !defined __CYGWIN32__)
-/* _isnan() */
-#include <float.h>
-#endif
-
#include "typedefs.h"
#include "smalloc.h"
#include "sysstuff.h"
#include "disre.h"
#include "orires.h"
#include "dihre.h"
-#include "pppm.h"
#include "pme.h"
#include "mdatoms.h"
#include "repl_ex.h"
#include "mtop_util.h"
#include "sighandler.h"
#include "string2.h"
+#include "membed.h"
#ifdef GMX_LIB_MPI
#include <mpi.h>
#endif
-#ifdef GMX_THREADS
+#ifdef GMX_THREAD_MPI
#include "tmpi.h"
#endif
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,
+ int repl_ex_nst,int repl_ex_seed,gmx_membed_t membed,
real cpt_period,real max_hours,
const char *deviceOptions,
unsigned long Flags,
double t,t0,lam0;
gmx_bool bGStatEveryStep,bGStat,bNstEner,bCalcEnerPres;
gmx_bool bNS,bNStList,bSimAnn,bStopCM,bRerunMD,bNotLastFrame=FALSE,
- bFirstStep,bStateFromTPX,bInitStep,bLastStep,
+ bFirstStep,bStateFromCP,bStateFromTPX,bInitStep,bLastStep,
bBornRadii,bStartingFromCpt;
gmx_bool bDoDHDL=FALSE;
gmx_bool do_ene,do_log,do_verbose,bRerunWarnNoV=TRUE,
gmx_bool bAppend;
gmx_bool bResetCountersHalfMaxH=FALSE;
gmx_bool bVV,bIterations,bFirstIterate,bTemp,bPres,bTrotter;
- real temp0,mu_aver=0,dvdl;
+ real mu_aver=0,dvdl;
int a0,a1,gnx=0,ii;
atom_id *grpindex=NULL;
char *grpname;
if (DEFORM(*ir))
{
-#ifdef GMX_THREADS
+#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_THREADS
+#ifdef GMX_THREAD_MPI
tMPI_Thread_mutex_unlock(&deform_init_box_mutex);
#endif
}
update_mdatoms(mdatoms,state->lambda);
+ if (opt2bSet("-cpi",nfile,fnm))
+ {
+ bStateFromCP = gmx_fexist_master(opt2fn_master("-cpi",nfile,fnm,cr),cr);
+ }
+ else
+ {
+ bStateFromCP = FALSE;
+ }
+
if (MASTER(cr))
{
- if (opt2bSet("-cpi",nfile,fnm))
+ if (bStateFromCP)
{
/* Update mdebin with energy history if appending to output files */
if ( Flags & MD_APPENDFILES )
enerd->term[F_TEMP] *= 2; /* result of averages being done over previous and current step,
and there is no previous step */
}
- temp0 = enerd->term[F_TEMP];
/* if using an iterative algorithm, we need to create a working directory for the state. */
if (bIterations)
"RMS relative constraint deviation after constraining: %.2e\n",
constr_rmsd(constr,FALSE));
}
- fprintf(fplog,"Initial temperature: %g K\n",enerd->term[F_TEMP]);
+ 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"
/* 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);
+ bStateFromTPX = !bStateFromCP;
bInitStep = bFirstStep && (bStateFromTPX || bVV);
bStartingFromCpt = (Flags & MD_STARTFROMCPT) && bInitStep;
bLastStep = FALSE;
{
if (update_forcefield(fplog,
nfile,fnm,fr,
- mdatoms->nr,state->x,state->box)) {
- if (gmx_parallel_env_initialized())
- {
- gmx_finalize();
- }
+ mdatoms->nr,state->x,state->box))
+ {
+ gmx_finalize_par();
+
exit(0);
}
}
/* Check whether everything is still allright */
if (((int)gmx_get_stop_condition() > handled_stop_condition)
-#ifdef GMX_THREADS
+#ifdef GMX_THREAD_MPI
&& MASTER(cr)
#endif
)
f,NULL,xcopy,
&(top_global->mols),mdatoms->massT,pres))
{
- if (gmx_parallel_env_initialized())
- {
- gmx_finalize();
- }
+ gmx_finalize_par();
+
fprintf(stderr,"\n");
exit(0);
}
}
/* ####### END SET VARIABLES FOR NEXT ITERATION ###### */
-
+
+ if ( (membed!=NULL) && (!bLastStep) )
+ {
+ rescale_membed(step_rel,membed,state_global->x);
+ }
+
if (bRerunMD)
{
if (MASTER(cr))
# Files called xxx_test.c are test drivers with a main() function for
# module xxx.c, so they should not be included in the library
-file(GLOB_RECURSE NOT_MDLIB_SOURCES *_test.c *\#*)
-list(REMOVE_ITEM MDLIB_SOURCES ${NOT_MDLIB_SOURCES})
+ if(NOT GMX_FFT_FFTPACK)
+ list(REMOVE_ITEM MDLIB_SOURCES ${CMAKE_CURRENT_SOURCE_DIR}/fftpack.c)
+ endif()
+
add_library(md ${MDLIB_SOURCES})
target_link_libraries(md gmx ${GMX_EXTRA_LIBRARIES} ${FFT_LIBRARIES} ${XML_LIBRARIES})
set_target_properties(md PROPERTIES OUTPUT_NAME "md${GMX_LIBS_SUFFIX}" SOVERSION ${SOVERSION} INSTALL_NAME_DIR "${LIB_INSTALL_DIR}")
#define FFTWPREFIX(name) fftwf_ ## name
#endif
-#ifdef GMX_THREADS
++#ifdef GMX_THREAD_MPI
+ #include "thread_mpi/threads.h"
+ #endif
-#ifdef GMX_THREADS
+#ifdef GMX_THREAD_MPI
/* none of the fftw3 calls, except execute(), are thread-safe, so
we need to serialize them with this mutex. */
static tMPI_Thread_mutex_t big_fftw_mutex=TMPI_THREAD_MUTEX_INITIALIZER;
static gmx_bool gmx_fft_threads_initialized=FALSE;
-#define FFTW_LOCK tMPI_Thread_mutex_lock(&big_fftw_mutex);
-#define FFTW_UNLOCK tMPI_Thread_mutex_unlock(&big_fftw_mutex);
-#else /* GMX_THREADS */
+#define FFTW_LOCK tMPI_Thread_mutex_lock(&big_fftw_mutex)
+#define FFTW_UNLOCK tMPI_Thread_mutex_unlock(&big_fftw_mutex)
+#else /* GMX_THREAD_MPI */
#define FFTW_LOCK
#define FFTW_UNLOCK
-#endif /* GMX_THREADS */
+#endif /* GMX_THREAD_MPI */
/* We assume here that aligned memory starts at multiple of 16 bytes and unaligned memory starts at multiple of 8 bytes. The later is guranteed for all malloc implementation.
Consequesences:
if (top)
gpbc = gmx_rmpbc_init(&top->idef,ePBC,nat,box);
- if (top)
- gmx_rmpbc_done(gpbc);
for(i=0; i<nat; i++)
all_at[i]=i;
"computed based on the Quasiharmonic approach and based on",
"Schlitter's formula."
};
- static int first=1,last=8,skip=1,nextr=2,nskip=6;
+ static int first=1,last=-1,skip=1,nextr=2,nskip=6;
static real max=0.0,temp=298.15;
static gmx_bool bSplit=FALSE,bEntropy=FALSE;
t_pargs pa[] = {
#include <ctype.h>
#include "sysstuff.h"
-#include "string.h"
+#include <string.h>
#include "string2.h"
#include "typedefs.h"
#include "smalloc.h"
}
void calc_electron_density(const char *fn, atom_id **index, int gnx[],
- real ***slDensity, int *nslices, t_topology *top,
+ double ***slDensity, int *nslices, t_topology *top,
int ePBC,
int axis, int nr_grps, real *slWidth,
t_electron eltab[], int nr,gmx_bool bCenter,
}
void calc_density(const char *fn, atom_id **index, int gnx[],
- real ***slDensity, int *nslices, t_topology *top, int ePBC,
+ double ***slDensity, int *nslices, t_topology *top, int ePBC,
int axis, int nr_grps, real *slWidth, gmx_bool bCenter,
const output_env_t oenv)
{
sfree(x0); /* free memory used by coordinate array */
}
- void plot_density(real *slDensity[], const char *afile, int nslices,
+ void plot_density(double *slDensity[], const char *afile, int nslices,
int nr_grps, char *grpname[], real slWidth,
const char **dens_opt,
gmx_bool bSymmetrize, const output_env_t oenv)
"When calculating electron densities, atomnames are used instead of types. This is bad.",
};
- real **density; /* density per slice */
+ double **density; /* density per slice */
real slWidth; /* width of one slice */
char **grpname; /* groupnames */
int nr_electrons; /* nr. electrons */
int nsatm;
t_simat sat[3];
} t_simlist;
-static const char *pdbtp[epdbNR] =
- { "ATOM ", "HETATM" };
real calc_mass(t_atoms *atoms, gmx_bool bGetMass, gmx_atomprop_t aps)
{
"-rvdw", FALSE, etREAL,
{ &rvdw },
"Default Van der Waals radius (in nm) if one can not be found in the database or if no parameters are present in the topology file" },
- { "-sig56", FALSE, etREAL,
+ { "-sig56", FALSE, etBOOL,
{ &bSig56 },
"Use rmin/2 (minimum in the Van der Waals potential) rather than [GRK]sigma[grk]/2 " },
{
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff