if (lincsd != nullptr)
{
- bOK = constrain_lincs(bLog || bEner, ir, step, lincsd, md, cr, ms, x, xprime, min_proj, box,
- pbc_null, lambda, dvdlambda, invdt, v.unpaddedArrayRef(),
- vir != nullptr, vir_r_m_dr, econq, nrnb, maxwarn, &warncount_lincs);
+ bOK = constrain_lincs(bLog || bEner, ir, step, lincsd, md.invmass, cr, ms, x, xprime,
+ min_proj, box, pbc_null, md.nMassPerturbed != 0, lambda, dvdlambda,
+ invdt, v.unpaddedArrayRef(), vir != nullptr, vir_r_m_dr, econq, nrnb,
+ maxwarn, &warncount_lincs);
if (!bOK && maxwarn < INT_MAX)
{
if (log != nullptr)
}
return bOK;
-}
+} // namespace gmx
ArrayRef<real> Constraints::rmsdData() const
{
*/
if (ir.eConstrAlg == econtLINCS)
{
- set_lincs(*idef, md, EI_DYNAMICS(ir.eI), cr, lincsd);
+ set_lincs(*idef, md.nr, md.invmass, md.lambda, EI_DYNAMICS(ir.eI), cr, lincsd);
}
if (ir.eConstrAlg == econtSHAKE)
{
}
else
{
- make_shake_sblock_serial(shaked.get(), &top->idef, md);
+ make_shake_sblock_serial(shaked.get(), &top->idef, md.nr);
}
}
}
#include "gromacs/mdtypes/commrec.h"
#include "gromacs/mdtypes/inputrec.h"
#include "gromacs/mdtypes/md_enums.h"
-#include "gromacs/mdtypes/mdatom.h"
#include "gromacs/pbcutil/pbc.h"
#include "gromacs/pbcutil/pbc_simd.h"
#include "gromacs/simd/simd.h"
t_pbc* pbc,
Lincs* lincsd,
int th,
- real* invmass,
+ const real* invmass,
ConstraintVariable econq,
bool bCalcDHDL,
bool bCalcVir,
}
/*! \brief Sets the elements in the LINCS matrix. */
-static void set_lincs_matrix(Lincs* li, real* invmass, real lambda)
+static void set_lincs_matrix(Lincs* li, const real* invmass, real lambda)
{
const real invsqrt2 = 0.7071067811865475244;
}
}
-void set_lincs(const InteractionDefinitions& idef, const t_mdatoms& md, bool bDynamics, const t_commrec* cr, Lincs* li)
+void set_lincs(const InteractionDefinitions& idef,
+ const int numAtoms,
+ const real* invmass,
+ const real lambda,
+ bool bDynamics,
+ const t_commrec* cr,
+ Lincs* li)
{
li->nc_real = 0;
li->nc = 0;
}
else
{
- natoms = md.homenr;
+ natoms = numAtoms;
}
const ListOfLists<int> at2con =
if (li->ntask > 1)
{
- lincs_thread_setup(li, md.nr);
+ lincs_thread_setup(li, numAtoms);
}
- set_lincs_matrix(li, md.invmass, md.lambda);
+ set_lincs_matrix(li, invmass, lambda);
}
//! Issues a warning when LINCS constraints cannot be satisfied.
const t_inputrec& ir,
int64_t step,
Lincs* lincsd,
- const t_mdatoms& md,
+ const real* invmass,
const t_commrec* cr,
const gmx_multisim_t* ms,
ArrayRefWithPadding<const RVec> xPadded,
ArrayRef<RVec> min_proj,
const matrix box,
t_pbc* pbc,
+ const bool hasMassPerturbed,
real lambda,
real* dvdlambda,
real invdt,
*/
if (ir.efep != efepNO)
{
- if (md.nMassPerturbed && lincsd->matlam != md.lambda)
+ if (hasMassPerturbed && lincsd->matlam != lambda)
{
- set_lincs_matrix(lincsd, md.invmass, md.lambda);
+ set_lincs_matrix(lincsd, invmass, lambda);
}
for (int i = 0; i < lincsd->nc; i++)
clear_mat(lincsd->task[th].vir_r_m_dr);
- do_lincs(xPadded, xprimePadded, box, pbc, lincsd, th, md.invmass, cr, bCalcDHDL,
+ do_lincs(xPadded, xprimePadded, box, pbc, lincsd, th, invmass, cr, bCalcDHDL,
ir.LincsWarnAngle, &bWarn, invdt, v, bCalcVir,
th == 0 ? vir_r_m_dr : lincsd->task[th].vir_r_m_dr);
}
{
int th = gmx_omp_get_thread_num();
- do_lincsp(xPadded, xprimePadded, min_proj, pbc, lincsd, th, md.invmass, econq,
+ do_lincsp(xPadded, xprimePadded, min_proj, pbc, lincsd, th, invmass, econq,
bCalcDHDL, bCalcVir, th == 0 ? vir_r_m_dr : lincsd->task[th].vir_r_m_dr);
}
GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
class InteractionDefinitions;
struct t_commrec;
struct t_inputrec;
-struct t_mdatoms;
struct t_nrnb;
struct t_pbc;
/*! \brief Initialize lincs stuff */
void set_lincs(const InteractionDefinitions& idef,
- const t_mdatoms& md,
+ int numAtoms,
+ const real* invmass,
+ real lambda,
bool bDynamics,
const t_commrec* cr,
Lincs* li);
const t_inputrec& ir,
int64_t step,
Lincs* lincsd,
- const t_mdatoms& md,
+ const real* invmass,
const t_commrec* cr,
const gmx_multisim_t* ms,
ArrayRefWithPadding<const RVec> x,
ArrayRef<RVec> min_proj,
const matrix box,
t_pbc* pbc,
+ bool hasMassPerturbed,
real lambda,
real* dvdlambda,
real invdt,
return true;
}
-void LincsGpu::set(const InteractionDefinitions& idef, const t_mdatoms& md)
+void LincsGpu::set(const InteractionDefinitions& idef, const int numAtoms, const real* invmass)
{
- int numAtoms = md.nr;
// List of constrained atoms (CPU memory)
std::vector<int2> constraintsHost;
// Equilibrium distances for the constraints (CPU)
int center = c1a1;
- float sqrtmu1 = 1.0 / sqrt(md.invmass[c1a1] + md.invmass[c1a2]);
- float sqrtmu2 = 1.0 / sqrt(md.invmass[c2a1] + md.invmass[c2a2]);
+ float sqrtmu1 = 1.0 / sqrt(invmass[c1a1] + invmass[c1a2]);
+ float sqrtmu2 = 1.0 / sqrt(invmass[c2a1] + invmass[c2a2]);
- massFactorsHost.at(index) = -sign * md.invmass[center] * sqrtmu1 * sqrtmu2;
+ massFactorsHost.at(index) = -sign * invmass[center] * sqrtmu1 * sqrtmu2;
coupledConstraintsCountsHost.at(splitMap.at(c1))++;
}
int center = c1a2;
- float sqrtmu1 = 1.0 / sqrt(md.invmass[c1a1] + md.invmass[c1a2]);
- float sqrtmu2 = 1.0 / sqrt(md.invmass[c2a1] + md.invmass[c2a2]);
+ float sqrtmu1 = 1.0 / sqrt(invmass[c1a1] + invmass[c1a2]);
+ float sqrtmu2 = 1.0 / sqrt(invmass[c2a1] + invmass[c2a2]);
- massFactorsHost.at(index) = sign * md.invmass[center] * sqrtmu1 * sqrtmu2;
+ massFactorsHost.at(index) = sign * invmass[center] * sqrtmu1 * sqrtmu2;
coupledConstraintsCountsHost.at(splitMap.at(c1))++;
}
maxCoupledConstraints * kernelParams_.numConstraintsThreads, deviceStream_,
GpuApiCallBehavior::Sync, nullptr);
- GMX_RELEASE_ASSERT(md.invmass != nullptr, "Masses of atoms should be specified.\n");
- copyToDeviceBuffer(&kernelParams_.d_inverseMasses, md.invmass, 0, numAtoms, deviceStream_,
+ GMX_RELEASE_ASSERT(invmass != nullptr, "Masses of atoms should be specified.\n");
+ copyToDeviceBuffer(&kernelParams_.d_inverseMasses, invmass, 0, numAtoms, deviceStream_,
GpuApiCallBehavior::Sync, nullptr);
}
#include "gromacs/gpu_utils/device_context.h"
#include "gromacs/gpu_utils/gputraits.cuh"
#include "gromacs/mdlib/constr.h"
-#include "gromacs/mdtypes/mdatom.h"
#include "gromacs/pbcutil/pbc_aiuc.h"
#include "gromacs/utility/classhelpers.h"
* Information about constraints is taken from:
* idef.il[F_CONSTR].iatoms --- type (T) of constraint and two atom indexes (i1, i2)
* idef.iparams[T].constr.dA --- target length for constraint of type T
- * From t_mdatom, the code takes:
- * md.invmass --- array of inverse square root of masses for each atom in the system.
*
- * \param[in] idef Local topology data to get information on constraints from.
- * \param[in] md Atoms data to get atom masses from.
+ * \param[in] idef Local topology data to get information on constraints from.
+ * \param[in] numAtoms Number of atoms.
+ * \param[in] invmass Inverse masses of atoms.
*/
- void set(const InteractionDefinitions& idef, const t_mdatoms& md);
+ void set(const InteractionDefinitions& idef, int numAtoms, const real* invmass);
/*! \brief
* Returns whether the maximum number of coupled constraints is supported
#include "gromacs/mdlib/splitter.h"
#include "gromacs/mdtypes/inputrec.h"
#include "gromacs/mdtypes/md_enums.h"
-#include "gromacs/mdtypes/mdatom.h"
#include "gromacs/pbcutil/pbc.h"
#include "gromacs/topology/invblock.h"
#include "gromacs/utility/fatalerror.h"
shaked->scaled_lagrange_multiplier.resize(ncons);
}
-void make_shake_sblock_serial(shakedata* shaked, InteractionDefinitions* idef, const t_mdatoms& md)
+void make_shake_sblock_serial(shakedata* shaked, InteractionDefinitions* idef, const int numAtoms)
{
int i, m, ncons;
int bstart, bnr;
init_blocka(&sblocks);
sfree(sblocks.index); // To solve memory leak
- gen_sblocks(nullptr, md.homenr, *idef, &sblocks, FALSE);
+ gen_sblocks(nullptr, numAtoms, *idef, &sblocks, FALSE);
/*
bstart=(idef->nodeid > 0) ? blocks->multinr[idef->nodeid-1] : 0;
}
/* Calculate block number for each atom */
- inv_sblock = make_invblocka(&sblocks, md.nr);
+ inv_sblock = make_invblocka(&sblocks, numAtoms);
done_blocka(&sblocks);
struct InteractionList;
class InteractionDefinitions;
struct t_inputrec;
-struct t_mdatoms;
struct t_nrnb;
struct t_pbc;
};
//! Make SHAKE blocks when not using DD.
-void make_shake_sblock_serial(shakedata* shaked, InteractionDefinitions* idef, const t_mdatoms& md);
+void make_shake_sblock_serial(shakedata* shaked, InteractionDefinitions* idef, int numAtoms);
//! Make SHAKE blocks when using DD.
void make_shake_sblock_dd(shakedata* shaked, const InteractionList& ilcon);
ir_.delta_t = timestep;
ir_.eI = 0;
- // MD atoms data
- md_.nMassPerturbed = 0;
- md_.lambda = 0.0;
- md_.invmass = invmass_.data();
- md_.nr = numAtoms;
- md_.homenr = numAtoms;
-
// Virial evaluation
computeVirial_ = computeVirial;
if (computeVirial)
}
}
-
// Free energy evaluation
+ hasMassPerturbed_ = false;
+ lambda_ = 0.0;
compute_dHdLambda_ = compute_dHdLambda;
dHdLambda_ = 0;
if (compute_dHdLambda_)
#include "gromacs/mdlib/lincs.h"
#include "gromacs/mdlib/shake.h"
#include "gromacs/mdtypes/inputrec.h"
-#include "gromacs/mdtypes/mdatom.h"
#include "gromacs/pbcutil/pbc.h"
#include "gromacs/topology/idef.h"
#include "gromacs/topology/ifunc.h"
t_inputrec ir_;
//! Local topology
std::unique_ptr<InteractionDefinitions> idef_;
- //! MD atoms
- t_mdatoms md_;
//! Computational time array (normally used to benchmark performance)
t_nrnb nrnb_;
tensor virialScaledRef_;
//! If the free energy is computed
bool compute_dHdLambda_;
+ //! If there are atoms with perturbed mass
+ bool hasMassPerturbed_ = false;
+ //! Lambda value
+ real lambda_ = 0.0;
//! For free energy computation
real dHdLambda_;
//! For free energy computation (reference value)
#include "gromacs/mdrunutility/multisim.h"
#include "gromacs/mdtypes/commrec.h"
#include "gromacs/mdtypes/inputrec.h"
-#include "gromacs/mdtypes/mdatom.h"
#include "gromacs/pbcutil/pbc.h"
#include "gromacs/topology/idef.h"
#include "gromacs/topology/ifunc.h"
void applyShake(ConstraintsTestData* testData, t_pbc gmx_unused pbc)
{
shakedata shaked;
- make_shake_sblock_serial(&shaked, testData->idef_.get(), testData->md_);
+ make_shake_sblock_serial(&shaked, testData->idef_.get(), testData->numAtoms_);
bool success = constrain_shake(
nullptr, &shaked, testData->invmass_.data(), *testData->idef_, testData->ir_, testData->x_,
- testData->xPrime_, testData->xPrime2_, nullptr, &testData->nrnb_, testData->md_.lambda,
+ testData->xPrime_, testData->xPrime2_, nullptr, &testData->nrnb_, testData->lambda_,
&testData->dHdLambda_, testData->invdt_, testData->v_, testData->computeVirial_,
testData->virialScaled_, false, gmx::ConstraintVariable::Positions);
EXPECT_TRUE(success) << "Test failed with a false return value in SHAKE.";
// Initialize LINCS
lincsd = init_lincs(nullptr, testData->mtop_, testData->nflexcon_, at2con_mt, false,
testData->ir_.nLincsIter, testData->ir_.nProjOrder);
- set_lincs(*testData->idef_, testData->md_, EI_DYNAMICS(testData->ir_.eI), &cr, lincsd);
+ set_lincs(*testData->idef_, testData->numAtoms_, testData->invmass_.data(), testData->lambda_,
+ EI_DYNAMICS(testData->ir_.eI), &cr, lincsd);
// Evaluate constraints
bool success = constrain_lincs(
- false, testData->ir_, 0, lincsd, testData->md_, &cr, &ms,
+ false, testData->ir_, 0, lincsd, testData->invmass_.data(), &cr, &ms,
testData->x_.arrayRefWithPadding(), testData->xPrime_.arrayRefWithPadding(),
testData->xPrime2_.arrayRefWithPadding().unpaddedArrayRef(), pbc.box, &pbc,
- testData->md_.lambda, &testData->dHdLambda_, testData->invdt_,
+ testData->hasMassPerturbed_, testData->lambda_, &testData->dHdLambda_, testData->invdt_,
testData->v_.arrayRefWithPadding().unpaddedArrayRef(), testData->computeVirial_,
testData->virialScaled_, gmx::ConstraintVariable::Positions, &testData->nrnb_, maxwarn,
&warncount_lincs);
int numAtoms = testData->numAtoms_;
float3 *d_x, *d_xp, *d_v;
- lincsGpu->set(*testData->idef_, testData->md_);
+ lincsGpu->set(*testData->idef_, testData->numAtoms_, testData->invmass_.data());
PbcAiuc pbcAiuc;
setPbcAiuc(pbc.ndim_ePBC, pbc.box, &pbcAiuc);
// Integrator should also update something, but it does not even have a method yet
integrator_->set(md, numTempScaleValues, md.cTC);
- lincsGpu_->set(idef, md);
+ lincsGpu_->set(idef, md.nr, md.invmass);
settleGpu_->set(idef, md);
coordinateScalingKernelLaunchConfig_.gridSize[0] =