of motion. For constant NVE simulations started from
corresponding points in the same trajectory, the trajectories
are analytically, but not binary, identical to the
- :mdp-value:`integrator=md` leap-frog integrator. The the kinetic
+ :mdp-value:`integrator=md` leap-frog integrator. The kinetic
energy, which is determined from the whole step velocities and
is therefore slightly too high. The advantage of this integrator
is more accurate, reversible Nose-Hoover and Parrinello-Rahman
significant part of the simulation time. The temperature for one
or more groups of atoms (:mdp:`tc-grps`) is set with
:mdp:`ref-t`, the inverse friction constant for each group is
- set with :mdp:`tau-t`. The parameter :mdp:`tcoupl` is
- ignored. The random generator is initialized with
+ set with :mdp:`tau-t`. The parameters :mdp:`tcoupl` and :mdp:`nsttcouple`
+ are ignored. The random generator is initialized with
:mdp:`ld-seed`. When used as a thermostat, an appropriate value
for :mdp:`tau-t` is 2 ps, since this results in a friction that
is lower than the internal friction of water, while it is high
Generate a pair list with buffering. The buffer size is
automatically set based on :mdp:`verlet-buffer-tolerance`,
unless this is set to -1, in which case :mdp:`rlist` will be
- used. This option has an explicit, exact cut-off at :mdp:`rvdw`
- equal to :mdp:`rcoulomb`, unless PME or Ewald is used, in which
- case :mdp:`rcoulomb` > :mdp:`rvdw` is allowed. Currently only
- cut-off, reaction-field, PME or Ewald electrostatics and plain
- LJ are supported. Some :ref:`gmx mdrun` functionality is not yet
- supported with the :mdp-value:`cutoff-scheme=Verlet` scheme, but :ref:`gmx grompp`
- checks for this. Native GPU acceleration is only supported with
- :mdp-value:`cutoff-scheme=Verlet`. With GPU-accelerated PME or with separate PME
- ranks, :ref:`gmx mdrun` will automatically tune the CPU/GPU load
- balance by scaling :mdp:`rcoulomb` and the grid spacing. This
- can be turned off with ``mdrun -notunepme``. :mdp-value:`cutoff-scheme=Verlet` is
- faster than :mdp-value:`cutoff-scheme=group` when there is no water, or if
- :mdp-value:`cutoff-scheme=group` would use a pair-list buffer to conserve energy.
+ used.
.. mdp-value:: group
- Generate a pair list for groups of atoms. These groups
- correspond to the charge groups in the topology. This was the
- only cut-off treatment scheme before version 4.6, and is
- **deprecated since 5.1**. There is no explicit buffering of
- the pair list. This enables efficient force calculations for
- water, but energy is only conserved when a buffer is explicitly
- added.
+ Generate a pair list for groups of atoms, corresponding
+ to the charge groups in the topology. This option is no longer
+ supported.
.. mdp:: nstlist
.. mdp-value:: >0
- Frequency to update the neighbor list. When this is 0, the
- neighbor list is made only once. With energy minimization the
- pair list will be updated for every energy evaluation when
- :mdp:`nstlist` is greater than 0. With :mdp-value:`cutoff-scheme=Verlet` and
+ Frequency to update the neighbor list. When dynamics and
:mdp:`verlet-buffer-tolerance` set, :mdp:`nstlist` is actually
a minimum value and :ref:`gmx mdrun` might increase it, unless
it is set to 1. With parallel simulations and/or non-bonded
force calculation on the GPU, a value of 20 or 40 often gives
- the best performance. With :mdp-value:`cutoff-scheme=group` and non-exact
- cut-off's, :mdp:`nstlist` will affect the accuracy of your
- simulation and it can not be chosen freely.
+ the best performance.
.. mdp-value:: 0
The neighbor list is only constructed once and never
updated. This is mainly useful for vacuum simulations in which
- all particles see each other.
+ all particles see each other. But vacuum simulations are
+ (temporarily) not supported.
.. mdp-value:: <0
Unused.
-.. mdp:: ns-type
-
- .. mdp-value:: grid
-
- Make a grid in the box and only check atoms in neighboring grid
- cells when constructing a new neighbor list every
- :mdp:`nstlist` steps. In large systems grid search is much
- faster than simple search.
-
- .. mdp-value:: simple
-
- Check every atom in the box when constructing a new neighbor
- list every :mdp:`nstlist` steps (only with :mdp-value:`cutoff-scheme=group`
- cut-off scheme).
-
.. mdp:: pbc
.. mdp-value:: xyz
(0.005) [kJ mol\ :sup:`-1` ps\ :sup:`-1`]
- Useful only with the :mdp-value:`cutoff-scheme=Verlet` :mdp:`cutoff-scheme`. This sets
+ Used when performing a simulation with dynamics. This sets
the maximum allowed error for pair interactions per particle caused
by the Verlet buffer, which indirectly sets :mdp:`rlist`. As both
:mdp:`nstlist` and the Verlet buffer size are fixed (for
.. mdp:: rlist
(1) [nm]
- Cut-off distance for the short-range neighbor list. With the
- :mdp-value:`cutoff-scheme=Verlet` :mdp:`cutoff-scheme`, this is by default set by the
- :mdp:`verlet-buffer-tolerance` option and the value of
- :mdp:`rlist` is ignored.
+ Cut-off distance for the short-range neighbor list. With dynamics,
+ this is by default set by the :mdp:`verlet-buffer-tolerance` option
+ and the value of :mdp:`rlist` is ignored. Without dynamics, this
+ is by default set to the maximum cut-off plus 5% buffer, except
+ for test particle insertion, where the buffer is managed exactly
+ and automatically. For NVE simulations, where the automated
+ setting is not possible, the advised procedure is to run :ref:`gmx grompp`
+ with an NVT setup with the expected temperature and copy the resulting
+ value of :mdp:`rlist` to the NVE setup.
Electrostatics
:mdp:`epsilon-rf`. The dielectric constant can be set to
infinity by setting :mdp:`epsilon-rf` =0.
- .. mdp-value:: Generalized-Reaction-Field
-
- Generalized reaction field with Coulomb cut-off
- :mdp:`rcoulomb`, where :mdp:`rlist` >= :mdp:`rcoulomb`. The
- dielectric constant beyond the cut-off is
- :mdp:`epsilon-rf`. The ionic strength is computed from the
- number of charged (*i.e.* with non zero charge) charge
- groups. The temperature for the GRF potential is set with
- :mdp:`ref-t`.
-
- .. mdp-value:: Reaction-Field-zero
-
- In |Gromacs|, normal reaction-field electrostatics with
- :mdp-value:`cutoff-scheme=group` leads to bad energy
- conservation. :mdp-value:`coulombtype=Reaction-Field-zero` solves this by making
- the potential zero beyond the cut-off. It can only be used with
- an infinite dielectric constant (:mdp:`epsilon-rf` =0), because
- only for that value the force vanishes at the
- cut-off. :mdp:`rlist` should be 0.1 to 0.3 nm larger than
- :mdp:`rcoulomb` to accommodate the size of charge groups
- and diffusion between neighbor list updates. This, and the fact
- that table lookups are used instead of analytical functions make
- reaction-field-zero computationally more expensive than
- normal reaction-field.
-
- .. mdp-value:: Shift
-
- Analogous to :mdp-value:`vdwtype=Shift` for :mdp:`vdwtype`. You
- might want to use :mdp-value:`coulombtype=Reaction-Field-zero` instead, which has
- a similar potential shape, but has a physical interpretation and
- has better energies due to the exclusion correction terms.
-
- .. mdp-value:: Encad-Shift
-
- The Coulomb potential is decreased over the whole range, using
- the definition from the Encad simulation package.
-
- .. mdp-value:: Switch
-
- Analogous to :mdp-value:`vdwtype=Switch` for
- :mdp:`vdwtype`. Switching the Coulomb potential can lead to
- serious artifacts, advice: use :mdp-value:`coulombtype=Reaction-Field-zero`
- instead.
-
.. mdp-value:: User
+ Currently unsupported.
:ref:`gmx mdrun` will now expect to find a file ``table.xvg``
with user-defined potential functions for repulsion, dispersion
and Coulomb. When pair interactions are present, :ref:`gmx
.. mdp-value:: PME-Switch
+ Currently unsupported.
A combination of PME and a switch function for the direct-space
part (see above). :mdp:`rcoulomb` is allowed to be smaller than
- :mdp:`rlist`. This is mainly useful constant energy simulations
- (note that using PME with :mdp-value:`cutoff-scheme=Verlet`
- will be more efficient).
+ :mdp:`rlist`.
.. mdp-value:: PME-User
+ Currently unsupported.
A combination of PME and user tables (see
above). :mdp:`rcoulomb` is allowed to be smaller than
:mdp:`rlist`. The PME mesh contribution is subtracted from the
.. mdp-value:: PME-User-Switch
+ Currently unsupported.
A combination of PME-User and a switching function (see
above). The switching function is applied to final
particle-particle interaction, *i.e.* both to the user supplied
.. mdp:: coulomb-modifier
- .. mdp-value:: Potential-shift-Verlet
-
- Selects Potential-shift with the Verlet cutoff-scheme, as it is
- (nearly) free; selects None with the group cutoff-scheme.
-
.. mdp-value:: Potential-shift
Shift the Coulomb potential by a constant such that it is zero
.. mdp-value:: None
- Use an unmodified Coulomb potential. With the group scheme this
- means no exact cut-off is used, energies and forces are
- calculated for all pairs in the pair list.
+ Use an unmodified Coulomb potential. This can be useful
+ when comparing energies with those computed with other software.
.. mdp:: rcoulomb-switch
:mdp-value:`vdwtype=Cut-off` with :mdp-value:`vdw-modifier=Force-switch`.
The LJ (not Buckingham) potential is decreased over the whole range and
the forces decay smoothly to zero between :mdp:`rvdw-switch` and
- :mdp:`rvdw`. The neighbor search cut-off :mdp:`rlist` should
- be 0.1 to 0.3 nm larger than :mdp:`rvdw` to accommodate the
- size of charge groups and diffusion between neighbor list
- updates.
+ :mdp:`rvdw`.
.. mdp-value:: Switch
potential and force functions are continuously smooth, but be
aware that all switch functions will give rise to a bulge
(increase) in the force (since we are switching the
- potential). The neighbor search cut-off :mdp:`rlist` should be
- 0.1 to 0.3 nm larger than :mdp:`rvdw` to accommodate the
- size of charge groups and diffusion between neighbor list
- updates.
-
- .. mdp-value:: Encad-Shift
-
- The LJ (not Buckingham) potential is decreased over the whole
- range, using the definition from the Encad simulation package.
+ potential).
.. mdp-value:: User
+ Currently unsupported.
See user for :mdp:`coulombtype`. The function value at zero is
not important. When you want to use LJ correction, make sure
that :mdp:`rvdw` corresponds to the cut-off in the user-defined
.. mdp:: vdw-modifier
- .. mdp-value:: Potential-shift-Verlet
-
- Selects Potential-shift with the Verlet cutoff-scheme, as it is
- (nearly) free; selects None with the group cutoff-scheme.
-
.. mdp-value:: Potential-shift
Shift the Van der Waals potential by a constant such that it is
.. mdp-value:: None
- Use an unmodified Van der Waals potential. With the group scheme
- this means no exact cut-off is used, energies and forces are
- calculated for all pairs in the pair list.
+ Use an unmodified Van der Waals potential. This can be useful
+ when comparing energies with those computed with other software.
.. mdp-value:: Force-switch
(1) [nm]
Extension of the non-bonded potential lookup tables beyond the
- largest cut-off distance. The value should be large enough to
- account for charge group sizes and the diffusion between
- neighbor-list updates. Without user defined potential the same
- table length is used for the lookup tables for the 1-4
- interactions, which are always tabulated irrespective of the use of
- tables for the non-bonded interactions. The value of
- :mdp:`table-extension` in no way affects the values of
- :mdp:`rlist`, :mdp:`rcoulomb`, or :mdp:`rvdw`.
+ largest cut-off distance. With actual non-bonded interactions
+ the tables are never accessed beyond the cut-off. But a longer
+ table length might be needed for the 1-4 interactions, which
+ are always tabulated irrespective of the use of tables for
+ the non-bonded interactions.
.. mdp:: energygrp-table
+ Currently unsupported.
When user tables are used for electrostatics and/or VdW, here one
can give pairs of energy groups for which seperate user tables
should be used. The two energy groups will be appended to the table
integrators, the leap-frog :mdp-value:`integrator=md` integrator
only supports 1. Data for the NH chain variables is not printed
to the :ref:`edr` file by default, but can be turned on with the
- :mdp:`print-nose-hoover-chains` option.
+ :mdp:`print-nose-hoover-chain-variables` option.
.. mdp:: print-nose-hoover-chain-variables
SHAKE is slightly slower and less stable than LINCS, but does
work with angle constraints. The relative tolerance is set with
:mdp:`shake-tol`, 0.0001 is a good value for "normal" MD. SHAKE
- does not support constraints between atoms on different nodes,
- thus it can not be used with domain decompositon when inter
- charge-group constraints are present. SHAKE can not be used with
- energy minimization.
+ does not support constraints between atoms on different
+ decomposition domains, so it can only be used with domain
+ decomposition when so-called update-groups are used, which is
+ usally the case when only bonds involving hydrogens are
+ constrained. SHAKE can not be used with energy minimization.
.. mdp:: continuation
.. mdp-value:: direction-periodic
- As :mdp-value:`pull-coord1-geometry=direction`, but allows the distance to be larger
- than half the box size. With this geometry the box should not be
+ As :mdp-value:`pull-coord1-geometry=direction`, but does not apply
+ periodic box vector corrections to keep the distance within half
+ the box length. This is (only) useful for pushing groups apart
+ by more than half the box length by continuously changing the reference
+ location using a pull rate. With this geometry the box should not be
dynamic (*e.g.* no pressure scaling) in the pull dimensions and
- the pull force is not added to virial.
+ the pull force is not added to the virial.
.. mdp-value:: direction-relative
.. mdp:: sc-r-power
(6)
- the power of the radial term in the soft-core equation. Possible
- values are 6 and 48. 6 is more standard, and is the default. When
- 48 is used, then sc-alpha should generally be much lower (between
- 0.001 and 0.003).
+ power 6 for the radial term in the soft-core equation.
+
+ (48)
+ (deprecated) power 48 for the radial term in the soft-core equation.
+ Note that sc-alpha should generally be much lower (between 0.001 and 0.003).
.. mdp:: sc-coul
(1.0) [nm] Lower extension of the split cylinder #1.
+Density-guided simulations
+^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+These options enable and control the calculation and application of additional
+forces that are derived from three-dimensional densities, e.g., from cryo
+electron-microscopy experiments. (See the `reference manual`_ for details)
+
+.. mdp:: density-guided-simulation-active
+
+ (no) Activate density-guided simulations.
+
+.. mdp:: density-guided-simulation-group
+
+ (protein) The atoms that are subject to the forces from the density-guided
+ simulation and contribute to the simulated density.
+
+.. mdp:: density-guided-simulation-similarity-measure
+
+ (inner-product) Similarity measure between the density that is calculated
+ from the atom positions and the reference density.
+
+ .. mdp-value:: inner-product
+
+ Takes the sum of the product of reference density and simulated density
+ voxel values.
+
+ .. mdp-value:: relative-entropy
+
+ Uses the negative relative entropy (or Kullback-Leibler divergence)
+ between reference density and simulated density as similarity measure.
+ Negative density values are ignored.
+
+.. mdp:: density-guided-simulation-atom-spreading-weight
+
+ (unity) Determines the multiplication factor for the Gaussian kernel when
+ spreading atoms on the grid.
+
+ .. mdp-value:: unity
+
+ Every atom in the density fitting group is assigned the same unit factor.
+
+ .. mdp-value:: mass
+
+ Atoms contribute to the simulated density proportional to their mass.
+
+ .. mdp-value:: charge
+
+ Atoms contribute to the simulated density proportional to their charge.
+
+.. mdp:: density-guided-simulation-force-constant
+
+ (1e+09) [kJ mol\ :sup:`-1`] The scaling factor for density-guided simulation
+ forces. May also be negative.
+
+.. mdp:: density-guided-simulation-gaussian-transform-spreading-width
+
+ (0.2) [nm] The Gaussian RMS width for the spread kernel for the simulated
+ density.
+
+.. mdp:: density-guided-simulation-gaussian-transform-spreading-range-in-multiples-of-width
+
+ (4) The range after which the gaussian is cut off in multiples of the Gaussian
+ RMS width described above.
+
+.. mdp:: density-guided-simulation-reference-density-filename
+
+ (reference.mrc) Reference density file name using an absolute path or a path
+ relative to the to the folder from which :ref:`gmx mdrun` is called.
+
+.. mdp:: density-guided-simulation-nst
+
+ (1) Interval in steps at which the density fitting forces are evaluated
+ and applied. The forces are scaled by this number when applied (See the
+ `reference manual`_ for details).
+
+.. mdp:: density-guided-simulation-normalize-densities
+
+ (true) Normalize the sum of density voxel values to one for the reference
+ density as well as the simulated density.
+
+.. mdp:: density-guided-simulation-adaptive-force-scaling
+
+ (false) Adapt the force constant to ensure a steady increase in similarity
+ between simulated and reference density.
+
+ .. mdp-value: false
+
+ Do not use adaptive force scaling.
+
+ .. mdp-value:: true
+
+ Use adaptive force scaling.
+
+.. mdp:: density-guided-simulation-adaptive-force-scaling-time-constant
+
+ (4) [ps] Couple force constant to increase in similarity with reference density
+ with this time constant. Larger times result in looser coupling.
User defined thingies
^^^^^^^^^^^^^^^^^^^^^
*/
class BiasParams
{
- public:
- /*! \brief Switch to turn off update skips, useful for testing.
- */
- enum class DisableUpdateSkips
- {
- no, /**< Allow update skips (when supported by the method) */
- yes /**< Disable update skips */
- };
-
- /*! \brief
- * Check if the parameters permit skipping updates.
- *
- * Generally, we can skip updates of points that are non-local
- * at the time of the update if we for later times, when the points
- * with skipped updates have become local, know exactly how to apply
- * the previous updates. The free energy updates only depend
- * on local sampling, but the histogram rescaling factors
- * generally depend on the histogram size (all samples).
- * If the histogram size is kept constant or the scaling factors
- * are trivial, this is not a problem. However, if the histogram growth
- * is scaled down by some factor the size at the time of the update
- * needs to be known. It would be fairly simple to, for a deterministically
- * growing histogram, backtrack and calculate this value, but currently
- * we just disallow this case. This is not a restriction because it
- * only affects the local Boltzmann target type for which every update
- * is currently anyway global because the target is always updated globally.
- *
- * \returns true when we can skip updates.
- */
- inline bool skipUpdates() const
- {
- return (!disableUpdateSkips_ && localWeightScaling == 1);
- }
-
- /*! \brief
- * Returns the the radius that needs to be sampled around a point before it is considered covered.
- */
- inline const awh_ivec &coverRadius() const
- {
- return coverRadius_;
- }
-
- /*! \brief
- * Returns whether we should sample the coordinate.
- *
- * \param[in] step The MD step number.
- */
- inline bool isSampleCoordStep(int64_t step) const
- {
- return (step > 0 && step % numStepsSampleCoord_ == 0);
- }
-
- /*! \brief
- * Returns whether we should update the free energy.
- *
- * \param[in] step The MD step number.
- */
- inline bool isUpdateFreeEnergyStep(int64_t step) const
- {
- int stepIntervalUpdateFreeEnergy = numSamplesUpdateFreeEnergy_*numStepsSampleCoord_;
- return (step > 0 && step % stepIntervalUpdateFreeEnergy == 0);
- }
-
- /*! \brief
- * Returns whether we should update the target distribution.
- *
- * \param[in] step The MD step number.
- */
- inline bool isUpdateTargetStep(int64_t step) const
- {
- return step % numStepsUpdateTarget_ == 0;
- }
-
- /*! \brief
- * Returns if to do checks for covering in the initial stage.
- *
- * To avoid overhead due to expensive checks, we do not check
- * at every free energy update. However, if checks are
- * performed too rarely the detection of coverings will be
- * delayed, ultimately affecting free energy convergence.
- *
- * \param[in] step Time step.
- * \returns true at steps where checks should be performed.
- * \note Only returns true at free energy update steps.
- */
- bool isCheckCoveringStep(int64_t step) const
- {
- return step > 0 && (step % numStepsCheckCovering_ == 0);
- }
-
- /*! \brief
- * Returns if to perform checks for anomalies in the histogram.
- *
- * To avoid overhead due to expensive checks, we do not check
- * at every free energy update. These checks are only used for
- * warning the user and can be made as infrequently as
- * neccessary without affecting the algorithm itself.
- *
- * \param[in] step Time step.
- * \returns true at steps where checks should be performed.
- * \note Only returns true at free energy update steps.
- * \todo Currently this function just calls isCheckCoveringStep but the checks could be done less frequently.
- */
- bool isCheckHistogramForAnomaliesStep(int64_t step) const
- {
- return isCheckCoveringStep(step);
- }
-
- /*! \brief Constructor.
- *
- * The local Boltzmann target distibution is defined by
- * 1) Adding the sampled weights instead of the target weights to the reference weight histogram.
- * 2) Scaling the weights of these samples by the beta scaling factor.
- * 3) Setting the target distribution equal the reference weight histogram.
- * This requires the following special update settings:
- * localWeightScaling = targetParam
- * idealWeighthistUpdate = false
- * Note: these variables could in principle be set to something else also for other target distribution types.
- * However, localWeightScaling < 1 is in general expected to give lower efficiency and, except for local Boltzmann,
- * idealWeightHistUpdate = false gives (in my experience) unstable, non-converging results.
- *
- * \param[in] awhParams AWH parameters.
- * \param[in] awhBiasParams Bias parameters.
- * \param[in] dimParams Bias dimension parameters.
- * \param[in] beta 1/(k_B T) in units of 1/(kJ/mol), should be > 0.
- * \param[in] mdTimeStep The MD time step.
- * \param[in] numSharingSimulations The number of simulations to share the bias across.
- * \param[in] gridAxis The grid axes.
- * \param[in] disableUpdateSkips If to disable update skips, useful for testing.
- * \param[in] biasIndex Index of the bias.
- */
- BiasParams(const AwhParams &awhParams,
- const AwhBiasParams &awhBiasParams,
- const std::vector<DimParams> &dimParams,
- double beta,
- double mdTimeStep,
- DisableUpdateSkips disableUpdateSkips,
- int numSharingSimulations,
- const std::vector<GridAxis> &gridAxis,
- int biasIndex);
-
- /* Data members */
- const double invBeta; /**< 1/beta = kT in kJ/mol */
- private:
- const int64_t numStepsSampleCoord_; /**< Number of steps per coordinate value sample. */
- public:
- const int numSamplesUpdateFreeEnergy_; /**< Number of samples per free energy update. */
- private:
- const int64_t numStepsUpdateTarget_; /**< Number of steps per updating the target distribution. */
- const int64_t numStepsCheckCovering_; /**< Number of steps per checking for covering. */
- public:
- const int eTarget; /**< Type of target distribution. */
- const double freeEnergyCutoffInKT; /**< Free energy cut-off in kT for cut-off target distribution. */
- const double temperatureScaleFactor; /**< Temperature scaling factor for temperature scaled targed distributions. */
- const bool idealWeighthistUpdate; /**< Update reference weighthistogram using the target distribution? Otherwise use the realized distribution. */
- const int numSharedUpdate; /**< The number of (multi-)simulations sharing the bias update */
- const double updateWeight; /**< The probability weight accumulated for each update. */
- const double localWeightScaling; /**< Scaling factor applied to a sample before adding it to the reference weight histogram (= 1, usually). */
- const double initialErrorInKT; /**< Estimated initial free energy error in kT. */
- const double initialHistogramSize; /**< Initial reference weight histogram size. */
- private:
- awh_ivec coverRadius_; /**< The radius (in points) that needs to be sampled around a point before it is considered covered. */
- public:
- const bool convolveForce; /**< True if we convolve the force, false means use MC between umbrellas. */
- const int biasIndex; /**< Index of the bias, used as a second random seed and for priting. */
- private:
- const bool disableUpdateSkips_; /**< If true, we disallow update skips, even when the method supports it. */
+public:
+ /*! \brief Switch to turn off update skips, useful for testing.
+ */
+ enum class DisableUpdateSkips
+ {
+ no, /**< Allow update skips (when supported by the method) */
+ yes /**< Disable update skips */
+ };
+
+ /*! \brief
+ * Check if the parameters permit skipping updates.
+ *
+ * Generally, we can skip updates of points that are non-local
+ * at the time of the update if we for later times, when the points
+ * with skipped updates have become local, know exactly how to apply
+ * the previous updates. The free energy updates only depend
+ * on local sampling, but the histogram rescaling factors
+ * generally depend on the histogram size (all samples).
+ * If the histogram size is kept constant or the scaling factors
+ * are trivial, this is not a problem. However, if the histogram growth
+ * is scaled down by some factor the size at the time of the update
+ * needs to be known. It would be fairly simple to, for a deterministically
+ * growing histogram, backtrack and calculate this value, but currently
+ * we just disallow this case. This is not a restriction because it
+ * only affects the local Boltzmann target type for which every update
+ * is currently anyway global because the target is always updated globally.
+ *
+ * \returns true when we can skip updates.
+ */
+ inline bool skipUpdates() const { return (!disableUpdateSkips_ && localWeightScaling == 1); }
+
+ /*! \brief
+ * Returns the radius that needs to be sampled around a point before it is considered covered.
+ */
+ inline const awh_ivec& coverRadius() const { return coverRadius_; }
+
+ /*! \brief
+ * Returns whether we should sample the coordinate.
+ *
+ * \param[in] step The MD step number.
+ */
+ inline bool isSampleCoordStep(int64_t step) const
+ {
+ return (step > 0 && step % numStepsSampleCoord_ == 0);
+ }
+
+ /*! \brief
+ * Returns whether we should update the free energy.
+ *
+ * \param[in] step The MD step number.
+ */
+ inline bool isUpdateFreeEnergyStep(int64_t step) const
+ {
+ int stepIntervalUpdateFreeEnergy = numSamplesUpdateFreeEnergy_ * numStepsSampleCoord_;
+ return (step > 0 && step % stepIntervalUpdateFreeEnergy == 0);
+ }
+
+ /*! \brief
+ * Returns whether we should update the target distribution.
+ *
+ * \param[in] step The MD step number.
+ */
+ inline bool isUpdateTargetStep(int64_t step) const { return step % numStepsUpdateTarget_ == 0; }
+
+ /*! \brief
+ * Returns if to do checks for covering in the initial stage.
+ *
+ * To avoid overhead due to expensive checks, we do not check
+ * at every free energy update. However, if checks are
+ * performed too rarely the detection of coverings will be
+ * delayed, ultimately affecting free energy convergence.
+ *
+ * \param[in] step Time step.
+ * \returns true at steps where checks should be performed.
+ * \note Only returns true at free energy update steps.
+ */
- bool isCheckCoveringStep(int64_t step) const { return step % numStepsCheckCovering_ == 0; }
++ bool isCheckCoveringStep(int64_t step) const
++ {
++ return step > 0 && (step % numStepsCheckCovering_ == 0);
++ }
+
+ /*! \brief
+ * Returns if to perform checks for anomalies in the histogram.
+ *
+ * To avoid overhead due to expensive checks, we do not check
+ * at every free energy update. These checks are only used for
+ * warning the user and can be made as infrequently as
+ * neccessary without affecting the algorithm itself.
+ *
+ * \param[in] step Time step.
+ * \returns true at steps where checks should be performed.
+ * \note Only returns true at free energy update steps.
+ * \todo Currently this function just calls isCheckCoveringStep but the checks could be done less frequently.
+ */
+ bool isCheckHistogramForAnomaliesStep(int64_t step) const { return isCheckCoveringStep(step); }
+
+ /*! \brief Constructor.
+ *
+ * The local Boltzmann target distibution is defined by
+ * 1) Adding the sampled weights instead of the target weights to the reference weight histogram.
+ * 2) Scaling the weights of these samples by the beta scaling factor.
+ * 3) Setting the target distribution equal the reference weight histogram.
+ * This requires the following special update settings:
+ * localWeightScaling = targetParam
+ * idealWeighthistUpdate = false
+ * Note: these variables could in principle be set to something else also for other target distribution types.
+ * However, localWeightScaling < 1 is in general expected to give lower efficiency and, except for local Boltzmann,
+ * idealWeightHistUpdate = false gives (in my experience) unstable, non-converging results.
+ *
+ * \param[in] awhParams AWH parameters.
+ * \param[in] awhBiasParams Bias parameters.
+ * \param[in] dimParams Bias dimension parameters.
+ * \param[in] beta 1/(k_B T) in units of 1/(kJ/mol), should be > 0.
+ * \param[in] mdTimeStep The MD time step.
+ * \param[in] numSharingSimulations The number of simulations to share the bias across.
+ * \param[in] gridAxis The grid axes.
+ * \param[in] disableUpdateSkips If to disable update skips, useful for testing.
+ * \param[in] biasIndex Index of the bias.
+ */
+ BiasParams(const AwhParams& awhParams,
+ const AwhBiasParams& awhBiasParams,
+ const std::vector<DimParams>& dimParams,
+ double beta,
+ double mdTimeStep,
+ DisableUpdateSkips disableUpdateSkips,
+ int numSharingSimulations,
+ const std::vector<GridAxis>& gridAxis,
+ int biasIndex);
+
+ /* Data members */
+ const double invBeta; /**< 1/beta = kT in kJ/mol */
+private:
+ const int64_t numStepsSampleCoord_; /**< Number of steps per coordinate value sample. */
+public:
+ const int numSamplesUpdateFreeEnergy_; /**< Number of samples per free energy update. */
+private:
+ const int64_t numStepsUpdateTarget_; /**< Number of steps per updating the target distribution. */
+ const int64_t numStepsCheckCovering_; /**< Number of steps per checking for covering. */
+public:
+ const int eTarget; /**< Type of target distribution. */
+ const double freeEnergyCutoffInKT; /**< Free energy cut-off in kT for cut-off target distribution. */
+ const double temperatureScaleFactor; /**< Temperature scaling factor for temperature scaled targed distributions. */
+ const bool idealWeighthistUpdate; /**< Update reference weighthistogram using the target distribution? Otherwise use the realized distribution. */
+ const int numSharedUpdate; /**< The number of (multi-)simulations sharing the bias update */
+ const double updateWeight; /**< The probability weight accumulated for each update. */
+ const double localWeightScaling; /**< Scaling factor applied to a sample before adding it to the reference weight histogram (= 1, usually). */
+ const double initialErrorInKT; /**< Estimated initial free energy error in kT. */
+ const double initialHistogramSize; /**< Initial reference weight histogram size. */
+private:
+ awh_ivec coverRadius_; /**< The radius (in points) that needs to be sampled around a point before it is considered covered. */
+public:
+ const bool convolveForce; /**< True if we convolve the force, false means use MC between umbrellas. */
+ const int biasIndex; /**< Index of the bias, used as a second random seed and for priting. */
+private:
+ const bool disableUpdateSkips_; /**< If true, we disallow update skips, even when the method supports it. */
};
-} // namespace gmx
+} // namespace gmx
#endif /* GMX_AWH_BIASPARAMS_H */
--- /dev/null
- GMX_RELEASE_ASSERT(intervalLength < (1 + margin) * boxLength,
- "We have checked before that interval <= period");
+/*
+ * 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,2018,2019, 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.
+ */
+#include "gmxpre.h"
+
+#include "read_params.h"
+
+#include "gromacs/awh/awh.h"
+#include "gromacs/fileio/readinp.h"
+#include "gromacs/fileio/warninp.h"
+#include "gromacs/math/units.h"
+#include "gromacs/math/utilities.h"
+#include "gromacs/math/vec.h"
+#include "gromacs/mdtypes/awh_params.h"
+#include "gromacs/mdtypes/inputrec.h"
+#include "gromacs/mdtypes/md_enums.h"
+#include "gromacs/mdtypes/pull_params.h"
+#include "gromacs/pbcutil/pbc.h"
+#include "gromacs/pulling/pull.h"
+#include "gromacs/random/seed.h"
+#include "gromacs/utility/cstringutil.h"
+#include "gromacs/utility/fatalerror.h"
+#include "gromacs/utility/smalloc.h"
+#include "gromacs/utility/stringutil.h"
+
+#include "biasparams.h"
+#include "biassharing.h"
+
+namespace gmx
+{
+
+const char* eawhtarget_names[eawhtargetNR + 1] = { "constant", "cutoff", "boltzmann",
+ "local-boltzmann", nullptr };
+
+const char* eawhgrowth_names[eawhgrowthNR + 1] = { "exp-linear", "linear", nullptr };
+
+const char* eawhpotential_names[eawhpotentialNR + 1] = { "convolved", "umbrella", nullptr };
+
+const char* eawhcoordprovider_names[eawhcoordproviderNR + 1] = { "pull", nullptr };
+
+/*! \brief
+ * Read parameters of an AWH bias dimension.
+ *
+ * \param[in,out] inp Input file entries.
+ * \param[in] prefix Prefix for dimension parameters.
+ * \param[in,out] dimParams AWH dimensional parameters.
+ * \param[in] pull_params Pull parameters.
+ * \param[in,out] wi Struct for bookeeping warnings.
+ * \param[in] bComment True if comments should be printed.
+ */
+static void readDimParams(std::vector<t_inpfile>* inp,
+ const std::string& prefix,
+ AwhDimParams* dimParams,
+ const pull_params_t* pull_params,
+ warninp_t wi,
+ bool bComment)
+{
+ std::string opt;
+ if (bComment)
+ {
+ printStringNoNewline(
+ inp, "The provider of the reaction coordinate, currently only pull is supported");
+ }
+
+ opt = prefix + "-coord-provider";
+ dimParams->eCoordProvider = get_eeenum(inp, opt, eawhcoordprovider_names, wi);
+
+ if (bComment)
+ {
+ printStringNoNewline(inp, "The coordinate index for this dimension");
+ }
+ opt = prefix + "-coord-index";
+ int coordIndexInput;
+ coordIndexInput = get_eint(inp, opt, 1, wi);
+ if (coordIndexInput < 1)
+ {
+ gmx_fatal(FARGS,
+ "Failed to read a valid coordinate index for %s. "
+ "Note that the pull coordinate indexing starts at 1.",
+ opt.c_str());
+ }
+
+ /* The pull coordinate indices start at 1 in the input file, at 0 internally */
+ dimParams->coordIndex = coordIndexInput - 1;
+
+ /* The pull settings need to be consistent with the AWH settings */
+ if (!(pull_params->coord[dimParams->coordIndex].eType == epullEXTERNAL))
+ {
+ gmx_fatal(FARGS, "AWH biasing can only be applied to pull type %s", EPULLTYPE(epullEXTERNAL));
+ }
+
+ if (dimParams->coordIndex >= pull_params->ncoord)
+ {
+ gmx_fatal(FARGS,
+ "The given AWH coordinate index (%d) is larger than the number of pull "
+ "coordinates (%d)",
+ coordIndexInput, pull_params->ncoord);
+ }
+ if (pull_params->coord[dimParams->coordIndex].rate != 0)
+ {
+ auto message = formatString(
+ "Setting pull-coord%d-rate (%g) is incompatible with AWH biasing this coordinate",
+ coordIndexInput, pull_params->coord[dimParams->coordIndex].rate);
+ warning_error(wi, message);
+ }
+
+ /* Grid params for each axis */
+ int eGeom = pull_params->coord[dimParams->coordIndex].eGeom;
+
+ if (bComment)
+ {
+ printStringNoNewline(inp, "Start and end values for each coordinate dimension");
+ }
+
+ opt = prefix + "-start";
+ dimParams->origin = get_ereal(inp, opt, 0., wi);
+
+ opt = prefix + "-end";
+ dimParams->end = get_ereal(inp, opt, 0., wi);
+
+ if (gmx_within_tol(dimParams->end - dimParams->origin, 0, GMX_REAL_EPS))
+ {
+ auto message = formatString(
+ "The given interval length given by %s-start (%g) and %s-end (%g) is zero. "
+ "This will result in only one point along this axis in the coordinate value grid.",
+ prefix.c_str(), dimParams->origin, prefix.c_str(), dimParams->end);
+ warning(wi, message);
+ }
+ /* Check that the requested interval is in allowed range */
+ if (eGeom == epullgDIST)
+ {
+ if (dimParams->origin < 0 || dimParams->end < 0)
+ {
+ gmx_fatal(FARGS,
+ "%s-start (%g) or %s-end (%g) set to a negative value. With pull geometry "
+ "distance coordinate values are non-negative. "
+ "Perhaps you want to use geometry %s instead?",
+ prefix.c_str(), dimParams->origin, prefix.c_str(), dimParams->end,
+ EPULLGEOM(epullgDIR));
+ }
+ }
+ else if (eGeom == epullgANGLE || eGeom == epullgANGLEAXIS)
+ {
+ if (dimParams->origin < 0 || dimParams->end > 180)
+ {
+ gmx_fatal(FARGS,
+ "%s-start (%g) and %s-end (%g) are outside of the allowed range 0 to 180 deg "
+ "for pull geometries %s and %s ",
+ prefix.c_str(), dimParams->origin, prefix.c_str(), dimParams->end,
+ EPULLGEOM(epullgANGLE), EPULLGEOM(epullgANGLEAXIS));
+ }
+ }
+ else if (eGeom == epullgDIHEDRAL)
+ {
+ if (dimParams->origin < -180 || dimParams->end > 180)
+ {
+ gmx_fatal(FARGS,
+ "%s-start (%g) and %s-end (%g) are outside of the allowed range -180 to 180 "
+ "deg for pull geometry %s. ",
+ prefix.c_str(), dimParams->origin, prefix.c_str(), dimParams->end,
+ EPULLGEOM(epullgDIHEDRAL));
+ }
+ }
+
+ if (bComment)
+ {
+ printStringNoNewline(
+ inp, "The force constant for this coordinate (kJ/mol/nm^2 or kJ/mol/rad^2)");
+ }
+ opt = prefix + "-force-constant";
+ dimParams->forceConstant = get_ereal(inp, opt, 0, wi);
+ if (dimParams->forceConstant <= 0)
+ {
+ warning_error(wi, "The force AWH bias force constant should be > 0");
+ }
+
+ if (bComment)
+ {
+ printStringNoNewline(inp, "Estimated diffusion constant (nm^2/ps or rad^2/ps)");
+ }
+ opt = prefix + "-diffusion";
+ dimParams->diffusion = get_ereal(inp, opt, 0, wi);
+
+ if (dimParams->diffusion <= 0)
+ {
+ const double diffusion_default = 1e-5;
+ auto message = formatString(
+ "%s not explicitly set by user. You can choose to use a default "
+ "value (%g nm^2/ps or rad^2/ps) but this may very well be "
+ "non-optimal for your system!",
+ opt.c_str(), diffusion_default);
+ warning(wi, message);
+ dimParams->diffusion = diffusion_default;
+ }
+
+ if (bComment)
+ {
+ printStringNoNewline(inp,
+ "Diameter that needs to be sampled around a point before it is "
+ "considered covered.");
+ }
+ opt = prefix + "-cover-diameter";
+ dimParams->coverDiameter = get_ereal(inp, opt, 0, wi);
+
+ if (dimParams->coverDiameter < 0)
+ {
+ gmx_fatal(FARGS, "%s (%g) cannot be negative.", opt.c_str(), dimParams->coverDiameter);
+ }
+}
+
+/*! \brief
+ * Check consistency of input at the AWH bias level.
+ *
+ * \param[in] awhBiasParams AWH bias parameters.
+ * \param[in,out] wi Struct for bookkeeping warnings.
+ */
+static void checkInputConsistencyAwhBias(const AwhBiasParams& awhBiasParams, warninp_t wi)
+{
+ /* Covering diameter and sharing warning. */
+ for (int d = 0; d < awhBiasParams.ndim; d++)
+ {
+ double coverDiameter = awhBiasParams.dimParams[d].coverDiameter;
+ if (awhBiasParams.shareGroup <= 0 && coverDiameter > 0)
+ {
+ warning(wi,
+ "The covering diameter is only relevant to set for bias sharing simulations.");
+ }
+ }
+}
+
+/*! \brief
+ * Read parameters of an AWH bias.
+ *
+ * \param[in,out] inp Input file entries.
+ * \param[in,out] awhBiasParams AWH dimensional parameters.
+ * \param[in] prefix Prefix for bias parameters.
+ * \param[in] ir Input parameter struct.
+ * \param[in,out] wi Struct for bookeeping warnings.
+ * \param[in] bComment True if comments should be printed.
+ */
+static void read_bias_params(std::vector<t_inpfile>* inp,
+ AwhBiasParams* awhBiasParams,
+ const std::string& prefix,
+ const t_inputrec* ir,
+ warninp_t wi,
+ bool bComment)
+{
+ if (bComment)
+ {
+ printStringNoNewline(inp, "Estimated initial PMF error (kJ/mol)");
+ }
+
+ std::string opt = prefix + "-error-init";
+ /* We allow using a default value here without warning (but warn the user if the diffusion constant is not set). */
+ awhBiasParams->errorInitial = get_ereal(inp, opt, 10, wi);
+ if (awhBiasParams->errorInitial <= 0)
+ {
+ gmx_fatal(FARGS, "%s needs to be > 0.", opt.c_str());
+ }
+
+ if (bComment)
+ {
+ printStringNoNewline(inp,
+ "Growth rate of the reference histogram determining the bias update "
+ "size: exp-linear or linear");
+ }
+ opt = prefix + "-growth";
+ awhBiasParams->eGrowth = get_eeenum(inp, opt, eawhgrowth_names, wi);
+
+ if (bComment)
+ {
+ printStringNoNewline(inp,
+ "Start the simulation by equilibrating histogram towards the target "
+ "distribution: no or yes");
+ }
+ opt = prefix + "-equilibrate-histogram";
+ awhBiasParams->equilibrateHistogram = (get_eeenum(inp, opt, yesno_names, wi) != 0);
+ if (awhBiasParams->equilibrateHistogram && awhBiasParams->eGrowth != eawhgrowthEXP_LINEAR)
+ {
+ auto message =
+ formatString("Option %s will only have an effect for histogram growth type '%s'.",
+ opt.c_str(), EAWHGROWTH(eawhgrowthEXP_LINEAR));
+ warning(wi, message);
+ }
+
+ if (bComment)
+ {
+ printStringNoNewline(
+ inp, "Target distribution type: constant, cutoff, boltzmann or local-boltzmann");
+ }
+ opt = prefix + "-target";
+ awhBiasParams->eTarget = get_eeenum(inp, opt, eawhtarget_names, wi);
+
+ if ((awhBiasParams->eTarget == eawhtargetLOCALBOLTZMANN)
+ && (awhBiasParams->eGrowth == eawhgrowthEXP_LINEAR))
+ {
+ auto message = formatString(
+ "Target type '%s' combined with histogram growth type '%s' is not "
+ "expected to give stable bias updates. You probably want to use growth type "
+ "'%s' instead.",
+ EAWHTARGET(eawhtargetLOCALBOLTZMANN), EAWHGROWTH(eawhgrowthEXP_LINEAR),
+ EAWHGROWTH(eawhgrowthLINEAR));
+ warning(wi, message);
+ }
+
+ if (bComment)
+ {
+ printStringNoNewline(inp,
+ "Boltzmann beta scaling factor for target distribution types "
+ "'boltzmann' and 'boltzmann-local'");
+ }
+ opt = prefix + "-target-beta-scaling";
+ awhBiasParams->targetBetaScaling = get_ereal(inp, opt, 0, wi);
+
+ switch (awhBiasParams->eTarget)
+ {
+ case eawhtargetBOLTZMANN:
+ case eawhtargetLOCALBOLTZMANN:
+ if (awhBiasParams->targetBetaScaling < 0 || awhBiasParams->targetBetaScaling > 1)
+ {
+ gmx_fatal(FARGS, "%s = %g is not useful for target type %s.", opt.c_str(),
+ awhBiasParams->targetBetaScaling, EAWHTARGET(awhBiasParams->eTarget));
+ }
+ break;
+ default:
+ if (awhBiasParams->targetBetaScaling != 0)
+ {
+ gmx_fatal(
+ FARGS,
+ "Value for %s (%g) set explicitly but will not be used for target type %s.",
+ opt.c_str(), awhBiasParams->targetBetaScaling,
+ EAWHTARGET(awhBiasParams->eTarget));
+ }
+ break;
+ }
+
+ if (bComment)
+ {
+ printStringNoNewline(inp, "Free energy cutoff value for target distribution type 'cutoff'");
+ }
+ opt = prefix + "-target-cutoff";
+ awhBiasParams->targetCutoff = get_ereal(inp, opt, 0, wi);
+
+ switch (awhBiasParams->eTarget)
+ {
+ case eawhtargetCUTOFF:
+ if (awhBiasParams->targetCutoff <= 0)
+ {
+ gmx_fatal(FARGS, "%s = %g is not useful for target type %s.", opt.c_str(),
+ awhBiasParams->targetCutoff, EAWHTARGET(awhBiasParams->eTarget));
+ }
+ break;
+ default:
+ if (awhBiasParams->targetCutoff != 0)
+ {
+ gmx_fatal(
+ FARGS,
+ "Value for %s (%g) set explicitly but will not be used for target type %s.",
+ opt.c_str(), awhBiasParams->targetCutoff, EAWHTARGET(awhBiasParams->eTarget));
+ }
+ break;
+ }
+
+ if (bComment)
+ {
+ printStringNoNewline(inp, "Initialize PMF and target with user data: no or yes");
+ }
+ opt = prefix + "-user-data";
+ awhBiasParams->bUserData = get_eeenum(inp, opt, yesno_names, wi);
+
+ if (bComment)
+ {
+ printStringNoNewline(inp, "Group index to share the bias with, 0 means not shared");
+ }
+ opt = prefix + "-share-group";
+ awhBiasParams->shareGroup = get_eint(inp, opt, 0, wi);
+ if (awhBiasParams->shareGroup < 0)
+ {
+ warning_error(wi, "AWH bias share-group should be >= 0");
+ }
+
+ if (bComment)
+ {
+ printStringNoNewline(inp, "Dimensionality of the coordinate");
+ }
+ opt = prefix + "-ndim";
+ awhBiasParams->ndim = get_eint(inp, opt, 0, wi);
+
+ if (awhBiasParams->ndim <= 0 || awhBiasParams->ndim > c_biasMaxNumDim)
+ {
+ gmx_fatal(FARGS, "%s (%d) needs to be > 0 and at most %d\n", opt.c_str(),
+ awhBiasParams->ndim, c_biasMaxNumDim);
+ }
+ if (awhBiasParams->ndim > 2)
+ {
+ warning_note(wi,
+ "For awh-dim > 2 the estimate based on the diffusion and the initial error is "
+ "currently only a rough guideline."
+ " You should verify its usefulness for your system before production runs!");
+ }
+ snew(awhBiasParams->dimParams, awhBiasParams->ndim);
+ for (int d = 0; d < awhBiasParams->ndim; d++)
+ {
+ bComment = bComment && d == 0;
+ std::string prefixdim = prefix + formatString("-dim%d", d + 1);
+ readDimParams(inp, prefixdim, &awhBiasParams->dimParams[d], ir->pull, wi, bComment);
+ }
+
+ /* Check consistencies here that cannot be checked at read time at a lower level. */
+ checkInputConsistencyAwhBias(*awhBiasParams, wi);
+}
+
+/*! \brief
+ * Check consistency of input at the AWH level.
+ *
+ * \param[in] awhParams AWH parameters.
+ * \param[in,out] wi Struct for bookkeeping warnings.
+ */
+static void checkInputConsistencyAwh(const AwhParams& awhParams, warninp_t wi)
+{
+ /* Each pull coord can map to at most 1 AWH coord.
+ * Check that we have a shared bias when requesting multisim sharing.
+ */
+ bool haveSharedBias = false;
+ for (int k1 = 0; k1 < awhParams.numBias; k1++)
+ {
+ const AwhBiasParams& awhBiasParams1 = awhParams.awhBiasParams[k1];
+
+ if (awhBiasParams1.shareGroup > 0)
+ {
+ haveSharedBias = true;
+ }
+
+ /* k1 is the reference AWH, k2 is the AWH we compare with (can be equal to k1) */
+ for (int k2 = k1; k2 < awhParams.numBias; k2++)
+ {
+ for (int d1 = 0; d1 < awhBiasParams1.ndim; d1++)
+ {
+ const AwhBiasParams& awhBiasParams2 = awhParams.awhBiasParams[k2];
+
+ /* d1 is the reference dimension of the reference AWH. d2 is the dim index of the AWH to compare with. */
+ for (int d2 = 0; d2 < awhBiasParams2.ndim; d2++)
+ {
+ /* Give an error if (d1, k1) is different from (d2, k2) but the pull coordinate is the same */
+ if ((d1 != d2 || k1 != k2)
+ && (awhBiasParams1.dimParams[d1].coordIndex == awhBiasParams2.dimParams[d2].coordIndex))
+ {
+ char errormsg[STRLEN];
+ sprintf(errormsg,
+ "One pull coordinate (%d) cannot be mapped to two separate AWH "
+ "dimensions (awh%d-dim%d and awh%d-dim%d). "
+ "If this is really what you want to do you will have to duplicate "
+ "this pull coordinate.",
+ awhBiasParams1.dimParams[d1].coordIndex + 1, k1 + 1, d1 + 1, k2 + 1,
+ d2 + 1);
+ gmx_fatal(FARGS, "%s", errormsg);
+ }
+ }
+ }
+ }
+ }
+
+ if (awhParams.shareBiasMultisim && !haveSharedBias)
+ {
+ warning(wi,
+ "Sharing of biases over multiple simulations is requested, but no bias is marked "
+ "as shared (share-group > 0)");
+ }
+
+ /* mdrun does not support this (yet), but will check again */
+ if (haveBiasSharingWithinSimulation(awhParams))
+ {
+ warning(wi,
+ "You have shared biases within a single simulation, but mdrun does not support "
+ "this (yet)");
+ }
+}
+
+AwhParams* readAndCheckAwhParams(std::vector<t_inpfile>* inp, const t_inputrec* ir, warninp_t wi)
+{
+ AwhParams* awhParams;
+ snew(awhParams, 1);
+ std::string opt;
+
+ /* Parameters common for all biases */
+
+ printStringNoNewline(inp, "The way to apply the biasing potential: convolved or umbrella");
+ opt = "awh-potential";
+ awhParams->ePotential = get_eeenum(inp, opt, eawhpotential_names, wi);
+
+ printStringNoNewline(inp,
+ "The random seed used for sampling the umbrella center in the case of "
+ "umbrella type potential");
+ opt = "awh-seed";
+ awhParams->seed = get_eint(inp, opt, -1, wi);
+ if (awhParams->seed == -1)
+ {
+ awhParams->seed = static_cast<int>(gmx::makeRandomSeed());
+ fprintf(stderr, "Setting the AWH bias MC random seed to %" PRId64 "\n", awhParams->seed);
+ }
+
+ printStringNoNewline(inp, "Data output interval in number of steps");
+ opt = "awh-nstout";
+ awhParams->nstOut = get_eint(inp, opt, 100000, wi);
+ if (awhParams->nstOut <= 0)
+ {
+ auto message = formatString("Not writing AWH output with AWH (%s = %d) does not make sense",
+ opt.c_str(), awhParams->nstOut);
+ warning_error(wi, message);
+ }
+ /* This restriction can be removed by changing a flag of print_ebin() */
+ if (ir->nstenergy == 0 || awhParams->nstOut % ir->nstenergy != 0)
+ {
+ auto message = formatString("%s (%d) should be a multiple of nstenergy (%d)", opt.c_str(),
+ awhParams->nstOut, ir->nstenergy);
+ warning_error(wi, message);
+ }
+
+ printStringNoNewline(inp, "Coordinate sampling interval in number of steps");
+ opt = "awh-nstsample";
+ awhParams->nstSampleCoord = get_eint(inp, opt, 10, wi);
+
+ printStringNoNewline(inp, "Free energy and bias update interval in number of samples");
+ opt = "awh-nsamples-update";
+ awhParams->numSamplesUpdateFreeEnergy = get_eint(inp, opt, 10, wi);
+ if (awhParams->numSamplesUpdateFreeEnergy <= 0)
+ {
+ warning_error(wi, opt + " needs to be an integer > 0");
+ }
+
+ printStringNoNewline(
+ inp, "When true, biases with share-group>0 are shared between multiple simulations");
+ opt = "awh-share-multisim";
+ awhParams->shareBiasMultisim = (get_eeenum(inp, opt, yesno_names, wi) != 0);
+
+ printStringNoNewline(inp, "The number of independent AWH biases");
+ opt = "awh-nbias";
+ awhParams->numBias = get_eint(inp, opt, 1, wi);
+ if (awhParams->numBias <= 0)
+ {
+ gmx_fatal(FARGS, "%s needs to be an integer > 0", opt.c_str());
+ }
+
+ /* Read the parameters specific to each AWH bias */
+ snew(awhParams->awhBiasParams, awhParams->numBias);
+
+ for (int k = 0; k < awhParams->numBias; k++)
+ {
+ bool bComment = (k == 0);
+ std::string prefixawh = formatString("awh%d", k + 1);
+ read_bias_params(inp, &awhParams->awhBiasParams[k], prefixawh, ir, wi, bComment);
+ }
+
+ /* Do a final consistency check before returning */
+ checkInputConsistencyAwh(*awhParams, wi);
+
+ if (ir->init_step != 0)
+ {
+ warning_error(wi, "With AWH init-step should be 0");
+ }
+
+ return awhParams;
+}
+
+/*! \brief
+ * Gets the period of a pull coordinate.
+ *
+ * \param[in] pullCoordParams The parameters for the pull coordinate.
+ * \param[in] pbc The PBC setup
+ * \param[in] intervalLength The length of the AWH interval for this pull coordinate
+ * \returns the period (or 0 if not periodic).
+ */
+static double get_pull_coord_period(const t_pull_coord& pullCoordParams, const t_pbc& pbc, const real intervalLength)
+{
+ double period = 0;
+
+ if (pullCoordParams.eGeom == epullgDIR)
+ {
+ const real margin = 0.001;
+ // Make dims periodic when the interval covers > 95%
+ const real periodicFraction = 0.95;
+
+ // Check if the pull direction is along a box vector
+ for (int dim = 0; dim < pbc.ndim_ePBC; dim++)
+ {
+ const real boxLength = norm(pbc.box[dim]);
+ const real innerProduct = iprod(pullCoordParams.vec, pbc.box[dim]);
+ if (innerProduct >= (1 - margin) * boxLength && innerProduct <= (1 + margin) * boxLength)
+ {
++ if (intervalLength > (1 + margin) * boxLength)
++ {
++ gmx_fatal(FARGS,
++ "The AWH interval (%f nm) for a pull coordinate is larger than the "
++ "box size (%f nm)",
++ intervalLength, boxLength);
++ }
++
+ if (intervalLength > periodicFraction * boxLength)
+ {
+ period = boxLength;
+ }
+ }
+ }
+ }
+ else if (pullCoordParams.eGeom == epullgDIHEDRAL)
+ {
+ /* The dihedral angle is periodic in -180 to 180 deg */
+ period = 360;
+ }
+
+ return period;
+}
+
+/*! \brief
+ * Checks if the given interval is defined in the correct periodic interval.
+ *
+ * \param[in] origin Start value of interval.
+ * \param[in] end End value of interval.
+ * \param[in] period Period (or 0 if not periodic).
+ * \returns true if the end point values are in the correct periodic interval.
+ */
+static bool intervalIsInPeriodicInterval(double origin, double end, double period)
+{
+ return (period == 0) || (std::fabs(origin) <= 0.5 * period && std::fabs(end) <= 0.5 * period);
+}
+
+/*! \brief
+ * Checks if a value is within an interval.
+ *
+ * \param[in] origin Start value of interval.
+ * \param[in] end End value of interval.
+ * \param[in] period Period (or 0 if not periodic).
+ * \param[in] value Value to check.
+ * \returns true if the value is within the interval.
+ */
+static bool valueIsInInterval(double origin, double end, double period, double value)
+{
+ bool bIn_interval;
+
+ if (period > 0)
+ {
+ if (origin < end)
+ {
+ /* The interval closes within the periodic interval */
+ bIn_interval = (value >= origin) && (value <= end);
+ }
+ else
+ {
+ /* The interval wraps around the periodic boundary */
+ bIn_interval = ((value >= origin) && (value <= 0.5 * period))
+ || ((value >= -0.5 * period) && (value <= end));
+ }
+ }
+ else
+ {
+ bIn_interval = (value >= origin) && (value <= end);
+ }
+
+ return bIn_interval;
+}
+
+/*! \brief
+ * Check if the starting configuration is consistent with the given interval.
+ *
+ * \param[in] awhParams AWH parameters.
+ * \param[in,out] wi Struct for bookeeping warnings.
+ */
+static void checkInputConsistencyInterval(const AwhParams* awhParams, warninp_t wi)
+{
+ for (int k = 0; k < awhParams->numBias; k++)
+ {
+ AwhBiasParams* awhBiasParams = &awhParams->awhBiasParams[k];
+ for (int d = 0; d < awhBiasParams->ndim; d++)
+ {
+ AwhDimParams* dimParams = &awhBiasParams->dimParams[d];
+ int coordIndex = dimParams->coordIndex;
+ double origin = dimParams->origin, end = dimParams->end, period = dimParams->period;
+ double coordValueInit = dimParams->coordValueInit;
+
+ if ((period == 0) && (origin > end))
+ {
+ gmx_fatal(FARGS,
+ "For the non-periodic pull coordinates awh%d-dim%d-start (%f) cannot be "
+ "larger than awh%d-dim%d-end (%f)",
+ k + 1, d + 1, origin, k + 1, d + 1, end);
+ }
+
+ /* Currently we assume symmetric periodic intervals, meaning we use [-period/2, period/2] as the reference interval.
+ Make sure the AWH interval is within this reference interval.
+
+ Note: we could fairly simply allow using a more general interval (e.g. [x, x + period]) but it complicates
+ things slightly and I don't see that there is a great need for it. It would also mean that the interval would
+ depend on AWH input. Also, for dihedral angles you would always want the reference interval to be -180, +180,
+ independent of AWH parameters.
+ */
+ if (!intervalIsInPeriodicInterval(origin, end, period))
+ {
+ gmx_fatal(FARGS,
+ "When using AWH with periodic pull coordinate geometries "
+ "awh%d-dim%d-start (%.8g) and "
+ "awh%d-dim%d-end (%.8g) should cover at most one period (%.8g) and take "
+ "values in between "
+ "minus half a period and plus half a period, i.e. in the interval [%.8g, "
+ "%.8g].",
+ k + 1, d + 1, origin, k + 1, d + 1, end, period, -0.5 * period, 0.5 * period);
+ }
+
+ /* Warn if the pull initial coordinate value is not in the grid */
+ if (!valueIsInInterval(origin, end, period, coordValueInit))
+ {
+ auto message = formatString(
+ "The initial coordinate value (%.8g) for pull coordinate index %d falls "
+ "outside "
+ "of the sampling nterval awh%d-dim%d-start (%.8g) to awh%d-dim%d-end "
+ "(%.8g). "
+ "This can lead to large initial forces pulling the coordinate towards the "
+ "sampling interval.",
+ coordValueInit, coordIndex + 1, k + 1, d + 1, origin, k + 1, d + 1, end);
+ warning(wi, message);
+ }
+ }
+ }
+}
+
+void setStateDependentAwhParams(AwhParams* awhParams,
+ const pull_params_t* pull_params,
+ pull_t* pull_work,
+ const matrix box,
+ int ePBC,
+ const tensor& compressibility,
+ const t_grpopts* inputrecGroupOptions,
+ warninp_t wi)
+{
+ /* The temperature is not really state depenendent but is not known
+ * when read_awhParams is called (in get ir).
+ * It is known first after do_index has been called in grompp.cpp.
+ */
+ if (inputrecGroupOptions->ref_t == nullptr || inputrecGroupOptions->ref_t[0] <= 0)
+ {
+ gmx_fatal(FARGS, "AWH biasing is only supported for temperatures > 0");
+ }
+ for (int i = 1; i < inputrecGroupOptions->ngtc; i++)
+ {
+ if (inputrecGroupOptions->ref_t[i] != inputrecGroupOptions->ref_t[0])
+ {
+ gmx_fatal(FARGS,
+ "AWH biasing is currently only supported for identical temperatures for all "
+ "temperature coupling groups");
+ }
+ }
+
+ t_pbc pbc;
+ set_pbc(&pbc, ePBC, box);
+
+ for (int k = 0; k < awhParams->numBias; k++)
+ {
+ AwhBiasParams* awhBiasParams = &awhParams->awhBiasParams[k];
+ for (int d = 0; d < awhBiasParams->ndim; d++)
+ {
+ AwhDimParams* dimParams = &awhBiasParams->dimParams[d];
+ const t_pull_coord& pullCoordParams = pull_params->coord[dimParams->coordIndex];
+
+ if (pullCoordParams.eGeom == epullgDIRPBC)
+ {
+ gmx_fatal(FARGS,
+ "AWH does not support pull geometry '%s'. "
+ "If the maximum distance between the groups is always less than half the "
+ "box size, "
+ "you can use geometry '%s' instead.",
+ EPULLGEOM(epullgDIRPBC), EPULLGEOM(epullgDIR));
+ }
+
+ dimParams->period =
+ get_pull_coord_period(pullCoordParams, pbc, dimParams->end - dimParams->origin);
+ // We would like to check for scaling, but we don't have the full inputrec available here
+ if (dimParams->period > 0
+ && !(pullCoordParams.eGeom == epullgANGLE || pullCoordParams.eGeom == epullgDIHEDRAL))
+ {
+ bool coordIsScaled = false;
+ for (int d2 = 0; d2 < DIM; d2++)
+ {
+ if (pullCoordParams.vec[d2] != 0 && norm2(compressibility[d2]) != 0)
+ {
+ coordIsScaled = true;
+ }
+ }
+ if (coordIsScaled)
+ {
+ std::string mesg = gmx::formatString(
+ "AWH dimension %d of bias %d is periodic with pull geometry '%s', "
+ "while you should are applying pressure scaling to the corresponding "
+ "box vector, this is not supported.",
+ d + 1, k + 1, EPULLGEOM(pullCoordParams.eGeom));
+ warning(wi, mesg.c_str());
+ }
+ }
+
+ /* The initial coordinate value, converted to external user units. */
+ dimParams->coordValueInit = get_pull_coord_value(pull_work, dimParams->coordIndex, &pbc);
+
+ dimParams->coordValueInit *= pull_conversion_factor_internal2userinput(&pullCoordParams);
+ }
+ }
+ checkInputConsistencyInterval(awhParams, wi);
+
+ /* Register AWH as external potential with pull to check consistency. */
+ Awh::registerAwhWithPull(*awhParams, pull_work);
+}
+
+} // namespace gmx
#include "gromacs/utility/smalloc.h"
#include "gromacs/utility/snprintf.h"
-typedef struct {
+typedef struct
+{
int ai, aj;
} gmx_conection_t;
-typedef struct gmx_conect_t {
+typedef struct gmx_conect_t
+{
int nconect;
gmx_bool bSorted;
- gmx_conection_t *conect;
+ gmx_conection_t* conect;
} gmx_conect_t;
-static const char *pdbtp[epdbNR] = {
- "ATOM ", "HETATM", "ANISOU", "CRYST1",
- "COMPND", "MODEL", "ENDMDL", "TER", "HEADER", "TITLE", "REMARK",
- "CONECT"
-};
-
+static const char* pdbtp[epdbNR] = { "ATOM ", "HETATM", "ANISOU", "CRYST1", "COMPND", "MODEL",
+ "ENDMDL", "TER", "HEADER", "TITLE", "REMARK", "CONECT" };
#define REMARK_SIM_BOX "REMARK THIS IS A SIMULATION BOX"
-static void xlate_atomname_pdb2gmx(char *name)
+static void xlate_atomname_pdb2gmx(char* name)
{
int i, length;
char temp;
temp = name[0];
for (i = 1; i < length; i++)
{
- name[i-1] = name[i];
+ name[i - 1] = name[i];
}
- name[length-1] = temp;
+ name[length - 1] = temp;
}
}
static std::string xlate_atomname_gmx2pdb(std::string name)
{
size_t length = name.size();
- if (length > 3 && std::isdigit(name[length-1]))
+ if (length > 3 && std::isdigit(name[length - 1]))
{
- char temp = name[length-1];
- for (size_t i = length-1; i > 0; --i)
+ char temp = name[length - 1];
+ for (size_t i = length - 1; i > 0; --i)
{
- name[i] = name[i-1];
+ name[i] = name[i - 1];
}
name[0] = temp;
}
}
-void gmx_write_pdb_box(FILE *out, int ePBC, const matrix box)
+void gmx_write_pdb_box(FILE* out, int ePBC, const matrix box)
{
real alpha, beta, gamma;
return;
}
- if (norm2(box[YY])*norm2(box[ZZ]) != 0)
+ if (norm2(box[YY]) * norm2(box[ZZ]) != 0)
{
- alpha = RAD2DEG*gmx_angle(box[YY], box[ZZ]);
+ alpha = RAD2DEG * gmx_angle(box[YY], box[ZZ]);
}
else
{
alpha = 90;
}
- if (norm2(box[XX])*norm2(box[ZZ]) != 0)
+ if (norm2(box[XX]) * norm2(box[ZZ]) != 0)
{
- beta = RAD2DEG*gmx_angle(box[XX], box[ZZ]);
+ beta = RAD2DEG * gmx_angle(box[XX], box[ZZ]);
}
else
{
- beta = 90;
+ beta = 90;
}
- if (norm2(box[XX])*norm2(box[YY]) != 0)
+ if (norm2(box[XX]) * norm2(box[YY]) != 0)
{
- gamma = RAD2DEG*gmx_angle(box[XX], box[YY]);
+ gamma = RAD2DEG * gmx_angle(box[XX], box[YY]);
}
else
{
fprintf(out, "REMARK THIS IS A SIMULATION BOX\n");
if (ePBC != epbcSCREW)
{
- fprintf(out, "CRYST1%9.3f%9.3f%9.3f%7.2f%7.2f%7.2f %-11s%4d\n",
- 10*norm(box[XX]), 10*norm(box[YY]), 10*norm(box[ZZ]),
- alpha, beta, gamma, "P 1", 1);
+ fprintf(out, "CRYST1%9.3f%9.3f%9.3f%7.2f%7.2f%7.2f %-11s%4d\n", 10 * norm(box[XX]),
+ 10 * norm(box[YY]), 10 * norm(box[ZZ]), alpha, beta, gamma, "P 1", 1);
}
else
{
/* Double the a-vector length and write the correct space group */
- fprintf(out, "CRYST1%9.3f%9.3f%9.3f%7.2f%7.2f%7.2f %-11s%4d\n",
- 20*norm(box[XX]), 10*norm(box[YY]), 10*norm(box[ZZ]),
- alpha, beta, gamma, "P 21 1 1", 1);
-
+ fprintf(out, "CRYST1%9.3f%9.3f%9.3f%7.2f%7.2f%7.2f %-11s%4d\n", 20 * norm(box[XX]),
+ 10 * norm(box[YY]), 10 * norm(box[ZZ]), alpha, beta, gamma, "P 21 1 1", 1);
}
}
-static void read_cryst1(char *line, int *ePBC, matrix box)
+static void read_cryst1(char* line, int* ePBC, matrix box)
{
#define SG_SIZE 11
- char sa[12], sb[12], sc[12], sg[SG_SIZE+1], ident;
+ char sa[12], sb[12], sc[12], sg[SG_SIZE + 1], ident;
double fa, fb, fc, alpha, beta, gamma, cosa, cosb, cosg, sing;
int syma, symb, symc;
int ePBC_file;
ePBC_file = -1;
if (strlen(line) >= 55)
{
- strncpy(sg, line+55, SG_SIZE);
+ strncpy(sg, line + 55, SG_SIZE);
sg[SG_SIZE] = '\0';
ident = ' ';
syma = 0;
symb = 0;
symc = 0;
sscanf(sg, "%c %d %d %d", &ident, &syma, &symb, &symc);
- if (ident == 'P' && syma == 1 && symb <= 1 && symc <= 1)
+ if (ident == 'P' && syma == 1 && symb <= 1 && symc <= 1)
{
- fc = strtod(sc, nullptr)*0.1;
+ fc = strtod(sc, nullptr) * 0.1;
ePBC_file = (fc > 0 ? epbcXYZ : epbcXY);
}
if (ident == 'P' && syma == 21 && symb == 1 && symc == 1)
if (box)
{
- fa = strtod(sa, nullptr)*0.1;
- fb = strtod(sb, nullptr)*0.1;
- fc = strtod(sc, nullptr)*0.1;
+ fa = strtod(sa, nullptr) * 0.1;
+ fb = strtod(sb, nullptr) * 0.1;
+ fc = strtod(sc, nullptr) * 0.1;
if (ePBC_file == epbcSCREW)
{
fa *= 0.5;
{
if (alpha != 90.0)
{
- cosa = std::cos(alpha*DEG2RAD);
+ cosa = std::cos(alpha * DEG2RAD);
}
else
{
}
if (beta != 90.0)
{
- cosb = std::cos(beta*DEG2RAD);
+ cosb = std::cos(beta * DEG2RAD);
}
else
{
}
if (gamma != 90.0)
{
- cosg = std::cos(gamma*DEG2RAD);
- sing = std::sin(gamma*DEG2RAD);
+ cosg = std::cos(gamma * DEG2RAD);
+ sing = std::sin(gamma * DEG2RAD);
}
else
{
cosg = 0;
sing = 1;
}
- box[YY][XX] = fb*cosg;
- box[YY][YY] = fb*sing;
- box[ZZ][XX] = fc*cosb;
- box[ZZ][YY] = fc*(cosa - cosb*cosg)/sing;
- box[ZZ][ZZ] = std::sqrt(fc*fc
- - box[ZZ][XX]*box[ZZ][XX] - box[ZZ][YY]*box[ZZ][YY]);
+ box[YY][XX] = fb * cosg;
+ box[YY][YY] = fb * sing;
+ box[ZZ][XX] = fc * cosb;
+ box[ZZ][YY] = fc * (cosa - cosb * cosg) / sing;
+ box[ZZ][ZZ] = std::sqrt(fc * fc - box[ZZ][XX] * box[ZZ][XX] - box[ZZ][YY] * box[ZZ][YY]);
}
else
{
}
}
-static int
-gmx_fprintf_pqr_atomline(FILE * fp,
- enum PDB_record record,
- int atom_seq_number,
- const char * atom_name,
- const char * res_name,
- char chain_id,
- int res_seq_number,
- real x,
- real y,
- real z,
- real occupancy,
- real b_factor)
+static int gmx_fprintf_pqr_atomline(FILE* fp,
+ enum PDB_record record,
+ int atom_seq_number,
+ const char* atom_name,
+ const char* res_name,
+ char chain_id,
+ int res_seq_number,
+ real x,
+ real y,
+ real z,
+ real occupancy,
+ real b_factor)
{
GMX_RELEASE_ASSERT(record == epdbATOM || record == epdbHETATM,
"Can only print PQR atom lines as ATOM or HETATM records");
/* Check atom name */
- GMX_RELEASE_ASSERT(atom_name != nullptr,
- "Need atom information to print pqr");
+ GMX_RELEASE_ASSERT(atom_name != nullptr, "Need atom information to print pqr");
/* Check residue name */
- GMX_RELEASE_ASSERT(res_name != nullptr,
- "Need residue information to print pqr");
+ GMX_RELEASE_ASSERT(res_name != nullptr, "Need residue information to print pqr");
/* Truncate integers so they fit */
atom_seq_number = atom_seq_number % 100000;
res_seq_number = res_seq_number % 10000;
- int n = fprintf(fp,
- "%-6s%5d %-4.4s%4.4s%c%4d %8.3f %8.3f %8.3f %6.2f %6.2f\n",
- pdbtp[record],
- atom_seq_number,
- atom_name,
- res_name,
- chain_id,
- res_seq_number,
- x, y, z,
- occupancy,
- b_factor);
+ int n = fprintf(fp, "%-6s%5d %-4.4s%4.4s%c%4d %8.3f %8.3f %8.3f %6.2f %6.2f\n", pdbtp[record],
+ atom_seq_number, atom_name, res_name, chain_id, res_seq_number, x, y, z,
+ occupancy, b_factor);
return n;
}
-void write_pdbfile_indexed(FILE *out, const char *title,
- const t_atoms *atoms, const rvec x[],
- int ePBC, const matrix box, char chainid,
- int model_nr, int nindex, const int index[],
- gmx_conect conect, gmx_bool bTerSepChains,
- bool usePqrFormat)
+void write_pdbfile_indexed(FILE* out,
+ const char* title,
+ const t_atoms* atoms,
+ const rvec x[],
+ int ePBC,
+ const matrix box,
+ char chainid,
+ int model_nr,
+ int nindex,
+ const int index[],
+ gmx_conect conect,
+ bool usePqrFormat)
{
- gmx_conect_t *gc = static_cast<gmx_conect_t *>(conect);
- int i, ii;
- int resind, resnr;
- enum PDB_record type;
- unsigned char resic, ch;
- char altloc;
- real occup, bfac;
- gmx_bool bOccup;
- int chainnum, lastchainnum;
- gmx_residuetype_t*rt;
- const char *p_restype;
- const char *p_lastrestype;
-
- gmx_residuetype_init(&rt);
+ gmx_conect_t* gc = static_cast<gmx_conect_t*>(conect);
+ enum PDB_record type;
+ char altloc;
+ real occup, bfac;
+ gmx_bool bOccup;
+
fprintf(out, "TITLE %s\n", (title && title[0]) ? title : gmx::bromacs().c_str());
- if (box && ( (norm2(box[XX]) != 0.0f) || (norm2(box[YY]) != 0.0f) || (norm2(box[ZZ]) != 0.0f) ) )
+ if (box && ((norm2(box[XX]) != 0.0F) || (norm2(box[YY]) != 0.0F) || (norm2(box[ZZ]) != 0.0F)))
{
gmx_write_pdb_box(out, ePBC, box);
}
* otherwise set them all to one
*/
bOccup = TRUE;
- for (ii = 0; (ii < nindex) && bOccup; ii++)
+ for (int ii = 0; (ii < nindex) && bOccup; ii++)
{
- i = index[ii];
+ int i = index[ii];
bOccup = bOccup && (atoms->pdbinfo[i].occup == 0.0);
}
}
fprintf(out, "MODEL %8d\n", model_nr > 0 ? model_nr : 1);
- lastchainnum = -1;
- p_restype = nullptr;
-
- for (ii = 0; ii < nindex; ii++)
+ ResidueType rt;
+ for (int ii = 0; ii < nindex; ii++)
{
- i = index[ii];
- resind = atoms->atom[i].resind;
- chainnum = atoms->resinfo[resind].chainnum;
- p_lastrestype = p_restype;
- gmx_residuetype_get_type(rt, *atoms->resinfo[resind].name, &p_restype);
-
- /* Add a TER record if we changed chain, and if either the previous or this chain is protein/DNA/RNA. */
- if (bTerSepChains && ii > 0 && chainnum != lastchainnum)
- {
- /* Only add TER if the previous chain contained protein/DNA/RNA. */
- if (gmx_residuetype_is_protein(rt, p_lastrestype) || gmx_residuetype_is_dna(rt, p_lastrestype) || gmx_residuetype_is_rna(rt, p_lastrestype))
- {
- fprintf(out, "TER\n");
- }
- lastchainnum = chainnum;
- }
-
- std::string resnm = *atoms->resinfo[resind].name;
- std::string nm = *atoms->atomname[i];
+ int i = index[ii];
+ int resind = atoms->atom[i].resind;
+ std::string resnm = *atoms->resinfo[resind].name;
+ std::string nm = *atoms->atomname[i];
/* rename HG12 to 2HG1, etc. */
- nm = xlate_atomname_gmx2pdb(nm);
- resnr = atoms->resinfo[resind].nr;
- resic = atoms->resinfo[resind].ic;
+ nm = xlate_atomname_gmx2pdb(nm);
+ int resnr = atoms->resinfo[resind].nr;
+ unsigned char resic = atoms->resinfo[resind].ic;
+ unsigned char ch;
if (chainid != ' ')
{
ch = chainid;
bfac = pdbinfo.bfac;
if (!usePqrFormat)
{
- gmx_fprintf_pdb_atomline(out,
- type,
- i+1,
- nm.c_str(),
- altloc,
- resnm.c_str(),
- ch,
- resnr,
- resic,
- 10*x[i][XX], 10*x[i][YY], 10*x[i][ZZ],
- occup,
- bfac,
- atoms->atom[i].elem);
+ gmx_fprintf_pdb_atomline(out, type, i + 1, nm.c_str(), altloc, resnm.c_str(), ch, resnr,
+ resic, 10 * x[i][XX], 10 * x[i][YY], 10 * x[i][ZZ], occup,
+ bfac, atoms->atom[i].elem);
if (atoms->pdbinfo && atoms->pdbinfo[i].bAnisotropic)
{
- fprintf(out, "ANISOU%5d %-4.4s%4.4s%c%4d%c %7d%7d%7d%7d%7d%7d\n",
- (i+1)%100000, nm.c_str(), resnm.c_str(), ch, resnr,
- (resic == '\0') ? ' ' : resic,
- atoms->pdbinfo[i].uij[0], atoms->pdbinfo[i].uij[1],
- atoms->pdbinfo[i].uij[2], atoms->pdbinfo[i].uij[3],
- atoms->pdbinfo[i].uij[4], atoms->pdbinfo[i].uij[5]);
+ fprintf(out, "ANISOU%5d %-4.4s%4.4s%c%4d%c %7d%7d%7d%7d%7d%7d\n", (i + 1) % 100000,
+ nm.c_str(), resnm.c_str(), ch, resnr, (resic == '\0') ? ' ' : resic,
+ atoms->pdbinfo[i].uij[0], atoms->pdbinfo[i].uij[1], atoms->pdbinfo[i].uij[2],
+ atoms->pdbinfo[i].uij[3], atoms->pdbinfo[i].uij[4], atoms->pdbinfo[i].uij[5]);
}
}
else
{
- gmx_fprintf_pqr_atomline(out,
- type,
- i+1,
- nm.c_str(),
- resnm.c_str(),
- ch,
- resnr,
- 10*x[i][XX], 10*x[i][YY], 10*x[i][ZZ],
- occup,
- bfac);
+ gmx_fprintf_pqr_atomline(out, type, i + 1, nm.c_str(), resnm.c_str(), ch, resnr,
+ 10 * x[i][XX], 10 * x[i][YY], 10 * x[i][ZZ], occup, bfac);
}
}
if (nullptr != gc)
{
/* Write conect records */
- for (i = 0; (i < gc->nconect); i++)
+ for (int i = 0; (i < gc->nconect); i++)
{
- fprintf(out, "CONECT%5d%5d\n", gc->conect[i].ai+1, gc->conect[i].aj+1);
+ fprintf(out, "CONECT%5d%5d\n", gc->conect[i].ai + 1, gc->conect[i].aj + 1);
}
}
-
- gmx_residuetype_destroy(rt);
}
-void write_pdbfile(FILE *out, const char *title, const t_atoms *atoms, const rvec x[],
- int ePBC, const matrix box, char chainid, int model_nr, gmx_conect conect, gmx_bool bTerSepChains)
+void write_pdbfile(FILE* out,
+ const char* title,
+ const t_atoms* atoms,
+ const rvec x[],
+ int ePBC,
+ const matrix box,
+ char chainid,
+ int model_nr,
+ gmx_conect conect)
{
int i, *index;
{
index[i] = i;
}
- write_pdbfile_indexed(out, title, atoms, x, ePBC, box, chainid, model_nr,
- atoms->nr, index, conect, bTerSepChains, false);
+ write_pdbfile_indexed(out, title, atoms, x, ePBC, box, chainid, model_nr, atoms->nr, index,
+ conect, false);
sfree(index);
}
-static int line2type(const char *line)
+static int line2type(const char* line)
{
int k;
char type[8];
return k;
}
-static void read_anisou(char line[], int natom, t_atoms *atoms)
+static void read_anisou(char line[], int natom, t_atoms* atoms)
{
int i, j, k, atomnr;
char nc = '\0';
/* Search backwards for number and name only */
atomnr = std::strtol(anr, nullptr, 10);
- for (i = natom-1; (i >= 0); i--)
+ for (i = natom - 1; (i >= 0); i--)
{
- if ((std::strcmp(anm, *(atoms->atomname[i])) == 0) &&
- (atomnr == atoms->pdbinfo[i].atomnr))
+ if ((std::strcmp(anm, *(atoms->atomname[i])) == 0) && (atomnr == atoms->pdbinfo[i].atomnr))
{
break;
}
}
if (i < 0)
{
- fprintf(stderr, "Skipping ANISOU record (atom %s %d not found)\n",
- anm, atomnr);
+ fprintf(stderr, "Skipping ANISOU record (atom %s %d not found)\n", anm, atomnr);
}
else
{
- if (sscanf(line+29, "%d%d%d%d%d%d",
- &atoms->pdbinfo[i].uij[U11], &atoms->pdbinfo[i].uij[U22],
+ if (sscanf(line + 29, "%d%d%d%d%d%d", &atoms->pdbinfo[i].uij[U11], &atoms->pdbinfo[i].uij[U22],
&atoms->pdbinfo[i].uij[U33], &atoms->pdbinfo[i].uij[U12],
&atoms->pdbinfo[i].uij[U13], &atoms->pdbinfo[i].uij[U23])
== 6)
}
}
-void get_pdb_atomnumber(const t_atoms *atoms, gmx_atomprop_t aps)
+void get_pdb_atomnumber(const t_atoms* atoms, AtomProperties* aps)
{
int i, atomnumber, len;
size_t k;
- char anm[6], anm_copy[6], *ptr;
+ char anm[6], anm_copy[6];
char nc = '\0';
real eval;
std::strcpy(anm, atoms->pdbinfo[i].atomnm);
std::strcpy(anm_copy, atoms->pdbinfo[i].atomnm);
bool atomNumberSet = false;
- len = strlen(anm);
+ len = strlen(anm);
if ((anm[0] != ' ') && ((len <= 2) || !std::isdigit(anm[2])))
{
anm_copy[2] = nc;
- if (gmx_atomprop_query(aps, epropElement, "???", anm_copy, &eval))
+ if (aps->setAtomProperty(epropElement, "???", anm_copy, &eval))
{
atomnumber = gmx::roundToInt(eval);
atomNumberSet = true;
else
{
anm_copy[1] = nc;
- if (gmx_atomprop_query(aps, epropElement, "???", anm_copy, &eval))
+ if (aps->setAtomProperty(epropElement, "???", anm_copy, &eval))
{
atomnumber = gmx::roundToInt(eval);
atomNumberSet = true;
}
anm_copy[0] = anm[k];
anm_copy[1] = nc;
- if (gmx_atomprop_query(aps, epropElement, "???", anm_copy, &eval))
+ if (aps->setAtomProperty(epropElement, "???", anm_copy, &eval))
{
atomnumber = gmx::roundToInt(eval);
atomNumberSet = true;
}
}
+ std::string buf;
if (atomNumberSet)
{
atoms->atom[i].atomnumber = atomnumber;
- ptr = gmx_atomprop_element(aps, atomnumber);
+ buf = aps->elementFromAtomNumber(atomnumber);
if (debug)
{
- fprintf(debug, "Atomnumber for atom '%s' is %d\n",
- anm, atomnumber);
+ fprintf(debug, "Atomnumber for atom '%s' is %d\n", anm, atomnumber);
}
}
- else
- {
- ptr = nullptr;
- }
- std::strncpy(atoms->atom[i].elem, ptr == nullptr ? "" : ptr, 4);
+ buf.resize(3);
+ std::strncpy(atoms->atom[i].elem, buf.c_str(), 4);
}
}
-static int read_atom(t_symtab *symtab,
- const char line[], int type, int natom,
- t_atoms *atoms, rvec x[], int chainnum, gmx_bool bChange)
+static int
+read_atom(t_symtab* symtab, const char line[], int type, int natom, t_atoms* atoms, rvec x[], int chainnum, gmx_bool bChange)
{
- t_atom *atomn;
+ t_atom* atomn;
int j, k;
char nc = '\0';
char anr[12], anm[12], anm_copy[12], altloc, resnm[12], rnr[12], elem[3];
if (natom >= atoms->nr)
{
- gmx_fatal(FARGS, "\nFound more atoms (%d) in pdb file than expected (%d)",
- natom+1, atoms->nr);
+ gmx_fatal(FARGS, "\nFound more atoms (%d) in pdb file than expected (%d)", natom + 1, atoms->nr);
}
/* Skip over type */
trim(rnr);
resnr = std::strtol(rnr, nullptr, 10);
resic = line[j];
- j += 4;
+ j += 4;
/* X,Y,Z Coordinate */
for (k = 0; (k < 8); k++, j++)
if (atoms->atom)
{
atomn = &(atoms->atom[natom]);
- if ((natom == 0) ||
- atoms->resinfo[atoms->atom[natom-1].resind].nr != resnr ||
- atoms->resinfo[atoms->atom[natom-1].resind].ic != resic ||
- (strcmp(*atoms->resinfo[atoms->atom[natom-1].resind].name, resnm) != 0))
+ if ((natom == 0) || atoms->resinfo[atoms->atom[natom - 1].resind].nr != resnr
+ || atoms->resinfo[atoms->atom[natom - 1].resind].ic != resic
+ || (strcmp(*atoms->resinfo[atoms->atom[natom - 1].resind].name, resnm) != 0))
{
if (natom == 0)
{
}
else
{
- atomn->resind = atoms->atom[natom-1].resind + 1;
+ atomn->resind = atoms->atom[natom - 1].resind + 1;
}
atoms->nres = atomn->resind + 1;
t_atoms_set_resinfo(atoms, natom, symtab, resnm, resnr, resic, chainnum, chainid);
}
else
{
- atomn->resind = atoms->atom[natom-1].resind;
+ atomn->resind = atoms->atom[natom - 1].resind;
}
if (bChange)
{
atomn->atomnumber = atomnumber;
strncpy(atomn->elem, elem, 4);
}
- x[natom][XX] = strtod(xc, nullptr)*0.1;
- x[natom][YY] = strtod(yc, nullptr)*0.1;
- x[natom][ZZ] = strtod(zc, nullptr)*0.1;
+ x[natom][XX] = strtod(xc, nullptr) * 0.1;
+ x[natom][YY] = strtod(yc, nullptr) * 0.1;
+ x[natom][ZZ] = strtod(zc, nullptr) * 0.1;
if (atoms->pdbinfo)
{
atoms->pdbinfo[natom].type = type;
return natom;
}
-gmx_bool is_hydrogen(const char *nm)
+gmx_bool is_hydrogen(const char* nm)
{
char buf[30];
return FALSE;
}
-gmx_bool is_dummymass(const char *nm)
+gmx_bool is_dummymass(const char* nm)
{
char buf[30];
std::strcpy(buf, nm);
trim(buf);
- return (buf[0] == 'M') && (std::isdigit(buf[strlen(buf)-1]) != 0);
+ return (buf[0] == 'M') && (std::isdigit(buf[strlen(buf) - 1]) != 0);
}
-static void gmx_conect_addline(gmx_conect_t *con, char *line)
+static void gmx_conect_addline(gmx_conect_t* con, char* line)
{
- int n, ai, aj;
+ int n, ai, aj;
std::string form2 = "%%*s";
std::string format = form2 + "%%d";
n = sscanf(line, format.c_str(), &aj);
if (n == 1)
{
- srenew(con->conect, ++con->nconect);
- con->conect[con->nconect - 1].ai = ai - 1;
- con->conect[con->nconect - 1].aj = aj - 1;
+ gmx_conect_add(con, ai - 1, aj - 1); /* to prevent duplicated records */
}
- }
- while (n == 1);
+ } while (n == 1);
}
}
-void gmx_conect_dump(FILE *fp, gmx_conect conect)
+void gmx_conect_dump(FILE* fp, gmx_conect conect)
{
- gmx_conect_t *gc = static_cast<gmx_conect_t *>(conect);
+ gmx_conect_t* gc = static_cast<gmx_conect_t*>(conect);
int i;
for (i = 0; (i < gc->nconect); i++)
{
- fprintf(fp, "%6s%5d%5d\n", "CONECT",
- gc->conect[i].ai+1, gc->conect[i].aj+1);
+ fprintf(fp, "%6s%5d%5d\n", "CONECT", gc->conect[i].ai + 1, gc->conect[i].aj + 1);
}
}
gmx_conect gmx_conect_init()
{
- gmx_conect_t *gc;
+ gmx_conect_t* gc;
snew(gc, 1);
void gmx_conect_done(gmx_conect conect)
{
- gmx_conect_t *gc = conect;
+ gmx_conect_t* gc = conect;
sfree(gc->conect);
}
gmx_bool gmx_conect_exist(gmx_conect conect, int ai, int aj)
{
- gmx_conect_t *gc = conect;
+ gmx_conect_t* gc = conect;
int i;
/* if (!gc->bSorted)
for (i = 0; (i < gc->nconect); i++)
{
- if (((gc->conect[i].ai == ai) &&
- (gc->conect[i].aj == aj)) ||
- ((gc->conect[i].aj == ai) &&
- (gc->conect[i].ai == aj)))
+ if (((gc->conect[i].ai == ai) && (gc->conect[i].aj == aj))
+ || ((gc->conect[i].aj == ai) && (gc->conect[i].ai == aj)))
{
return TRUE;
}
void gmx_conect_add(gmx_conect conect, int ai, int aj)
{
- gmx_conect_t *gc = static_cast<gmx_conect_t *>(conect);
+ gmx_conect_t* gc = static_cast<gmx_conect_t*>(conect);
/* if (!gc->bSorted)
sort_conect(gc);*/
if (!gmx_conect_exist(conect, ai, aj))
{
srenew(gc->conect, ++gc->nconect);
- gc->conect[gc->nconect-1].ai = ai;
- gc->conect[gc->nconect-1].aj = aj;
+ gc->conect[gc->nconect - 1].ai = ai;
+ gc->conect[gc->nconect - 1].aj = aj;
}
}
-int read_pdbfile(FILE *in, char *title, int *model_nr,
- t_atoms *atoms, t_symtab *symtab, rvec x[], int *ePBC,
- matrix box, gmx_bool bChange, gmx_conect conect)
+int read_pdbfile(FILE* in,
+ char* title,
+ int* model_nr,
+ t_atoms* atoms,
+ t_symtab* symtab,
+ rvec x[],
+ int* ePBC,
+ matrix box,
+ gmx_bool bChange,
+ gmx_conect conect)
{
- gmx_conect_t *gc = conect;
+ gmx_conect_t* gc = conect;
gmx_bool bCOMPND;
gmx_bool bConnWarn = FALSE;
- char line[STRLEN+1];
+ char line[STRLEN + 1];
int line_type;
- char *c, *d;
+ char * c, *d;
int natom, chainnum;
- gmx_bool bStop = FALSE;
+ gmx_bool bStop = FALSE;
if (ePBC)
{
}
break;
- case epdbCRYST1:
- read_cryst1(line, ePBC, box);
- break;
+ case epdbCRYST1: read_cryst1(line, ePBC, box); break;
case epdbTITLE:
case epdbHEADER:
if (std::strlen(line) > 6)
{
- c = line+6;
+ c = line + 6;
/* skip HEADER or TITLE and spaces */
while (c[0] != ' ')
{
break;
case epdbCOMPND:
- if ((!std::strstr(line, ": ")) || (std::strstr(line+6, "MOLECULE:")))
+ if ((!std::strstr(line, ": ")) || (std::strstr(line + 6, "MOLECULE:")))
{
- if (!(c = std::strstr(line+6, "MOLECULE:")) )
+ if (!(c = std::strstr(line + 6, "MOLECULE:")))
{
c = line;
}
d = strstr(c, " ");
if (d)
{
- while ( (d[-1] == ';') && d > c)
+ while ((d[-1] == ';') && d > c)
{
d--;
}
}
break;
- case epdbTER:
- chainnum++;
- break;
+ case epdbTER: chainnum++; break;
case epdbMODEL:
if (model_nr)
}
break;
- case epdbENDMDL:
- bStop = TRUE;
- break;
+ case epdbENDMDL: bStop = TRUE; break;
case epdbCONECT:
if (gc)
{
}
break;
- default:
- break;
+ default: break;
}
}
return natom;
}
-void get_pdb_coordnum(FILE *in, int *natoms)
+void get_pdb_coordnum(FILE* in, int* natoms)
{
char line[STRLEN];
}
}
-void gmx_pdb_read_conf(const char *infile,
- t_symtab *symtab, char **name, t_atoms *atoms,
- rvec x[], int *ePBC, matrix box)
+void gmx_pdb_read_conf(const char* infile, t_symtab* symtab, char** name, t_atoms* atoms, rvec x[], int* ePBC, matrix box)
{
- FILE *in = gmx_fio_fopen(infile, "r");
+ FILE* in = gmx_fio_fopen(infile, "r");
char title[STRLEN];
read_pdbfile(in, title, nullptr, atoms, symtab, x, ePBC, box, TRUE, nullptr);
if (name != nullptr)
gmx_fio_fclose(in);
}
-gmx_conect gmx_conect_generate(const t_topology *top)
+gmx_conect gmx_conect_generate(const t_topology* top)
{
int f, i;
gmx_conect gc;
/* Fill the conect records */
- gc = gmx_conect_init();
+ gc = gmx_conect_init();
for (f = 0; (f < F_NRE); f++)
{
if (IS_CHEMBOND(f))
{
- for (i = 0; (i < top->idef.il[f].nr); i += interaction_function[f].nratoms+1)
+ for (i = 0; (i < top->idef.il[f].nr); i += interaction_function[f].nratoms + 1)
{
- gmx_conect_add(gc, top->idef.il[f].iatoms[i+1],
- top->idef.il[f].iatoms[i+2]);
+ gmx_conect_add(gc, top->idef.il[f].iatoms[i + 1], top->idef.il[f].iatoms[i + 2]);
}
}
}
return gc;
}
-int
-gmx_fprintf_pdb_atomline(FILE * fp,
- enum PDB_record record,
- int atom_seq_number,
- const char * atom_name,
- char alternate_location,
- const char * res_name,
- char chain_id,
- int res_seq_number,
- char res_insertion_code,
- real x,
- real y,
- real z,
- real occupancy,
- real b_factor,
- const char * element)
+int gmx_fprintf_pdb_atomline(FILE* fp,
+ enum PDB_record record,
+ int atom_seq_number,
+ const char* atom_name,
+ char alternate_location,
+ const char* res_name,
+ char chain_id,
+ int res_seq_number,
+ char res_insertion_code,
+ real x,
+ real y,
+ real z,
+ real occupancy,
+ real b_factor,
+ const char* element)
{
char tmp_atomname[6], tmp_resname[6];
gmx_bool start_name_in_col13;
/* If the atom name is an element name with two chars, it should start already in column 13.
* Otherwise it should start in column 14, unless the name length is 4 chars.
*/
- if ( (element != nullptr) && (std::strlen(element) >= 2) && (gmx_strncasecmp(atom_name, element, 2) == 0) )
+ if ((element != nullptr) && (std::strlen(element) >= 2)
+ && (gmx_strncasecmp(atom_name, element, 2) == 0))
{
start_name_in_col13 = TRUE;
}
atom_seq_number = atom_seq_number % 100000;
res_seq_number = res_seq_number % 10000;
- n = fprintf(fp,
- "%-6s%5d %-4.4s%c%4.4s%c%4d%c %8.3f%8.3f%8.3f%6.2f%6.2f %2s\n",
- pdbtp[record],
- atom_seq_number,
- tmp_atomname,
- alternate_location,
- tmp_resname,
- chain_id,
- res_seq_number,
- res_insertion_code,
- x, y, z,
- occupancy,
- b_factor,
+ n = fprintf(fp, "%-6s%5d %-4.4s%c%4.4s%c%4d%c %8.3f%8.3f%8.3f%6.2f%6.2f %2s\n",
+ pdbtp[record], atom_seq_number, tmp_atomname, alternate_location, tmp_resname,
+ chain_id, res_seq_number, res_insertion_code, x, y, z, occupancy, b_factor,
(element != nullptr) ? element : "");
return n;
#include <cstdlib>
#include <algorithm>
+#include <memory>
#include "gromacs/commandline/filenm.h"
-#include "gromacs/compat/make_unique.h"
#include "gromacs/domdec/domdec_struct.h"
#include "gromacs/domdec/localatomset.h"
#include "gromacs/domdec/localatomsetmanager.h"
#include "gromacs/math/vec.h"
#include "gromacs/math/vectypes.h"
#include "gromacs/mdlib/gmx_omp_nthreads.h"
-#include "gromacs/mdlib/mdrun.h"
#include "gromacs/mdtypes/commrec.h"
#include "gromacs/mdtypes/forceoutput.h"
#include "gromacs/mdtypes/inputrec.h"
#include "pull_internal.h"
-static int groupPbcFromParams(const t_pull_group ¶ms, bool setPbcRefToPrevStepCOM)
+namespace gmx
+{
+extern template LocalAtomSet LocalAtomSetManager::add<void, void>(ArrayRef<const int> globalAtomIndex);
+} // namespace gmx
+
+static int groupPbcFromParams(const t_pull_group& params, bool setPbcRefToPrevStepCOM)
{
if (params.nat <= 1)
{
* should not end before that of this object.
* This will be fixed when the pointers are replacted by std::vector.
*/
-pull_group_work_t::pull_group_work_t(const t_pull_group ¶ms,
+pull_group_work_t::pull_group_work_t(const t_pull_group& params,
gmx::LocalAtomSet atomSet,
bool bSetPbcRefToPrevStepCOM) :
params(params),
clear_dvec(xp);
};
-bool pull_coordinate_is_angletype(const t_pull_coord *pcrd)
+bool pull_coordinate_is_angletype(const t_pull_coord* pcrd)
{
- return (pcrd->eGeom == epullgANGLE ||
- pcrd->eGeom == epullgDIHEDRAL ||
- pcrd->eGeom == epullgANGLEAXIS);
+ return (pcrd->eGeom == epullgANGLE || pcrd->eGeom == epullgDIHEDRAL || pcrd->eGeom == epullgANGLEAXIS);
}
-static bool pull_coordinate_is_directional(const t_pull_coord *pcrd)
+static bool pull_coordinate_is_directional(const t_pull_coord* pcrd)
{
- return (pcrd->eGeom == epullgDIR ||
- pcrd->eGeom == epullgDIRPBC ||
- pcrd->eGeom == epullgDIRRELATIVE ||
- pcrd->eGeom == epullgCYL);
+ return (pcrd->eGeom == epullgDIR || pcrd->eGeom == epullgDIRPBC
+ || pcrd->eGeom == epullgDIRRELATIVE || pcrd->eGeom == epullgCYL);
}
-const char *pull_coordinate_units(const t_pull_coord *pcrd)
+const char* pull_coordinate_units(const t_pull_coord* pcrd)
{
- return pull_coordinate_is_angletype(pcrd) ? "deg" : "nm";
+ return pull_coordinate_is_angletype(pcrd) ? "deg" : "nm";
}
-double pull_conversion_factor_userinput2internal(const t_pull_coord *pcrd)
+double pull_conversion_factor_userinput2internal(const t_pull_coord* pcrd)
{
if (pull_coordinate_is_angletype(pcrd))
{
}
}
-double pull_conversion_factor_internal2userinput(const t_pull_coord *pcrd)
+double pull_conversion_factor_internal2userinput(const t_pull_coord* pcrd)
{
if (pull_coordinate_is_angletype(pcrd))
{
}
/* Apply forces in a mass weighted fashion for part of the pull group */
-static void apply_forces_grp_part(const pull_group_work_t *pgrp,
- int ind_start, int ind_end,
- const t_mdatoms *md,
- const dvec f_pull, int sign, rvec *f)
+static void apply_forces_grp_part(const pull_group_work_t* pgrp,
+ int ind_start,
+ int ind_end,
+ const t_mdatoms* md,
+ const dvec f_pull,
+ int sign,
+ rvec* f)
{
double inv_wm = pgrp->mwscale;
- auto localAtomIndices = pgrp->atomSet.localIndex();
+ auto localAtomIndices = pgrp->atomSet.localIndex();
for (int i = ind_start; i < ind_end; i++)
{
int ii = localAtomIndices[i];
}
/* Apply forces in a mass weighted fashion */
-static void apply_forces_grp(const pull_group_work_t *pgrp,
- const t_mdatoms *md,
- const dvec f_pull, int sign, rvec *f,
- int nthreads)
+static void apply_forces_grp(const pull_group_work_t* pgrp,
+ const t_mdatoms* md,
+ const dvec f_pull,
+ int sign,
+ rvec* f,
+ int nthreads)
{
auto localAtomIndices = pgrp->atomSet.localIndex();
*/
for (int d = 0; d < DIM; d++)
{
- f[localAtomIndices[0]][d] += sign*f_pull[d];
+ f[localAtomIndices[0]][d] += sign * f_pull[d];
}
}
else
#pragma omp parallel for num_threads(nthreads) schedule(static)
for (int th = 0; th < nthreads; th++)
{
- int ind_start = (localAtomIndices.size()*(th + 0))/nthreads;
- int ind_end = (localAtomIndices.size()*(th + 1))/nthreads;
- apply_forces_grp_part(pgrp, ind_start, ind_end,
- md, f_pull, sign, f);
+ int ind_start = (localAtomIndices.size() * (th + 0)) / nthreads;
+ int ind_end = (localAtomIndices.size() * (th + 1)) / nthreads;
+ apply_forces_grp_part(pgrp, ind_start, ind_end, md, f_pull, sign, f);
}
}
}
}
/* Apply forces in a mass weighted fashion to a cylinder group */
-static void apply_forces_cyl_grp(const pull_group_work_t *pgrp,
- const double dv_corr,
- const t_mdatoms *md,
- const dvec f_pull, double f_scal,
- int sign, rvec *f,
+static void apply_forces_cyl_grp(const pull_group_work_t* pgrp,
+ const double dv_corr,
+ const t_mdatoms* md,
+ const dvec f_pull,
+ double f_scal,
+ int sign,
+ rvec* f,
int gmx_unused nthreads)
{
double inv_wm = pgrp->mwscale;
- auto localAtomIndices = pgrp->atomSet.localIndex();
+ auto localAtomIndices = pgrp->atomSet.localIndex();
/* The cylinder group is always a slab in the system, thus large.
* Therefore we always thread-parallelize this group.
{
continue;
}
- int ii = localAtomIndices[i];
- double mass = md->massT[ii];
+ int ii = localAtomIndices[i];
+ double mass = md->massT[ii];
/* The stored axial distance from the cylinder center (dv) needs
* to be corrected for an offset (dv_corr), which was unknown when
* we calculated dv.
*/
for (int m = 0; m < DIM; m++)
{
- f[ii][m] += sign*inv_wm*(mass*weight*f_pull[m] +
- pgrp->mdw[i][m]*dv_com*f_scal);
+ f[ii][m] += sign * inv_wm * (mass * weight * f_pull[m] + pgrp->mdw[i][m] * dv_com * f_scal);
}
}
}
/* Apply torque forces in a mass weighted fashion to the groups that define
* the pull vector direction for pull coordinate pcrd.
*/
-static void apply_forces_vec_torque(const struct pull_t *pull,
- const pull_coord_work_t *pcrd,
- const t_mdatoms *md,
- rvec *f)
+static void apply_forces_vec_torque(const struct pull_t* pull,
+ const pull_coord_work_t* pcrd,
+ const t_mdatoms* md,
+ rvec* f)
{
- const PullCoordSpatialData &spatialData = pcrd->spatialData;
+ const PullCoordSpatialData& spatialData = pcrd->spatialData;
/* The component inpr along the pull vector is accounted for in the usual
* way. Here we account for the component perpendicular to vec.
double inpr = 0;
for (int m = 0; m < DIM; m++)
{
- inpr += spatialData.dr01[m]*spatialData.vec[m];
+ inpr += spatialData.dr01[m] * spatialData.vec[m];
}
/* The torque force works along the component of the distance vector
* of between the two "usual" pull groups that is perpendicular to
dvec f_perp;
for (int m = 0; m < DIM; m++)
{
- f_perp[m] = (spatialData.dr01[m] - inpr*spatialData.vec[m])/spatialData.vec_len*pcrd->scalarForce;
+ f_perp[m] = (spatialData.dr01[m] - inpr * spatialData.vec[m]) / spatialData.vec_len * pcrd->scalarForce;
}
/* Apply the force to the groups defining the vector using opposite signs */
- apply_forces_grp(&pull->group[pcrd->params.group[2]], md,
- f_perp, -1, f, pull->nthreads);
- apply_forces_grp(&pull->group[pcrd->params.group[3]], md,
- f_perp, 1, f, pull->nthreads);
+ apply_forces_grp(&pull->group[pcrd->params.group[2]], md, f_perp, -1, f, pull->nthreads);
+ apply_forces_grp(&pull->group[pcrd->params.group[3]], md, f_perp, 1, f, pull->nthreads);
}
/* Apply forces in a mass weighted fashion */
-static void apply_forces_coord(struct pull_t * pull, int coord,
- const PullCoordVectorForces &forces,
- const t_mdatoms * md,
- rvec *f)
+static void apply_forces_coord(struct pull_t* pull,
+ int coord,
+ const PullCoordVectorForces& forces,
+ const t_mdatoms* md,
+ rvec* f)
{
/* Here it would be more efficient to use one large thread-parallel
* region instead of potential parallel regions within apply_forces_grp.
* to data races.
*/
- const pull_coord_work_t &pcrd = pull->coord[coord];
+ const pull_coord_work_t& pcrd = pull->coord[coord];
if (pcrd.params.eGeom == epullgCYL)
{
- apply_forces_cyl_grp(&pull->dyna[coord], pcrd.spatialData.cyl_dev, md,
- forces.force01, pcrd.scalarForce, -1, f,
- pull->nthreads);
+ apply_forces_cyl_grp(&pull->dyna[coord], pcrd.spatialData.cyl_dev, md, forces.force01,
+ pcrd.scalarForce, -1, f, pull->nthreads);
/* Sum the force along the vector and the radial force */
dvec f_tot;
for (int m = 0; m < DIM; m++)
{
- f_tot[m] = forces.force01[m] + pcrd.scalarForce*pcrd.spatialData.ffrad[m];
+ f_tot[m] = forces.force01[m] + pcrd.scalarForce * pcrd.spatialData.ffrad[m];
}
- apply_forces_grp(&pull->group[pcrd.params.group[1]], md,
- f_tot, 1, f, pull->nthreads);
+ apply_forces_grp(&pull->group[pcrd.params.group[1]], md, f_tot, 1, f, pull->nthreads);
}
else
{
if (pull->group[pcrd.params.group[0]].params.nat > 0)
{
- apply_forces_grp(&pull->group[pcrd.params.group[0]], md,
- forces.force01, -1, f, pull->nthreads);
+ apply_forces_grp(&pull->group[pcrd.params.group[0]], md, forces.force01, -1, f, pull->nthreads);
}
- apply_forces_grp(&pull->group[pcrd.params.group[1]], md,
- forces.force01, 1, f, pull->nthreads);
+ apply_forces_grp(&pull->group[pcrd.params.group[1]], md, forces.force01, 1, f, pull->nthreads);
if (pcrd.params.ngroup >= 4)
{
- apply_forces_grp(&pull->group[pcrd.params.group[2]], md,
- forces.force23, -1, f, pull->nthreads);
- apply_forces_grp(&pull->group[pcrd.params.group[3]], md,
- forces.force23, 1, f, pull->nthreads);
+ apply_forces_grp(&pull->group[pcrd.params.group[2]], md, forces.force23, -1, f, pull->nthreads);
+ apply_forces_grp(&pull->group[pcrd.params.group[3]], md, forces.force23, 1, f, pull->nthreads);
}
if (pcrd.params.ngroup >= 6)
{
- apply_forces_grp(&pull->group[pcrd.params.group[4]], md,
- forces.force45, -1, f, pull->nthreads);
- apply_forces_grp(&pull->group[pcrd.params.group[5]], md,
- forces.force45, 1, f, pull->nthreads);
+ apply_forces_grp(&pull->group[pcrd.params.group[4]], md, forces.force45, -1, f, pull->nthreads);
+ apply_forces_grp(&pull->group[pcrd.params.group[5]], md, forces.force45, 1, f, pull->nthreads);
}
}
}
-real max_pull_distance2(const pull_coord_work_t *pcrd,
- const t_pbc *pbc)
+real max_pull_distance2(const pull_coord_work_t* pcrd, const t_pbc* pbc)
{
/* Note that this maximum distance calculation is more complex than
* most other cases in GROMACS, since here we have to take care of
}
}
- return 0.25*max_d2;
+ return 0.25 * max_d2;
}
/* This function returns the distance based on coordinates xg and xref.
* Note that the pull coordinate struct pcrd is not modified.
*/
-static void low_get_pull_coord_dr(const struct pull_t *pull,
- const pull_coord_work_t *pcrd,
- const t_pbc *pbc,
- dvec xg, dvec xref, double max_dist2,
- dvec dr)
+static void low_get_pull_coord_dr(const struct pull_t* pull,
+ const pull_coord_work_t* pcrd,
+ const t_pbc* pbc,
+ dvec xg,
+ dvec xref,
+ double max_dist2,
+ dvec dr)
{
- const pull_group_work_t *pgrp0 = &pull->group[pcrd->params.group[0]];
+ const pull_group_work_t* pgrp0 = &pull->group[pcrd->params.group[0]];
/* Only the first group can be an absolute reference, in that case nat=0 */
if (pgrp0->params.nat == 0)
dvec xrefr;
copy_dvec(xref, xrefr);
- dvec dref = {0, 0, 0};
+ dvec dref = { 0, 0, 0 };
if (pcrd->params.eGeom == epullgDIRPBC)
{
for (int m = 0; m < DIM; m++)
{
- dref[m] = pcrd->value_ref*pcrd->spatialData.vec[m];
+ dref[m] = pcrd->value_ref * pcrd->spatialData.vec[m];
}
/* Add the reference position, so we use the correct periodic image */
dvec_inc(xrefr, dref);
dr[m] *= pcrd->params.dim[m];
if (pcrd->params.dim[m] && !(directional && pcrd->spatialData.vec[m] == 0))
{
- dr2 += dr[m]*dr[m];
+ dr2 += dr[m] * dr[m];
}
}
/* Check if we are close to switching to another periodic image */
- if (max_dist2 > 0 && dr2 > 0.98*0.98*max_dist2)
+ if (max_dist2 > 0 && dr2 > 0.98 * 0.98 * max_dist2)
{
/* Note that technically there is no issue with switching periodic
* image, as pbc_dx_d returns the distance to the closest periodic
* image. However in all cases where periodic image switches occur,
* the pull results will be useless in practice.
*/
- gmx_fatal(FARGS, "Distance between pull groups %d and %d (%f nm) is larger than 0.49 times the box size (%f).\n%s",
- pcrd->params.group[0], pcrd->params.group[1],
- sqrt(dr2), sqrt(0.98*0.98*max_dist2),
- pcrd->params.eGeom == epullgDIR ? "You might want to consider using \"pull-geometry = direction-periodic\" instead.\n" : "");
+ gmx_fatal(FARGS,
+ "Distance between pull groups %d and %d (%f nm) is larger than 0.49 times the "
+ "box size (%f).\n%s",
+ pcrd->params.group[0], pcrd->params.group[1], sqrt(dr2),
+ sqrt(0.98 * 0.98 * max_dist2), pcrd->params.eGeom == epullgDIR ? "You might want to consider using \"pull-geometry = direction-periodic\" instead.\n" : "");
}
if (pcrd->params.eGeom == epullgDIRPBC)
/* This function returns the distance based on the contents of the pull struct.
* pull->coord[coord_ind].dr, and potentially vector, are updated.
*/
-static void get_pull_coord_dr(struct pull_t *pull,
- int coord_ind,
- const t_pbc *pbc)
+static void get_pull_coord_dr(struct pull_t* pull, int coord_ind, const t_pbc* pbc)
{
- pull_coord_work_t *pcrd = &pull->coord[coord_ind];
- PullCoordSpatialData &spatialData = pcrd->spatialData;
+ pull_coord_work_t* pcrd = &pull->coord[coord_ind];
+ PullCoordSpatialData& spatialData = pcrd->spatialData;
- double md2;
+ double md2;
- if (pcrd->params.eGeom == epullgDIRPBC)
+ /* With AWH pulling we allow for periodic pulling with geometry=direction.
+ * TODO: Store a periodicity flag instead of checking for external pull provider.
+ */
- if (pcrd->params.eGeom == epullgDIRPBC || (pcrd->params.eGeom == epullgDIR &&
- pcrd->params.eType == epullEXTERNAL))
++ if (pcrd->params.eGeom == epullgDIRPBC
++ || (pcrd->params.eGeom == epullgDIR && pcrd->params.eType == epullEXTERNAL))
{
md2 = -1;
}
if (pcrd->params.eGeom == epullgDIRRELATIVE)
{
/* We need to determine the pull vector */
- dvec vec;
- int m;
+ dvec vec;
+ int m;
- const pull_group_work_t &pgrp2 = pull->group[pcrd->params.group[2]];
- const pull_group_work_t &pgrp3 = pull->group[pcrd->params.group[3]];
+ const pull_group_work_t& pgrp2 = pull->group[pcrd->params.group[2]];
+ const pull_group_work_t& pgrp3 = pull->group[pcrd->params.group[3]];
pbc_dx_d(pbc, pgrp3.x, pgrp2.x, vec);
spatialData.vec_len = dnorm(vec);
for (m = 0; m < DIM; m++)
{
- spatialData.vec[m] = vec[m]/spatialData.vec_len;
+ spatialData.vec[m] = vec[m] / spatialData.vec_len;
}
if (debug)
{
fprintf(debug, "pull coord %d vector: %6.3f %6.3f %6.3f normalized: %6.3f %6.3f %6.3f\n",
- coord_ind,
- vec[XX], vec[YY], vec[ZZ],
- spatialData.vec[XX], spatialData.vec[YY], spatialData.vec[ZZ]);
+ coord_ind, vec[XX], vec[YY], vec[ZZ], spatialData.vec[XX], spatialData.vec[YY],
+ spatialData.vec[ZZ]);
}
}
- pull_group_work_t *pgrp0 = &pull->group[pcrd->params.group[0]];
- pull_group_work_t *pgrp1 = &pull->group[pcrd->params.group[1]];
+ pull_group_work_t* pgrp0 = &pull->group[pcrd->params.group[0]];
+ pull_group_work_t* pgrp1 = &pull->group[pcrd->params.group[1]];
- low_get_pull_coord_dr(pull, pcrd, pbc,
- pgrp1->x,
- pcrd->params.eGeom == epullgCYL ? pull->dyna[coord_ind].x : pgrp0->x,
- md2,
+ low_get_pull_coord_dr(pull, pcrd, pbc, pgrp1->x,
+ pcrd->params.eGeom == epullgCYL ? pull->dyna[coord_ind].x : pgrp0->x, md2,
spatialData.dr01);
if (pcrd->params.ngroup >= 4)
pgrp2 = &pull->group[pcrd->params.group[2]];
pgrp3 = &pull->group[pcrd->params.group[3]];
- low_get_pull_coord_dr(pull, pcrd, pbc,
- pgrp3->x,
- pgrp2->x,
- md2,
- spatialData.dr23);
+ low_get_pull_coord_dr(pull, pcrd, pbc, pgrp3->x, pgrp2->x, md2, spatialData.dr23);
}
if (pcrd->params.ngroup >= 6)
{
pgrp4 = &pull->group[pcrd->params.group[4]];
pgrp5 = &pull->group[pcrd->params.group[5]];
- low_get_pull_coord_dr(pull, pcrd, pbc,
- pgrp5->x,
- pgrp4->x,
- md2,
- spatialData.dr45);
+ low_get_pull_coord_dr(pull, pcrd, pbc, pgrp5->x, pgrp4->x, md2, spatialData.dr45);
}
}
* It is assumed that x is in [-3pi, 3pi) so that x
* needs to be shifted by at most one period.
*/
-static void make_periodic_2pi(double *x)
+static void make_periodic_2pi(double* x)
{
if (*x >= M_PI)
{
}
/* This function should always be used to modify pcrd->value_ref */
-static void low_set_pull_coord_reference_value(pull_coord_work_t *pcrd,
- int coord_ind,
- double value_ref)
+static void low_set_pull_coord_reference_value(pull_coord_work_t* pcrd, int coord_ind, double value_ref)
{
- GMX_ASSERT(pcrd->params.eType != epullEXTERNAL, "The pull coord reference value should not be used with type external-potential");
+ GMX_ASSERT(pcrd->params.eType != epullEXTERNAL,
+ "The pull coord reference value should not be used with type external-potential");
if (pcrd->params.eGeom == epullgDIST)
{
if (value_ref < 0)
{
- gmx_fatal(FARGS, "Pull reference distance for coordinate %d (%f) needs to be non-negative", coord_ind + 1, value_ref);
+ gmx_fatal(FARGS,
+ "Pull reference distance for coordinate %d (%f) needs to be non-negative",
+ coord_ind + 1, value_ref);
}
}
else if (pcrd->params.eGeom == epullgANGLE || pcrd->params.eGeom == epullgANGLEAXIS)
{
if (value_ref < 0 || value_ref > M_PI)
{
- gmx_fatal(FARGS, "Pull reference angle for coordinate %d (%f) needs to be in the allowed interval [0,180] deg", coord_ind + 1, value_ref*pull_conversion_factor_internal2userinput(&pcrd->params));
+ gmx_fatal(FARGS,
+ "Pull reference angle for coordinate %d (%f) needs to be in the allowed "
+ "interval [0,180] deg",
+ coord_ind + 1,
+ value_ref * pull_conversion_factor_internal2userinput(&pcrd->params));
}
}
else if (pcrd->params.eGeom == epullgDIHEDRAL)
pcrd->value_ref = value_ref;
}
-static void update_pull_coord_reference_value(pull_coord_work_t *pcrd, int coord_ind, double t)
+static void update_pull_coord_reference_value(pull_coord_work_t* pcrd, int coord_ind, double t)
{
/* With zero rate the reference value is set initially and doesn't change */
if (pcrd->params.rate != 0)
{
- double value_ref = (pcrd->params.init + pcrd->params.rate*t)*pull_conversion_factor_userinput2internal(&pcrd->params);
+ double value_ref = (pcrd->params.init + pcrd->params.rate * t)
+ * pull_conversion_factor_userinput2internal(&pcrd->params);
low_set_pull_coord_reference_value(pcrd, coord_ind, value_ref);
}
}
/* Returns the dihedral angle. Updates the plane normal vectors m, n. */
-static double get_dihedral_angle_coord(PullCoordSpatialData *spatialData)
+static double get_dihedral_angle_coord(PullCoordSpatialData* spatialData)
{
double phi, sign;
dvec dr32; /* store instead of dr23? */
dsvmul(-1, spatialData->dr23, dr32);
- dcprod(spatialData->dr01, dr32, spatialData->planevec_m); /* Normal of first plane */
- dcprod(dr32, spatialData->dr45, spatialData->planevec_n); /* Normal of second plane */
+ dcprod(spatialData->dr01, dr32, spatialData->planevec_m); /* Normal of first plane */
+ dcprod(dr32, spatialData->dr45, spatialData->planevec_n); /* Normal of second plane */
phi = gmx_angle_between_dvecs(spatialData->planevec_m, spatialData->planevec_n);
/* Note 1: the sign below is opposite of that in the bondeds or Bekker 1994
* we get a positive sign below. Thus, the range of the dihedral angle is (-180, 180].
*/
sign = (diprod(spatialData->dr01, spatialData->planevec_n) < 0.0) ? 1.0 : -1.0;
- return sign*phi;
+ return sign * phi;
}
/* Calculates pull->coord[coord_ind].value.
* This function also updates pull->coord[coord_ind].dr.
*/
-static void get_pull_coord_distance(struct pull_t *pull,
- int coord_ind,
- const t_pbc *pbc)
+static void get_pull_coord_distance(struct pull_t* pull, int coord_ind, const t_pbc* pbc)
{
- pull_coord_work_t *pcrd;
+ pull_coord_work_t* pcrd;
int m;
pcrd = &pull->coord[coord_ind];
get_pull_coord_dr(pull, coord_ind, pbc);
- PullCoordSpatialData &spatialData = pcrd->spatialData;
+ PullCoordSpatialData& spatialData = pcrd->spatialData;
switch (pcrd->params.eGeom)
{
spatialData.value = 0;
for (m = 0; m < DIM; m++)
{
- spatialData.value += spatialData.vec[m]*spatialData.dr01[m];
+ spatialData.value += spatialData.vec[m] * spatialData.dr01[m];
}
break;
case epullgANGLE:
spatialData.value = gmx_angle_between_dvecs(spatialData.dr01, spatialData.dr23);
break;
- case epullgDIHEDRAL:
- spatialData.value = get_dihedral_angle_coord(&spatialData);
- break;
+ case epullgDIHEDRAL: spatialData.value = get_dihedral_angle_coord(&spatialData); break;
case epullgANGLEAXIS:
spatialData.value = gmx_angle_between_dvecs(spatialData.dr01, spatialData.vec);
break;
- default:
- gmx_incons("Unsupported pull type in get_pull_coord_distance");
+ default: gmx_incons("Unsupported pull type in get_pull_coord_distance");
}
}
/* Returns the deviation from the reference value.
* Updates pull->coord[coord_ind].dr, .value and .value_ref.
*/
-static double get_pull_coord_deviation(struct pull_t *pull,
- int coord_ind,
- const t_pbc *pbc,
- double t)
+static double get_pull_coord_deviation(struct pull_t* pull, int coord_ind, const t_pbc* pbc, double t)
{
- pull_coord_work_t *pcrd;
+ pull_coord_work_t* pcrd;
double dev = 0;
pcrd = &pull->coord[coord_ind];
return dev;
}
-double get_pull_coord_value(struct pull_t *pull,
- int coord_ind,
- const t_pbc *pbc)
+double get_pull_coord_value(struct pull_t* pull, int coord_ind, const t_pbc* pbc)
{
get_pull_coord_distance(pull, coord_ind, pbc);
return pull->coord[coord_ind].spatialData.value;
}
-void clear_pull_forces(pull_t *pull)
+void clear_pull_forces(pull_t* pull)
{
/* Zeroing the forces is only required for constraint pulling.
* It can happen that multiple constraint steps need to be applied
* and therefore the constraint forces need to be accumulated.
*/
- for (pull_coord_work_t &coord : pull->coord)
+ for (pull_coord_work_t& coord : pull->coord)
{
coord.scalarForce = 0;
}
}
/* Apply constraint using SHAKE */
-static void do_constraint(struct pull_t *pull, t_pbc *pbc,
- rvec *x, rvec *v,
- gmx_bool bMaster, tensor vir,
- double dt, double t)
+static void
+do_constraint(struct pull_t* pull, t_pbc* pbc, rvec* x, rvec* v, gmx_bool bMaster, tensor vir, double dt, double t)
{
- dvec *r_ij; /* x[i] com of i in prev. step. Obeys constr. -> r_ij[i] */
- dvec unc_ij; /* xp[i] com of i this step, before constr. -> unc_ij */
- dvec *rnew; /* current 'new' positions of the groups */
- double *dr_tot; /* the total update of the coords */
- dvec vec;
- double inpr;
- double lambda, rm, invdt = 0;
- gmx_bool bConverged_all, bConverged = FALSE;
- int niter = 0, ii, j, m, max_iter = 100;
- double a;
- dvec tmp, tmp3;
-
- snew(r_ij, pull->coord.size());
+ dvec* r_ij; /* x[i] com of i in prev. step. Obeys constr. -> r_ij[i] */
+ dvec unc_ij; /* xp[i] com of i this step, before constr. -> unc_ij */
+ dvec* rnew; /* current 'new' positions of the groups */
+ double* dr_tot; /* the total update of the coords */
+ dvec vec;
+ double inpr;
+ double lambda, rm, invdt = 0;
+ gmx_bool bConverged_all, bConverged = FALSE;
+ int niter = 0, ii, j, m, max_iter = 100;
+ double a;
+ dvec tmp, tmp3;
+
+ snew(r_ij, pull->coord.size());
snew(dr_tot, pull->coord.size());
snew(rnew, pull->group.size());
/* Determine the constraint directions from the old positions */
for (size_t c = 0; c < pull->coord.size(); c++)
{
- pull_coord_work_t *pcrd;
+ pull_coord_work_t* pcrd;
pcrd = &pull->coord[c];
*/
get_pull_coord_distance(pull, c, pbc);
- const PullCoordSpatialData &spatialData = pcrd->spatialData;
+ const PullCoordSpatialData& spatialData = pcrd->spatialData;
if (debug)
{
- fprintf(debug, "Pull coord %zu dr %f %f %f\n",
- c, spatialData.dr01[XX], spatialData.dr01[YY], spatialData.dr01[ZZ]);
+ fprintf(debug, "Pull coord %zu dr %f %f %f\n", c, spatialData.dr01[XX],
+ spatialData.dr01[YY], spatialData.dr01[ZZ]);
}
- if (pcrd->params.eGeom == epullgDIR ||
- pcrd->params.eGeom == epullgDIRPBC)
+ if (pcrd->params.eGeom == epullgDIR || pcrd->params.eGeom == epullgDIRPBC)
{
/* Select the component along vec */
a = 0;
for (m = 0; m < DIM; m++)
{
- a += spatialData.vec[m]*spatialData.dr01[m];
+ a += spatialData.vec[m] * spatialData.dr01[m];
}
for (m = 0; m < DIM; m++)
{
- r_ij[c][m] = a*spatialData.vec[m];
+ r_ij[c][m] = a * spatialData.vec[m];
}
}
else
if (dnorm2(r_ij[c]) == 0)
{
- gmx_fatal(FARGS, "Distance for pull coordinate %zu is zero with constraint pulling, which is not allowed.", c + 1);
+ gmx_fatal(FARGS,
+ "Distance for pull coordinate %zu is zero with constraint pulling, which is "
+ "not allowed.",
+ c + 1);
}
}
/* loop over all constraints */
for (size_t c = 0; c < pull->coord.size(); c++)
{
- pull_coord_work_t *pcrd;
+ pull_coord_work_t* pcrd;
pull_group_work_t *pgrp0, *pgrp1;
dvec dr0, dr1;
- pcrd = &pull->coord[c];
+ pcrd = &pull->coord[c];
if (pcrd->params.eType != epullCONSTRAINT)
{
pgrp1 = &pull->group[pcrd->params.group[1]];
/* Get the current difference vector */
- low_get_pull_coord_dr(pull, pcrd, pbc,
- rnew[pcrd->params.group[1]],
- rnew[pcrd->params.group[0]],
- -1, unc_ij);
+ low_get_pull_coord_dr(pull, pcrd, pbc, rnew[pcrd->params.group[1]],
+ rnew[pcrd->params.group[0]], -1, unc_ij);
if (debug)
{
fprintf(debug, "Pull coord %zu, iteration %d\n", c, niter);
}
- rm = 1.0/(pgrp0->invtm + pgrp1->invtm);
+ rm = 1.0 / (pgrp0->invtm + pgrp1->invtm);
switch (pcrd->params.eGeom)
{
case epullgDIST:
if (pcrd->value_ref <= 0)
{
- gmx_fatal(FARGS, "The pull constraint reference distance for group %zu is <= 0 (%f)", c, pcrd->value_ref);
+ gmx_fatal(
+ FARGS,
+ "The pull constraint reference distance for group %zu is <= 0 (%f)",
+ c, pcrd->value_ref);
}
{
double q, c_a, c_b, c_c;
c_a = diprod(r_ij[c], r_ij[c]);
- c_b = diprod(unc_ij, r_ij[c])*2;
+ c_b = diprod(unc_ij, r_ij[c]) * 2;
c_c = diprod(unc_ij, unc_ij) - gmx::square(pcrd->value_ref);
if (c_b < 0)
{
- q = -0.5*(c_b - std::sqrt(c_b*c_b - 4*c_a*c_c));
- lambda = -q/c_a;
+ q = -0.5 * (c_b - std::sqrt(c_b * c_b - 4 * c_a * c_c));
+ lambda = -q / c_a;
}
else
{
- q = -0.5*(c_b + std::sqrt(c_b*c_b - 4*c_a*c_c));
- lambda = -c_c/q;
+ q = -0.5 * (c_b + std::sqrt(c_b * c_b - 4 * c_a * c_c));
+ lambda = -c_c / q;
}
if (debug)
{
- fprintf(debug,
- "Pull ax^2+bx+c=0: a=%e b=%e c=%e lambda=%e\n",
- c_a, c_b, c_c, lambda);
+ fprintf(debug, "Pull ax^2+bx+c=0: a=%e b=%e c=%e lambda=%e\n", c_a, c_b,
+ c_c, lambda);
}
}
/* The position corrections dr due to the constraints */
- dsvmul(-lambda*rm*pgrp1->invtm, r_ij[c], dr1);
- dsvmul( lambda*rm*pgrp0->invtm, r_ij[c], dr0);
- dr_tot[c] += -lambda*dnorm(r_ij[c]);
+ dsvmul(-lambda * rm * pgrp1->invtm, r_ij[c], dr1);
+ dsvmul(lambda * rm * pgrp0->invtm, r_ij[c], dr0);
+ dr_tot[c] += -lambda * dnorm(r_ij[c]);
break;
case epullgDIR:
case epullgDIRPBC:
for (m = 0; m < DIM; m++)
{
vec[m] = pcrd->spatialData.vec[m];
- a += unc_ij[m]*vec[m];
+ a += unc_ij[m] * vec[m];
}
/* Select only the component along the vector */
dsvmul(a, vec, unc_ij);
}
/* The position corrections dr due to the constraints */
- dsvmul(-lambda*rm*pgrp1->invtm, vec, dr1);
- dsvmul( lambda*rm*pgrp0->invtm, vec, dr0);
+ dsvmul(-lambda * rm * pgrp1->invtm, vec, dr1);
+ dsvmul(lambda * rm * pgrp0->invtm, vec, dr0);
dr_tot[c] += -lambda;
break;
- default:
- gmx_incons("Invalid enumeration value for eGeom");
+ default: gmx_incons("Invalid enumeration value for eGeom");
}
/* DEBUG */
g1 = pcrd->params.group[1];
low_get_pull_coord_dr(pull, pcrd, pbc, rnew[g1], rnew[g0], -1, tmp);
low_get_pull_coord_dr(pull, pcrd, pbc, dr1, dr0, -1, tmp3);
- fprintf(debug,
- "Pull cur %8.5f %8.5f %8.5f j:%8.5f %8.5f %8.5f d: %8.5f\n",
- rnew[g0][0], rnew[g0][1], rnew[g0][2],
- rnew[g1][0], rnew[g1][1], rnew[g1][2], dnorm(tmp));
- fprintf(debug,
- "Pull ref %8s %8s %8s %8s %8s %8s d: %8.5f\n",
- "", "", "", "", "", "", pcrd->value_ref);
- fprintf(debug,
- "Pull cor %8.5f %8.5f %8.5f j:%8.5f %8.5f %8.5f d: %8.5f\n",
- dr0[0], dr0[1], dr0[2],
- dr1[0], dr1[1], dr1[2],
- dnorm(tmp3));
+ fprintf(debug, "Pull cur %8.5f %8.5f %8.5f j:%8.5f %8.5f %8.5f d: %8.5f\n", rnew[g0][0],
+ rnew[g0][1], rnew[g0][2], rnew[g1][0], rnew[g1][1], rnew[g1][2], dnorm(tmp));
+ fprintf(debug, "Pull ref %8s %8s %8s %8s %8s %8s d: %8.5f\n", "", "", "", "", "",
+ "", pcrd->value_ref);
+ fprintf(debug, "Pull cor %8.5f %8.5f %8.5f j:%8.5f %8.5f %8.5f d: %8.5f\n", dr0[0],
+ dr0[1], dr0[2], dr1[0], dr1[1], dr1[2], dnorm(tmp3));
} /* END DEBUG */
/* Update the COMs with dr */
}
/* Check if all constraints are fullfilled now */
- for (pull_coord_work_t &coord : pull->coord)
+ for (pull_coord_work_t& coord : pull->coord)
{
if (coord.params.eType != epullCONSTRAINT)
{
continue;
}
- low_get_pull_coord_dr(pull, &coord, pbc,
- rnew[coord.params.group[1]],
- rnew[coord.params.group[0]],
- -1, unc_ij);
+ low_get_pull_coord_dr(pull, &coord, pbc, rnew[coord.params.group[1]],
+ rnew[coord.params.group[0]], -1, unc_ij);
switch (coord.params.eGeom)
{
case epullgDIST:
- bConverged =
- fabs(dnorm(unc_ij) - coord.value_ref) < pull->params.constr_tol;
+ bConverged = fabs(dnorm(unc_ij) - coord.value_ref) < pull->params.constr_tol;
break;
case epullgDIR:
case epullgDIRPBC:
}
inpr = diprod(unc_ij, vec);
dsvmul(inpr, vec, unc_ij);
- bConverged =
- fabs(diprod(unc_ij, vec) - coord.value_ref) < pull->params.constr_tol;
+ bConverged = fabs(diprod(unc_ij, vec) - coord.value_ref) < pull->params.constr_tol;
break;
}
{
if (debug)
{
- fprintf(debug, "Pull constraint not converged: "
+ fprintf(debug,
+ "Pull constraint not converged: "
"groups %d %d,"
"d_ref = %f, current d = %f\n",
- coord.params.group[0], coord.params.group[1],
- coord.value_ref, dnorm(unc_ij));
+ coord.params.group[0], coord.params.group[1], coord.value_ref, dnorm(unc_ij));
}
bConverged_all = FALSE;
if (v)
{
- invdt = 1/dt;
+ invdt = 1 / dt;
}
/* update atoms in the groups */
for (size_t g = 0; g < pull->group.size(); g++)
{
- const pull_group_work_t *pgrp;
+ const pull_group_work_t* pgrp;
dvec dr;
pgrp = &pull->group[g];
/* update the atom positions */
auto localAtomIndices = pgrp->atomSet.localIndex();
copy_dvec(dr, tmp);
- for (gmx::index j = 0; j < localAtomIndices.size(); j++)
+ for (gmx::index j = 0; j < localAtomIndices.ssize(); j++)
{
ii = localAtomIndices[j];
if (!pgrp->localWeights.empty())
{
- dsvmul(pgrp->wscale*pgrp->localWeights[j], dr, tmp);
+ dsvmul(pgrp->wscale * pgrp->localWeights[j], dr, tmp);
}
for (m = 0; m < DIM; m++)
{
{
for (m = 0; m < DIM; m++)
{
- v[ii][m] += invdt*tmp[m];
+ v[ii][m] += invdt * tmp[m];
}
}
}
/* calculate the constraint forces, used for output and virial only */
for (size_t c = 0; c < pull->coord.size(); c++)
{
- pull_coord_work_t *pcrd;
+ pull_coord_work_t* pcrd;
pcrd = &pull->coord[c];
}
/* Accumulate the forces, in case we have multiple constraint steps */
- double force = dr_tot[c]/((pull->group[pcrd->params.group[0]].invtm + pull->group[pcrd->params.group[1]].invtm)*dt*dt);
+ double force =
+ dr_tot[c]
+ / ((pull->group[pcrd->params.group[0]].invtm + pull->group[pcrd->params.group[1]].invtm)
+ * dt * dt);
pcrd->scalarForce += force;
if (vir != nullptr && pcrd->params.eGeom != epullgDIRPBC && bMaster)
/* Add the pull contribution to the virial */
/* We have already checked above that r_ij[c] != 0 */
- f_invr = pcrd->scalarForce/dnorm(r_ij[c]);
+ f_invr = pcrd->scalarForce / dnorm(r_ij[c]);
for (j = 0; j < DIM; j++)
{
for (m = 0; m < DIM; m++)
{
- vir[j][m] -= 0.5*f_invr*r_ij[c][j]*r_ij[c][m];
+ vir[j][m] -= 0.5 * f_invr * r_ij[c][j] * r_ij[c][m];
}
}
}
{
for (int m = 0; m < DIM; m++)
{
- vir[j][m] -= 0.5*f[j]*dr[m];
+ vir[j][m] -= 0.5 * f[j] * dr[m];
}
}
}
/* Adds the pull contribution to the virial */
-static void add_virial_coord(tensor vir,
- const pull_coord_work_t &pcrd,
- const PullCoordVectorForces &forces)
+static void add_virial_coord(tensor vir, const pull_coord_work_t& pcrd, const PullCoordVectorForces& forces)
{
if (vir != nullptr && pcrd.params.eGeom != epullgDIRPBC)
{
}
}
-static void calc_pull_coord_scalar_force_and_potential(pull_coord_work_t *pcrd,
- double dev, real lambda,
- real *V, real *dVdl)
+static void calc_pull_coord_scalar_force_and_potential(pull_coord_work_t* pcrd,
+ double dev,
+ real lambda,
+ real* V,
+ real* dVdl)
{
- real k, dkdl;
+ real k, dkdl;
- k = (1.0 - lambda)*pcrd->params.k + lambda*pcrd->params.kB;
+ k = (1.0 - lambda) * pcrd->params.k + lambda * pcrd->params.kB;
dkdl = pcrd->params.kB - pcrd->params.k;
switch (pcrd->params.eType)
* potential is that a flat-bottom is zero above or below
the reference value.
*/
- if ((pcrd->params.eType == epullFLATBOTTOM && dev < 0) ||
- (pcrd->params.eType == epullFLATBOTTOMHIGH && dev > 0))
+ if ((pcrd->params.eType == epullFLATBOTTOM && dev < 0)
+ || (pcrd->params.eType == epullFLATBOTTOMHIGH && dev > 0))
{
dev = 0;
}
- pcrd->scalarForce = -k*dev;
- *V += 0.5* k*gmx::square(dev);
- *dVdl += 0.5*dkdl*gmx::square(dev);
+ pcrd->scalarForce = -k * dev;
+ *V += 0.5 * k * gmx::square(dev);
+ *dVdl += 0.5 * dkdl * gmx::square(dev);
break;
case epullCONST_F:
- pcrd->scalarForce = -k;
- *V += k*pcrd->spatialData.value;
- *dVdl += dkdl*pcrd->spatialData.value;
+ pcrd->scalarForce = -k;
+ *V += k * pcrd->spatialData.value;
+ *dVdl += dkdl * pcrd->spatialData.value;
break;
case epullEXTERNAL:
- gmx_incons("the scalar pull force should not be calculated internally for pull type external");
- default:
- gmx_incons("Unsupported pull type in do_pull_pot");
+ gmx_incons(
+ "the scalar pull force should not be calculated internally for pull type "
+ "external");
+ default: gmx_incons("Unsupported pull type in do_pull_pot");
}
}
-static PullCoordVectorForces
-calculateVectorForces(const pull_coord_work_t &pcrd)
+static PullCoordVectorForces calculateVectorForces(const pull_coord_work_t& pcrd)
{
- const t_pull_coord ¶ms = pcrd.params;
- const PullCoordSpatialData &spatialData = pcrd.spatialData;
+ const t_pull_coord& params = pcrd.params;
+ const PullCoordSpatialData& spatialData = pcrd.spatialData;
/* The geometry of the coordinate determines how the scalar force relates to the force on each group */
- PullCoordVectorForces forces;
+ PullCoordVectorForces forces;
if (params.eGeom == epullgDIST)
{
- double invdr01 = spatialData.value > 0 ? 1./spatialData.value : 0.;
+ double invdr01 = spatialData.value > 0 ? 1. / spatialData.value : 0.;
for (int m = 0; m < DIM; m++)
{
- forces.force01[m] = pcrd.scalarForce*spatialData.dr01[m]*invdr01;
+ forces.force01[m] = pcrd.scalarForce * spatialData.dr01[m] * invdr01;
}
}
else if (params.eGeom == epullgANGLE)
* f23 = -dV/dtheta*(a*dr01/|dr01| - b*dr23*|dr23|)/|dr23|,
*/
double a = -gmx::invsqrt(1 - cos_theta2); /* comes from d/dx acos(x) */
- double b = a*cos_theta;
- double invdr01 = 1./dnorm(spatialData.dr01);
- double invdr23 = 1./dnorm(spatialData.dr23);
+ double b = a * cos_theta;
+ double invdr01 = 1. / dnorm(spatialData.dr01);
+ double invdr23 = 1. / dnorm(spatialData.dr23);
dvec normalized_dr01, normalized_dr23;
dsvmul(invdr01, spatialData.dr01, normalized_dr01);
dsvmul(invdr23, spatialData.dr23, normalized_dr23);
for (int m = 0; m < DIM; m++)
{
/* Here, f_scal is -dV/dtheta */
- forces.force01[m] = pcrd.scalarForce*invdr01*(a*normalized_dr23[m] - b*normalized_dr01[m]);
- forces.force23[m] = pcrd.scalarForce*invdr23*(a*normalized_dr01[m] - b*normalized_dr23[m]);
+ forces.force01[m] =
+ pcrd.scalarForce * invdr01 * (a * normalized_dr23[m] - b * normalized_dr01[m]);
+ forces.force23[m] =
+ pcrd.scalarForce * invdr23 * (a * normalized_dr01[m] - b * normalized_dr23[m]);
}
}
else
double invdr01;
dvec normalized_dr01;
- invdr01 = 1./dnorm(spatialData.dr01);
+ invdr01 = 1. / dnorm(spatialData.dr01);
dsvmul(invdr01, spatialData.dr01, normalized_dr01);
- a = -gmx::invsqrt(1 - cos_theta2); /* comes from d/dx acos(x) */
- b = a*cos_theta;
+ a = -gmx::invsqrt(1 - cos_theta2); /* comes from d/dx acos(x) */
+ b = a * cos_theta;
for (int m = 0; m < DIM; m++)
{
- forces.force01[m] = pcrd.scalarForce*invdr01*(a*spatialData.vec[m] - b*normalized_dr01[m]);
+ forces.force01[m] =
+ pcrd.scalarForce * invdr01 * (a * spatialData.vec[m] - b * normalized_dr01[m]);
}
}
else
n2 = diprod(spatialData.planevec_n, spatialData.planevec_n);
dsvmul(-1, spatialData.dr23, dr32);
sqrdist_32 = diprod(dr32, dr32);
- tol = sqrdist_32*GMX_REAL_EPS; /* Avoid tiny angles */
+ tol = sqrdist_32 * GMX_REAL_EPS; /* Avoid tiny angles */
if ((m2 > tol) && (n2 > tol))
{
double a_01, a_23_01, a_23_45, a_45;
double inv_dist_32, inv_sqrdist_32, dist_32;
dvec u, v;
inv_dist_32 = gmx::invsqrt(sqrdist_32);
- inv_sqrdist_32 = inv_dist_32*inv_dist_32;
- dist_32 = sqrdist_32*inv_dist_32;
+ inv_sqrdist_32 = inv_dist_32 * inv_dist_32;
+ dist_32 = sqrdist_32 * inv_dist_32;
/* Forces on groups 0, 1 */
- a_01 = pcrd.scalarForce*dist_32/m2; /* scalarForce is -dV/dphi */
- dsvmul(-a_01, spatialData.planevec_m, forces.force01); /* added sign to get force on group 1, not 0 */
+ a_01 = pcrd.scalarForce * dist_32 / m2; /* scalarForce is -dV/dphi */
+ dsvmul(-a_01, spatialData.planevec_m,
+ forces.force01); /* added sign to get force on group 1, not 0 */
/* Forces on groups 4, 5 */
- a_45 = -pcrd.scalarForce*dist_32/n2;
- dsvmul(a_45, spatialData.planevec_n, forces.force45); /* force on group 5 */
+ a_45 = -pcrd.scalarForce * dist_32 / n2;
+ dsvmul(a_45, spatialData.planevec_n, forces.force45); /* force on group 5 */
/* Force on groups 2, 3 (defining the axis) */
- a_23_01 = -diprod(spatialData.dr01, dr32)*inv_sqrdist_32;
- a_23_45 = -diprod(spatialData.dr45, dr32)*inv_sqrdist_32;
+ a_23_01 = -diprod(spatialData.dr01, dr32) * inv_sqrdist_32;
+ a_23_45 = -diprod(spatialData.dr45, dr32) * inv_sqrdist_32;
dsvmul(-a_23_01, forces.force01, u); /* added sign to get force from group 0, not 1 */
dsvmul(a_23_45, forces.force45, v);
- dvec_sub(u, v, forces.force23); /* force on group 3 */
+ dvec_sub(u, v, forces.force23); /* force on group 3 */
}
else
{
{
for (int m = 0; m < DIM; m++)
{
- forces.force01[m] = pcrd.scalarForce*spatialData.vec[m];
+ forces.force01[m] = pcrd.scalarForce * spatialData.vec[m];
}
}
using Lock = gmx::lock_guard<gmx::Mutex>;
-void register_external_pull_potential(struct pull_t *pull,
- int coord_index,
- const char *provider)
+void register_external_pull_potential(struct pull_t* pull, int coord_index, const char* provider)
{
GMX_RELEASE_ASSERT(pull != nullptr, "register_external_pull_potential called before init_pull");
- GMX_RELEASE_ASSERT(provider != nullptr, "register_external_pull_potential called with NULL as provider name");
+ GMX_RELEASE_ASSERT(provider != nullptr,
+ "register_external_pull_potential called with NULL as provider name");
- if (coord_index < 0 || coord_index >= static_cast<int>(pull->coord.size()))
+ if (coord_index < 0 || coord_index >= gmx::ssize(pull->coord))
{
- gmx_fatal(FARGS, "Module '%s' attempted to register an external potential for pull coordinate %d which is out of the pull coordinate range %d - %zu\n",
+ gmx_fatal(FARGS,
+ "Module '%s' attempted to register an external potential for pull coordinate %d "
+ "which is out of the pull coordinate range %d - %zu\n",
provider, coord_index + 1, 1, pull->coord.size());
}
- pull_coord_work_t *pcrd = &pull->coord[coord_index];
+ pull_coord_work_t* pcrd = &pull->coord[coord_index];
if (pcrd->params.eType != epullEXTERNAL)
{
- gmx_fatal(FARGS, "Module '%s' attempted to register an external potential for pull coordinate %d which of type '%s', whereas external potentials are only supported with type '%s'",
- provider, coord_index + 1, epull_names[pcrd->params.eType], epull_names[epullEXTERNAL]);
+ gmx_fatal(
+ FARGS,
+ "Module '%s' attempted to register an external potential for pull coordinate %d "
+ "which of type '%s', whereas external potentials are only supported with type '%s'",
+ provider, coord_index + 1, epull_names[pcrd->params.eType], epull_names[epullEXTERNAL]);
}
- GMX_RELEASE_ASSERT(pcrd->params.externalPotentialProvider != nullptr, "The external potential provider string for a pull coordinate is NULL");
+ GMX_RELEASE_ASSERT(pcrd->params.externalPotentialProvider != nullptr,
+ "The external potential provider string for a pull coordinate is NULL");
if (gmx_strcasecmp(provider, pcrd->params.externalPotentialProvider) != 0)
{
- gmx_fatal(FARGS, "Module '%s' attempted to register an external potential for pull coordinate %d which expects the external potential to be provided by a module named '%s'",
+ gmx_fatal(FARGS,
+ "Module '%s' attempted to register an external potential for pull coordinate %d "
+ "which expects the external potential to be provided by a module named '%s'",
provider, coord_index + 1, pcrd->params.externalPotentialProvider);
}
if (pcrd->bExternalPotentialProviderHasBeenRegistered)
{
- gmx_fatal(FARGS, "Module '%s' attempted to register an external potential for pull coordinate %d more than once",
+ gmx_fatal(FARGS,
+ "Module '%s' attempted to register an external potential for pull coordinate %d "
+ "more than once",
provider, coord_index + 1);
}
pcrd->bExternalPotentialProviderHasBeenRegistered = true;
pull->numUnregisteredExternalPotentials--;
- GMX_RELEASE_ASSERT(pull->numUnregisteredExternalPotentials >= 0, "Negative unregistered potentials, the pull code is inconsistent");
+ GMX_RELEASE_ASSERT(pull->numUnregisteredExternalPotentials >= 0,
+ "Negative unregistered potentials, the pull code is inconsistent");
}
-static void check_external_potential_registration(const struct pull_t *pull)
+static void check_external_potential_registration(const struct pull_t* pull)
{
if (pull->numUnregisteredExternalPotentials > 0)
{
size_t c;
for (c = 0; c < pull->coord.size(); c++)
{
- if (pull->coord[c].params.eType == epullEXTERNAL &&
- !pull->coord[c].bExternalPotentialProviderHasBeenRegistered)
+ if (pull->coord[c].params.eType == epullEXTERNAL
+ && !pull->coord[c].bExternalPotentialProviderHasBeenRegistered)
{
break;
}
}
- GMX_RELEASE_ASSERT(c < pull->coord.size(), "Internal inconsistency in the pull potential provider counting");
+ GMX_RELEASE_ASSERT(c < pull->coord.size(),
+ "Internal inconsistency in the pull potential provider counting");
- gmx_fatal(FARGS, "No external provider for external pull potentials have been provided for %d pull coordinates. The first coordinate without provider is number %zu, which expects a module named '%s' to provide the external potential.",
- pull->numUnregisteredExternalPotentials,
- c + 1, pull->coord[c].params.externalPotentialProvider);
+ gmx_fatal(FARGS,
+ "No external provider for external pull potentials have been provided for %d "
+ "pull coordinates. The first coordinate without provider is number %zu, which "
+ "expects a module named '%s' to provide the external potential.",
+ pull->numUnregisteredExternalPotentials, c + 1,
+ pull->coord[c].params.externalPotentialProvider);
}
}
* potential energy is added either to the pull term or to a term
* specific to the external module.
*/
-void apply_external_pull_coord_force(struct pull_t *pull,
+void apply_external_pull_coord_force(struct pull_t* pull,
int coord_index,
double coord_force,
- const t_mdatoms *mdatoms,
- gmx::ForceWithVirial *forceWithVirial)
+ const t_mdatoms* mdatoms,
+ gmx::ForceWithVirial* forceWithVirial)
{
- pull_coord_work_t *pcrd;
+ pull_coord_work_t* pcrd;
- GMX_ASSERT(coord_index >= 0 && coord_index < static_cast<int>(pull->coord.size()), "apply_external_pull_coord_force called with coord_index out of range");
+ GMX_ASSERT(coord_index >= 0 && coord_index < gmx::ssize(pull->coord),
+ "apply_external_pull_coord_force called with coord_index out of range");
if (pull->comm.bParticipate)
{
pcrd = &pull->coord[coord_index];
- GMX_RELEASE_ASSERT(pcrd->params.eType == epullEXTERNAL, "The pull force can only be set externally on pull coordinates of external type");
+ GMX_RELEASE_ASSERT(
+ pcrd->params.eType == epullEXTERNAL,
+ "The pull force can only be set externally on pull coordinates of external type");
- GMX_ASSERT(pcrd->bExternalPotentialProviderHasBeenRegistered, "apply_external_pull_coord_force called for an unregistered pull coordinate");
+ GMX_ASSERT(pcrd->bExternalPotentialProviderHasBeenRegistered,
+ "apply_external_pull_coord_force called for an unregistered pull coordinate");
/* Set the force */
pcrd->scalarForce = coord_force;
/* Calculate the pull potential and scalar force for a pull coordinate.
* Returns the vector forces for the pull coordinate.
*/
-static PullCoordVectorForces
-do_pull_pot_coord(struct pull_t *pull, int coord_ind, t_pbc *pbc,
- double t, real lambda,
- real *V, tensor vir, real *dVdl)
+static PullCoordVectorForces do_pull_pot_coord(struct pull_t* pull,
+ int coord_ind,
+ t_pbc* pbc,
+ double t,
+ real lambda,
+ real* V,
+ tensor vir,
+ real* dVdl)
{
- pull_coord_work_t &pcrd = pull->coord[coord_ind];
+ pull_coord_work_t& pcrd = pull->coord[coord_ind];
assert(pcrd.params.eType != epullCONSTRAINT);
return pullCoordForces;
}
-real pull_potential(struct pull_t *pull, const t_mdatoms *md, t_pbc *pbc,
- const t_commrec *cr, double t, real lambda,
- const rvec *x, gmx::ForceWithVirial *force, real *dvdlambda)
+real pull_potential(struct pull_t* pull,
+ const t_mdatoms* md,
+ t_pbc* pbc,
+ const t_commrec* cr,
+ double t,
+ real lambda,
+ const rvec* x,
+ gmx::ForceWithVirial* force,
+ real* dvdlambda)
{
real V = 0;
pull_calc_coms(cr, pull, md, pbc, t, x, nullptr);
- rvec *f = as_rvec_array(force->force_.data());
- matrix virial = { { 0 } };
- const bool computeVirial = (force->computeVirial_ && MASTER(cr));
+ rvec* f = as_rvec_array(force->force_.data());
+ matrix virial = { { 0 } };
+ const bool computeVirial = (force->computeVirial_ && MASTER(cr));
for (size_t c = 0; c < pull->coord.size(); c++)
{
- pull_coord_work_t *pcrd;
+ pull_coord_work_t* pcrd;
pcrd = &pull->coord[c];
- /* For external potential the force is assumed to be given by an external module by a call to
- apply_pull_coord_external_force */
+ /* For external potential the force is assumed to be given by an external module by a
+ call to apply_pull_coord_external_force */
if (pcrd->params.eType == epullCONSTRAINT || pcrd->params.eType == epullEXTERNAL)
{
continue;
}
- PullCoordVectorForces pullCoordForces =
- do_pull_pot_coord(pull, c, pbc, t, lambda,
- &V,
- computeVirial ? virial : nullptr,
- &dVdl);
+ PullCoordVectorForces pullCoordForces = do_pull_pot_coord(
+ pull, c, pbc, t, lambda, &V, computeVirial ? virial : nullptr, &dVdl);
/* Distribute the force over the atoms in the pulled groups */
apply_forces_coord(pull, c, pullCoordForces, md, f);
}
}
- GMX_ASSERT(pull->numExternalPotentialsStillToBeAppliedThisStep == 0, "Too few or too many external pull potentials have been applied the previous step");
+ GMX_ASSERT(pull->numExternalPotentialsStillToBeAppliedThisStep == 0,
+ "Too few or too many external pull potentials have been applied the previous step");
/* All external pull potentials still need to be applied */
- pull->numExternalPotentialsStillToBeAppliedThisStep =
- pull->numCoordinatesWithExternalPotential;
+ pull->numExternalPotentialsStillToBeAppliedThisStep = pull->numCoordinatesWithExternalPotential;
return (MASTER(cr) ? V : 0.0);
}
-void pull_constraint(struct pull_t *pull, const t_mdatoms *md, t_pbc *pbc,
- const t_commrec *cr, double dt, double t,
- rvec *x, rvec *xp, rvec *v, tensor vir)
+void pull_constraint(struct pull_t* pull,
+ const t_mdatoms* md,
+ t_pbc* pbc,
+ const t_commrec* cr,
+ double dt,
+ double t,
+ rvec* x,
+ rvec* xp,
+ rvec* v,
+ tensor vir)
{
assert(pull != nullptr);
}
}
-void dd_make_local_pull_groups(const t_commrec *cr, struct pull_t *pull)
+void dd_make_local_pull_groups(const t_commrec* cr, struct pull_t* pull)
{
- gmx_domdec_t *dd;
- pull_comm_t *comm;
- gmx_bool bMustParticipate;
+ gmx_domdec_t* dd;
+ pull_comm_t* comm;
+ gmx_bool bMustParticipate;
dd = cr->dd;
*/
bMustParticipate = (comm->bParticipateAll || comm->isMasterRank);
- for (pull_group_work_t &group : pull->group)
+ for (pull_group_work_t& group : pull->group)
{
if (!group.globalWeights.empty())
{
}
}
- GMX_ASSERT(bMustParticipate || dd != nullptr, "Either all ranks (including this rank) participate, or we use DD and need to have access to dd here");
+ GMX_ASSERT(bMustParticipate || dd != nullptr,
+ "Either all ranks (including this rank) participate, or we use DD and need to "
+ "have access to dd here");
/* We should participate if we have pull or pbc atoms */
- if (!bMustParticipate &&
- (group.atomSet.numAtomsLocal() > 0 || (group.epgrppbc == epgrppbcREFAT &&
- group.pbcAtomSet->numAtomsLocal() > 0)))
+ if (!bMustParticipate
+ && (group.atomSet.numAtomsLocal() > 0
+ || (group.epgrppbc == epgrppbcREFAT && group.pbcAtomSet->numAtomsLocal() > 0)))
{
bMustParticipate = TRUE;
}
* if they needed to participate up to 20 decompositions ago.
* This avoids frequent rebuilds due to atoms jumping back and forth.
*/
- const int64_t history_count = 20;
- gmx_bool bWillParticipate;
- int count[2];
+ const int64_t history_count = 20;
+ gmx_bool bWillParticipate;
+ int count[2];
/* Increase the decomposition counter for the current call */
comm->setup_count++;
}
bWillParticipate =
- bMustParticipate ||
- (comm->bParticipate &&
- comm->must_count >= comm->setup_count - history_count);
+ bMustParticipate
+ || (comm->bParticipate && comm->must_count >= comm->setup_count - history_count);
if (debug && dd != nullptr)
{
- fprintf(debug, "Our DD rank (%3d) pull #atoms>0 or master: %s, will be part %s\n",
- dd->rank, gmx::boolToString(bMustParticipate), gmx::boolToString(bWillParticipate));
+ fprintf(debug, "Our DD rank (%3d) pull #atoms>0 or master: %s, will be part %s\n", dd->rank,
+ gmx::boolToString(bMustParticipate), gmx::boolToString(bWillParticipate));
}
if (bWillParticipate)
/* If we are missing ranks or if we have 20% more ranks than needed
* we make a new sub-communicator.
*/
- if (count[1] > 0 || 6*count[0] < 5*comm->nparticipate)
+ if (count[1] > 0 || 6 * count[0] < 5 * comm->nparticipate)
{
if (debug)
{
- fprintf(debug, "Creating new pull subcommunicator of size %d\n",
- count[0]);
+ fprintf(debug, "Creating new pull subcommunicator of size %d\n", count[0]);
}
#if GMX_MPI
* to avoid this splitting as much as possible.
*/
assert(dd != nullptr);
- MPI_Comm_split(dd->mpi_comm_all, bWillParticipate ? 0 : 1, dd->rank,
- &comm->mpi_comm_com);
+ MPI_Comm_split(dd->mpi_comm_all, bWillParticipate ? 0 : 1, dd->rank, &comm->mpi_comm_com);
#endif
comm->bParticipate = bWillParticipate;
/* When we use the previous COM for PBC, we need to broadcast
* the previous COM to ranks that have joined the communicator.
*/
- for (pull_group_work_t &group : pull->group)
+ for (pull_group_work_t& group : pull->group)
{
if (group.epgrppbc == epgrppbcPREVSTEPCOM)
{
GMX_ASSERT(comm->bParticipate || !MASTER(cr),
- "The master rank has to participate, as it should pass an up to date prev. COM "
+ "The master rank has to participate, as it should pass an up to "
+ "date prev. COM "
"to bcast here as well as to e.g. checkpointing");
gmx_bcast(sizeof(dvec), group.x_prev_step, cr);
pull->bSetPBCatoms = TRUE;
}
-static void init_pull_group_index(FILE *fplog, const t_commrec *cr,
- int g, pull_group_work_t *pg,
- gmx_bool bConstraint, const ivec pulldim_con,
- const gmx_mtop_t *mtop,
- const t_inputrec *ir, real lambda)
+static void init_pull_group_index(FILE* fplog,
+ const t_commrec* cr,
+ int g,
+ pull_group_work_t* pg,
+ gmx_bool bConstraint,
+ const ivec pulldim_con,
+ const gmx_mtop_t* mtop,
+ const t_inputrec* ir,
+ real lambda)
{
/* With EM and BD there are no masses in the integrator.
* But we still want to have the correct mass-weighted COMs.
* So we store the real masses in the weights.
*/
- const bool setWeights = (pg->params.nweight > 0 ||
- EI_ENERGY_MINIMIZATION(ir->eI) ||
- ir->eI == eiBD);
+ const bool setWeights = (pg->params.nweight > 0 || EI_ENERGY_MINIMIZATION(ir->eI) || ir->eI == eiBD);
/* In parallel, store we need to extract localWeights from weights at DD time */
- std::vector<real> &weights = ((cr && PAR(cr)) ? pg->globalWeights : pg->localWeights);
+ std::vector<real>& weights = ((cr && PAR(cr)) ? pg->globalWeights : pg->localWeights);
- const gmx_groups_t *groups = &mtop->groups;
+ const SimulationGroups& groups = mtop->groups;
/* Count frozen dimensions and (weighted) mass */
int nfrozen = 0;
{
for (int d = 0; d < DIM; d++)
{
- if (pulldim_con[d] == 1 &&
- ir->opts.nFreeze[getGroupType(groups, egcFREEZE, ii)][d])
+ if (pulldim_con[d] == 1
+ && ir->opts.nFreeze[getGroupType(groups, SimulationAtomGroupType::Freeze, ii)][d])
{
nfrozen++;
}
}
}
- const t_atom &atom = mtopGetAtomParameters(mtop, ii, &molb);
+ const t_atom& atom = mtopGetAtomParameters(mtop, ii, &molb);
real m;
if (ir->efep == efepNO)
{
}
else
{
- m = (1 - lambda)*atom.m + lambda*atom.mB;
+ m = (1 - lambda) * atom.m + lambda * atom.mB;
}
real w;
if (pg->params.nweight > 0)
if (EI_ENERGY_MINIMIZATION(ir->eI))
{
/* Move the mass to the weight */
- w *= m;
- m = 1;
+ w *= m;
+ m = 1;
}
else if (ir->eI == eiBD)
{
real mbd;
- if (ir->bd_fric != 0.0f)
+ if (ir->bd_fric != 0.0F)
{
- mbd = ir->bd_fric*ir->delta_t;
+ mbd = ir->bd_fric * ir->delta_t;
}
else
{
- if (groups->grpnr[egcTC] == nullptr)
+ if (groups.groupNumbers[SimulationAtomGroupType::TemperatureCoupling].empty())
{
- mbd = ir->delta_t/ir->opts.tau_t[0];
+ mbd = ir->delta_t / ir->opts.tau_t[0];
}
else
{
- mbd = ir->delta_t/ir->opts.tau_t[groups->grpnr[egcTC][ii]];
+ mbd = ir->delta_t
+ / ir->opts.tau_t[groups.groupNumbers[SimulationAtomGroupType::TemperatureCoupling][ii]];
}
}
- w *= m/mbd;
- m = mbd;
+ w *= m / mbd;
+ m = mbd;
}
if (setWeights)
{
weights.push_back(w);
}
- tmass += m;
- wmass += m*w;
- wwmass += m*w*w;
+ tmass += m;
+ wmass += m * w;
+ wwmass += m * w * w;
}
if (wmass == 0)
}
if (fplog)
{
- fprintf(fplog,
- "Pull group %d: %5d atoms, mass %9.3f",
- g, pg->params.nat, tmass);
- if (pg->params.weight ||
- EI_ENERGY_MINIMIZATION(ir->eI) || ir->eI == eiBD)
+ fprintf(fplog, "Pull group %d: %5d atoms, mass %9.3f", g, pg->params.nat, tmass);
+ if (pg->params.weight || EI_ENERGY_MINIMIZATION(ir->eI) || ir->eI == eiBD)
{
- fprintf(fplog, ", weighted mass %9.3f", wmass*wmass/wwmass);
+ fprintf(fplog, ", weighted mass %9.3f", wmass * wmass / wwmass);
}
if (pg->epgrppbc == epgrppbcCOS)
{
if (nfrozen == 0)
{
/* A value != 0 signals not frozen, it is updated later */
- pg->invtm = -1.0;
+ pg->invtm = -1.0;
}
else
{
int ndim = 0;
for (int d = 0; d < DIM; d++)
{
- ndim += pulldim_con[d]*pg->params.nat;
+ ndim += pulldim_con[d] * pg->params.nat;
}
if (fplog && nfrozen > 0 && nfrozen < ndim)
{
}
}
-struct pull_t *
-init_pull(FILE *fplog, const pull_params_t *pull_params, const t_inputrec *ir,
- const gmx_mtop_t *mtop, const t_commrec *cr, gmx::LocalAtomSetManager *atomSets,
- real lambda)
+struct pull_t* init_pull(FILE* fplog,
+ const pull_params_t* pull_params,
+ const t_inputrec* ir,
+ const gmx_mtop_t* mtop,
+ const t_commrec* cr,
+ gmx::LocalAtomSetManager* atomSets,
+ real lambda)
{
- struct pull_t *pull;
- pull_comm_t *comm;
+ struct pull_t* pull;
+ pull_comm_t* comm;
pull = new pull_t;
/* Copy the pull parameters */
- pull->params = *pull_params;
+ pull->params = *pull_params;
/* Avoid pointer copies */
pull->params.group = nullptr;
pull->params.coord = nullptr;
for (int i = 0; i < pull_params->ngroup; ++i)
{
- pull->group.emplace_back(pull_params->group[i], atomSets->add({pull_params->group[i].ind, pull_params->group[i].ind+pull_params->group[i].nat}),
+ pull->group.emplace_back(pull_params->group[i],
+ atomSets->add({ pull_params->group[i].ind,
+ pull_params->group[i].ind + pull_params->group[i].nat }),
pull_params->bSetPbcRefToPrevStepCOM);
}
if (cr != nullptr && DOMAINDECOMP(cr))
{
/* Set up the global to local atom mapping for PBC atoms */
- for (pull_group_work_t &group : pull->group)
+ for (pull_group_work_t& group : pull->group)
{
if (group.epgrppbc == epgrppbcREFAT || group.epgrppbc == epgrppbcPREVSTEPCOM)
{
/* pbcAtomSet consists of a single atom */
- group.pbcAtomSet = gmx::compat::make_unique<gmx::LocalAtomSet>(atomSets->add({&group.params.pbcatom, &group.params.pbcatom + 1}));
+ group.pbcAtomSet = std::make_unique<gmx::LocalAtomSet>(
+ atomSets->add({ &group.params.pbcatom, &group.params.pbcatom + 1 }));
}
}
}
pull->bXOutAverage = pull_params->bXOutAverage;
pull->bFOutAverage = pull_params->bFOutAverage;
- GMX_RELEASE_ASSERT(pull->group[0].params.nat == 0, "pull group 0 is an absolute reference group and should not contain atoms");
+ GMX_RELEASE_ASSERT(pull->group[0].params.nat == 0,
+ "pull group 0 is an absolute reference group and should not contain atoms");
pull->numCoordinatesWithExternalPotential = 0;
/* Construct a pull coordinate, copying all coordinate parameters */
pull->coord.emplace_back(pull_params->coord[c]);
- pull_coord_work_t *pcrd = &pull->coord.back();
+ pull_coord_work_t* pcrd = &pull->coord.back();
switch (pcrd->params.eGeom)
{
case epullgDIST:
- case epullgDIRRELATIVE: /* Direction vector is determined at each step */
+ case epullgDIRRELATIVE: /* Direction vector is determined at each step */
case epullgANGLE:
- case epullgDIHEDRAL:
- break;
+ case epullgDIHEDRAL: break;
case epullgDIR:
case epullgDIRPBC:
case epullgCYL:
* the string corresponding to the geometry, which will result
* in printing "UNDEFINED".
*/
- gmx_fatal(FARGS, "Pull geometry not supported for pull coordinate %d. The geometry enum %d in the input is larger than that supported by the code (up to %d). You are probably reading a tpr file generated with a newer version of Gromacs with an binary from an older version of Gromacs.",
+ gmx_fatal(FARGS,
+ "Pull geometry not supported for pull coordinate %d. The geometry enum "
+ "%d in the input is larger than that supported by the code (up to %d). "
+ "You are probably reading a tpr file generated with a newer version of "
+ "Gromacs with an binary from an older version of Gromacs.",
c + 1, pcrd->params.eGeom, epullgNR - 1);
}
if (pcrd->params.eType == epullCONSTRAINT)
{
/* Check restrictions of the constraint pull code */
- if (pcrd->params.eGeom == epullgCYL ||
- pcrd->params.eGeom == epullgDIRRELATIVE ||
- pcrd->params.eGeom == epullgANGLE ||
- pcrd->params.eGeom == epullgDIHEDRAL ||
- pcrd->params.eGeom == epullgANGLEAXIS)
+ if (pcrd->params.eGeom == epullgCYL || pcrd->params.eGeom == epullgDIRRELATIVE
+ || pcrd->params.eGeom == epullgANGLE || pcrd->params.eGeom == epullgDIHEDRAL
+ || pcrd->params.eGeom == epullgANGLEAXIS)
{
- gmx_fatal(FARGS, "Pulling of type %s can not be combined with geometry %s. Consider using pull type %s.",
- epull_names[pcrd->params.eType],
- epullg_names[pcrd->params.eGeom],
+ gmx_fatal(FARGS,
+ "Pulling of type %s can not be combined with geometry %s. Consider using "
+ "pull type %s.",
+ epull_names[pcrd->params.eType], epullg_names[pcrd->params.eGeom],
epull_names[epullUMBRELLA]);
}
{
pull->bCylinder = TRUE;
}
- else if (pcrd->params.eGeom == epullgANGLE || pcrd->params.eGeom == epullgDIHEDRAL || pcrd->params.eGeom == epullgANGLEAXIS)
+ else if (pcrd->params.eGeom == epullgANGLE || pcrd->params.eGeom == epullgDIHEDRAL
+ || pcrd->params.eGeom == epullgANGLEAXIS)
{
pull->bAngle = TRUE;
}
for (int i = 0; i < pcrd->params.ngroup; i++)
{
int groupIndex = pcrd->params.group[i];
- if (groupIndex > 0 &&
- !(pcrd->params.eGeom == epullgCYL && i == 0))
+ if (groupIndex > 0 && !(pcrd->params.eGeom == epullgCYL && i == 0))
{
pull->group[groupIndex].needToCalcCom = true;
}
/* Initialize the constant reference value */
if (pcrd->params.eType != epullEXTERNAL)
{
- low_set_pull_coord_reference_value(pcrd, c, pcrd->params.init*pull_conversion_factor_userinput2internal(&pcrd->params));
+ low_set_pull_coord_reference_value(
+ pcrd, c,
+ pcrd->params.init * pull_conversion_factor_userinput2internal(&pcrd->params));
}
else
{
if (pcrd->params.eType == epullEXTERNAL)
{
- GMX_RELEASE_ASSERT(pcrd->params.rate == 0, "With an external potential, a pull coordinate should have rate = 0");
+ GMX_RELEASE_ASSERT(
+ pcrd->params.rate == 0,
+ "With an external potential, a pull coordinate should have rate = 0");
/* This external potential needs to be registered later */
pull->numCoordinatesWithExternalPotential++;
pcrd->bExternalPotentialProviderHasBeenRegistered = false;
}
- pull->numUnregisteredExternalPotentials =
- pull->numCoordinatesWithExternalPotential;
+ pull->numUnregisteredExternalPotentials = pull->numCoordinatesWithExternalPotential;
pull->numExternalPotentialsStillToBeAppliedThisStep = 0;
pull->ePBC = ir->ePBC;
switch (pull->ePBC)
{
case epbcNONE: pull->npbcdim = 0; break;
- case epbcXY: pull->npbcdim = 2; break;
- default: pull->npbcdim = 3; break;
+ case epbcXY: pull->npbcdim = 2; break;
+ default: pull->npbcdim = 3; break;
}
if (fplog)
gmx_bool bAbs, bCos;
bAbs = FALSE;
- for (const pull_coord_work_t &coord : pull->coord)
+ for (const pull_coord_work_t& coord : pull->coord)
{
- if (pull->group[coord.params.group[0]].params.nat == 0 ||
- pull->group[coord.params.group[1]].params.nat == 0)
+ if (pull->group[coord.params.group[0]].params.nat == 0
+ || pull->group[coord.params.group[1]].params.nat == 0)
{
bAbs = TRUE;
}
}
// Don't include the reference group 0 in output, so we report ngroup-1
int numRealGroups = pull->group.size() - 1;
- GMX_RELEASE_ASSERT(numRealGroups > 0, "The reference absolute position pull group should always be present");
- fprintf(fplog, "with %zu pull coordinate%s and %d group%s\n",
- pull->coord.size(), pull->coord.size() == 1 ? "" : "s",
- numRealGroups, numRealGroups == 1 ? "" : "s");
+ GMX_RELEASE_ASSERT(numRealGroups > 0,
+ "The reference absolute position pull group should always be present");
+ fprintf(fplog, "with %zu pull coordinate%s and %d group%s\n", pull->coord.size(),
+ pull->coord.size() == 1 ? "" : "s", numRealGroups, numRealGroups == 1 ? "" : "s");
if (bAbs)
{
fprintf(fplog, "with an absolute reference\n");
// Don't include the reference group 0 in loop
for (size_t g = 1; g < pull->group.size(); g++)
{
- if (pull->group[g].params.nat > 1 &&
- pull->group[g].params.pbcatom < 0)
+ if (pull->group[g].params.nat > 1 && pull->group[g].params.pbcatom < 0)
{
/* We are using cosine weighting */
fprintf(fplog, "Cosine weighting is used for group %zu\n", g);
}
}
- pull->bRefAt = FALSE;
- pull->cosdim = -1;
+ pull->bRefAt = FALSE;
+ pull->cosdim = -1;
for (size_t g = 0; g < pull->group.size(); g++)
{
- pull_group_work_t *pgrp;
+ pull_group_work_t* pgrp;
- pgrp = &pull->group[g];
+ pgrp = &pull->group[g];
if (pgrp->params.nat > 0)
{
/* There is an issue when a group is used in multiple coordinates
bConstraint = FALSE;
clear_ivec(pulldim);
clear_ivec(pulldim_con);
- for (const pull_coord_work_t &coord : pull->coord)
+ for (const pull_coord_work_t& coord : pull->coord)
{
gmx_bool bGroupUsed = FALSE;
for (int gi = 0; gi < coord.params.ngroup; gi++)
*/
switch (pgrp->epgrppbc)
{
- case epgrppbcREFAT:
- pull->bRefAt = TRUE;
- break;
+ case epgrppbcREFAT: pull->bRefAt = TRUE; break;
case epgrppbcCOS:
if (pgrp->params.weight != nullptr)
{
- gmx_fatal(FARGS, "Pull groups can not have relative weights and cosine weighting at same time");
+ gmx_fatal(FARGS,
+ "Pull groups can not have relative weights and cosine weighting "
+ "at same time");
}
for (int m = 0; m < DIM; m++)
{
{
if (pull->cosdim >= 0 && pull->cosdim != m)
{
- gmx_fatal(FARGS, "Can only use cosine weighting with pulling in one dimension (use mdp option pull-coord?-dim)");
+ gmx_fatal(FARGS,
+ "Can only use cosine weighting with pulling in one "
+ "dimension (use mdp option pull-coord?-dim)");
}
pull->cosdim = m;
}
}
break;
- case epgrppbcNONE:
- break;
+ case epgrppbcNONE: break;
}
/* Set the indices */
- init_pull_group_index(fplog, cr, g, pgrp,
- bConstraint, pulldim_con,
- mtop, ir, lambda);
+ init_pull_group_index(fplog, cr, g, pgrp, bConstraint, pulldim_con, mtop, ir, lambda);
}
else
{
/* If we use cylinder coordinates, do some initialising for them */
if (pull->bCylinder)
{
- for (const pull_coord_work_t &coord : pull->coord)
+ for (const pull_coord_work_t& coord : pull->coord)
{
if (coord.params.eGeom == epullgCYL)
{
if (pull->group[coord.params.group[0]].params.nat == 0)
{
- gmx_fatal(FARGS, "A cylinder pull group is not supported when using absolute reference!\n");
+ gmx_fatal(FARGS,
+ "A cylinder pull group is not supported when using absolute "
+ "reference!\n");
}
}
- const auto &referenceGroup = pull->group[coord.params.group[0]];
- pull->dyna.emplace_back(referenceGroup.params, referenceGroup.atomSet, pull->params.bSetPbcRefToPrevStepCOM);
+ const auto& referenceGroup = pull->group[coord.params.group[0]];
+ pull->dyna.emplace_back(referenceGroup.params, referenceGroup.atomSet,
+ pull->params.bSetPbcRefToPrevStepCOM);
}
}
* when we have an external pull potential, since then the external
* potential provider expects each rank to have the coordinate.
*/
- comm->bParticipateAll = (cr == nullptr || !DOMAINDECOMP(cr) ||
- cr->dd->nnodes <= 32 ||
- pull->numCoordinatesWithExternalPotential > 0 ||
- getenv("GMX_PULL_PARTICIPATE_ALL") != nullptr);
+ comm->bParticipateAll = (cr == nullptr || !DOMAINDECOMP(cr) || cr->dd->nnodes <= 32
+ || pull->numCoordinatesWithExternalPotential > 0
+ || getenv("GMX_PULL_PARTICIPATE_ALL") != nullptr);
/* This sub-commicator is not used with comm->bParticipateAll,
* so we can always initialize it to NULL.
*/
- comm->mpi_comm_com = MPI_COMM_NULL;
- comm->nparticipate = 0;
- comm->isMasterRank = (cr == nullptr || MASTER(cr));
+ comm->mpi_comm_com = MPI_COMM_NULL;
+ comm->nparticipate = 0;
+ comm->isMasterRank = (cr == nullptr || MASTER(cr));
#else
/* No MPI: 1 rank: all ranks pull */
comm->bParticipateAll = TRUE;
comm->isMasterRank = true;
#endif
- comm->bParticipate = comm->bParticipateAll;
- comm->setup_count = 0;
- comm->must_count = 0;
+ comm->bParticipate = comm->bParticipateAll;
+ comm->setup_count = 0;
+ comm->must_count = 0;
if (!comm->bParticipateAll && fplog != nullptr)
{
}
comm->pbcAtomBuffer.resize(pull->group.size());
- comm->comBuffer.resize(pull->group.size()*c_comBufferStride);
+ comm->comBuffer.resize(pull->group.size() * c_comBufferStride);
if (pull->bCylinder)
{
- comm->cylinderBuffer.resize(pull->coord.size()*c_cylinderBufferStride);
+ comm->cylinderBuffer.resize(pull->coord.size() * c_cylinderBufferStride);
}
/* We still need to initialize the PBC reference coordinates */
return pull;
}
-static void destroy_pull(struct pull_t *pull)
+static void destroy_pull(struct pull_t* pull)
{
#if GMX_MPI
if (pull->comm.mpi_comm_com != MPI_COMM_NULL)
delete pull;
}
-void preparePrevStepPullCom(const t_inputrec *ir, const t_mdatoms *md, t_state *state, const t_state *state_global, const t_commrec *cr, bool startingFromCheckpoint)
+void preparePrevStepPullCom(const t_inputrec* ir,
+ pull_t* pull_work,
+ const t_mdatoms* md,
+ t_state* state,
+ const t_state* state_global,
+ const t_commrec* cr,
+ bool startingFromCheckpoint)
{
if (!ir->pull || !ir->pull->bSetPbcRefToPrevStepCOM)
{
return;
}
- allocStatePrevStepPullCom(state, ir->pull_work);
+ allocStatePrevStepPullCom(state, pull_work);
if (startingFromCheckpoint)
{
if (MASTER(cr))
/* Only the master rank has the checkpointed COM from the previous step */
gmx_bcast(sizeof(double) * state->pull_com_prev_step.size(), &state->pull_com_prev_step[0], cr);
}
- setPrevStepPullComFromState(ir->pull_work, state);
+ setPrevStepPullComFromState(pull_work, state);
}
else
{
t_pbc pbc;
set_pbc(&pbc, ir->ePBC, state->box);
- initPullComFromPrevStep(cr, ir->pull_work, md, &pbc, state->x.rvec_array());
- updatePrevStepPullCom(ir->pull_work, state);
+ initPullComFromPrevStep(cr, pull_work, md, &pbc, state->x.rvec_array());
+ updatePrevStepPullCom(pull_work, state);
}
}
-void finish_pull(struct pull_t *pull)
+void finish_pull(struct pull_t* pull)
{
check_external_potential_registration(pull);
destroy_pull(pull);
}
-gmx_bool pull_have_potential(const struct pull_t *pull)
+gmx_bool pull_have_potential(const struct pull_t* pull)
{
return pull->bPotential;
}
-gmx_bool pull_have_constraint(const struct pull_t *pull)
+gmx_bool pull_have_constraint(const struct pull_t* pull)
{
return pull->bConstraint;
}
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff