From: Kevin Boyd Date: Sun, 23 Dec 2018 01:13:26 +0000 (-0500) Subject: Convert t_vcm to C++, plugging a memory leak X-Git-Url: http://biod.pnpi.spb.ru/gitweb/?a=commitdiff_plain;h=d76bcf15bdadfc3ec6f2faf003fed7b35e94cfa9;p=alexxy%2Fgromacs.git Convert t_vcm to C++, plugging a memory leak Change init function to constructor Extract COM log file writing from init to a function Make do_md owner of t_vcm with unique_ptr Move initialization of t_vcm out of init_md Change pointer arrays to vectors where possible. Where not possible, we now free dynamically allocated resources in destructor. Change-Id: I711a8d81ec0343969dcb437d9a3b8a474d291e87 --- diff --git a/src/gromacs/mdlib/sim_util.cpp b/src/gromacs/mdlib/sim_util.cpp index 469e7b7458..22b13e8825 100644 --- a/src/gromacs/mdlib/sim_util.cpp +++ b/src/gromacs/mdlib/sim_util.cpp @@ -2944,7 +2944,7 @@ void init_md(FILE *fplog, 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; @@ -2988,11 +2988,6 @@ void init_md(FILE *fplog, 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) diff --git a/src/gromacs/mdlib/sim_util.h b/src/gromacs/mdlib/sim_util.h index a9bd066d75..4ca378b515 100644 --- a/src/gromacs/mdlib/sim_util.h +++ b/src/gromacs/mdlib/sim_util.h @@ -3,7 +3,7 @@ * * 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. @@ -156,7 +156,7 @@ void init_md(FILE *fplog, 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, diff --git a/src/gromacs/mdlib/stat.cpp b/src/gromacs/mdlib/stat.cpp index ba9a135ecf..a2159835c7 100644 --- a/src/gromacs/mdlib/stat.cpp +++ b/src/gromacs/mdlib/stat.cpp @@ -3,7 +3,7 @@ * * 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. @@ -261,7 +261,7 @@ void global_stat(const gmx_global_stat *gs, 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]); @@ -359,7 +359,7 @@ void global_stat(const gmx_global_stat *gs, 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]); diff --git a/src/gromacs/mdlib/vcm.cpp b/src/gromacs/mdlib/vcm.cpp index d7a56714be..b0e8881668 100644 --- a/src/gromacs/mdlib/vcm.cpp +++ b/src/gromacs/mdlib/vcm.cpp @@ -3,7 +3,7 @@ * * 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. @@ -54,74 +54,83 @@ #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; diff --git a/src/gromacs/mdlib/vcm.h b/src/gromacs/mdlib/vcm.h index 716cf07f70..8351d97e4e 100644 --- a/src/gromacs/mdlib/vcm.h +++ b/src/gromacs/mdlib/vcm.h @@ -3,7 +3,7 @@ * * 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. @@ -39,6 +39,8 @@ #include +#include + #include "gromacs/math/vectypes.h" #include "gromacs/mdtypes/mdatom.h" #include "gromacs/utility/basedefinitions.h" @@ -48,35 +50,41 @@ struct gmx_groups_t; 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 group_ndf; /* Number of degrees of freedom */ + std::vector group_mass; /* Mass per group */ + std::vector group_p; /* Linear momentum per group */ + std::vector group_v; /* Linear velocity per group */ + std::vector group_x; /* Center of mass per group */ + std::vector group_j; /* Angular momentum per group */ + std::vector group_w; /* Angular velocity (omega) */ + tensor *group_i; /* Moment of inertia per group */ + std::vector group_name; /* These two are copies to pointers in */ + ivec *nFreeze; /* Tells whether dimensions are frozen per freeze group */ + std::vector 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, diff --git a/src/gromacs/mdrun/md.cpp b/src/gromacs/mdrun/md.cpp index 34c5e0b08b..1d108d726f 100644 --- a/src/gromacs/mdrun/md.cpp +++ b/src/gromacs/mdrun/md.cpp @@ -166,7 +166,6 @@ void gmx::Integrator::do_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; @@ -256,7 +255,7 @@ void gmx::Integrator::do_md() &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); @@ -486,6 +485,10 @@ void gmx::Integrator::do_md() | (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 @@ -499,7 +502,7 @@ void gmx::Integrator::do_md() 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 @@ -516,7 +519,7 @@ void gmx::Integrator::do_md() 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, @@ -787,7 +790,7 @@ void gmx::Integrator::do_md() /* 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, @@ -935,7 +938,7 @@ void gmx::Integrator::do_md() 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, @@ -990,7 +993,7 @@ void gmx::Integrator::do_md() /* 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, @@ -1171,7 +1174,7 @@ void gmx::Integrator::do_md() { /* 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, @@ -1258,7 +1261,7 @@ void gmx::Integrator::do_md() 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,