#include <climits>
#include <cmath>
+#include "gromacs/domdec/ga2la.h"
+#include "gromacs/domdec/localatomsetmanager.h"
#include "gromacs/gmxlib/network.h"
#include "gromacs/linearalgebra/nrjac.h"
#include "gromacs/math/do_fit.h"
// TODO This implementation of ensemble orientation restraints is nasty because
// a user can't just do multi-sim with single-sim orientation restraints.
-t_oriresdata::t_oriresdata(FILE* fplog,
- const gmx_mtop_t& mtop,
- const t_inputrec& ir,
- const t_commrec* cr,
- const gmx_multisim_t* ms,
- t_state* globalState) :
- numRestraints(gmx_mtop_ftype_count(mtop, F_ORIRES))
+void extendStateWithOriresHistory(const gmx_mtop_t& mtop, const t_inputrec& ir, t_state* globalState)
+{
+ GMX_RELEASE_ASSERT(globalState != nullptr,
+ "We need a valid global state in extendStateWithOriresHistory()");
+
+ const int numRestraints = gmx_mtop_ftype_count(mtop, F_ORIRES);
+ if (numRestraints > 0 && ir.orires_tau > 0)
+ {
+ /* Extend the state with the orires history */
+ globalState->flags |= enumValueToBitMask(StateEntry::OrireInitF);
+ globalState->hist.orire_initf = 1;
+ globalState->flags |= enumValueToBitMask(StateEntry::OrireDtav);
+ globalState->hist.orire_Dtav.resize(numRestraints * 5);
+ }
+}
+
+namespace
+{
+
+//! Creates and returns a list of global atom indices of the orientation restraint fit group
+std::vector<gmx::index> fitGlobalAtomIndices(const gmx_mtop_t& mtop)
+{
+ std::vector<gmx::index> indices;
+
+ for (int i = 0; i < mtop.natoms; i++)
+ {
+ if (getGroupType(mtop.groups, SimulationAtomGroupType::OrientationRestraintsFit, i) == 0)
+ {
+ indices.push_back(i);
+ }
+ }
+
+ return indices;
+}
+
+} // namespace
+
+t_oriresdata::t_oriresdata(FILE* fplog,
+ const gmx_mtop_t& mtop,
+ const t_inputrec& ir,
+ const gmx_multisim_t* ms,
+ t_state* globalState,
+ gmx::LocalAtomSetManager* localAtomSetManager) :
+ numRestraints(gmx_mtop_ftype_count(mtop, F_ORIRES)),
+ fitLocalAtomSet_(localAtomSetManager->add(fitGlobalAtomIndices(mtop)))
{
GMX_RELEASE_ASSERT(numRestraints > 0,
"orires() should only be called with orientation restraints present");
"in the system"));
}
- if (cr && PAR(cr))
- {
- GMX_THROW(gmx::InvalidInputError(
- "Orientation restraints do not work with MPI parallelization. Choose 1 MPI rank, "
- "if possible."));
- }
-
GMX_RELEASE_ASSERT(globalState != nullptr, "We need a valid global state in t_oriresdata()");
fc = ir.orires_fc;
edt = std::exp(-ir.delta_t / ir.orires_tau);
edt_1 = 1.0 - edt;
- /* Extend the state with the orires history */
- globalState->flags |= enumValueToBitMask(StateEntry::OrireInitF);
- globalState->hist.orire_initf = 1;
- globalState->flags |= enumValueToBitMask(StateEntry::OrireDtav);
- globalState->hist.orire_Dtav.resize(numRestraints * 5);
+ timeAveragingInitFactor_ = std::reference_wrapper<real>(globalState->hist.orire_initf);
+ DTensorsTimeAveragedHistory_ = globalState->hist.orire_Dtav;
}
orientations.resize(numRestraints);
}
tmpEq.resize(numExperiments);
- numReferenceAtoms = 0;
- for (int i = 0; i < mtop.natoms; i++)
- {
- if (getGroupType(mtop.groups, SimulationAtomGroupType::OrientationRestraintsFit, i) == 0)
- {
- numReferenceAtoms++;
- }
- }
- mref.resize(numReferenceAtoms);
- xref.resize(numReferenceAtoms);
- xtmp.resize(numReferenceAtoms);
-
eigenOutput.resize(numExperiments * c_numEigenRealsPerExperiment);
/* Determine the reference structure on the master node.
* Copy it to the other nodes after checking multi compatibility,
* so we are sure the subsystems match before copying.
*/
- auto x = makeArrayRef(globalState->x);
+ auto x = makeConstArrayRef(globalState->x);
rvec com = { 0, 0, 0 };
double mtot = 0.0;
int j = 0;
{
const t_atom& local = atomP.atom();
int i = atomP.globalAtomNumber();
- if (mtop.groups.groupNumbers[SimulationAtomGroupType::OrientationRestraintsFit].empty()
- || mtop.groups.groupNumbers[SimulationAtomGroupType::OrientationRestraintsFit][i] == 0)
+ if (getGroupType(mtop.groups, SimulationAtomGroupType::OrientationRestraintsFit, i) == 0)
{
- /* Not correct for free-energy with changing masses */
- mref[j] = local.m;
+ // Not correct for free-energy with changing masses
+ const real mass = local.m;
// Note that only one rank per sim is supported.
if (isMasterSim(ms))
{
- copy_rvec(x[i], xref[j]);
+ referenceCoordinates_.push_back(x[i]);
for (int d = 0; d < DIM; d++)
{
- com[d] += mref[j] * x[i][d];
+ com[d] += mass * x[i][d];
}
}
- mtot += mref[j];
+ fitMasses_.push_back(mass);
+ mtot += mass;
j++;
}
}
+
svmul(1.0 / mtot, com, com);
if (isMasterSim(ms))
{
- for (int j = 0; j < numReferenceAtoms; j++)
+ for (auto& refCoord : referenceCoordinates_)
{
- rvec_dec(xref[j], com);
+ refCoord -= com;
}
}
+ const size_t numFitAtoms = referenceCoordinates_.size();
+ fitLocalAtomIndices_.resize(numFitAtoms);
+ xTmp_.resize(numFitAtoms);
+
if (fplog)
{
fprintf(fplog, "Found %d orientation experiments\n", numExperiments);
fprintf(fplog, " experiment %d has %d restraints\n", i + 1, nr_ex[i]);
}
- fprintf(fplog, " the fit group consists of %d atoms and has total mass %g\n", numReferenceAtoms, mtot);
+ fprintf(fplog, " the fit group consists of %zu atoms and has total mass %g\n", numFitAtoms, mtot);
}
if (ms)
check_multi_int(fplog, ms, numRestraints, "the number of orientation restraints", FALSE);
check_multi_int(
- fplog, ms, numReferenceAtoms, "the number of fit atoms for orientation restraining", FALSE);
+ fplog, ms, numFitAtoms, "the number of fit atoms for orientation restraining", FALSE);
check_multi_int(fplog, ms, ir.nsteps, "nsteps", FALSE);
/* Copy the reference coordinates from the master to the other nodes */
- gmx_sum_sim(DIM * numReferenceAtoms, xref[0], ms);
+ gmx_sum_sim(DIM * referenceCoordinates_.size(), as_rvec_array(referenceCoordinates_.data())[0], ms);
}
please_cite(fplog, "Hess2003");
int nfa,
const t_iatom forceatoms[],
const t_iparams ip[],
- const t_mdatoms* md,
ArrayRef<const RVec> xWholeMolecules,
const rvec x[],
const t_pbc* pbc,
- t_oriresdata* od,
- const history_t* hist)
+ t_oriresdata* od)
{
real invn, pfac, r2, invr, corrfac, wsv2, sw, dev;
- double mtot;
rvec com, r_unrot, r;
const real two_thr = 2.0 / 3.0;
- const bool bTAV = (od->edt != 0);
- const real edt = od->edt;
- const real edt_1 = od->edt_1;
- gmx::ArrayRef<OriresMatEq> matEq = od->tmpEq;
- const int nref = od->numReferenceAtoms;
- gmx::ArrayRef<real> mref = od->mref;
- gmx::ArrayRef<const gmx::RVec> xref = od->xref;
- gmx::ArrayRef<gmx::RVec> xtmp = od->xtmp;
+ const bool bTAV = (od->edt != 0);
+ const real edt = od->edt;
+ const real edt_1 = od->edt_1;
+ gmx::ArrayRef<OriresMatEq> matEq = od->tmpEq;
+ gmx::ArrayRef<gmx::RVec> xFit = od->xTmp();
if (bTAV)
{
- od->exp_min_t_tau = hist->orire_initf * edt;
+ od->exp_min_t_tau = od->timeAveragingInitFactor() * edt;
/* Correction factor to correct for the lack of history
* at short times.
invn = 1.0;
}
+ // Extract the local atom indices involved in the fit group
+ const auto fitLocalAtomIndices = od->fitLocalAtomSet().localIndex();
+ // We need all atoms in the fit group to be local available. This means that
+ // orientation restraining is limited to one PP-rank. This should be ensured
+ // by the mdrun setup code. We assert here to catch incorrect setup code.
+ GMX_RELEASE_ASSERT(fitLocalAtomIndices.size() == od->referenceCoordinates().size(),
+ "All fit atoms should be locally available");
+
clear_rvec(com);
- mtot = 0;
- int j = 0;
- auto* massT = md->massT;
- auto* cORF = md->cORF;
- for (int i = 0; i < md->nr; i++)
- {
- if (cORF[i] == 0)
+ double mtot = 0.0;
+ gmx::index referenceListIndex = 0;
+ for (const int fitLocalAtomIndex : fitLocalAtomIndices)
+ {
+ const gmx::RVec& x = xWholeMolecules[fitLocalAtomIndex];
+ const real mass = od->fitMasses()[referenceListIndex];
+ xFit[referenceListIndex] = x;
+ for (int d = 0; d < DIM; d++)
{
- copy_rvec(xWholeMolecules[i], xtmp[j]);
- mref[j] = massT[i];
- for (int d = 0; d < DIM; d++)
- {
- com[d] += mref[j] * xtmp[j][d];
- }
- mtot += mref[j];
- j++;
+ com[d] += mass * x[d];
}
+ mtot += mass;
+ referenceListIndex++;
}
svmul(1.0 / mtot, com, com);
- for (int j = 0; j < nref; j++)
+ for (auto& xFitCoord : xFit)
{
- rvec_dec(xtmp[j], com);
+ xFitCoord -= com;
}
/* Calculate the rotation matrix to rotate x to the reference orientation */
- calc_fit_R(DIM, nref, mref.data(), as_rvec_array(xref.data()), as_rvec_array(xtmp.data()), od->rotationMatrix);
+ calc_fit_R(DIM,
+ xFit.size(),
+ od->fitMasses().data(),
+ as_rvec_array(od->referenceCoordinates().data()),
+ as_rvec_array(xFit.data()),
+ od->rotationMatrix);
for (int fa = 0; fa < nfa; fa += 3)
{
*/
for (int i = 0; i < 5; i++)
{
- Dtav[i] = edt * hist->orire_Dtav[restraintIndex * 5 + i]
+ Dtav[i] = edt * od->DTensorsTimeAveragedHistory()[restraintIndex * 5 + i]
+ edt_1 * od->DTensorsEnsembleAv[restraintIndex][i];
}
}
/* Approx. 80*nfa/3 flops */
}
-void update_orires_history(const t_oriresdata& od, history_t* hist)
+void t_oriresdata::updateHistory()
{
- if (od.edt != 0)
+ if (edt != 0)
{
/* Copy the new time averages that have been calculated
* in calc_orires_dev.
*/
- hist->orire_initf = od.exp_min_t_tau;
- for (int pair = 0; pair < od.numRestraints; pair++)
+ *timeAveragingInitFactor_ = exp_min_t_tau;
+ for (int pair = 0; pair < numRestraints; pair++)
{
for (int i = 0; i < 5; i++)
{
- hist->orire_Dtav[pair * 5 + i] = od.DTensorsTimeAndEnsembleAv[pair][i];
+ DTensorsTimeAveragedHistory_[pair * 5 + i] = DTensorsTimeAndEnsembleAv[pair][i];
}
}
}
#ifndef GMX_MDTYPES_FCDATA_H
#define GMX_MDTYPES_FCDATA_H
+#include <functional>
#include <memory>
+#include <optional>
#include <vector>
+#include "gromacs/domdec/localatomset.h"
#include "gromacs/math/vectypes.h"
#include "gromacs/topology/idef.h"
#include "gromacs/utility/arrayref.h"
#include "gromacs/utility/real.h"
enum class DistanceRestraintWeighting : int;
+class gmx_ga2la_t;
struct gmx_mtop_t;
struct gmx_multisim_t;
-struct t_commrec;
struct t_inputrec;
class t_state;
+namespace gmx
+{
+class LocalAtomSetManager;
+}
+
typedef real rvec5[5];
/* Distance restraining stuff */
* \param[in] fplog Log file, can be nullptr
* \param[in] mtop The global topology
* \param[in] ir The input record
- * \param[in] cr The commrec, can be nullptr when not running in parallel
* \param[in] ms The multisim communicator, pass nullptr to avoid ensemble averaging
- * \param[in,out] globalState The global state, orientation restraint entires are added
+ * \param[in] globalState The global state, references are set to members
+ * \param[in,out] localAtomSetManager The local atom set manager
*
* \throws InvalidInputError when there is domain decomposition, fewer than 5 restraints,
* periodic molecules or more than 1 molecule for a moleculetype with restraints.
*/
- t_oriresdata(FILE* fplog,
- const gmx_mtop_t& mtop,
- const t_inputrec& ir,
- const t_commrec* cr,
- const gmx_multisim_t* ms,
- t_state* globalState);
+ t_oriresdata(FILE* fplog,
+ const gmx_mtop_t& mtop,
+ const t_inputrec& ir,
+ const gmx_multisim_t* ms,
+ t_state* globalState,
+ gmx::LocalAtomSetManager* localAtomSetManager);
//! Destructor
~t_oriresdata();
+ //! Returns the local atom set for fitting
+ const gmx::LocalAtomSet& fitLocalAtomSet() const { return fitLocalAtomSet_; }
+
+ //! Returns the list of reference coordinates
+ gmx::ArrayRef<const gmx::RVec> referenceCoordinates() const { return referenceCoordinates_; }
+
+ //! Returns the list of masses for fitting
+ gmx::ArrayRef<const real> fitMasses() const { return fitMasses_; }
+
+ //! Returns the list of local atoms for fitting, matching the order of referenceCoordinates
+ gmx::ArrayRef<const int> fitLocalAtomIndices() const { return fitLocalAtomIndices_; }
+
+ //! Returns the list of coordinates for temporary use, size matches referenceCoordinates
+ gmx::ArrayRef<gmx::RVec> xTmp() { return xTmp_; }
+
+ //! Returns the factor for initializing the time averaging
+ real timeAveragingInitFactor() const { return *timeAveragingInitFactor_; }
+
+ //! Returns a const view on the time averaged D tensor history
+ gmx::ArrayRef<const real> DTensorsTimeAveragedHistory() const
+ {
+ return DTensorsTimeAveragedHistory_;
+ }
+
+ //! Updates the history with the current values
+ void updateHistory();
+
//! Force constant for the restraints
real fc;
//! Multiplication factor for time averaging
int numExperiments;
//! The minimum iparam type index for restraints
int typeMin;
- //! The number of atoms for the fit
- int numReferenceAtoms;
- //! The masses of the reference atoms
- std::vector<real> mref;
+
+private:
+ //! List of local atom corresponding to the fit group
+ gmx::LocalAtomSet fitLocalAtomSet_;
//! The reference coordinates for the fit
- std::vector<gmx::RVec> xref;
- //! Temporary array for fitting
- std::vector<gmx::RVec> xtmp;
+ std::vector<gmx::RVec> referenceCoordinates_;
+ //! The masses for fitting
+ std::vector<real> fitMasses_;
+ //! List of reference atoms for fitting
+ std::vector<int> fitLocalAtomIndices_;
+ //! Temporary array, used for fitting
+ std::vector<gmx::RVec> xTmp_;
+ //! The factor for initializing the time averaging, only present when time averaging is used
+ //! This references the value stored in the global state, which depends on time.
+ std::optional<std::reference_wrapper<real>> timeAveragingInitFactor_;
+ //! View on the time averaged history of the orientation tensors
+ gmx::ArrayRef<real> DTensorsTimeAveragedHistory_;
+
+public:
//! Rotation matrix to rotate to the reference coordinates
matrix rotationMatrix;
//! Array of order tensors, one for each experiment
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff