From 05e9a6d19035cd45034bcbf74e5a951b7bfc5a9e Mon Sep 17 00:00:00 2001 From: Artem Zhmurov Date: Tue, 9 Jun 2020 14:49:28 +0000 Subject: [PATCH] Make SD stuff private for update.cpp 1. Make update_sd_second_half(...) into a method of Update class 2. Make update_temperature_constants(...) into a method of Update class 3. Unexpose gmx_stochd_t TODO: Need beter solution for Andersen stuff. --- src/gromacs/mdlib/coupling.cpp | 225 +++++++++++++++++ src/gromacs/mdlib/coupling.h | 231 ++++++++++++++++++ src/gromacs/mdlib/md_support.cpp | 1 + src/gromacs/mdlib/tgroup.cpp | 2 +- src/gromacs/mdlib/update.cpp | 221 ----------------- src/gromacs/mdlib/update.h | 167 ------------- src/gromacs/mdrun/md.cpp | 1 + src/gromacs/mdrun/minimize.cpp | 1 + .../modularsimulator/energyelement.cpp | 1 + .../modularsimulator/modularsimulator.cpp | 1 + .../parrinellorahmanbarostat.cpp | 2 +- .../modularsimulator/vrescalethermostat.cpp | 4 +- 12 files changed, 465 insertions(+), 392 deletions(-) create mode 100644 src/gromacs/mdlib/coupling.h diff --git a/src/gromacs/mdlib/coupling.cpp b/src/gromacs/mdlib/coupling.cpp index 6ea616fe94..5545af97c8 100644 --- a/src/gromacs/mdlib/coupling.cpp +++ b/src/gromacs/mdlib/coupling.cpp @@ -37,6 +37,8 @@ */ #include "gmxpre.h" +#include "coupling.h" + #include #include @@ -49,6 +51,7 @@ #include "gromacs/math/units.h" #include "gromacs/math/vec.h" #include "gromacs/math/vecdump.h" +#include "gromacs/mdlib/boxdeformation.h" #include "gromacs/mdlib/expanded.h" #include "gromacs/mdlib/gmx_omp_nthreads.h" #include "gromacs/mdlib/stat.h" @@ -89,6 +92,228 @@ static const double sy_const_5[] = { 0.2967324292201065, 0.2967324292201065, -0. static const double* sy_const[] = { nullptr, sy_const_1, nullptr, sy_const_3, nullptr, sy_const_5 }; + +void update_tcouple(int64_t step, + const t_inputrec* inputrec, + t_state* state, + gmx_ekindata_t* ekind, + const t_extmass* MassQ, + const t_mdatoms* md) + +{ + // This condition was explicitly checked in previous version, but should have never been satisfied + GMX_ASSERT(!(EI_VV(inputrec->eI) + && (inputrecNvtTrotter(inputrec) || inputrecNptTrotter(inputrec) + || inputrecNphTrotter(inputrec))), + "Temperature coupling was requested with velocity verlet and trotter"); + + bool doTemperatureCoupling = false; + + // For VV temperature coupling parameters are updated on the current + // step, for the others - one step before. + if (inputrec->etc == etcNO) + { + doTemperatureCoupling = false; + } + else if (EI_VV(inputrec->eI)) + { + doTemperatureCoupling = do_per_step(step, inputrec->nsttcouple); + } + else + { + doTemperatureCoupling = do_per_step(step + inputrec->nsttcouple - 1, inputrec->nsttcouple); + } + + if (doTemperatureCoupling) + { + real dttc = inputrec->nsttcouple * inputrec->delta_t; + + // TODO: berendsen_tcoupl(...), nosehoover_tcoupl(...) and vrescale_tcoupl(...) update + // temperature coupling parameters, which should be reflected in the name of these + // subroutines + switch (inputrec->etc) + { + case etcNO: break; + case etcBERENDSEN: + berendsen_tcoupl(inputrec, ekind, dttc, state->therm_integral); + break; + case etcNOSEHOOVER: + nosehoover_tcoupl(&(inputrec->opts), ekind, dttc, state->nosehoover_xi.data(), + state->nosehoover_vxi.data(), MassQ); + break; + case etcVRESCALE: + vrescale_tcoupl(inputrec, step, ekind, dttc, state->therm_integral.data()); + break; + } + /* rescale in place here */ + if (EI_VV(inputrec->eI)) + { + rescale_velocities(ekind, md, 0, md->homenr, state->v.rvec_array()); + } + } + else + { + // Set the T scaling lambda to 1 to have no scaling + // TODO: Do we have to do it on every non-t-couple step? + for (int i = 0; (i < inputrec->opts.ngtc); i++) + { + ekind->tcstat[i].lambda = 1.0; + } + } +} + +void update_pcouple_before_coordinates(FILE* fplog, + int64_t step, + const t_inputrec* inputrec, + t_state* state, + matrix parrinellorahmanMu, + matrix M, + gmx_bool bInitStep) +{ + /* Berendsen P-coupling is completely handled after the coordinate update. + * Trotter P-coupling is handled by separate calls to trotter_update(). + */ + if (inputrec->epc == epcPARRINELLORAHMAN + && do_per_step(step + inputrec->nstpcouple - 1, inputrec->nstpcouple)) + { + real dtpc = inputrec->nstpcouple * inputrec->delta_t; + + parrinellorahman_pcoupl(fplog, step, inputrec, dtpc, state->pres_prev, state->box, + state->box_rel, state->boxv, M, parrinellorahmanMu, bInitStep); + } +} + +void update_pcouple_after_coordinates(FILE* fplog, + int64_t step, + const t_inputrec* inputrec, + const t_mdatoms* md, + const matrix pressure, + const matrix forceVirial, + const matrix constraintVirial, + matrix pressureCouplingMu, + t_state* state, + t_nrnb* nrnb, + gmx::BoxDeformation* boxDeformation, + const bool scaleCoordinates) +{ + int start = 0; + int homenr = md->homenr; + + /* Cast to real for faster code, no loss in precision (see comment above) */ + real dt = inputrec->delta_t; + + + /* now update boxes */ + switch (inputrec->epc) + { + case (epcNO): break; + case (epcBERENDSEN): + if (do_per_step(step, inputrec->nstpcouple)) + { + real dtpc = inputrec->nstpcouple * dt; + berendsen_pcoupl(fplog, step, inputrec, dtpc, pressure, state->box, forceVirial, + constraintVirial, pressureCouplingMu, &state->baros_integral); + berendsen_pscale(inputrec, pressureCouplingMu, state->box, state->box_rel, start, + homenr, state->x.rvec_array(), md->cFREEZE, nrnb, scaleCoordinates); + } + break; + case (epcPARRINELLORAHMAN): + if (do_per_step(step + inputrec->nstpcouple - 1, inputrec->nstpcouple)) + { + /* The box velocities were updated in do_pr_pcoupl, + * but we dont change the box vectors until we get here + * since we need to be able to shift/unshift above. + */ + real dtpc = inputrec->nstpcouple * dt; + for (int i = 0; i < DIM; i++) + { + for (int m = 0; m <= i; m++) + { + state->box[i][m] += dtpc * state->boxv[i][m]; + } + } + preserve_box_shape(inputrec, state->box_rel, state->box); + + /* Scale the coordinates */ + if (scaleCoordinates) + { + auto x = state->x.rvec_array(); + for (int n = start; n < start + homenr; n++) + { + tmvmul_ur0(pressureCouplingMu, x[n], x[n]); + } + } + } + break; + case (epcMTTK): + switch (inputrec->epct) + { + case (epctISOTROPIC): + /* DIM * eta = ln V. so DIM*eta_new = DIM*eta_old + DIM*dt*veta => + ln V_new = ln V_old + 3*dt*veta => V_new = V_old*exp(3*dt*veta) => + Side length scales as exp(veta*dt) */ + + msmul(state->box, std::exp(state->veta * dt), state->box); + + /* Relate veta to boxv. veta = d(eta)/dT = (1/DIM)*1/V dV/dT. + o If we assume isotropic scaling, and box length scaling + factor L, then V = L^DIM (det(M)). So dV/dt = DIM + L^(DIM-1) dL/dt det(M), and veta = (1/L) dL/dt. The + determinant of B is L^DIM det(M), and the determinant + of dB/dt is (dL/dT)^DIM det (M). veta will be + (det(dB/dT)/det(B))^(1/3). Then since M = + B_new*(vol_new)^(1/3), dB/dT_new = (veta_new)*B(new). */ + + msmul(state->box, state->veta, state->boxv); + break; + default: break; + } + break; + default: break; + } + + if (boxDeformation) + { + auto localX = makeArrayRef(state->x).subArray(start, homenr); + boxDeformation->apply(localX, state->box, step); + } +} + +extern gmx_bool update_randomize_velocities(const t_inputrec* ir, + int64_t step, + const t_commrec* cr, + const t_mdatoms* md, + gmx::ArrayRef v, + const gmx::Update* upd, + const gmx::Constraints* constr) +{ + + real rate = (ir->delta_t) / ir->opts.tau_t[0]; + + if (ir->etc == etcANDERSEN && constr != nullptr) + { + /* Currently, Andersen thermostat does not support constrained + systems. Functionality exists in the andersen_tcoupl + function in GROMACS 4.5.7 to allow this combination. That + code could be ported to the current random-number + generation approach, but has not yet been done because of + lack of time and resources. */ + 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, gmx::roundToInt(1.0 / rate))) + { + andersen_tcoupl(ir, step, cr, md, v, rate, upd->getAndersenRandomizeGroup(), + upd->getBoltzmanFactor()); + return TRUE; + } + return FALSE; +} + /* static const double sy_const[MAX_SUZUKI_YOSHIDA_NUM+1][MAX_SUZUKI_YOSHIDA_NUM+1] = { {}, diff --git a/src/gromacs/mdlib/coupling.h b/src/gromacs/mdlib/coupling.h new file mode 100644 index 0000000000..37c0166554 --- /dev/null +++ b/src/gromacs/mdlib/coupling.h @@ -0,0 +1,231 @@ +/* + * This file is part of the GROMACS molecular simulation package. + * + * Copyright (c) 1991-2000, University of Groningen, The Netherlands. + * Copyright (c) 2001-2004, The GROMACS development team. + * Copyright (c) 2013,2014,2015,2016,2017 by the GROMACS development team. + * Copyright (c) 2018,2019,2020, 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. + * + * 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. + * + * 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. + * + * 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. + * + * To help us fund GROMACS development, we humbly ask that you cite + * the research papers on the package. Check out http://www.gromacs.org. + */ +#ifndef GMX_MDLIB_COUPLING_H +#define GMX_MDLIB_COUPLING_H + +#include "gromacs/math/paddedvector.h" +#include "gromacs/math/vectypes.h" +#include "gromacs/mdtypes/md_enums.h" +#include "gromacs/utility/arrayref.h" +#include "gromacs/utility/basedefinitions.h" +#include "gromacs/utility/real.h" + +struct gmx_ekindata_t; +struct gmx_enerdata_t; +struct t_commrec; +struct t_extmass; +struct t_grpopts; +struct t_inputrec; +struct t_mdatoms; +struct t_nrnb; +class t_state; + +enum class PbcType; + +namespace gmx +{ +class BoxDeformation; +class Constraints; +class Update; +}; // namespace gmx + +/* Update the size of per-atom arrays (e.g. after DD re-partitioning, + which might increase the number of home atoms). */ + +void update_tcouple(int64_t step, + const t_inputrec* inputrec, + t_state* state, + gmx_ekindata_t* ekind, + const t_extmass* MassQ, + const t_mdatoms* md); + +/* Update Parrinello-Rahman, to be called before the coordinate update */ +void update_pcouple_before_coordinates(FILE* fplog, + int64_t step, + const t_inputrec* inputrec, + t_state* state, + matrix parrinellorahmanMu, + matrix M, + gmx_bool bInitStep); + +/* Update the box, to be called after the coordinate update. + * For Berendsen P-coupling, also calculates the scaling factor + * and scales the coordinates. + * When the deform option is used, scales coordinates and box here. + */ +void update_pcouple_after_coordinates(FILE* fplog, + int64_t step, + const t_inputrec* inputrec, + const t_mdatoms* md, + const matrix pressure, + const matrix forceVirial, + const matrix constraintVirial, + matrix pressureCouplingMu, + t_state* state, + t_nrnb* nrnb, + gmx::BoxDeformation* boxDeformation, + bool scaleCoordinates); + +/* Return TRUE if OK, FALSE in case of Shake Error */ + +extern gmx_bool update_randomize_velocities(const t_inputrec* ir, + int64_t step, + const t_commrec* cr, + const t_mdatoms* md, + gmx::ArrayRef v, + const gmx::Update* upd, + const gmx::Constraints* constr); + +void berendsen_tcoupl(const t_inputrec* ir, + gmx_ekindata_t* ekind, + real dt, + std::vector& therm_integral); //NOLINT(google-runtime-references) + +void andersen_tcoupl(const t_inputrec* ir, + int64_t step, + const t_commrec* cr, + const t_mdatoms* md, + gmx::ArrayRef v, + real rate, + const std::vector& randomize, + gmx::ArrayRef boltzfac); + +void nosehoover_tcoupl(const t_grpopts* opts, + const gmx_ekindata_t* ekind, + real dt, + double xi[], + double vxi[], + const t_extmass* MassQ); + +void trotter_update(const t_inputrec* ir, + int64_t step, + gmx_ekindata_t* ekind, + const gmx_enerdata_t* enerd, + t_state* state, + const tensor vir, + const t_mdatoms* md, + const t_extmass* MassQ, + gmx::ArrayRef> trotter_seqlist, + int trotter_seqno); + +std::array, ettTSEQMAX> +init_npt_vars(const t_inputrec* ir, t_state* state, t_extmass* Mass, gmx_bool bTrotter); + +real NPT_energy(const t_inputrec* ir, const t_state* state, const t_extmass* MassQ); +/* computes all the pressure/tempertature control energy terms to get a conserved energy */ + +void vrescale_tcoupl(const t_inputrec* ir, int64_t step, gmx_ekindata_t* ekind, real dt, double therm_integral[]); +/* Compute temperature scaling. For V-rescale it is done in update. */ + +void rescale_velocities(const gmx_ekindata_t* ekind, const t_mdatoms* mdatoms, int start, int end, rvec v[]); +/* Rescale the velocities with the scaling factor in ekind */ + +//! Check whether we do simulated annealing. +bool doSimulatedAnnealing(const t_inputrec* ir); + +//! Initialize simulated annealing. +bool initSimulatedAnnealing(t_inputrec* ir, gmx::Update* upd); + +// TODO: This is the only function in update.h altering the inputrec +void update_annealing_target_temp(t_inputrec* ir, real t, gmx::Update* upd); +/* Set reference temp for simulated annealing at time t*/ + +real calc_temp(real ekin, real nrdf); +/* Calculate the temperature */ + +real calc_pres(PbcType pbcType, int nwall, const matrix box, const tensor ekin, const tensor vir, tensor pres); +/* Calculate the pressure tensor, returns the scalar pressure. + * The unit of pressure is bar. + */ + +void parrinellorahman_pcoupl(FILE* fplog, + int64_t step, + const t_inputrec* ir, + real dt, + const tensor pres, + const tensor box, + tensor box_rel, + tensor boxv, + tensor M, + matrix mu, + gmx_bool bFirstStep); + +void berendsen_pcoupl(FILE* fplog, + int64_t step, + const t_inputrec* ir, + real dt, + const tensor pres, + const matrix box, + const matrix force_vir, + const matrix constraint_vir, + matrix mu, + double* baros_integral); + +void berendsen_pscale(const t_inputrec* ir, + const matrix mu, + matrix box, + matrix box_rel, + int start, + int nr_atoms, + rvec x[], + const unsigned short cFREEZE[], + t_nrnb* nrnb, + bool scaleCoordinates); + +void pleaseCiteCouplingAlgorithms(FILE* fplog, const t_inputrec& ir); + +/*! \brief Generate a new kinetic energy for the v-rescale thermostat + * + * Generates a new value for the kinetic energy, according to + * Bussi et al JCP (2007), Eq. (A7) + * + * This is used by update_tcoupl(), and by the VRescaleThermostat of the modular + * simulator. + * \todo Move this to the VRescaleThermostat once the modular simulator becomes + * the default code path. + * + * \param[in] kk present value of the kinetic energy of the atoms to be thermalized (in + * arbitrary units) \param[in] sigma target average value of the kinetic energy (ndeg k_b T/2) (in + * the same units as kk) \param[in] ndeg number of degrees of freedom of the atoms to be + * thermalized \param[in] taut relaxation time of the thermostat, in units of 'how often this + * routine is called' \param[in] step the time step this routine is called on \param[in] seed the + * random number generator seed \return the new kinetic energy + */ +real vrescale_resamplekin(real kk, real sigma, real ndeg, real taut, int64_t step, int64_t seed); + +#endif // GMX_MDLIB_COUPLING_H diff --git a/src/gromacs/mdlib/md_support.cpp b/src/gromacs/mdlib/md_support.cpp index 23b08394db..8803c5a99b 100644 --- a/src/gromacs/mdlib/md_support.cpp +++ b/src/gromacs/mdlib/md_support.cpp @@ -49,6 +49,7 @@ #include "gromacs/gmxlib/network.h" #include "gromacs/gmxlib/nrnb.h" #include "gromacs/math/vec.h" +#include "gromacs/mdlib/coupling.h" #include "gromacs/mdlib/dispersioncorrection.h" #include "gromacs/mdlib/gmx_omp_nthreads.h" #include "gromacs/mdlib/simulationsignal.h" diff --git a/src/gromacs/mdlib/tgroup.cpp b/src/gromacs/mdlib/tgroup.cpp index 1591ef9861..3505a1a1ea 100644 --- a/src/gromacs/mdlib/tgroup.cpp +++ b/src/gromacs/mdlib/tgroup.cpp @@ -44,9 +44,9 @@ #include "gromacs/gmxlib/network.h" #include "gromacs/math/vec.h" +#include "gromacs/mdlib/coupling.h" #include "gromacs/mdlib/gmx_omp_nthreads.h" #include "gromacs/mdlib/rbin.h" -#include "gromacs/mdlib/update.h" #include "gromacs/mdtypes/group.h" #include "gromacs/mdtypes/inputrec.h" #include "gromacs/mdtypes/mdatom.h" diff --git a/src/gromacs/mdlib/update.cpp b/src/gromacs/mdlib/update.cpp index 9c0c659533..9b47cf8ca3 100644 --- a/src/gromacs/mdlib/update.cpp +++ b/src/gromacs/mdlib/update.cpp @@ -1262,75 +1262,6 @@ void restore_ekinstate_from_state(const t_commrec* cr, gmx_ekindata_t* ekind, co } } -void update_tcouple(int64_t step, - const t_inputrec* inputrec, - t_state* state, - gmx_ekindata_t* ekind, - const t_extmass* MassQ, - const t_mdatoms* md) - -{ - // This condition was explicitly checked in previous version, but should have never been satisfied - GMX_ASSERT(!(EI_VV(inputrec->eI) - && (inputrecNvtTrotter(inputrec) || inputrecNptTrotter(inputrec) - || inputrecNphTrotter(inputrec))), - "Temperature coupling was requested with velocity verlet and trotter"); - - bool doTemperatureCoupling = false; - - // For VV temperature coupling parameters are updated on the current - // step, for the others - one step before. - if (inputrec->etc == etcNO) - { - doTemperatureCoupling = false; - } - else if (EI_VV(inputrec->eI)) - { - doTemperatureCoupling = do_per_step(step, inputrec->nsttcouple); - } - else - { - doTemperatureCoupling = do_per_step(step + inputrec->nsttcouple - 1, inputrec->nsttcouple); - } - - if (doTemperatureCoupling) - { - real dttc = inputrec->nsttcouple * inputrec->delta_t; - - // TODO: berendsen_tcoupl(...), nosehoover_tcoupl(...) and vrescale_tcoupl(...) update - // temperature coupling parameters, which should be reflected in the name of these - // subroutines - switch (inputrec->etc) - { - case etcNO: break; - case etcBERENDSEN: - berendsen_tcoupl(inputrec, ekind, dttc, state->therm_integral); - break; - case etcNOSEHOOVER: - nosehoover_tcoupl(&(inputrec->opts), ekind, dttc, state->nosehoover_xi.data(), - state->nosehoover_vxi.data(), MassQ); - break; - case etcVRESCALE: - vrescale_tcoupl(inputrec, step, ekind, dttc, state->therm_integral.data()); - break; - } - /* rescale in place here */ - if (EI_VV(inputrec->eI)) - { - rescale_velocities(ekind, md, 0, md->homenr, state->v.rvec_array()); - } - } - else - { - // Set the T scaling lambda to 1 to have no scaling - // TODO: Do we have to do it on every non-t-couple step? - for (int i = 0; (i < inputrec->opts.ngtc); i++) - { - ekind->tcstat[i].lambda = 1.0; - } - } -} - void getThreadAtomRange(int numThreads, int threadIndex, int numAtoms, int* startAtom, int* endAtom) { #if GMX_HAVE_SIMD_UPDATE @@ -1348,27 +1279,6 @@ void getThreadAtomRange(int numThreads, int threadIndex, int numAtoms, int* star } } -void update_pcouple_before_coordinates(FILE* fplog, - int64_t step, - const t_inputrec* inputrec, - t_state* state, - matrix parrinellorahmanMu, - matrix M, - gmx_bool bInitStep) -{ - /* Berendsen P-coupling is completely handled after the coordinate update. - * Trotter P-coupling is handled by separate calls to trotter_update(). - */ - if (inputrec->epc == epcPARRINELLORAHMAN - && do_per_step(step + inputrec->nstpcouple - 1, inputrec->nstpcouple)) - { - real dtpc = inputrec->nstpcouple * inputrec->delta_t; - - parrinellorahman_pcoupl(fplog, step, inputrec, dtpc, state->pres_prev, state->box, - state->box_rel, state->boxv, M, parrinellorahmanMu, bInitStep); - } -} - void Update::Impl::update_sd_second_half(const t_inputrec& inputRecord, int64_t step, real* dvdlambda, @@ -1484,102 +1394,6 @@ void Update::Impl::finish_update(const t_inputrec& inputRecord, wallcycle_stop(wcycle, ewcUPDATE); } -void update_pcouple_after_coordinates(FILE* fplog, - int64_t step, - const t_inputrec* inputrec, - const t_mdatoms* md, - const matrix pressure, - const matrix forceVirial, - const matrix constraintVirial, - matrix pressureCouplingMu, - t_state* state, - t_nrnb* nrnb, - gmx::BoxDeformation* boxDeformation, - const bool scaleCoordinates) -{ - int start = 0; - int homenr = md->homenr; - - /* Cast to real for faster code, no loss in precision (see comment above) */ - real dt = inputrec->delta_t; - - - /* now update boxes */ - switch (inputrec->epc) - { - case (epcNO): break; - case (epcBERENDSEN): - if (do_per_step(step, inputrec->nstpcouple)) - { - real dtpc = inputrec->nstpcouple * dt; - berendsen_pcoupl(fplog, step, inputrec, dtpc, pressure, state->box, forceVirial, - constraintVirial, pressureCouplingMu, &state->baros_integral); - berendsen_pscale(inputrec, pressureCouplingMu, state->box, state->box_rel, start, - homenr, state->x.rvec_array(), md->cFREEZE, nrnb, scaleCoordinates); - } - break; - case (epcPARRINELLORAHMAN): - if (do_per_step(step + inputrec->nstpcouple - 1, inputrec->nstpcouple)) - { - /* The box velocities were updated in do_pr_pcoupl, - * but we dont change the box vectors until we get here - * since we need to be able to shift/unshift above. - */ - real dtpc = inputrec->nstpcouple * dt; - for (int i = 0; i < DIM; i++) - { - for (int m = 0; m <= i; m++) - { - state->box[i][m] += dtpc * state->boxv[i][m]; - } - } - preserve_box_shape(inputrec, state->box_rel, state->box); - - /* Scale the coordinates */ - if (scaleCoordinates) - { - auto x = state->x.rvec_array(); - for (int n = start; n < start + homenr; n++) - { - tmvmul_ur0(pressureCouplingMu, x[n], x[n]); - } - } - } - break; - case (epcMTTK): - switch (inputrec->epct) - { - case (epctISOTROPIC): - /* DIM * eta = ln V. so DIM*eta_new = DIM*eta_old + DIM*dt*veta => - ln V_new = ln V_old + 3*dt*veta => V_new = V_old*exp(3*dt*veta) => - Side length scales as exp(veta*dt) */ - - msmul(state->box, std::exp(state->veta * dt), state->box); - - /* Relate veta to boxv. veta = d(eta)/dT = (1/DIM)*1/V dV/dT. - o If we assume isotropic scaling, and box length scaling - factor L, then V = L^DIM (det(M)). So dV/dt = DIM - L^(DIM-1) dL/dt det(M), and veta = (1/L) dL/dt. The - determinant of B is L^DIM det(M), and the determinant - of dB/dt is (dL/dT)^DIM det (M). veta will be - (det(dB/dT)/det(B))^(1/3). Then since M = - B_new*(vol_new)^(1/3), dB/dT_new = (veta_new)*B(new). */ - - msmul(state->box, state->veta, state->boxv); - break; - default: break; - } - break; - default: break; - } - - if (boxDeformation) - { - auto localX = makeArrayRef(state->x).subArray(start, homenr); - boxDeformation->apply(localX, state->box, step); - } -} - void Update::Impl::update_coords(const t_inputrec& inputRecord, int64_t step, const t_mdatoms* md, @@ -1679,38 +1493,3 @@ void Update::Impl::update_coords(const t_inputrec& GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR } } - -extern gmx_bool update_randomize_velocities(const t_inputrec* ir, - int64_t step, - const t_commrec* cr, - const t_mdatoms* md, - gmx::ArrayRef v, - const Update* upd, - const gmx::Constraints* constr) -{ - - real rate = (ir->delta_t) / ir->opts.tau_t[0]; - - if (ir->etc == etcANDERSEN && constr != nullptr) - { - /* Currently, Andersen thermostat does not support constrained - systems. Functionality exists in the andersen_tcoupl - function in GROMACS 4.5.7 to allow this combination. That - code could be ported to the current random-number - generation approach, but has not yet been done because of - lack of time and resources. */ - 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, roundToInt(1.0 / rate))) - { - andersen_tcoupl(ir, step, cr, md, v, rate, upd->getAndersenRandomizeGroup(), - upd->getBoltzmanFactor()); - return TRUE; - } - return FALSE; -} diff --git a/src/gromacs/mdlib/update.h b/src/gromacs/mdlib/update.h index 7cc6207625..b7853f3902 100644 --- a/src/gromacs/mdlib/update.h +++ b/src/gromacs/mdlib/update.h @@ -51,7 +51,6 @@ class ekinstate_t; struct gmx_ekindata_t; struct gmx_enerdata_t; enum class PbcType; -struct t_extmass; struct t_fcdata; struct t_graph; struct t_grpopts; @@ -199,54 +198,6 @@ private: }; // namespace gmx - -/* Update the size of per-atom arrays (e.g. after DD re-partitioning, - which might increase the number of home atoms). */ - -void update_tcouple(int64_t step, - const t_inputrec* inputrec, - t_state* state, - gmx_ekindata_t* ekind, - const t_extmass* MassQ, - const t_mdatoms* md); - -/* Update Parrinello-Rahman, to be called before the coordinate update */ -void update_pcouple_before_coordinates(FILE* fplog, - int64_t step, - const t_inputrec* inputrec, - t_state* state, - matrix parrinellorahmanMu, - matrix M, - gmx_bool bInitStep); - -/* Update the box, to be called after the coordinate update. - * For Berendsen P-coupling, also calculates the scaling factor - * and scales the coordinates. - * When the deform option is used, scales coordinates and box here. - */ -void update_pcouple_after_coordinates(FILE* fplog, - int64_t step, - const t_inputrec* inputrec, - const t_mdatoms* md, - const matrix pressure, - const matrix forceVirial, - const matrix constraintVirial, - matrix pressureCouplingMu, - t_state* state, - t_nrnb* nrnb, - gmx::BoxDeformation* boxDeformation, - bool scaleCoordinates); - -/* Return TRUE if OK, FALSE in case of Shake Error */ - -extern gmx_bool update_randomize_velocities(const t_inputrec* ir, - int64_t step, - const t_commrec* cr, - const t_mdatoms* md, - gmx::ArrayRef v, - const gmx::Update* upd, - const gmx::Constraints* constr); - /* * Compute the partial kinetic energy for home particles; * will be accumulated in the calling routine. @@ -272,104 +223,6 @@ void update_ekinstate(ekinstate_t* ekinstate, const gmx_ekindata_t* ekind); to the rest of the simulation */ void restore_ekinstate_from_state(const t_commrec* cr, gmx_ekindata_t* ekind, const ekinstate_t* ekinstate); -void berendsen_tcoupl(const t_inputrec* ir, - gmx_ekindata_t* ekind, - real dt, - std::vector& therm_integral); //NOLINT(google-runtime-references) - -void andersen_tcoupl(const t_inputrec* ir, - int64_t step, - const t_commrec* cr, - const t_mdatoms* md, - gmx::ArrayRef v, - real rate, - const std::vector& randomize, - gmx::ArrayRef boltzfac); - -void nosehoover_tcoupl(const t_grpopts* opts, - const gmx_ekindata_t* ekind, - real dt, - double xi[], - double vxi[], - const t_extmass* MassQ); - -void trotter_update(const t_inputrec* ir, - int64_t step, - gmx_ekindata_t* ekind, - const gmx_enerdata_t* enerd, - t_state* state, - const tensor vir, - const t_mdatoms* md, - const t_extmass* MassQ, - gmx::ArrayRef> trotter_seqlist, - int trotter_seqno); - -std::array, ettTSEQMAX> -init_npt_vars(const t_inputrec* ir, t_state* state, t_extmass* Mass, gmx_bool bTrotter); - -real NPT_energy(const t_inputrec* ir, const t_state* state, const t_extmass* MassQ); -/* computes all the pressure/tempertature control energy terms to get a conserved energy */ - -void vrescale_tcoupl(const t_inputrec* ir, int64_t step, gmx_ekindata_t* ekind, real dt, double therm_integral[]); -/* Compute temperature scaling. For V-rescale it is done in update. */ - -void rescale_velocities(const gmx_ekindata_t* ekind, const t_mdatoms* mdatoms, int start, int end, rvec v[]); -/* Rescale the velocities with the scaling factor in ekind */ - -//! Check whether we do simulated annealing. -bool doSimulatedAnnealing(const t_inputrec* ir); - -//! Initialize simulated annealing. -bool initSimulatedAnnealing(t_inputrec* ir, gmx::Update* upd); - -// TODO: This is the only function in update.h altering the inputrec -void update_annealing_target_temp(t_inputrec* ir, real t, gmx::Update* upd); -/* Set reference temp for simulated annealing at time t*/ - -real calc_temp(real ekin, real nrdf); -/* Calculate the temperature */ - -real calc_pres(PbcType pbcType, int nwall, const matrix box, const tensor ekin, const tensor vir, tensor pres); -/* Calculate the pressure tensor, returns the scalar pressure. - * The unit of pressure is bar. - */ - -void parrinellorahman_pcoupl(FILE* fplog, - int64_t step, - const t_inputrec* ir, - real dt, - const tensor pres, - const tensor box, - tensor box_rel, - tensor boxv, - tensor M, - matrix mu, - gmx_bool bFirstStep); - -void berendsen_pcoupl(FILE* fplog, - int64_t step, - const t_inputrec* ir, - real dt, - const tensor pres, - const matrix box, - const matrix force_vir, - const matrix constraint_vir, - matrix mu, - double* baros_integral); - -void berendsen_pscale(const t_inputrec* ir, - const matrix mu, - matrix box, - matrix box_rel, - int start, - int nr_atoms, - rvec x[], - const unsigned short cFREEZE[], - t_nrnb* nrnb, - bool scaleCoordinates); - -void pleaseCiteCouplingAlgorithms(FILE* fplog, const t_inputrec& ir); - /*! \brief Computes the atom range for a thread to operate on, ensuring SIMD aligned ranges * * \param[in] numThreads The number of threads to divide atoms over @@ -380,24 +233,4 @@ void pleaseCiteCouplingAlgorithms(FILE* fplog, const t_inputrec& ir); */ void getThreadAtomRange(int numThreads, int threadIndex, int numAtoms, int* startAtom, int* endAtom); -/*! \brief Generate a new kinetic energy for the v-rescale thermostat - * - * Generates a new value for the kinetic energy, according to - * Bussi et al JCP (2007), Eq. (A7) - * - * This is used by update_tcoupl(), and by the VRescaleThermostat of the modular - * simulator. - * TODO: Move this to the VRescaleThermostat once the modular simulator becomes - * the default code path. - * - * @param kk present value of the kinetic energy of the atoms to be thermalized (in arbitrary units) - * @param sigma target average value of the kinetic energy (ndeg k_b T/2) (in the same units as kk) - * @param ndeg number of degrees of freedom of the atoms to be thermalized - * @param taut relaxation time of the thermostat, in units of 'how often this routine is called' - * @param step the time step this routine is called on - * @param seed the random number generator seed - * @return the new kinetic energy - */ -real vrescale_resamplekin(real kk, real sigma, real ndeg, real taut, int64_t step, int64_t seed); - #endif diff --git a/src/gromacs/mdrun/md.cpp b/src/gromacs/mdrun/md.cpp index 9969a9720e..b2eb8fa378 100644 --- a/src/gromacs/mdrun/md.cpp +++ b/src/gromacs/mdrun/md.cpp @@ -78,6 +78,7 @@ #include "gromacs/mdlib/checkpointhandler.h" #include "gromacs/mdlib/compute_io.h" #include "gromacs/mdlib/constr.h" +#include "gromacs/mdlib/coupling.h" #include "gromacs/mdlib/ebin.h" #include "gromacs/mdlib/enerdata_utils.h" #include "gromacs/mdlib/energyoutput.h" diff --git a/src/gromacs/mdrun/minimize.cpp b/src/gromacs/mdrun/minimize.cpp index b13d54070f..ddd0870ba6 100644 --- a/src/gromacs/mdrun/minimize.cpp +++ b/src/gromacs/mdrun/minimize.cpp @@ -72,6 +72,7 @@ #include "gromacs/math/functions.h" #include "gromacs/math/vec.h" #include "gromacs/mdlib/constr.h" +#include "gromacs/mdlib/coupling.h" #include "gromacs/mdlib/dispersioncorrection.h" #include "gromacs/mdlib/ebin.h" #include "gromacs/mdlib/enerdata_utils.h" diff --git a/src/gromacs/modularsimulator/energyelement.cpp b/src/gromacs/modularsimulator/energyelement.cpp index bec2d3eb89..5e8abff78e 100644 --- a/src/gromacs/modularsimulator/energyelement.cpp +++ b/src/gromacs/modularsimulator/energyelement.cpp @@ -45,6 +45,7 @@ #include "gromacs/math/vec.h" #include "gromacs/mdlib/compute_io.h" +#include "gromacs/mdlib/coupling.h" #include "gromacs/mdlib/enerdata_utils.h" #include "gromacs/mdlib/energyoutput.h" #include "gromacs/mdlib/mdatoms.h" diff --git a/src/gromacs/modularsimulator/modularsimulator.cpp b/src/gromacs/modularsimulator/modularsimulator.cpp index d4bf961841..fa05c5f550 100644 --- a/src/gromacs/modularsimulator/modularsimulator.cpp +++ b/src/gromacs/modularsimulator/modularsimulator.cpp @@ -53,6 +53,7 @@ #include "gromacs/math/vec.h" #include "gromacs/mdlib/checkpointhandler.h" #include "gromacs/mdlib/constr.h" +#include "gromacs/mdlib/coupling.h" #include "gromacs/mdlib/energyoutput.h" #include "gromacs/mdlib/forcerec.h" #include "gromacs/mdlib/mdatoms.h" diff --git a/src/gromacs/modularsimulator/parrinellorahmanbarostat.cpp b/src/gromacs/modularsimulator/parrinellorahmanbarostat.cpp index 21ffd812bc..e1dda4b99c 100644 --- a/src/gromacs/modularsimulator/parrinellorahmanbarostat.cpp +++ b/src/gromacs/modularsimulator/parrinellorahmanbarostat.cpp @@ -45,9 +45,9 @@ #include "gromacs/domdec/domdec_network.h" #include "gromacs/math/vec.h" +#include "gromacs/mdlib/coupling.h" #include "gromacs/mdlib/mdatoms.h" #include "gromacs/mdlib/stat.h" -#include "gromacs/mdlib/update.h" #include "gromacs/mdtypes/commrec.h" #include "gromacs/mdtypes/inputrec.h" #include "gromacs/mdtypes/mdatom.h" diff --git a/src/gromacs/modularsimulator/vrescalethermostat.cpp b/src/gromacs/modularsimulator/vrescalethermostat.cpp index be1c4bcff8..9aa43be58d 100644 --- a/src/gromacs/modularsimulator/vrescalethermostat.cpp +++ b/src/gromacs/modularsimulator/vrescalethermostat.cpp @@ -1,7 +1,7 @@ /* * This file is part of the GROMACS molecular simulation package. * - * Copyright (c) 2019, by the GROMACS development team, led by + * Copyright (c) 2019,2020, 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. @@ -46,8 +46,8 @@ #include "gromacs/domdec/domdec_network.h" #include "gromacs/math/units.h" #include "gromacs/math/vec.h" +#include "gromacs/mdlib/coupling.h" #include "gromacs/mdlib/stat.h" -#include "gromacs/mdlib/update.h" #include "gromacs/mdtypes/commrec.h" #include "gromacs/mdtypes/group.h" #include "gromacs/mdtypes/state.h" -- 2.22.0