*/
#include "gmxpre.h"
+#include "coupling.h"
+
#include <cassert>
#include <cmath>
#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"
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<gmx::RVec> 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] = {
{},
--- /dev/null
+/*
+ * 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<gmx::RVec> v,
+ const gmx::Update* upd,
+ const gmx::Constraints* constr);
+
+void berendsen_tcoupl(const t_inputrec* ir,
+ gmx_ekindata_t* ekind,
+ real dt,
+ std::vector<double>& 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<gmx::RVec> v,
+ real rate,
+ const std::vector<bool>& randomize,
+ gmx::ArrayRef<const real> 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<std::vector<int>> trotter_seqlist,
+ int trotter_seqno);
+
+std::array<std::vector<int>, 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
#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"
#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"
}
}
-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
}
}
-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,
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,
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<gmx::RVec> 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;
-}
struct gmx_ekindata_t;
struct gmx_enerdata_t;
enum class PbcType;
-struct t_extmass;
struct t_fcdata;
struct t_graph;
struct t_grpopts;
}; // 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<gmx::RVec> v,
- const gmx::Update* upd,
- const gmx::Constraints* constr);
-
/*
* Compute the partial kinetic energy for home particles;
* will be accumulated in the calling routine.
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<double>& 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<gmx::RVec> v,
- real rate,
- const std::vector<bool>& randomize,
- gmx::ArrayRef<const real> 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<std::vector<int>> trotter_seqlist,
- int trotter_seqno);
-
-std::array<std::vector<int>, 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
*/
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
#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"
#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"
#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"
#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"
#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"
/*
* 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.
#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"
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff