From a75db633a2e363d48c66a440d2890592f4431d52 Mon Sep 17 00:00:00 2001 From: Berk Hess Date: Wed, 14 May 2014 22:36:29 +0200 Subject: [PATCH] Fix constraint virial with multiple time stepping With multiple time stepping the additional nstcalclr-1 force contribution was constrained to remove it from the virial. This procedure neglected the non-linear contribution due to the rotation of constraints. Now the contribution of this force component to the coordinate update is constrained instead and the corresponding virial contribution is subtracted from the constraint virial. Fixes #1400 Change-Id: If3217f52808bf7491998324f8dc3161bc003ec1b --- include/types/forcerec.h | 2 + include/update.h | 1 + src/kernel/md.c | 20 +++++-- src/mdlib/update.c | 118 ++++++++++++++++++++++++--------------- 4 files changed, 92 insertions(+), 49 deletions(-) diff --git a/include/types/forcerec.h b/include/types/forcerec.h index 920a81db54..7e14b7d5d4 100644 --- a/include/types/forcerec.h +++ b/include/types/forcerec.h @@ -332,6 +332,8 @@ typedef struct { 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 diff --git a/include/update.h b/include/update.h index 8fde948258..1d53f4a23a 100644 --- a/include/update.h +++ b/include/update.h @@ -106,6 +106,7 @@ void update_coords(FILE *fplog, 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, diff --git a/src/kernel/md.c b/src/kernel/md.c index 6b864b7a9e..da68da10f2 100644 --- a/src/kernel/md.c +++ b/src/kernel/md.c @@ -1237,7 +1237,7 @@ double do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[], 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, wcycle, upd, bInitStep, etrtVELOCITY1, cr, nrnb, constr, &top->idef); @@ -1288,6 +1288,12 @@ double do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[], cr, nrnb, wcycle, upd, constr, bInitStep, 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 && !bFFscan) { gmx_fatal(FARGS, "Constraint error: Shake, Lincs or Settle could not solve the constrains"); @@ -1739,7 +1745,7 @@ double do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[], /* 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, wcycle, upd, FALSE, etrtVELOCITY2, cr, nrnb, constr, &top->idef); } @@ -1756,7 +1762,7 @@ double do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[], 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, wcycle, upd, bInitStep, etrtPOSITION, cr, nrnb, constr, &top->idef); wallcycle_stop(wcycle, ewcUPDATE); @@ -1766,6 +1772,12 @@ double do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[], cr, nrnb, wcycle, upd, constr, bInitStep, 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? */ @@ -1784,7 +1796,7 @@ double do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[], 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, wcycle, upd, bInitStep, etrtPOSITION, cr, nrnb, constr, &top->idef); wallcycle_stop(wcycle, ewcUPDATE); diff --git a/src/mdlib/update.c b/src/mdlib/update.c index d1f851fbfc..d19d6015bb 100644 --- a/src/mdlib/update.c +++ b/src/mdlib/update.c @@ -1333,49 +1333,6 @@ static void deform(gmx_update_t upd, } } -static void combine_forces(int nstcalclr, - gmx_constr_t constr, - t_inputrec *ir, t_mdatoms *md, t_idef *idef, - t_commrec *cr, - gmx_large_int_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(FILE *fplog, gmx_large_int_t step, t_inputrec *inputrec, @@ -1519,6 +1476,75 @@ static rvec *get_xprime(const t_state *state, gmx_update_t upd) 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_large_int_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 = md->start; i < md->start + 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_large_int_t step, real *dvdlambda, /* the contribution to be added to the bonded interactions */ @@ -1837,6 +1863,7 @@ void update_coords(FILE *fplog, 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, @@ -1898,9 +1925,10 @@ void update_coords(FILE *fplog, * 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 -- 2.22.0