-/* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
+/*
+ * This file is part of the GROMACS molecular simulation package.
*
- *
- * 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
+ * Copyright (c) 2001-2004, The GROMACS development team.
+ * Copyright (c) 2013,2014, by the GROMACS development team, led by
+ * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+ * and including many others, as listed in the AUTHORS file in the
+ * top-level source directory and at http://www.gromacs.org.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
* 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.
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
*
- * 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.
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
*
- * For more info, check our website at http://www.gromacs.org
+ * If you want to redistribute modifications to GROMACS, 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 http://www.gromacs.org.
*
- * And Hey:
- * GROwing Monsters And Cloning Shrimps
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
*/
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
+#include "gmxpre.h"
+#include "gromacs/legacyheaders/update.h"
-#include <stdio.h>
#include <math.h>
-
-#include "types/commrec.h"
-#include "sysstuff.h"
-#include "smalloc.h"
-#include "typedefs.h"
-#include "nrnb.h"
-#include "physics.h"
-#include "macros.h"
-#include "vec.h"
-#include "main.h"
-#include "update.h"
-#include "gmx_random.h"
-#include "mshift.h"
-#include "tgroup.h"
-#include "force.h"
-#include "names.h"
-#include "txtdump.h"
-#include "mdrun.h"
-#include "constr.h"
-#include "edsam.h"
-#include "pull.h"
-#include "disre.h"
-#include "orires.h"
-#include "gmx_omp_nthreads.h"
+#include <stdio.h>
#include "gromacs/fileio/confio.h"
-#include "gromacs/fileio/futil.h"
+#include "gromacs/legacyheaders/constr.h"
+#include "gromacs/legacyheaders/disre.h"
+#include "gromacs/legacyheaders/force.h"
+#include "gromacs/legacyheaders/gmx_omp_nthreads.h"
+#include "gromacs/legacyheaders/macros.h"
+#include "gromacs/legacyheaders/mdrun.h"
+#include "gromacs/legacyheaders/names.h"
+#include "gromacs/legacyheaders/nrnb.h"
+#include "gromacs/legacyheaders/orires.h"
+#include "gromacs/legacyheaders/tgroup.h"
+#include "gromacs/legacyheaders/txtdump.h"
+#include "gromacs/legacyheaders/typedefs.h"
+#include "gromacs/legacyheaders/types/commrec.h"
+#include "gromacs/math/units.h"
+#include "gromacs/math/vec.h"
+#include "gromacs/pbcutil/mshift.h"
+#include "gromacs/pbcutil/pbc.h"
+#include "gromacs/pulling/pull.h"
+#include "gromacs/random/random.h"
#include "gromacs/timing/wallcycle.h"
+#include "gromacs/utility/futil.h"
#include "gromacs/utility/gmxomp.h"
+#include "gromacs/utility/smalloc.h"
/*For debugging, start at v(-dt/2) for velolcity verlet -- uncomment next line */
/*#define STARTFROMDT2*/
} gmx_sd_sigma_t;
typedef struct {
- /* The random state for ngaussrand threads.
- * Normal thermostats need just 1 random number generator,
- * but SD and BD with OpenMP parallelization need 1 for each thread.
- */
- int ngaussrand;
- gmx_rng_t *gaussrand;
/* BD stuff */
real *bd_rf;
/* SD stuff */
rvec *xp;
int xp_nalloc;
- /* variable size arrays for andersen */
- gmx_bool *randatom;
- int *randatom_list;
- gmx_bool randatom_list_init;
-
/* Variables for the deform algorithm */
- gmx_large_int_t deformref_step;
+ gmx_int64_t deformref_step;
matrix deformref_box;
} t_gmx_update;
}
}
-/* Allocates and initializes sd->gaussrand[i] for i=1, i<sd->ngaussrand,
- * Using seeds generated from sd->gaussrand[0].
- */
-static void init_multiple_gaussrand(gmx_stochd_t *sd)
-{
- int ngr, i;
- unsigned int *seed;
-
- ngr = sd->ngaussrand;
- snew(seed, ngr);
-
- for (i = 1; i < ngr; i++)
- {
- seed[i] = gmx_rng_uniform_uint32(sd->gaussrand[0]);
- }
-
- if (ngr != gmx_omp_nthreads_get(emntUpdate))
- {
- gmx_incons("The number of Gaussian number generators should be equal to gmx_omp_nthreads_get(emntUpdate)");
- }
-
-#pragma omp parallel num_threads(gmx_omp_nthreads_get(emntUpdate))
- {
- int th;
-
- th = gmx_omp_get_thread_num();
- if (th > 0)
- {
- /* Initialize on each thread to get memory allocated thread-local */
- sd->gaussrand[th] = gmx_rng_init(seed[th]);
- }
- }
-
- sfree(seed);
-}
-
-static gmx_stochd_t *init_stochd(t_inputrec *ir, int nthreads)
+static gmx_stochd_t *init_stochd(t_inputrec *ir)
{
gmx_stochd_t *sd;
gmx_sd_const_t *sdc;
snew(sd, 1);
- /* Initiate random number generator for langevin type dynamics,
- * for BD, SD or velocity rescaling temperature coupling.
- */
- if (ir->eI == eiBD || EI_SD(ir->eI))
- {
- sd->ngaussrand = nthreads;
- }
- else
- {
- sd->ngaussrand = 1;
- }
- snew(sd->gaussrand, sd->ngaussrand);
-
- /* Initialize the first random generator */
- sd->gaussrand[0] = gmx_rng_init(ir->ld_seed);
-
- if (sd->ngaussrand > 1)
- {
- /* Initialize the rest of the random number generators,
- * using the first one to generate seeds.
- */
- init_multiple_gaussrand(sd);
- }
-
ngtc = ir->opts.ngtc;
if (ir->eI == eiBD)
return sd;
}
-void get_stochd_state(gmx_update_t upd, t_state *state)
-{
- /* Note that we only get the state of the first random generator,
- * even if there are multiple. This avoids repetition.
- */
- gmx_rng_get_state(upd->sd->gaussrand[0], state->ld_rng, state->ld_rngi);
-}
-
-void set_stochd_state(gmx_update_t upd, t_state *state)
-{
- gmx_stochd_t *sd;
- int i;
-
- sd = upd->sd;
-
- gmx_rng_set_state(sd->gaussrand[0], state->ld_rng, state->ld_rngi[0]);
-
- if (sd->ngaussrand > 1)
- {
- /* We only end up here with SD or BD with OpenMP.
- * Destroy and reinitialize the rest of the random number generators,
- * using seeds generated from the first one.
- * Although this doesn't recover the previous state,
- * it at least avoids repetition, which is most important.
- * Exaclty restoring states with all MPI+OpenMP setups is difficult
- * and as the integrator is random to start with, doesn't gain us much.
- */
- for (i = 1; i < sd->ngaussrand; i++)
- {
- gmx_rng_destroy(sd->gaussrand[i]);
- }
-
- init_multiple_gaussrand(sd);
- }
-}
-
gmx_update_t init_update(t_inputrec *ir)
{
t_gmx_update *upd;
if (ir->eI == eiBD || EI_SD(ir->eI) || ir->etc == etcVRESCALE || ETC_ANDERSEN(ir->etc))
{
- upd->sd = init_stochd(ir, gmx_omp_nthreads_get(emntUpdate));
+ upd->sd = init_stochd(ir);
}
- upd->xp = NULL;
- upd->xp_nalloc = 0;
- upd->randatom = NULL;
- upd->randatom_list = NULL;
- upd->randatom_list_init = FALSE; /* we have not yet cleared the data structure at this point */
+ upd->xp = NULL;
+ upd->xp_nalloc = 0;
return upd;
}
static void do_update_sd1(gmx_stochd_t *sd,
- gmx_rng_t gaussrand,
int start, int nrend, double dt,
rvec accel[], ivec nFreeze[],
real invmass[], unsigned short ptype[],
unsigned short cFREEZE[], unsigned short cACC[],
unsigned short cTC[],
rvec x[], rvec xprime[], rvec v[], rvec f[],
- int ngtc, real tau_t[], real ref_t[])
+ int ngtc, real ref_t[],
+ gmx_bool bDoConstr,
+ gmx_bool bFirstHalfConstr,
+ gmx_int64_t step, int seed, int* gatindex)
{
gmx_sd_const_t *sdc;
gmx_sd_sigma_t *sig;
real kT;
int gf = 0, ga = 0, gt = 0;
- real ism, sd_V;
+ real ism;
int n, d;
sdc = sd->sdc;
sig[n].V = sqrt(kT*(1 - sdc[n].em*sdc[n].em));
}
- for (n = start; n < nrend; n++)
+ if (!bDoConstr)
{
- ism = sqrt(invmass[n]);
- if (cFREEZE)
- {
- gf = cFREEZE[n];
- }
- if (cACC)
+ for (n = start; n < nrend; n++)
{
- ga = cACC[n];
+ real rnd[3];
+ int ng = gatindex ? gatindex[n] : n;
+
+ ism = sqrt(invmass[n]);
+ if (cFREEZE)
+ {
+ gf = cFREEZE[n];
+ }
+ if (cACC)
+ {
+ ga = cACC[n];
+ }
+ if (cTC)
+ {
+ gt = cTC[n];
+ }
+
+ gmx_rng_cycle_3gaussian_table(step, ng, seed, RND_SEED_UPDATE, rnd);
+
+ for (d = 0; d < DIM; d++)
+ {
+ if ((ptype[n] != eptVSite) && (ptype[n] != eptShell) && !nFreeze[gf][d])
+ {
+ real sd_V, vn;
+
+ sd_V = ism*sig[gt].V*rnd[d];
+ vn = v[n][d] + (invmass[n]*f[n][d] + accel[ga][d])*dt;
+ v[n][d] = vn*sdc[gt].em + sd_V;
+ /* Here we include half of the friction+noise
+ * update of v into the integration of x.
+ */
+ xprime[n][d] = x[n][d] + 0.5*(vn + v[n][d])*dt;
+ }
+ else
+ {
+ v[n][d] = 0.0;
+ xprime[n][d] = x[n][d];
+ }
+ }
}
- if (cTC)
+ }
+ else
+ {
+ /* We do have constraints */
+ if (bFirstHalfConstr)
{
- gt = cTC[n];
- }
+ /* First update without friction and noise */
+ real im;
- for (d = 0; d < DIM; d++)
- {
- if ((ptype[n] != eptVSite) && (ptype[n] != eptShell) && !nFreeze[gf][d])
+ for (n = start; n < nrend; n++)
{
- sd_V = ism*sig[gt].V*gmx_rng_gaussian_table(gaussrand);
+ im = invmass[n];
- v[n][d] = v[n][d]*sdc[gt].em
- + (invmass[n]*f[n][d] + accel[ga][d])*tau_t[gt]*(1 - sdc[gt].em)
- + sd_V;
+ if (cFREEZE)
+ {
+ gf = cFREEZE[n];
+ }
+ if (cACC)
+ {
+ ga = cACC[n];
+ }
+ if (cTC)
+ {
+ gt = cTC[n];
+ }
- xprime[n][d] = x[n][d] + v[n][d]*dt;
+ for (d = 0; d < DIM; d++)
+ {
+ if ((ptype[n] != eptVSite) && (ptype[n] != eptShell) && !nFreeze[gf][d])
+ {
+ v[n][d] = v[n][d] + (im*f[n][d] + accel[ga][d])*dt;
+ xprime[n][d] = x[n][d] + v[n][d]*dt;
+ }
+ else
+ {
+ v[n][d] = 0.0;
+ xprime[n][d] = x[n][d];
+ }
+ }
}
- else
+ }
+ else
+ {
+ /* Update friction and noise only */
+ for (n = start; n < nrend; n++)
{
- v[n][d] = 0.0;
- xprime[n][d] = x[n][d];
+ real rnd[3];
+ int ng = gatindex ? gatindex[n] : n;
+
+ ism = sqrt(invmass[n]);
+ if (cFREEZE)
+ {
+ gf = cFREEZE[n];
+ }
+ if (cACC)
+ {
+ ga = cACC[n];
+ }
+ if (cTC)
+ {
+ gt = cTC[n];
+ }
+
+ gmx_rng_cycle_3gaussian_table(step, ng, seed, RND_SEED_UPDATE, rnd);
+
+ for (d = 0; d < DIM; d++)
+ {
+ if ((ptype[n] != eptVSite) && (ptype[n] != eptShell) && !nFreeze[gf][d])
+ {
+ real sd_V, vn;
+
+ sd_V = ism*sig[gt].V*rnd[d];
+ vn = v[n][d];
+ v[n][d] = vn*sdc[gt].em + sd_V;
+ /* Add the friction and noise contribution only */
+ xprime[n][d] = xprime[n][d] + 0.5*(v[n][d] - vn)*dt;
+ }
+ }
}
}
}
}
static void do_update_sd2(gmx_stochd_t *sd,
- gmx_rng_t gaussrand,
gmx_bool bInitStep,
int start, int nrend,
rvec accel[], ivec nFreeze[],
rvec x[], rvec xprime[], rvec v[], rvec f[],
rvec sd_X[],
const real tau_t[],
- gmx_bool bFirstHalf)
+ gmx_bool bFirstHalf, gmx_int64_t step, int seed,
+ int* gatindex)
{
gmx_sd_const_t *sdc;
gmx_sd_sigma_t *sig;
int gf = 0, ga = 0, gt = 0;
real vn = 0, Vmh, Xmh;
real ism;
- int n, d;
+ int n, d, ng;
sdc = sd->sdc;
sig = sd->sdsig;
for (n = start; n < nrend; n++)
{
+ real rnd[6], rndi[3];
+ ng = gatindex ? gatindex[n] : n;
ism = sqrt(invmass[n]);
if (cFREEZE)
{
gt = cTC[n];
}
+ gmx_rng_cycle_6gaussian_table(step*2+(bFirstHalf ? 1 : 2), ng, seed, RND_SEED_UPDATE, rnd);
+ if (bInitStep)
+ {
+ gmx_rng_cycle_3gaussian_table(step*2, ng, seed, RND_SEED_UPDATE, rndi);
+ }
for (d = 0; d < DIM; d++)
{
if (bFirstHalf)
{
if (bInitStep)
{
- sd_X[n][d] = ism*sig[gt].X*gmx_rng_gaussian_table(gaussrand);
+ sd_X[n][d] = ism*sig[gt].X*rndi[d];
}
Vmh = sd_X[n][d]*sdc[gt].d/(tau_t[gt]*sdc[gt].c)
- + ism*sig[gt].Yv*gmx_rng_gaussian_table(gaussrand);
- sd_V[n][d] = ism*sig[gt].V*gmx_rng_gaussian_table(gaussrand);
+ + ism*sig[gt].Yv*rnd[d*2];
+ sd_V[n][d] = ism*sig[gt].V*rnd[d*2+1];
v[n][d] = vn*sdc[gt].em
+ (invmass[n]*f[n][d] + accel[ga][d])*tau_t[gt]*(1 - sdc[gt].em)
}
else
{
-
/* Correct the velocities for the constraints.
* This operation introduces some inaccuracy,
* since the velocity is determined from differences in coordinates.
(xprime[n][d] - x[n][d])/(tau_t[gt]*(sdc[gt].eph - sdc[gt].emh));
Xmh = sd_V[n][d]*tau_t[gt]*sdc[gt].d/(sdc[gt].em-1)
- + ism*sig[gt].Yx*gmx_rng_gaussian_table(gaussrand);
- sd_X[n][d] = ism*sig[gt].X*gmx_rng_gaussian_table(gaussrand);
+ + ism*sig[gt].Yx*rnd[d*2];
+ sd_X[n][d] = ism*sig[gt].X*rnd[d*2+1];
xprime[n][d] += sd_X[n][d] - Xmh;
unsigned short cFREEZE[], unsigned short cTC[],
rvec x[], rvec xprime[], rvec v[],
rvec f[], real friction_coefficient,
- real *rf, gmx_rng_t gaussrand)
+ real *rf, gmx_int64_t step, int seed,
+ int* gatindex)
{
/* note -- these appear to be full step velocities . . . */
int gf = 0, gt = 0;
for (n = start; (n < nrend); n++)
{
+ real rnd[3];
+ int ng = gatindex ? gatindex[n] : n;
+
if (cFREEZE)
{
gf = cFREEZE[n];
{
gt = cTC[n];
}
+ gmx_rng_cycle_3gaussian_table(step, ng, seed, RND_SEED_UPDATE, rnd);
for (d = 0; (d < DIM); d++)
{
if ((ptype[n] != eptVSite) && (ptype[n] != eptShell) && !nFreeze[gf][d])
{
if (friction_coefficient != 0)
{
- vn = invfr*f[n][d] + rf[gt]*gmx_rng_gaussian_table(gaussrand);
+ vn = invfr*f[n][d] + rf[gt]*rnd[d];
}
else
{
/* NOTE: invmass = 2/(mass*friction_constant*dt) */
vn = 0.5*invmass[n]*f[n][d]*dt
- + sqrt(0.5*invmass[n])*rf[gt]*gmx_rng_gaussian_table(gaussrand);
+ + sqrt(0.5*invmass[n])*rf[gt]*rnd[d];
}
v[n][d] = vn;
matrix *ekin_sum;
real *dekindl_sum;
- start_t = md->start + ((thread+0)*md->homenr)/nthread;
- end_t = md->start + ((thread+1)*md->homenr)/nthread;
+ start_t = ((thread+0)*md->homenr)/nthread;
+ end_t = ((thread+1)*md->homenr)/nthread;
ekin_sum = ekind->ekin_work[thread];
dekindl_sum = ekind->dekindl_work[thread];
gmx_ekindata_t *ekind,
t_nrnb *nrnb, gmx_bool bEkinAveVel)
{
- int start = md->start, homenr = md->homenr;
+ int start = 0, homenr = md->homenr;
int g, d, n, m, gt = 0;
rvec v_corrt;
real hm;
}
if (md->nPerturbed && md->bPerturbed[n])
{
- dekindl += 0.5*(md->massB[n] - md->massA[n])*iprod(v_corrt, v_corrt);
+ /* The minus sign here might be confusing.
+ * The kinetic contribution from dH/dl doesn't come from
+ * d m(l)/2 v^2 / dl, but rather from d p^2/2m(l) / dl,
+ * where p are the momenta. The difference is only a minus sign.
+ */
+ dekindl -= 0.5*(md->massB[n] - md->massA[n])*iprod(v_corrt, v_corrt);
}
}
ekind->dekindl = dekindl;
}
}
-void set_deform_reference_box(gmx_update_t upd, gmx_large_int_t step, matrix box)
+void set_deform_reference_box(gmx_update_t upd, gmx_int64_t step, matrix box)
{
upd->deformref_step = step;
copy_mat(box, upd->deformref_box);
static void deform(gmx_update_t upd,
int start, int homenr, rvec x[], matrix box, matrix *scale_tot,
- const t_inputrec *ir, gmx_large_int_t step)
+ const t_inputrec *ir, gmx_int64_t step)
{
matrix bnew, invbox, mu;
real elapsed_time;
x[i][YY] = mu[YY][YY]*x[i][YY]+mu[ZZ][YY]*x[i][ZZ];
x[i][ZZ] = mu[ZZ][ZZ]*x[i][ZZ];
}
- if (*scale_tot)
+ if (scale_tot != NULL)
{
/* The transposes of the scaling matrices are stored,
* so we need to do matrix multiplication in the inverse order.
}
}
-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(gmx_large_int_t step,
+void update_tcouple(gmx_int64_t step,
t_inputrec *inputrec,
t_state *state,
gmx_ekindata_t *ekind,
- gmx_update_t upd,
t_extmass *MassQ,
t_mdatoms *md)
state->nosehoover_xi, state->nosehoover_vxi, MassQ);
break;
case etcVRESCALE:
- vrescale_tcoupl(inputrec, ekind, dttc,
- state->therm_integral, upd->sd->gaussrand[0]);
+ vrescale_tcoupl(inputrec, step, ekind, dttc,
+ state->therm_integral);
break;
}
/* rescale in place here */
if (EI_VV(inputrec->eI))
{
- rescale_velocities(ekind, md, md->start, md->start+md->homenr, state->v);
+ rescale_velocities(ekind, md, 0, md->homenr, state->v);
}
}
else
}
void update_pcouple(FILE *fplog,
- gmx_large_int_t step,
+ gmx_int64_t step,
t_inputrec *inputrec,
t_state *state,
matrix pcoupl_mu,
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, 1.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,
+ gmx_int64_t step,
real *dvdlambda, /* the contribution to be added to the bonded interactions */
t_inputrec *inputrec, /* input record and box stuff */
gmx_ekindata_t *ekind,
/* for now, SD update is here -- though it really seems like it
should be reformulated as a velocity verlet method, since it has two parts */
- start = md->start;
+ start = 0;
homenr = md->homenr;
nrend = start+homenr;
if (EI_VV(inputrec->eI) && bFirstHalf)
{
constrain(NULL, bLog, bEner, constr, idef,
- inputrec, ekind, cr, step, 1, md,
+ inputrec, ekind, cr, step, 1, 1.0, md,
state->x, state->v, state->v,
bMolPBC, state->box,
state->lambda[efptBONDED], dvdlambda,
else
{
constrain(NULL, bLog, bEner, constr, idef,
- inputrec, ekind, cr, step, 1, md,
+ inputrec, ekind, cr, step, 1, 1.0, md,
state->x, xprime, NULL,
bMolPBC, state->box,
state->lambda[efptBONDED], dvdlambda,
}
where();
+
+ if (inputrec->eI == eiSD1 && bDoConstr && !bFirstHalf)
+ {
+ wallcycle_start(wcycle, ewcUPDATE);
+ xprime = get_xprime(state, upd);
+
+ nth = gmx_omp_nthreads_get(emntUpdate);
+
+#pragma omp parallel for num_threads(nth) schedule(static)
+
+ for (th = 0; th < nth; th++)
+ {
+ int start_th, end_th;
+
+ start_th = start + ((nrend-start)* th )/nth;
+ end_th = start + ((nrend-start)*(th+1))/nth;
+
+ /* The second part of the SD integration */
+ do_update_sd1(upd->sd,
+ start_th, end_th, dt,
+ inputrec->opts.acc, inputrec->opts.nFreeze,
+ md->invmass, md->ptype,
+ md->cFREEZE, md->cACC, md->cTC,
+ state->x, xprime, state->v, force,
+ inputrec->opts.ngtc, inputrec->opts.ref_t,
+ bDoConstr, FALSE,
+ step, inputrec->ld_seed,
+ DOMAINDECOMP(cr) ? cr->dd->gatindex : NULL);
+ }
+ inc_nrnb(nrnb, eNR_UPDATE, homenr);
+ wallcycle_stop(wcycle, ewcUPDATE);
+
+ if (bDoConstr)
+ {
+ /* Constrain the coordinates xprime for half a time step */
+ wallcycle_start(wcycle, ewcCONSTR);
+
+ constrain(NULL, bLog, bEner, constr, idef,
+ inputrec, NULL, cr, step, 1, 0.5, md,
+ state->x, xprime, NULL,
+ bMolPBC, state->box,
+ state->lambda[efptBONDED], dvdlambda,
+ state->v, NULL, nrnb, econqCoord, FALSE, 0, 0);
+
+ wallcycle_stop(wcycle, ewcCONSTR);
+ }
+ }
+
if ((inputrec->eI == eiSD2) && !(bFirstHalf))
{
+ wallcycle_start(wcycle, ewcUPDATE);
xprime = get_xprime(state, upd);
nth = gmx_omp_nthreads_get(emntUpdate);
end_th = start + ((nrend-start)*(th+1))/nth;
/* The second part of the SD integration */
- do_update_sd2(upd->sd, upd->sd->gaussrand[th],
+ do_update_sd2(upd->sd,
FALSE, start_th, end_th,
inputrec->opts.acc, inputrec->opts.nFreeze,
md->invmass, md->ptype,
md->cFREEZE, md->cACC, md->cTC,
state->x, xprime, state->v, force, state->sd_X,
inputrec->opts.tau_t,
- FALSE);
+ FALSE, step, inputrec->ld_seed,
+ DOMAINDECOMP(cr) ? cr->dd->gatindex : NULL);
}
inc_nrnb(nrnb, eNR_UPDATE, homenr);
+ wallcycle_stop(wcycle, ewcUPDATE);
if (bDoConstr)
{
/* Constrain the coordinates xprime */
wallcycle_start(wcycle, ewcCONSTR);
constrain(NULL, bLog, bEner, constr, idef,
- inputrec, NULL, cr, step, 1, md,
+ inputrec, NULL, cr, step, 1, 1.0, md,
state->x, xprime, NULL,
bMolPBC, state->box,
state->lambda[efptBONDED], dvdlambda,
}
}
+
/* We must always unshift after updating coordinates; if we did not shake
x was shifted in do_force */
}
else
{
-#pragma omp parallel for num_threads(gmx_omp_nthreads_get(emntUpdate)) schedule(static)
+ nth = gmx_omp_nthreads_get(emntUpdate);
+
+#pragma omp parallel for num_threads(nth) schedule(static)
for (i = start; i < nrend; i++)
{
copy_rvec(upd->xp[i], state->x[i]);
}
void update_box(FILE *fplog,
- gmx_large_int_t step,
+ gmx_int64_t step,
t_inputrec *inputrec, /* input record and box stuff */
t_mdatoms *md,
t_state *state,
int start, homenr, nrend, i, n, m, g;
tensor vir_con;
- start = md->start;
+ start = 0;
homenr = md->homenr;
nrend = start+homenr;
}
void update_coords(FILE *fplog,
- gmx_large_int_t step,
+ gmx_int64_t step,
t_inputrec *inputrec, /* input record and box stuff */
t_mdatoms *md,
t_state *state,
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,
gmx_constr_t constr,
t_idef *idef)
{
- gmx_bool bNH, bPR, bLastStep, bLog = FALSE, bEner = FALSE;
+ gmx_bool bNH, bPR, bLastStep, bLog = FALSE, bEner = FALSE, bDoConstr = FALSE;
double dt, alpha;
real *imass, *imassin;
rvec *force;
rvec *vcom, *xcom, *vall, *xall, *xin, *vin, *forcein, *fall, *xpall, *xprimein, *xprime;
int nth, th;
+ bDoConstr = (NULL != constr);
+
/* Running the velocity half does nothing except for velocity verlet */
if ((UpdatePart == etrtVELOCITY1 || UpdatePart == etrtVELOCITY2) &&
!EI_VV(inputrec->eI))
gmx_incons("update_coords called for velocity without VV integrator");
}
- start = md->start;
+ start = 0;
homenr = md->homenr;
nrend = start+homenr;
* 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
}
break;
case (eiSD1):
- do_update_sd1(upd->sd, upd->sd->gaussrand[th],
+ /* With constraints, the SD1 update is done in 2 parts */
+ do_update_sd1(upd->sd,
start_th, end_th, dt,
inputrec->opts.acc, inputrec->opts.nFreeze,
md->invmass, md->ptype,
md->cFREEZE, md->cACC, md->cTC,
state->x, xprime, state->v, force,
- inputrec->opts.ngtc, inputrec->opts.tau_t, inputrec->opts.ref_t);
+ inputrec->opts.ngtc, inputrec->opts.ref_t,
+ bDoConstr, TRUE,
+ step, inputrec->ld_seed, DOMAINDECOMP(cr) ? cr->dd->gatindex : NULL);
break;
case (eiSD2):
- /* The SD update is done in 2 parts, because an extra constraint step
- * is needed
+ /* The SD2 update is always done in 2 parts,
+ * because an extra constraint step is needed
*/
- do_update_sd2(upd->sd, upd->sd->gaussrand[th],
+ do_update_sd2(upd->sd,
bInitStep, start_th, end_th,
inputrec->opts.acc, inputrec->opts.nFreeze,
md->invmass, md->ptype,
md->cFREEZE, md->cACC, md->cTC,
state->x, xprime, state->v, force, state->sd_X,
inputrec->opts.tau_t,
- TRUE);
+ TRUE, step, inputrec->ld_seed,
+ DOMAINDECOMP(cr) ? cr->dd->gatindex : NULL);
break;
case (eiBD):
do_update_bd(start_th, end_th, dt,
md->cFREEZE, md->cTC,
state->x, xprime, state->v, force,
inputrec->bd_fric,
- upd->sd->bd_rf, upd->sd->gaussrand[th]);
+ upd->sd->bd_rf,
+ step, inputrec->ld_seed, DOMAINDECOMP(cr) ? cr->dd->gatindex : NULL);
break;
case (eiVV):
case (eiVVAK):
mv[XX], mv[YY], mv[ZZ]);
}
-extern gmx_bool update_randomize_velocities(t_inputrec *ir, gmx_large_int_t step, t_mdatoms *md, t_state *state, gmx_update_t upd, t_idef *idef, gmx_constr_t constr)
+extern gmx_bool update_randomize_velocities(t_inputrec *ir, gmx_int64_t step, const t_commrec *cr,
+ t_mdatoms *md, t_state *state, gmx_update_t upd, gmx_constr_t constr)
{
int i;
real rate = (ir->delta_t)/ir->opts.tau_t[0];
+
+ if (ir->etc == etcANDERSEN && constr != NULL)
+ {
+ gmx_fatal(FARGS, "Normal Andersen is currently not supported with constraints, use massive Andersen instead");
+ }
+
/* proceed with andersen if 1) it's fixed probability per
particle andersen or 2) it's massive andersen and it's tau_t/dt */
if ((ir->etc == etcANDERSEN) || do_per_step(step, (int)(1.0/rate)))
{
- srenew(upd->randatom, state->nalloc);
- srenew(upd->randatom_list, state->nalloc);
- if (upd->randatom_list_init == FALSE)
- {
- for (i = 0; i < state->nalloc; i++)
- {
- upd->randatom[i] = FALSE;
- upd->randatom_list[i] = 0;
- }
- upd->randatom_list_init = TRUE;
- }
- andersen_tcoupl(ir, md, state, upd->sd->gaussrand[0], rate,
- (ir->etc == etcANDERSEN) ? idef : NULL,
- constr ? get_nblocks(constr) : 0,
- constr ? get_sblock(constr) : NULL,
- upd->randatom, upd->randatom_list,
+ andersen_tcoupl(ir, step, cr, md, state, rate,
upd->sd->randomize_group, upd->sd->boltzfac);
return TRUE;
}
biod.pnpi.spb.ru / alexxy / gromacs.git / blobdiff