gmx_mdoutf_t *outf, t_mdebin **mdebin,
tensor force_vir, tensor shake_vir,
tensor total_vir, tensor pres, rvec mu_tot,
- gmx_bool *bSimAnn, t_vcm **vcm,
+ gmx_bool *bSimAnn,
gmx_wallcycle_t wcycle)
{
int i;
update_annealing_target_temp(ir, ir->init_t, upd ? *upd : nullptr);
}
- if (vcm != nullptr)
- {
- *vcm = init_vcm(fplog, &mtop->groups, ir);
- }
-
if (EI_DYNAMICS(ir->eI) && !mdrunOptions.continuationOptions.appendFiles)
{
if (ir->etc == etcBERENDSEN)
*
* Copyright (c) 1991-2000, University of Groningen, The Netherlands.
* Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2013,2014,2015,2016,2017,2018, by the GROMACS development team, led by
+ * Copyright (c) 2013,2014,2015,2016,2017,2018,2019, 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.
tensor force_vir, tensor shake_vir,
tensor total_vir, tensor pres,
rvec mu_tot,
- gmx_bool *bSimAnn, t_vcm **vcm,
+ gmx_bool *bSimAnn,
gmx_wallcycle_t wcycle);
void init_rerun(FILE *fplog,
*
* Copyright (c) 1991-2000, University of Groningen, The Netherlands.
* Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2013,2014,2015,2016,2017,2018, by the GROMACS development team, led by
+ * Copyright (c) 2013,2014,2015,2016,2017,2018,2019, 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.
if (vcm)
{
icm = add_binr(rb, DIM*vcm->nr, vcm->group_p[0]);
- imass = add_binr(rb, vcm->nr, vcm->group_mass);
+ imass = add_binr(rb, vcm->nr, vcm->group_mass.data());
if (vcm->mode == ecmANGULAR)
{
icj = add_binr(rb, DIM*vcm->nr, vcm->group_j[0]);
if (vcm)
{
extract_binr(rb, icm, DIM*vcm->nr, vcm->group_p[0]);
- extract_binr(rb, imass, vcm->nr, vcm->group_mass);
+ extract_binr(rb, imass, vcm->nr, vcm->group_mass.data());
if (vcm->mode == ecmANGULAR)
{
extract_binr(rb, icj, DIM*vcm->nr, vcm->group_j[0]);
*
* Copyright (c) 1991-2000, University of Groningen, The Netherlands.
* Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2013,2014,2015,2016,2017,2018, by the GROMACS development team, led by
+ * Copyright (c) 2013,2014,2015,2016,2017,2018,2019, 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.
#include "gromacs/utility/gmxomp.h"
#include "gromacs/utility/smalloc.h"
-t_vcm *init_vcm(FILE *fp, const gmx_groups_t *groups, const t_inputrec *ir)
+t_vcm::t_vcm(const gmx_groups_t &groups, const t_inputrec &ir)
{
- t_vcm *vcm;
- int g;
+ mode = (ir.nstcomm > 0) ? ir.comm_mode : ecmNO;
+ ndim = ndof_com(&ir);
+ timeStep = ir.nstcomm*ir.delta_t;
- snew(vcm, 1);
-
- vcm->mode = (ir->nstcomm > 0) ? ir->comm_mode : ecmNO;
- vcm->ndim = ndof_com(ir);
- vcm->timeStep = ir->nstcomm*ir->delta_t;
-
- if (vcm->mode == ecmANGULAR && vcm->ndim < 3)
+ if (mode == ecmANGULAR && ndim < 3)
{
gmx_fatal(FARGS, "Can not have angular comm removal with pbc=%s",
- epbc_names[ir->ePBC]);
+ epbc_names[ir.ePBC]);
}
- if (vcm->mode != ecmNO)
+ if (mode != ecmNO)
{
- vcm->nr = groups->grps[egcVCM].nr;
+ nr = groups.grps[egcVCM].nr;
/* Allocate one extra for a possible rest group */
- vcm->size = vcm->nr + 1;
+ size = nr + 1;
/* We need vcm->nr+1 elements per thread, but to avoid cache
* invalidation we add 2 elements to get a 152 byte separation.
*/
- vcm->stride = vcm->nr + 3;
- if (vcm->mode == ecmANGULAR)
+ stride = nr + 3;
+ if (mode == ecmANGULAR)
{
- snew(vcm->group_j, vcm->size);
- snew(vcm->group_x, vcm->size);
- snew(vcm->group_i, vcm->size);
- snew(vcm->group_w, vcm->size);
- }
- snew(vcm->group_p, vcm->size);
- snew(vcm->group_v, vcm->size);
- snew(vcm->group_mass, vcm->size);
- snew(vcm->group_name, vcm->size);
- snew(vcm->group_ndf, vcm->size);
- for (g = 0; (g < vcm->nr); g++)
- {
- vcm->group_ndf[g] = ir->opts.nrdf[g];
- }
+ snew(group_i, size);
+
+ group_j.resize(size);
+ group_x.resize(size);
+ group_w.resize(size);
- /* Copy pointer to group names and print it. */
- if (fp)
- {
- fprintf(fp, "Center of mass motion removal mode is %s\n",
- ECOM(vcm->mode));
- fprintf(fp, "We have the following groups for center of"
- " mass motion removal:\n");
}
- for (g = 0; (g < vcm->nr); g++)
+
+ group_name.resize(size);
+ group_p.resize(size);
+ group_v.resize(size);
+ group_mass.resize(size);
+ group_ndf.resize(size);
+ for (int g = 0; (g < nr); g++)
{
- vcm->group_name[g] = *groups->grpname[groups->grps[egcVCM].nm_ind[g]];
- if (fp)
- {
- fprintf(fp, "%3d: %s\n", g, vcm->group_name[g]);
- }
+ group_ndf[g] = ir.opts.nrdf[g];
+ group_name[g] = *groups.grpname[groups.grps[egcVCM].nm_ind[g]];
+
}
- snew(vcm->thread_vcm, gmx_omp_nthreads_get(emntDefault) * vcm->stride);
+ thread_vcm.resize(gmx_omp_nthreads_get(emntDefault) * stride);
}
- vcm->nFreeze = ir->opts.nFreeze;
+ nFreeze = ir.opts.nFreeze;
+}
- return vcm;
+t_vcm::~t_vcm()
+{
+ if (mode == ecmANGULAR)
+ {
+ sfree(group_i);
+ }
}
+void reportComRemovalInfo(FILE * fp, const t_vcm &vcm)
+{
+
+ /* Copy pointer to group names and print it. */
+ if (fp && vcm.mode != ecmNO)
+ {
+ fprintf(fp, "Center of mass motion removal mode is %s\n",
+ ECOM(vcm.mode));
+ fprintf(fp, "We have the following groups for center of"
+ " mass motion removal:\n");
+
+ for (int g = 0; (g < vcm.nr); g++)
+ {
+
+ fprintf(fp, "%3d: %s\n", g, vcm.group_name[g]);
+ }
+ }
+}
+
+
static void update_tensor(const rvec x, real m0, tensor I)
{
real xy, xz, yz;
*
* Copyright (c) 1991-2000, University of Groningen, The Netherlands.
* Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2013,2014,2015,2016,2017,2018, by the GROMACS development team, led by
+ * Copyright (c) 2013,2014,2015,2016,2017,2018,2019, 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.
#include <stdio.h>
+#include <vector>
+
#include "gromacs/math/vectypes.h"
#include "gromacs/mdtypes/mdatom.h"
#include "gromacs/utility/basedefinitions.h"
struct t_inputrec;
struct t_mdatoms;
-typedef struct {
- rvec p; /* Linear momentum */
- rvec x; /* Center of mass */
- rvec j; /* Angular momentum */
- tensor i; /* Moment of inertia */
- real mass; /* Mass */
-} t_vcm_thread;
+struct t_vcm_thread{
+ rvec p = {0}; /* Linear momentum */
+ rvec x = {0}; /* Center of mass */
+ rvec j = {0}; /* Angular momentum */
+ tensor i = {{0}}; /* Moment of inertia */
+ real mass = 0; /* Mass */
+};
+
+struct t_vcm
+{
+ int nr; /* Number of groups */
+ int size; /* Size of group arrays */
+ int stride; /* Stride for thread data */
+ int mode; /* One of the enums above */
+ int ndim; /* The number of dimensions for corr. */
+ real timeStep; /* The time step for COMM removal */
+ std::vector<real> group_ndf; /* Number of degrees of freedom */
+ std::vector<real> group_mass; /* Mass per group */
+ std::vector<gmx::RVec> group_p; /* Linear momentum per group */
+ std::vector<gmx::RVec> group_v; /* Linear velocity per group */
+ std::vector<gmx::RVec> group_x; /* Center of mass per group */
+ std::vector<gmx::RVec> group_j; /* Angular momentum per group */
+ std::vector<gmx::RVec> group_w; /* Angular velocity (omega) */
+ tensor *group_i; /* Moment of inertia per group */
+ std::vector<char *> group_name; /* These two are copies to pointers in */
+ ivec *nFreeze; /* Tells whether dimensions are frozen per freeze group */
+ std::vector<t_vcm_thread> thread_vcm; /* Temporary data per thread and group */
+
+ t_vcm(const gmx_groups_t &groups, const t_inputrec &ir);
+ ~t_vcm();
+};
-typedef struct {
- int nr; /* Number of groups */
- int size; /* Size of group arrays */
- int stride; /* Stride for thread data */
- int mode; /* One of the enums above */
- int ndim; /* The number of dimensions for corr. */
- real timeStep; /* The time step for COMM removal */
- real *group_ndf; /* Number of degrees of freedom */
- rvec *group_p; /* Linear momentum per group */
- rvec *group_v; /* Linear velocity per group */
- rvec *group_x; /* Center of mass per group */
- rvec *group_j; /* Angular momentum per group */
- rvec *group_w; /* Angular velocity (omega) */
- tensor *group_i; /* Moment of inertia per group */
- real *group_mass; /* Mass per group */
- char **group_name; /* These two are copies to pointers in */
- ivec *nFreeze; /* Tells whether dimensions are frozen per freeze group */
- t_vcm_thread *thread_vcm; /* Temporary data per thread and group */
-} t_vcm;
+/* print COM removal info to log */
+void reportComRemovalInfo(FILE * fp, const t_vcm &vcm);
-t_vcm *init_vcm(FILE *fp, const gmx_groups_t *groups, const t_inputrec *ir);
/* Do a per group center of mass things */
void calc_vcm_grp(int start, int homenr, t_mdatoms *md,
tensor force_vir, shake_vir, total_vir, tmp_vir, pres;
int i, m;
rvec mu_tot;
- t_vcm *vcm;
matrix parrinellorahmanMu, M;
gmx_repl_ex_t repl_ex = nullptr;
gmx_localtop_t top;
&t, &t0, state_global, lam0,
nrnb, top_global, &upd, deform,
nfile, fnm, &outf, &mdebin,
- force_vir, shake_vir, total_vir, pres, mu_tot, &bSimAnn, &vcm, wcycle);
+ force_vir, shake_vir, total_vir, pres, mu_tot, &bSimAnn, wcycle);
/* Energy terms and groups */
snew(enerd, 1);
| (continuationOptions.haveReadEkin ? CGLO_READEKIN : 0));
bSumEkinhOld = FALSE;
+
+ t_vcm vcm(top_global->groups, *ir);
+ reportComRemovalInfo(fplog, vcm);
+
/* To minimize communication, compute_globals computes the COM velocity
* and the kinetic energy for the velocities without COM motion removed.
* Thus to get the kinetic energy without the COM contribution, we need
cglo_flags_iteration |= CGLO_STOPCM;
cglo_flags_iteration &= ~CGLO_TEMPERATURE;
}
- compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, vcm,
+ compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, &vcm,
nullptr, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
constr, &nullSignaller, state->box,
&totalNumberOfBondedInteractions, &bSumEkinhOld, cglo_flags_iteration
kinetic energy calculation. This minimized excess variables, but
perhaps loses some logic?*/
- compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, vcm,
+ compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, &vcm,
nullptr, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
constr, &nullSignaller, state->box,
nullptr, &bSumEkinhOld,
/* We need the kinetic energy at minus the half step for determining
* the full step kinetic energy and possibly for T-coupling.*/
/* This may not be quite working correctly yet . . . . */
- compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, vcm,
+ compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, &vcm,
wcycle, enerd, nullptr, nullptr, nullptr, nullptr, mu_tot,
constr, &nullSignaller, state->box,
&totalNumberOfBondedInteractions, &bSumEkinhOld,
if (bGStat || do_per_step(step-1, nstglobalcomm))
{
wallcycle_stop(wcycle, ewcUPDATE);
- compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, vcm,
+ compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, &vcm,
wcycle, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
constr, &nullSignaller, state->box,
&totalNumberOfBondedInteractions, &bSumEkinhOld,
/* We need the kinetic energy at minus the half step for determining
* the full step kinetic energy and possibly for T-coupling.*/
/* This may not be quite working correctly yet . . . . */
- compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, vcm,
+ compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, &vcm,
wcycle, enerd, nullptr, nullptr, nullptr, nullptr, mu_tot,
constr, &nullSignaller, state->box,
nullptr, &bSumEkinhOld,
{
/* erase F_EKIN and F_TEMP here? */
/* just compute the kinetic energy at the half step to perform a trotter step */
- compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, vcm,
+ compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, &vcm,
wcycle, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
constr, &nullSignaller, lastbox,
nullptr, &bSumEkinhOld,
bool doIntraSimSignal = true;
SimulationSignaller signaller(&signals, cr, ms, doInterSimSignal, doIntraSimSignal);
- compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, vcm,
+ compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, &vcm,
wcycle, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
constr, &signaller,
lastbox,
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff