From 0860b3cc8f617063dc97fdbea85d7167f8d01e09 Mon Sep 17 00:00:00 2001 From: Michael Shirts Date: Sat, 14 Nov 2009 01:26:42 -0500 Subject: [PATCH] changed setup of trotter decomposition to make it easier to change the decomposition. Also fixed a variable initialization problem that just manifested itself now. --- include/types/enums.h | 5 + include/types/group.h | 1 + include/update.h | 13 ++- src/kernel/md.c | 20 ++-- src/mdlib/coupling.c | 246 +++++++++++++++++++++++++++++++++++++++++++++++++- src/mdlib/tgroup.c | 37 ++++---- src/mdlib/update.c | 6 +- 7 files changed, 297 insertions(+), 31 deletions(-) diff --git a/include/types/enums.h b/include/types/enums.h index a69ab6f617..25fa7a43b7 100644 --- a/include/types/enums.h +++ b/include/types/enums.h @@ -50,6 +50,11 @@ enum { epcNO, epcBERENDSEN, epcPARRINELLORAHMAN, epcISOTROPIC, epcMTTK, epcNR }; /* isotropic is an alias for berendsen */ +/* trotter decomposition extended variable parts */ +enum { + etrtNONE, etrtNHC, etrtBAROSTATV, etrtBAROSTATNHC, etrtNR +}; + enum { epctISOTROPIC, epctSEMIISOTROPIC, epctANISOTROPIC, epctSURFACETENSION, epctNR diff --git a/include/types/group.h b/include/types/group.h index 181466c40b..f95574ec6e 100644 --- a/include/types/group.h +++ b/include/types/group.h @@ -49,6 +49,7 @@ typedef struct { typedef struct { int nat; /* Number of atoms in this group */ rvec u; /* Mean velocities of home particles */ + rvec uold; /* Previous mean velocities of home particles */ double mA; /* Mass for topology A */ double mB; /* Mass for topology B */ } t_grp_acc; diff --git a/include/update.h b/include/update.h index 9e9e81840d..fc522eb2df 100644 --- a/include/update.h +++ b/include/update.h @@ -178,9 +178,20 @@ extern void nosehoover_tcoupl(t_grpopts *opts,gmx_ekindata_t *ekind,real dt, extern t_state *init_bufstate(int size, int ntc); -extern void trotter_update(t_inputrec *ir,gmx_ekindata_t *ekind, gmx_enerdata_t *enerd, + +/* +extern trotter_update(t_inputrec *ir,gmx_ekindata_t *ekind, gmx_enerdata_t *enerd, t_state *state, tensor vir, t_mdatoms *md, t_extmass *MassQ, bool bFirstHalf, bool bThermo, bool bBaro, bool bInitStep); +*/ + + +extern void trotter_update(t_inputrec *ir,gmx_ekindata_t *ekind, gmx_enerdata_t *enerd, + t_state *state, tensor vir, t_mdatoms *md, + t_extmass *MassQ, int *trotter_seq); + +extern int **init_trotter(t_inputrec *ir, t_state *state, t_extmass *Mass); + extern real NPT_energy(t_inputrec *ir, double *xi, double *vxi, real veta, tensor box, t_extmass *MassQ); /* computes all the pressure/tempertature control energy terms to get a conserved energy */ diff --git a/src/kernel/md.c b/src/kernel/md.c index dbfdcebca7..6f506d640b 100644 --- a/src/kernel/md.c +++ b/src/kernel/md.c @@ -1047,7 +1047,7 @@ static bool done_iterating(const t_commrec *cr,FILE *fplog, bool *bFirstIterate, /* Definitions for convergence. Could be optimized . . . */ #ifdef GMX_DOUBLE -#define CONVERGEITER 0.00000001 +#define CONVERGEITER 0.000000001 #else #define CONVERGEITER 0.0001 #endif @@ -1496,7 +1496,7 @@ double do_md(FILE *fplog,t_commrec *cr,int nfile,const t_filenm fnm[], real last_ekin = 0; int iter_i; t_extmass MassQ; - + int **trotter_seq; char sbuf[22],sbuf2[22]; bool bHandledSignal=FALSE; #ifdef GMX_FAHCORE @@ -1781,7 +1781,8 @@ double do_md(FILE *fplog,t_commrec *cr,int nfile,const t_filenm fnm[], if (bTrotter) { /* need to make an initial call to get the Trotter variables set. */ - trotter_update(ir,ekind,enerd,state,total_vir,mdatoms,&MassQ,FALSE,FALSE,FALSE,FALSE); + //trotter_update(ir,ekind,enerd,state,total_vir,mdatoms,&MassQ,FALSE,FALSE,FALSE,FALSE); + trotter_seq = init_trotter(ir,state,&MassQ); } if (MASTER(cr)) @@ -2271,7 +2272,8 @@ double do_md(FILE *fplog,t_commrec *cr,int nfile,const t_filenm fnm[], of veta. */ veta_save = state->veta; - trotter_update(ir,ekind,enerd,state,total_vir,mdatoms,&MassQ,TRUE,FALSE,TRUE,bInitStep); + //trotter_update(ir,ekind,enerd,state,total_vir,mdatoms,&MassQ,TRUE,FALSE,TRUE,bInitStep); + trotter_update(ir,ekind,enerd,state,total_vir,mdatoms,&MassQ,trotter_seq[0]); vetanew = state->veta; state->veta = veta_save; } @@ -2319,9 +2321,10 @@ double do_md(FILE *fplog,t_commrec *cr,int nfile,const t_filenm fnm[], ); /* temperature scaling and pressure scaling to produce the extended variables at t+dt */ - if (bTrotter) + if (bTrotter && !bInitStep) { - trotter_update(ir,ekind,enerd,state,total_vir,mdatoms,&MassQ,TRUE,TRUE,TRUE,bInitStep); + //trotter_update(ir,ekind,enerd,state,total_vir,mdatoms,&MassQ,TRUE,TRUE,TRUE,bInitStep); + trotter_update(ir,ekind,enerd,state,total_vir,mdatoms,&MassQ,trotter_seq[1]); } if (done_iterating(cr,fplog,&bFirstIterate,&bIterate,state->veta,&vetanew,0)) break; @@ -2539,7 +2542,8 @@ double do_md(FILE *fplog,t_commrec *cr,int nfile,const t_filenm fnm[], if ((!bFirstStep) || tracevir != 0) { /* don't apply the trotter update if it's the first step. */ - trotter_update(ir,ekind,enerd,state,total_vir,mdatoms,&MassQ,FALSE,TRUE,TRUE,bInitStep); + //trotter_update(ir,ekind,enerd,state,total_vir,mdatoms,&MassQ,FALSE,TRUE,TRUE,bInitStep); + trotter_update(ir,ekind,enerd,state,total_vir,mdatoms,&MassQ,trotter_seq[2]); } } /* We can only do Berendsen coupling after we have summed the kinetic @@ -2613,7 +2617,7 @@ double do_md(FILE *fplog,t_commrec *cr,int nfile,const t_filenm fnm[], } update_box(fplog,step,ir,mdatoms,state,graph,f, - scale_tot,pcoupl_mu,nrnb,wcycle,upd,bInitStep,FALSE); + ir->nstlist==-1 ? &nlh.scale_tot : NULL,pcoupl_mu,nrnb,wcycle,upd,bInitStep,FALSE); /* ################# END UPDATE STEP 2 ################# */ diff --git a/src/mdlib/coupling.c b/src/mdlib/coupling.c index 707f31ef0c..7d701e8a2f 100644 --- a/src/mdlib/coupling.c +++ b/src/mdlib/coupling.c @@ -48,6 +48,9 @@ #include "nrnb.h" #include "gmx_random.h" +#define NTROTTERCALLS 5 +#define NTROTTERPARTS 3 + /* these integration routines are only references inside this file */ static void NHC_trotter(t_grpopts *opts,gmx_ekindata_t *ekind,real dtfull, double xi[],double vxi[], double scalefac[], real *veta, t_extmass *MassQ) @@ -631,6 +634,7 @@ t_state *init_bufstate(int size, int ntc) return state; } +#if 0 void trotter_update(t_inputrec *ir,gmx_ekindata_t *ekind, gmx_enerdata_t *enerd, t_state *state, tensor vir, t_mdatoms *md, t_extmass *MassQ, bool bFirstHalf, @@ -749,9 +753,9 @@ void trotter_update(t_inputrec *ir,gmx_ekindata_t *ekind, } return; /* first time, we just initialize, not do anything else */ } - + scalefac = ir->opts.vscale_nhc; - + for (i=0; iopts); /* just for ease of referencing */ + ngtc = opts->ngtc; + scalefac = ir->opts.vscale_nhc; + for (i=0; iveta),ir->delta_t,state->box,ekind->ekin,vir, + enerd->term[F_PDISPCORR],enerd->term[F_DISPCORR],MassQ); + break; + case etrtBAROSTATNHC: + NHC_trotter(&(ir->opts),ekind,ir->delta_t,state->nosehoover_xi, + state->nosehoover_vxi,NULL,&(state->veta),MassQ); + break; + case etrtNHC: + NHC_trotter(opts,ekind,ir->delta_t,state->nosehoover_xi, + state->nosehoover_vxi,scalefac,NULL,MassQ); + /* 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]; + msmul(tcstat->ekinh, scalefac[i]*scalefac[i], tcstat->ekinh); + msmul(tcstat->ekin, scalefac[i]*scalefac[i], tcstat->ekin); + } + /* now that we've scaled the groupwise velocities, we can add them up to get the total */ + sum_ekin(FALSE,opts,ekind,&(enerd->term[F_DKDL]),TRUE); + + /* 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; + } + } +#if 0 + /* check for conserved momentum */ + if (debug) + { + if (bFirstHalf) + { + for (d=0;dnrdf[0])*((1.0/DIM)*log(det(state->box)/state->vol0)) + state->nosehoover_xi[0]); + } + fprintf(debug,"Consk: %15.8f %15.8f %15.8f\n",consk[0],consk[1],consk[2]); + } + } +#endif +} + +int **init_trotter(t_inputrec *ir, t_state *state, t_extmass *MassQ) +{ + int n,i,j,d,ntgrp,ngtc,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; + + snew(trotter_seq,NTROTTERCALLS); + /* first, initialize clear all the trotter calls */ + for (i=0;iopts); /* just for ease of referencing */ + ngtc = opts->ngtc; + + /* 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*/ + } + } + + /* now, allocate space for temperature variables */ + snew(MassQ->Qinv,(ngtc+1)*NNHCHAIN); + /* 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*NNHCHAIN+j] = 1.0/(sqr(opts->tau_t[i]/M_2PI)*ndj*kT); + } + } + else + { + reft=0.0; + for (j=0;jQinv[i*NNHCHAIN+j] = 0.0; + } + } + } + + switch (ir->epct) + { + case epctISOTROPIC: + default: + bmass = DIM*DIM; /* recommended mass parameters for isotropic barostat */ + } + + /* barostat temperature */ + i = ngtc; /* barostat temperature control variables are one after the last T-group */ + + if ((ir->tau_p > 0) && (opts->ref_t[0] > 0)) + { + reft = max(0.0,opts->ref_t[0]); + kT = BOLTZ*reft; + for (j=0;jQinv[i*NNHCHAIN+j] = 1.0/(sqr(opts->tau_t[0]/M_2PI)*qmass*kT); + } + } + else + { + for (j=0;jQinv[i*NNHCHAIN+j]=0.0; + } + } + return trotter_seq; +} real NPT_energy(t_inputrec *ir, double *xi, double *vxi, real veta, tensor box, t_extmass *MassQ) { diff --git a/src/mdlib/tgroup.c b/src/mdlib/tgroup.c index 615a8ce6c7..72372a38d2 100644 --- a/src/mdlib/tgroup.c +++ b/src/mdlib/tgroup.c @@ -64,24 +64,27 @@ static void init_grptcstat(int ngtc,t_grp_tcstat tcstat[]) static void init_grpstat(FILE *log, gmx_mtop_t *mtop,int ngacc,t_grp_acc gstat[]) { - gmx_groups_t *groups; - gmx_mtop_atomloop_all_t aloop; - int i,grp; - t_atom *atom; - - if (ngacc > 0) { - groups = &mtop->groups; - aloop = gmx_mtop_atomloop_all_init(mtop); - while (gmx_mtop_atomloop_all_next(aloop,&i,&atom)) { - grp = ggrpnr(groups,egcACC,i); - if ((grp < 0) && (grp >= ngacc)) - gmx_incons("Input for acceleration groups wrong"); - gstat[grp].nat++; - /* This will not work for integrator BD */ - gstat[grp].mA += atom->m; - gstat[grp].mB += atom->mB; + gmx_groups_t *groups; + gmx_mtop_atomloop_all_t aloop; + int i,grp; + t_atom *atom; + + if (ngacc > 0) { + groups = &mtop->groups; + aloop = gmx_mtop_atomloop_all_init(mtop); + while (gmx_mtop_atomloop_all_next(aloop,&i,&atom)) + { + grp = ggrpnr(groups,egcACC,i); + if ((grp < 0) && (grp >= ngacc)) + { + gmx_incons("Input for acceleration groups wrong"); + } + gstat[grp].nat++; + /* This will not work for integrator BD */ + gstat[grp].mA += atom->m; + gstat[grp].mB += atom->mB; + } } - } } void init_ekindata(FILE *log,gmx_mtop_t *mtop,t_grpopts *opts, diff --git a/src/mdlib/update.c b/src/mdlib/update.c index baac7592dd..14a2f1744b 100644 --- a/src/mdlib/update.c +++ b/src/mdlib/update.c @@ -787,8 +787,8 @@ static void dump_it_all(FILE *fp,const char *title, #endif } -static void calc_ke_part_normal(rvec v[],t_grpopts *opts,t_mdatoms *md, - gmx_ekindata_t *ekind,t_nrnb *nrnb,bool bFullStepV) +static void calc_ke_part_normal(rvec v[], t_grpopts *opts,t_mdatoms *md, + gmx_ekindata_t *ekind,t_nrnb *nrnb,bool bEkinAveVel) { int start=md->start,homenr=md->homenr; int g,d,n,m,ga=0,gt=0; @@ -832,7 +832,7 @@ static void calc_ke_part_normal(rvec v[],t_grpopts *opts,t_mdatoms *md, for (m=0;(m