- 402
-413
++414
If You Want Something Done You Have to Do It Yourself_(Highlander II)
I Live the Life They Wish They Did_(Tricky)
Jesus Built My Hotrod_(Ministry)
It was something to at least have a choice of nightmares_(Joseph Conrad)
You fight, work, sweat, nearly kill yourself, sometimes you do kill yourself, trying to accomplish something - and you can't._(Joseph Conrad)
And after some more talk we agreed that the wisdom of rats had been grossly overrated, being in fact no greater than that of men_(Joseph Conrad)
+It's an easy game, just don't let the ball past!_(Szilard Pall)
The soul? There's nothing but chemistry here_(Breaking Bad)
You got one part of that wrong. This is not meth._(Breaking Bad)
It's easy to remember: a half a kT is equal to five fourths of a kJ/mol._(Anders Gabrielsson)
-
+ Ubiquitin's just a rock_(Berk Hess)
+ ... an excellent man, almost worthy of such a wife ..._(Jane Eyre in Jane Eyre by Charlotte Bronte)
+ Humbug! Most things free-born will submit to anything for a salary_(Mr. Rochester in Jane Eyre by Charlotte Bronte)
+ Like other defaulters, I like to lay half the blame on ill-fortune and adverse circumstances_(Mr. Rochester in Jane Eyre by Charlotte Bronte)
+ Either you will be dashed to atoms on crag points, or lifted up and borne by some master-wave into a calmer current_(Charlotte Bronte)
+ I ought to warn you, I have no faith_(Jane Eyre in Jane Eyre by Charlotte Bronte)
+ ... yet the [economic] profession continued to churn out purely theoretical results without even knowing what facts needed to be explained._(Thomas Piketty)
+ Scientists think they are born with logic; God forbid they should study this discipline with a history of more than two and a half millennia._(Roald Hoffmann)
+ In the processing of models we must be especially cautious of the human weakness to think that models can be verified or validated. Especially one's own._(Roald Hoffmann)
+ ... and that dream of dreams, a computational model that predicts everything accurately._(Roald Hoffmann)
+ You see it through a charmed medium: you can not discern that the gilding is slime and the silk draperies cobwebs; that the marble is sordid slate, and the polished woods mere refuse chips and scale bark._(Mr. Rochester in Jane Eyre by Charlotte Bronte)
+ I know poetry is not dead, nor genius lost; nor has Mammon gained power over either, to bind or slay; they will both assert their existence, their presence, their liberty and strength again one day._(Jane Eyre in Jane Eyre by Charlotte Bronte)
#endif
#include <math.h>
+#include <stdlib.h>
#include <string.h>
#include "typedefs.h"
-#include "gmx_fatal.h"
-#include "vec.h"
+#include "gromacs/utility/fatalerror.h"
+#include "gromacs/math/vec.h"
#include "gromacs/utility/cstringutil.h"
#include "gromacs/utility/smalloc.h"
#include "gromacs/fileio/enxio.h"
#include "names.h"
#include "copyrite.h"
#include "macros.h"
-#include "xvgr.h"
+#include "gromacs/fileio/xvgr.h"
#include "gstat.h"
-#include "physics.h"
+#include "gromacs/math/units.h"
#include "gromacs/fileio/tpxio.h"
#include "gromacs/fileio/trxio.h"
#include "viewit.h"
-#include "mtop_util.h"
+#include "gromacs/topology/mtop_util.h"
#include "gmx_ana.h"
#include "mdebin.h"
}
static void einstein_visco(const char *fn, const char *fni, int nsets,
- int nframes, real **sum,
- real V, real T, int nsteps, double time[],
+ int nint, real **eneint,
+ real V, real T, double dt,
const output_env_t oenv)
{
FILE *fp0, *fp1;
double av[4], avold[4];
- double fac, dt, di;
+ double fac, di;
int i, j, m, nf4;
- if (nframes < 1)
- {
- return;
- }
-
- dt = (time[1]-time[0]);
- nf4 = nframes/4+1;
+ nf4 = nint/4 + 1;
for (i = 0; i <= nsets; i++)
{
"Time (ps)", "(kg m\\S-1\\N s\\S-1\\N ps)", oenv);
fp1 = xvgropen(fn, "Shear viscosity using Einstein relation",
"Time (ps)", "(kg m\\S-1\\N s\\S-1\\N)", oenv);
- for (i = 1; i < nf4; i++)
+ for (i = 0; i < nf4; i++)
{
- fac = dt*nframes/nsteps;
for (m = 0; m <= nsets; m++)
{
av[m] = 0;
}
- for (j = 0; j < nframes-i; j++)
+ for (j = 0; j < nint - i; j++)
{
for (m = 0; m < nsets; m++)
{
- di = sqr(fac*(sum[m][j+i]-sum[m][j]));
+ di = sqr(eneint[m][j+i] - eneint[m][j]);
av[m] += di;
av[nsets] += di/nsets;
}
}
/* Convert to SI for the viscosity */
- fac = (V*NANO*NANO*NANO*PICO*1e10)/(2*BOLTZMANN*T)/(nframes-i);
- fprintf(fp0, "%10g", time[i]-time[0]);
+ fac = (V*NANO*NANO*NANO*PICO*1e10)/(2*BOLTZMANN*T)/(nint - i);
+ fprintf(fp0, "%10g", i*dt);
for (m = 0; (m <= nsets); m++)
{
av[m] = fac*av[m];
fprintf(fp0, " %10g", av[m]);
}
fprintf(fp0, "\n");
- fprintf(fp1, "%10g", 0.5*(time[i]+time[i-1])-time[0]);
+ fprintf(fp1, "%10g", (i + 0.5)*dt);
for (m = 0; (m <= nsets); m++)
{
fprintf(fp1, " %10g", (av[m]-avold[m])/dt);
gmx_bool bVisco, const char *visfn, int nmol,
gmx_int64_t start_step, double start_t,
gmx_int64_t step, double t,
- double time[], real reftemp,
+ real reftemp,
enerdata_t *edat,
int nset, int set[], gmx_bool *bIsEner,
char **leg, gmx_enxnm_t *enm,
const char* leg[] = { "Shear", "Bulk" };
real factor;
real **eneset;
- real **enesum;
+ real **eneint;
/* Assume pressure tensor is in Pxx Pxy Pxz Pyx Pyy Pyz Pzx Pzy Pzz */
{
snew(eneset[i], edat->nframes);
}
- snew(enesum, 3);
- for (i = 0; i < 3; i++)
- {
- snew(enesum[i], edat->nframes);
- }
for (i = 0; (i < edat->nframes); i++)
{
eneset[0][i] = 0.5*(edat->s[1].ener[i]+edat->s[3].ener[i]);
eneset[j][i] = edat->s[j].ener[i];
}
eneset[11][i] -= Pres;
- enesum[0][i] = 0.5*(edat->s[1].es[i].sum+edat->s[3].es[i].sum);
- enesum[1][i] = 0.5*(edat->s[2].es[i].sum+edat->s[6].es[i].sum);
- enesum[2][i] = 0.5*(edat->s[5].es[i].sum+edat->s[7].es[i].sum);
+ }
+
+ /* Determine integrals of the off-diagonal pressure elements */
+ snew(eneint, 3);
+ for (i = 0; i < 3; i++)
+ {
+ snew(eneint[i], edat->nframes + 1);
+ }
+ eneint[0][0] = 0;
+ eneint[1][0] = 0;
+ eneint[2][0] = 0;
+ for (i = 0; i < edat->nframes; i++)
+ {
+ eneint[0][i+1] = eneint[0][i] + 0.5*(edat->s[1].es[i].sum + edat->s[3].es[i].sum)*Dt/edat->points[i];
+ eneint[1][i+1] = eneint[1][i] + 0.5*(edat->s[2].es[i].sum + edat->s[6].es[i].sum)*Dt/edat->points[i];
+ eneint[2][i+1] = eneint[2][i] + 0.5*(edat->s[5].es[i].sum + edat->s[7].es[i].sum)*Dt/edat->points[i];
}
einstein_visco("evisco.xvg", "eviscoi.xvg",
- 3, edat->nframes, enesum, Vaver, Temp, nsteps, time, oenv);
+ 3, edat->nframes+1, eneint, Vaver, Temp, Dt, oenv);
+
+ for (i = 0; i < 3; i++)
+ {
+ sfree(eneint[i]);
+ }
+ sfree(eneint);
/*do_autocorr(corrfn,buf,nenergy,3,eneset,Dt,eacNormal,TRUE);*/
/* Do it for shear viscosity */
fprintf(fp, "%10g %10g %10g\n", (i*Dt), integral, intBulk);
}
gmx_ffclose(fp);
+
+ for (i = 0; i < 12; i++)
+ {
+ sfree(eneset[i]);
+ }
+ sfree(eneset);
}
else if (bCorr)
{
bVisco, opt2fn("-vis", NFILE, fnm),
nmol,
start_step, start_t, frame[cur].step, frame[cur].t,
- time, reftemp, &edat,
+ reftemp, &edat,
nset, set, bIsEner, leg, enm, Vaver, ezero, nbmin, nbmax,
oenv);
if (bFluctProps)
#include "ns.h"
#include "genborn.h"
#include "qmmmrec.h"
-#include "idef.h"
+#include "../../topology/idef.h"
#include "nb_verlet.h"
#include "interaction_const.h"
#include "hw_info.h"
gmx_bool bTwinRange;
int nlr;
rvec *f_twin;
+ /* Constraint virial correction for multiple time stepping */
+ tensor vir_twin_constr;
/* Forces that should not enter into the virial summation:
* PPPM/PME/Ewald/posres
#define _update_h
#include "typedefs.h"
-#include "mshift.h"
#include "tgroup.h"
#include "network.h"
-#include "vec.h"
+#include "../timing/wallcycle.h"
#ifdef __cplusplus
extern "C" {
#endif
+struct t_graph;
+
/* Abstract type for stochastic dynamics */
typedef struct gmx_update *gmx_update_t;
rvec *f, /* forces on home particles */
gmx_bool bDoLR,
rvec *f_lr,
+ tensor *vir_lr_constr,
t_fcdata *fcd,
gmx_ekindata_t *ekind,
matrix M,
t_mdatoms *md,
t_state *state,
gmx_bool bMolPBC,
- t_graph *graph,
+ struct t_graph *graph,
rvec force[], /* forces on home particles */
t_idef *idef,
tensor vir_part,
#include <math.h>
#include "types/commrec.h"
-#include "sysstuff.h"
-#include "gromacs/utility/smalloc.h"
#include "typedefs.h"
#include "nrnb.h"
-#include "physics.h"
+#include "gromacs/math/units.h"
#include "macros.h"
-#include "vec.h"
-#include "main.h"
+#include "gromacs/math/vec.h"
#include "update.h"
#include "gromacs/random/random.h"
-#include "mshift.h"
#include "tgroup.h"
#include "force.h"
#include "names.h"
#include "gmx_omp_nthreads.h"
#include "gromacs/fileio/confio.h"
-#include "gromacs/fileio/futil.h"
+#include "gromacs/pbcutil/mshift.h"
+#include "gromacs/pbcutil/pbc.h"
+#include "gromacs/pulling/pull.h"
#include "gromacs/timing/wallcycle.h"
+#include "gromacs/utility/futil.h"
#include "gromacs/utility/gmxomp.h"
-#include "gromacs/pulling/pull.h"
+#include "gromacs/utility/smalloc.h"
/*For debugging, start at v(-dt/2) for velolcity verlet -- uncomment next line */
/*#define STARTFROMDT2*/
}
}
- static void combine_forces(int nstcalclr,
- gmx_constr_t constr,
- t_inputrec *ir, t_mdatoms *md, t_idef *idef,
- t_commrec *cr,
- gmx_int64_t step,
- t_state *state, gmx_bool bMolPBC,
- int start, int nrend,
- rvec f[], rvec f_lr[],
- t_nrnb *nrnb)
- {
- int i, d, nm1;
-
- /* f contains the short-range forces + the long range forces
- * which are stored separately in f_lr.
- */
-
- if (constr != NULL && !(ir->eConstrAlg == econtSHAKE && ir->epc == epcNO))
- {
- /* We need to constrain the LR forces separately,
- * because due to the different pre-factor for the SR and LR
- * forces in the update algorithm, we can not determine
- * the constraint force for the coordinate constraining.
- * Constrain only the additional LR part of the force.
- */
- /* MRS -- need to make sure this works with trotter integration -- the constraint calls may not be right.*/
- constrain(NULL, FALSE, FALSE, constr, idef, ir, NULL, cr, step, 0, md,
- state->x, f_lr, f_lr, bMolPBC, state->box, state->lambda[efptBONDED], NULL,
- NULL, NULL, nrnb, econqForce, ir->epc == epcMTTK, state->veta, state->veta);
- }
-
- /* Add nstcalclr-1 times the LR force to the sum of both forces
- * and store the result in forces_lr.
- */
- nm1 = nstcalclr - 1;
- for (i = start; i < nrend; i++)
- {
- for (d = 0; d < DIM; d++)
- {
- f_lr[i][d] = f[i][d] + nm1*f_lr[i][d];
- }
- }
- }
-
void update_tcouple(gmx_int64_t step,
t_inputrec *inputrec,
t_state *state,
return upd->xp;
}
+ static void combine_forces(gmx_update_t upd,
+ int nstcalclr,
+ gmx_constr_t constr,
+ t_inputrec *ir, t_mdatoms *md, t_idef *idef,
+ t_commrec *cr,
+ gmx_int64_t step,
+ t_state *state, gmx_bool bMolPBC,
+ int start, int nrend,
+ rvec f[], rvec f_lr[],
+ tensor *vir_lr_constr,
+ t_nrnb *nrnb)
+ {
+ int i, d;
+
+ /* f contains the short-range forces + the long range forces
+ * which are stored separately in f_lr.
+ */
+
+ if (constr != NULL && vir_lr_constr != NULL &&
+ !(ir->eConstrAlg == econtSHAKE && ir->epc == epcNO))
+ {
+ /* We need to constrain the LR forces separately,
+ * because due to the different pre-factor for the SR and LR
+ * forces in the update algorithm, we have to correct
+ * the constraint virial for the nstcalclr-1 extra f_lr.
+ * Constrain only the additional LR part of the force.
+ */
+ /* MRS -- need to make sure this works with trotter integration -- the constraint calls may not be right.*/
+ rvec *xp;
+ real fac;
+ int gf = 0;
+
+ xp = get_xprime(state, upd);
+
+ fac = (nstcalclr - 1)*ir->delta_t*ir->delta_t;
+
+ for (i = 0; i < md->homenr; i++)
+ {
+ if (md->cFREEZE != NULL)
+ {
+ gf = md->cFREEZE[i];
+ }
+ for (d = 0; d < DIM; d++)
+ {
+ if ((md->ptype[i] != eptVSite) &&
+ (md->ptype[i] != eptShell) &&
+ !ir->opts.nFreeze[gf][d])
+ {
+ xp[i][d] = state->x[i][d] + fac*f_lr[i][d]*md->invmass[i];
+ }
+ }
+ }
+ constrain(NULL, FALSE, FALSE, constr, idef, ir, NULL, cr, step, 0, md,
+ state->x, xp, xp, bMolPBC, state->box, state->lambda[efptBONDED], NULL,
+ NULL, vir_lr_constr, nrnb, econqCoord, ir->epc == epcMTTK, state->veta, state->veta);
+ }
+
+ /* Add nstcalclr-1 times the LR force to the sum of both forces
+ * and store the result in forces_lr.
+ */
+ for (i = start; i < nrend; i++)
+ {
+ for (d = 0; d < DIM; d++)
+ {
+ f_lr[i][d] = f[i][d] + (nstcalclr - 1)*f_lr[i][d];
+ }
+ }
+ }
+
void update_constraints(FILE *fplog,
gmx_int64_t step,
real *dvdlambda, /* the contribution to be added to the bonded interactions */
rvec *f, /* forces on home particles */
gmx_bool bDoLR,
rvec *f_lr,
+ tensor *vir_lr_constr,
t_fcdata *fcd,
gmx_ekindata_t *ekind,
matrix M,
* to produce twin time stepping.
*/
/* is this correct in the new construction? MRS */
- combine_forces(inputrec->nstcalclr, constr, inputrec, md, idef, cr,
+ combine_forces(upd,
+ inputrec->nstcalclr, constr, inputrec, md, idef, cr,
step, state, bMolPBC,
- start, nrend, f, f_lr, nrnb);
+ start, nrend, f, f_lr, vir_lr_constr, nrnb);
force = f_lr;
}
else
#include <config.h>
#endif
+#include <stdlib.h>
+
#include "typedefs.h"
-#include "gromacs/utility/smalloc.h"
-#include "sysstuff.h"
-#include "vec.h"
+#include "gromacs/math/vec.h"
#include "vcm.h"
#include "mdebin.h"
#include "nrnb.h"
#include "md_support.h"
#include "md_logging.h"
#include "network.h"
-#include "xvgr.h"
-#include "physics.h"
#include "names.h"
#include "force.h"
#include "disre.h"
#include "qmmm.h"
#include "domdec.h"
#include "domdec_network.h"
-#include "gromacs/gmxlib/topsort.h"
#include "coulomb.h"
#include "constr.h"
#include "shellfc.h"
#include "gromacs/gmxpreprocess/compute_io.h"
#include "checkpoint.h"
-#include "mtop_util.h"
+#include "gromacs/topology/mtop_util.h"
#include "sighandler.h"
#include "txtdump.h"
#include "gromacs/utility/cstringutil.h"
#include "types/iteratedconstraints.h"
#include "nbnxn_cuda_data_mgmt.h"
-#include "gromacs/utility/gmxmpi.h"
#include "gromacs/fileio/confio.h"
+#include "gromacs/fileio/mdoutf.h"
#include "gromacs/fileio/trajectory_writing.h"
#include "gromacs/fileio/trnio.h"
#include "gromacs/fileio/trxio.h"
#include "gromacs/fileio/xtcio.h"
-#include "gromacs/timing/wallcycle.h"
-#include "gromacs/timing/walltime_accounting.h"
+#include "gromacs/imd/imd.h"
+#include "gromacs/pbcutil/mshift.h"
+#include "gromacs/pbcutil/pbc.h"
#include "gromacs/pulling/pull.h"
#include "gromacs/swap/swapcoords.h"
-#include "gromacs/imd/imd.h"
+#include "gromacs/timing/wallcycle.h"
+#include "gromacs/timing/walltime_accounting.h"
+#include "gromacs/utility/gmxmpi.h"
+#include "gromacs/utility/smalloc.h"
#ifdef GMX_FAHCORE
#include "corewrap.h"
rvec mu_tot;
t_vcm *vcm;
t_state *bufstate = NULL;
- matrix *scale_tot, pcoupl_mu, M, ebox;
+ matrix pcoupl_mu, M;
gmx_nlheur_t nlh;
t_trxframe rerun_fr;
gmx_repl_ex_t repl_ex = NULL;
gmx_ekindata_t *ekind, *ekind_save;
gmx_shellfc_t shellfc;
int count, nconverged = 0;
- real timestep = 0;
- double tcount = 0;
- gmx_bool bConverged = TRUE, bOK, bSumEkinhOld, bExchanged, bNeedRepartition;
- gmx_bool bAppend;
+ double tcount = 0;
+ gmx_bool bConverged = TRUE, bOK, bSumEkinhOld, bExchanged, bNeedRepartition;
gmx_bool bResetCountersHalfMaxH = FALSE;
gmx_bool bVV, bIterativeCase, bFirstIterate, bTemp, bPres, bTrotter;
gmx_bool bUpdateDoLR;
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];
/* Check for special mdrun options */
bRerunMD = (Flags & MD_RERUN);
- bAppend = (Flags & MD_APPENDFILES);
if (Flags & MD_RESETCOUNTERSHALFWAY)
{
if (ir->nsteps > 0)
bStateFromTPX = !bStateFromCP;
bInitStep = bFirstStep && (bStateFromTPX || bVV);
bStartingFromCpt = (Flags & MD_STARTFROMCPT) && bInitStep;
- bLastStep = FALSE;
bSumEkinhOld = FALSE;
bExchanged = FALSE;
bNeedRepartition = FALSE;
step = ir->init_step;
step_rel = 0;
- if (ir->nstlist == -1)
- {
- init_nlistheuristics(&nlh, bGStatEveryStep, step);
- }
+ init_nlistheuristics(&nlh, bGStatEveryStep, step);
if (MULTISIM(cr) && (repl_ex_nst <= 0 ))
{
bUpdateDoLR = (fr->bTwinRange && do_per_step(step, ir->nstcalclr));
update_coords(fplog, step, ir, mdatoms, state, fr->bMolPBC,
- f, bUpdateDoLR, fr->f_twin, fcd,
+ f, bUpdateDoLR, fr->f_twin, bCalcVir ? &fr->vir_twin_constr : NULL, fcd,
ekind, M, upd, bInitStep, etrtVELOCITY1,
cr, nrnb, constr, &top->idef);
cr, nrnb, wcycle, upd, constr,
TRUE, bCalcVir, vetanew);
+ if (bCalcVir && bUpdateDoLR && ir->nstcalclr > 1)
+ {
+ /* Correct the virial for multiple time stepping */
+ m_sub(shake_vir, fr->vir_twin_constr, shake_vir);
+ }
+
if (!bOK)
{
gmx_fatal(FARGS, "Constraint error: Shake, Lincs or Settle could not solve the constrains");
/* velocity half-step update */
update_coords(fplog, step, ir, mdatoms, state, fr->bMolPBC, f,
- bUpdateDoLR, fr->f_twin, fcd,
+ bUpdateDoLR, fr->f_twin, bCalcVir ? &fr->vir_twin_constr : NULL, fcd,
ekind, M, upd, FALSE, etrtVELOCITY2,
cr, nrnb, constr, &top->idef);
}
bUpdateDoLR = (fr->bTwinRange && do_per_step(step, ir->nstcalclr));
update_coords(fplog, step, ir, mdatoms, state, fr->bMolPBC, f,
- bUpdateDoLR, fr->f_twin, fcd,
+ bUpdateDoLR, fr->f_twin, bCalcVir ? &fr->vir_twin_constr : NULL, fcd,
ekind, M, upd, bInitStep, etrtPOSITION, cr, nrnb, constr, &top->idef);
wallcycle_stop(wcycle, ewcUPDATE);
cr, nrnb, wcycle, upd, constr,
FALSE, bCalcVir, state->veta);
+ if (bCalcVir && bUpdateDoLR && ir->nstcalclr > 1)
+ {
+ /* Correct the virial for multiple time stepping */
+ m_sub(shake_vir, fr->vir_twin_constr, shake_vir);
+ }
+
if (ir->eI == eiVVAK)
{
/* erase F_EKIN and F_TEMP here? */
bUpdateDoLR = (fr->bTwinRange && do_per_step(step, ir->nstcalclr));
update_coords(fplog, step, ir, mdatoms, state, fr->bMolPBC, f,
- bUpdateDoLR, fr->f_twin, fcd,
+ bUpdateDoLR, fr->f_twin, bCalcVir ? &fr->vir_twin_constr : NULL, fcd,
ekind, M, upd, bInitStep, etrtPOSITION, cr, nrnb, constr, &top->idef);
wallcycle_stop(wcycle, ewcUPDATE);
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff