namespace gmx
{
-struct settleparam_t
-{
- real mO;
- real mH;
- real wh;
- real dOH;
- real dHH;
- real ra;
- real rb;
- real rc;
- real irc2;
- /* For projection */
- real imO;
- real imH;
- real invdOH;
- real invdHH;
- matrix invmat;
-};
-
struct settledata
{
- settleparam_t massw; /* Parameters for SETTLE for coordinates */
- settleparam_t mass1; /* Parameters with all masses 1, for forces */
+ SettleParameters massw; /* Parameters for SETTLE for coordinates */
+ SettleParameters mass1; /* Parameters with all masses 1, for forces */
int nsettle; /* The number of settles on our rank */
int* ow1; /* Index to OW1 atoms, size nsettle + SIMD padding */
bool bUseSimd; /* Use SIMD intrinsics code, if possible */
};
-
-//! Initializes a projection matrix.
-static void init_proj_matrix(real invmO, real invmH, real dOH, real dHH, matrix inverseCouplingMatrix)
+/*! \brief Initializes a projection matrix.
+ *
+ * \param[in] invmO Reciprocal oxygen mass
+ * \param[in] invmH Reciprocal hydrogen mass
+ * \param[in] dOH Target O-H bond length
+ * \param[in] dHH Target H-H bond length
+ * \param[out] inverseCouplingMatrix Inverse bond coupling matrix for the projection version of SETTLE
+ */
+static void initializeProjectionMatrix(const real invmO,
+ const real invmH,
+ const real dOH,
+ const real dHH,
+ matrix inverseCouplingMatrix)
{
- /* We normalize the inverse masses with invmO for the matrix inversion.
- * so we can keep using masses of almost zero for frozen particles,
- * without running out of the float range in invertMatrix.
- */
+ // We normalize the inverse masses with invmO for the matrix inversion.
+ // so we can keep using masses of almost zero for frozen particles,
+ // without running out of the float range in invertMatrix.
double invmORelative = 1.0;
double invmHRelative = invmH / static_cast<double>(invmO);
double distanceRatio = dHH / static_cast<double>(dOH);
msmul(inverseCouplingMatrix, 1 / invmO, inverseCouplingMatrix);
}
-//! Initializes settle parameters.
-static void settleparam_init(settleparam_t* p, real mO, real mH, real invmO, real invmH, real dOH, real dHH)
+void initSettleParameters(SettleParameters* p,
+ const real mO,
+ const real mH,
+ const real invmO,
+ const real invmH,
+ const real dOH,
+ const real dHH)
{
- /* We calculate parameters in double precision to minimize errors.
- * The velocity correction applied during SETTLE coordinate constraining
- * introduces a systematic error of approximately 1 bit per atom,
- * depending on what the compiler does with the code.
- */
+ // We calculate parameters in double precision to minimize errors.
+ // The velocity correction applied during SETTLE coordinate constraining
+ // introduces a systematic error of approximately 1 bit per atom,
+ // depending on what the compiler does with the code.
double wohh;
p->mO = mO;
p->ra = ra;
p->irc2 = 1.0 / dHH;
- /* For projection: inverse masses and coupling matrix inversion */
+ // For projection: inverse masses and coupling matrix inversion
p->imO = invmO;
p->imH = invmH;
p->invdOH = 1.0 / dOH;
p->invdHH = 1.0 / dHH;
- init_proj_matrix(invmO, invmH, dOH, dHH, p->invmat);
+ initializeProjectionMatrix(invmO, invmH, dOH, dHH, p->invmat);
if (debug)
{
real dOH = mtop.ffparams.iparams[settle_type].settle.doh;
real dHH = mtop.ffparams.iparams[settle_type].settle.dhh;
- settleparam_init(&settled->mass1, 1.0, 1.0, 1.0, 1.0, dOH, dHH);
+ initSettleParameters(&settled->mass1, 1.0, 1.0, 1.0, 1.0, dOH, dHH);
settled->ow1 = nullptr;
settled->hw2 = nullptr;
{
int firstO = iatoms[1];
int firstH = iatoms[2];
- settleparam_init(&settled->massw, masses[firstO], masses[firstH], inverseMasses[firstO],
- inverseMasses[firstH], settled->mass1.dOH, settled->mass1.dHH);
+ initSettleParameters(&settled->massw, masses[firstO], masses[firstH], inverseMasses[firstO],
+ inverseMasses[firstH], settled->mass1.dOH, settled->mass1.dHH);
}
if (nsettle + pack_size > settled->nalloc)
* Berk Hess 2008-1-10
*/
- settleparam_t* p;
- real imO, imH, dOH, dHH, invdOH, invdHH;
- matrix invmat;
- int i, m, m2, ow1, hw2, hw3;
- rvec roh2, roh3, rhh, dc, fc;
+ SettleParameters* p;
+ real imO, imH, dOH, dHH, invdOH, invdHH;
+ matrix invmat;
+ int i, m, m2, ow1, hw2, hw3;
+ rvec roh2, roh3, rhh, dc, fc;
calcvir_atom_end *= DIM;
TypeBool bError = TypeBool(false);
- const settleparam_t* p = &settled->massw;
- T wh = T(p->wh);
- T rc = T(p->rc);
- T ra = T(p->ra);
- T rb = T(p->rb);
- T irc2 = T(p->irc2);
- T mO = T(p->mO);
- T mH = T(p->mH);
+ const SettleParameters* p = &settled->massw;
+ T wh = T(p->wh);
+ T rc = T(p->rc);
+ T ra = T(p->ra);
+ T rb = T(p->rb);
+ T irc2 = T(p->irc2);
+ T mO = T(p->mO);
+ T mH = T(p->mH);
T almost_zero = T(1e-12);
/* Abstract type for SETTLE that is defined only in the file that uses it */
struct settledata;
+/* \brief Parameters for SETTLE algorithm.
+ *
+ * \todo Remove duplicates, check if recomputing makes more sense in some cases.
+ * \todo Move the projection parameters into separate structure.
+ */
+struct SettleParameters
+{
+ //! Mass of oxygen atom
+ real mO;
+ //! Mass of hydrogen atom
+ real mH;
+ //! Relative hydrogen mass (i.e. mH/(mO+2*mH))
+ real wh;
+ //! Target distance between oxygen and hydrogen
+ real dOH;
+ //! Target distance between hydrogen atoms
+ real dHH;
+ //! Double relative H mass times height of the H-O-H triangle.
+ real ra;
+ //! Height of the H-O-H triangle minus ra.
+ real rb;
+ //! Half the H-H distance.
+ real rc;
+ //! Reciprocal H-H distance
+ real irc2;
+
+ /* Parameters below are used for projection */
+ //! Reciprocal oxygen mass
+ real imO;
+ //! Reciprocal hydrogen mass
+ real imH;
+ //! Reciprocal O-H distance
+ real invdOH;
+ //! Reciprocal H-H distance (again)
+ real invdHH;
+ //! Reciprocal projection matrix
+ matrix invmat;
+};
+
+/*! \brief Initializes settle parameters.
+ *
+ * \param[out] p SettleParameters structure to initialize
+ * \param[in] mO Mass of oxygen atom
+ * \param[in] mH Mass of hydrogen atom
+ * \param[in] invmO Reciprocal mass of oxygen atom
+ * \param[in] invmH Reciprocal mass of hydrogen atom
+ * \param[in] dOH Target O-H bond length
+ * \param[in] dHH Target H-H bond length
+ */
+void initSettleParameters(SettleParameters* p, real mO, real mH, real dOH, real invmO, real invmH, real dHH);
+
/*! \brief Initializes and returns a structure with SETTLE parameters */
settledata* settle_init(const gmx_mtop_t& mtop);
real dOH = mtop.ffparams.iparams[settle_type].settle.doh;
real dHH = mtop.ffparams.iparams[settle_type].settle.dhh;
- initSettleParameters(&settleParameters_, mO, mH, dOH, dHH);
+ initSettleParameters(&settleParameters_, mO, mH, 1.0 / mO, 1.0 / mH, dOH, dHH);
allocateDeviceBuffer(&d_virialScaled_, 6, deviceContext_);
h_virialScaled_.resize(6);
namespace gmx
{
-/* \brief Parameters for SETTLE algorithm.
- *
- * Following structure and subroutine are copy-pasted from the CPU version of SETTLE.
- * This is a temporary solution and they will be recycled in more appropriate way.
- * \todo Remove duplicates, check if recomputing makes more sense in some cases.
- * \todo Move the projection parameters into separate structure.
- */
-struct SettleParameters
-{
- //! Mass of oxygen atom
- float mO;
- //! Mass of hydrogen atom
- float mH;
- //! Relative hydrogen mass (i.e. mH/(mO+2*mH))
- float wh;
- //! Target distance between oxygen and hydrogen
- float dOH;
- //! Target distance between hydrogen atoms
- float dHH;
- //! Double relative H mass times height of the H-O-H triangle.
- float ra;
- //! Height of the H-O-H triangle minus ra.
- float rb;
- //! Half the H-H distance.
- float rc;
- //! Reciprocal H-H distance
- float irc2;
-
- /* Parameters below are used for projection */
- //! Reciprocal oxygen mass
- float imO;
- //! Reciprocal hydrogen mass
- float imH;
- //! Reciprocal O-H distance
- float invdOH;
- //! Reciprocal H-H distance (again)
- float invdHH;
- //! Reciprocal projection matrix
- matrix invmat;
-};
-
-/*! \brief Initializes a projection matrix.
- *
- * \todo This is identical to to CPU code. Unification is needed.
- *
- * \param[in] invmO Reciprocal oxygen mass
- * \param[in] invmH Reciprocal hydrogen mass
- * \param[in] dOH Target O-H bond length
- * \param[in] dHH Target H-H bond length
- * \param[out] inverseCouplingMatrix Inverse bond coupling matrix for the projection version of SETTLE
- */
-static void initializeProjectionMatrix(const real invmO,
- const real invmH,
- const real dOH,
- const real dHH,
- matrix inverseCouplingMatrix)
-{
- // We normalize the inverse masses with invmO for the matrix inversion.
- // so we can keep using masses of almost zero for frozen particles,
- // without running out of the float range in invertMatrix.
- double invmORelative = 1.0;
- double invmHRelative = invmH / static_cast<double>(invmO);
- double distanceRatio = dHH / static_cast<double>(dOH);
-
- /* Construct the constraint coupling matrix */
- matrix mat;
- mat[0][0] = invmORelative + invmHRelative;
- mat[0][1] = invmORelative * (1.0 - 0.5 * gmx::square(distanceRatio));
- mat[0][2] = invmHRelative * 0.5 * distanceRatio;
- mat[1][1] = mat[0][0];
- mat[1][2] = mat[0][2];
- mat[2][2] = invmHRelative + invmHRelative;
- mat[1][0] = mat[0][1];
- mat[2][0] = mat[0][2];
- mat[2][1] = mat[1][2];
-
- gmx::invertMatrix(mat, inverseCouplingMatrix);
-
- msmul(inverseCouplingMatrix, 1 / invmO, inverseCouplingMatrix);
-}
-
-/*! \brief Initializes settle parameters.
- *
- * \todo This is identical to the CPU code. Unification is needed.
- *
- * \param[out] p SettleParameters structure to initialize
- * \param[in] mO Mass of oxygen atom
- * \param[in] mH Mass of hydrogen atom
- * \param[in] dOH Target O-H bond length
- * \param[in] dHH Target H-H bond length
- */
-gmx_unused // Temporary solution to keep clang happy
- static void
- initSettleParameters(SettleParameters* p, const real mO, const real mH, const real dOH, const real dHH)
-{
- // We calculate parameters in double precision to minimize errors.
- // The velocity correction applied during SETTLE coordinate constraining
- // introduces a systematic error of approximately 1 bit per atom,
- // depending on what the compiler does with the code.
- double wohh;
-
- real invmO = 1.0 / mO;
- real invmH = 1.0 / mH;
-
- p->mO = mO;
- p->mH = mH;
- wohh = mO + 2.0 * mH;
- p->wh = mH / wohh;
- p->dOH = dOH;
- p->dHH = dHH;
- double rc = dHH / 2.0;
- double ra = 2.0 * mH * std::sqrt(dOH * dOH - rc * rc) / wohh;
- p->rb = std::sqrt(dOH * dOH - rc * rc) - ra;
- p->rc = rc;
- p->ra = ra;
- p->irc2 = 1.0 / dHH;
-
- // For projection: inverse masses and coupling matrix inversion
- p->imO = invmO;
- p->imH = invmH;
-
- p->invdOH = 1.0 / dOH;
- p->invdHH = 1.0 / dHH;
-
- initializeProjectionMatrix(invmO, invmH, dOH, dHH, p->invmat);
-}
-
/*! \internal \brief Class with interfaces and data for GPU version of SETTLE. */
class SettleGpu
{
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff