2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team.
6 * Copyright (c) 2012,2014,2015,2016,2017 by the GROMACS development team.
7 * Copyright (c) 2018,2019,2020,2021, by the GROMACS development team, led by
8 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
9 * and including many others, as listed in the AUTHORS file in the
10 * top-level source directory and at http://www.gromacs.org.
12 * GROMACS is free software; you can redistribute it and/or
13 * modify it under the terms of the GNU Lesser General Public License
14 * as published by the Free Software Foundation; either version 2.1
15 * of the License, or (at your option) any later version.
17 * GROMACS is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 * Lesser General Public License for more details.
22 * You should have received a copy of the GNU Lesser General Public
23 * License along with GROMACS; if not, see
24 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
25 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
27 * If you want to redistribute modifications to GROMACS, please
28 * consider that scientific software is very special. Version
29 * control is crucial - bugs must be traceable. We will be happy to
30 * consider code for inclusion in the official distribution, but
31 * derived work must not be called official GROMACS. Details are found
32 * in the README & COPYING files - if they are missing, get the
33 * official version at http://www.gromacs.org.
35 * To help us fund GROMACS development, we humbly ask that you cite
36 * the research papers on the package. Check out http://www.gromacs.org.
38 #ifndef GMX_MDTYPES_FCDATA_H
39 #define GMX_MDTYPES_FCDATA_H
46 #include "gromacs/domdec/localatomset.h"
47 #include "gromacs/math/vectypes.h"
48 #include "gromacs/topology/idef.h"
49 #include "gromacs/utility/arrayref.h"
50 #include "gromacs/utility/classhelpers.h"
51 #include "gromacs/utility/real.h"
53 enum class DistanceRestraintWeighting : int;
56 struct gmx_multisim_t;
62 class LocalAtomSetManager;
65 typedef real rvec5[5];
67 /* Distance restraining stuff */
68 typedef struct t_disresdata
70 DistanceRestraintWeighting dr_weighting; /* Weighting of pairs in one restraint */
71 bool dr_bMixed; /* Use sqrt of the instantaneous times *
72 * the time averaged violation */
73 real dr_fc; /* Force constant for disres, *
74 * which is multiplied by a (possibly) *
75 * different factor for each restraint */
76 real dr_tau; /* Time constant for disres */
77 real ETerm; /* multiplication factor for time averaging */
78 real ETerm1; /* 1 - ETerm1 */
79 real exp_min_t_tau; /* Factor for slowly switching on the force */
80 int nres; /* The number of distance restraints */
81 int npair; /* The number of distance restraint pairs */
82 int type_min; /* The minimum iparam type index for restraints */
83 real sumviol; /* The sum of violations */
84 real* rt; /* The instantaneous distance (npair) */
85 real* rm3tav; /* The time averaged distance (npair) */
86 real* Rtl_6; /* The instantaneous r^-6 (nres) */
87 real* Rt_6; /* The instantaneous ensemble averaged r^-6 (nres) */
88 real* Rtav_6; /* The time and ensemble averaged r^-6 (nres) */
89 int nsystems; /* The number of systems for ensemble averaging */
91 /* TODO: Implement a proper solution for parallel disre indexing */
92 const t_iatom* forceatomsStart; /* Pointer to the start of the disre forceatoms */
95 /* All coefficients for the matrix equation for the orientation tensor */
98 real rhs[5]; /* The right hand side of the matrix equation */
99 real mat[5][5]; /* The matrix */
102 //! \brief Orientation restraining stuff
105 /*! \brief Constructor
107 * \param[in] fplog Log file, can be nullptr
108 * \param[in] mtop The global topology
109 * \param[in] ir The input record
110 * \param[in] ms The multisim communicator, pass nullptr to avoid ensemble averaging
111 * \param[in] globalState The global state, references are set to members
112 * \param[in,out] localAtomSetManager The local atom set manager
114 * \throws InvalidInputError when there is domain decomposition, fewer than 5 restraints,
115 * periodic molecules or more than 1 molecule for a moleculetype with restraints.
117 t_oriresdata(FILE* fplog,
118 const gmx_mtop_t& mtop,
119 const t_inputrec& ir,
120 const gmx_multisim_t* ms,
121 t_state* globalState,
122 gmx::LocalAtomSetManager* localAtomSetManager);
127 //! Returns the local atom set for fitting
128 const gmx::LocalAtomSet& fitLocalAtomSet() const { return fitLocalAtomSet_; }
130 //! Returns the list of reference coordinates
131 gmx::ArrayRef<const gmx::RVec> referenceCoordinates() const { return referenceCoordinates_; }
133 //! Returns the list of masses for fitting
134 gmx::ArrayRef<const real> fitMasses() const { return fitMasses_; }
136 //! Returns the list of local atoms for fitting, matching the order of referenceCoordinates
137 gmx::ArrayRef<const int> fitLocalAtomIndices() const { return fitLocalAtomIndices_; }
139 //! Returns the list of coordinates for temporary use, size matches referenceCoordinates
140 gmx::ArrayRef<gmx::RVec> xTmp() { return xTmp_; }
142 //! Returns the factor for initializing the time averaging
143 real timeAveragingInitFactor() const { return *timeAveragingInitFactor_; }
145 //! Returns a const view on the time averaged D tensor history
146 gmx::ArrayRef<const real> DTensorsTimeAveragedHistory() const
148 return DTensorsTimeAveragedHistory_;
151 //! Updates the history with the current values
152 void updateHistory();
154 //! Force constant for the restraints
156 //! Multiplication factor for time averaging
160 //! Factor for slowly switching on the force
162 //! The number of orientation restraints
163 const int numRestraints;
164 //! The number of experiments
166 //! The minimum iparam type index for restraints
170 //! List of local atom corresponding to the fit group
171 gmx::LocalAtomSet fitLocalAtomSet_;
172 //! The reference coordinates for the fit
173 std::vector<gmx::RVec> referenceCoordinates_;
174 //! The masses for fitting
175 std::vector<real> fitMasses_;
176 //! List of reference atoms for fitting
177 std::vector<int> fitLocalAtomIndices_;
178 //! Temporary array, used for fitting
179 std::vector<gmx::RVec> xTmp_;
180 //! The factor for initializing the time averaging, only present when time averaging is used
181 //! This references the value stored in the global state, which depends on time.
182 std::optional<std::reference_wrapper<real>> timeAveragingInitFactor_;
183 //! View on the time averaged history of the orientation tensors
184 gmx::ArrayRef<real> DTensorsTimeAveragedHistory_;
187 //! Rotation matrix to rotate to the reference coordinates
188 matrix rotationMatrix;
189 //! Array of order tensors, one for each experiment
190 tensor* orderTensors = nullptr;
191 //! The order tensor D for all restraints
192 rvec5* DTensors = nullptr;
193 //! The ensemble averaged D for all restraints
194 rvec5* DTensorsEnsembleAv = nullptr;
195 //! The time and ensemble averaged D restraints
196 rvec5* DTensorsTimeAndEnsembleAv = nullptr;
197 //! The calculated instantaneous orientations
198 std::vector<real> orientations;
199 //! The calculated emsemble averaged orientations
200 gmx::ArrayRef<real> orientationsEnsembleAv;
201 //! Buffer for the calculated emsemble averaged orientations, only used with ensemble averaging
202 std::vector<real> orientationsEnsembleAvBuffer;
203 //! The calculated time and ensemble averaged orientations
204 gmx::ArrayRef<real> orientationsTimeAndEnsembleAv;
205 //! The weighted (using kfac) RMS deviation
206 std::vector<real> orientationsTimeAndEnsembleAvBuffer;
207 //! Buffer for the weighted (using kfac) RMS deviation, only used with time averaging
209 //! An temporary array of matrix + rhs
210 std::vector<OriresMatEq> tmpEq;
211 //! The number of eigenvalues + eigenvectors per experiment
212 static constexpr int c_numEigenRealsPerExperiment = 12;
213 //! Eigenvalues/vectors, for output only (numExperiments x 12)
214 std::vector<real> eigenOutput;
216 // variables for diagonalization with diagonalize_orires_tensors()
217 //! Tensor to diagonalize
218 std::array<gmx::DVec, DIM> M;
220 std::array<double, DIM> eig_diag;
222 std::array<gmx::DVec, DIM> v;
224 // Default copy and assign would be incorrect and manual versions are not yet implemented.
225 GMX_DISALLOW_COPY_AND_ASSIGN(t_oriresdata);
228 /* Cubic spline table for tabulated bonded interactions */
231 int n; /* n+1 is the number of points */
232 real scale; /* distance between two points */
233 std::vector<real> data; /* the actual table data, per point there are 4 numbers */
237 * Data struct used in the force calculation routines
238 * for storing the tables for bonded interactions and
239 * for storing information which is needed in following steps
240 * (for instance for time averaging in distance retraints)
241 * or for storing output, since force routines only return the potential.
245 std::vector<bondedtable_t> bondtab;
246 std::vector<bondedtable_t> angletab;
247 std::vector<bondedtable_t> dihtab;
249 // TODO: Convert to C++ and unique_ptr (currently this data is not freed)
250 t_disresdata* disres = nullptr;
251 std::unique_ptr<t_oriresdata> orires;