2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2019, by the GROMACS development team, led by
5 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
6 * and including many others, as listed in the AUTHORS file in the
7 * top-level source directory and at http://www.gromacs.org.
9 * GROMACS is free software; you can redistribute it and/or
10 * modify it under the terms of the GNU Lesser General Public License
11 * as published by the Free Software Foundation; either version 2.1
12 * of the License, or (at your option) any later version.
14 * GROMACS is distributed in the hope that it will be useful,
15 * but WITHOUT ANY WARRANTY; without even the implied warranty of
16 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17 * Lesser General Public License for more details.
19 * You should have received a copy of the GNU Lesser General Public
20 * License along with GROMACS; if not, see
21 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
22 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
24 * If you want to redistribute modifications to GROMACS, please
25 * consider that scientific software is very special. Version
26 * control is crucial - bugs must be traceable. We will be happy to
27 * consider code for inclusion in the official distribution, but
28 * derived work must not be called official GROMACS. Details are found
29 * in the README & COPYING files - if they are missing, get the
30 * official version at http://www.gromacs.org.
32 * To help us fund GROMACS development, we humbly ask that you cite
33 * the research papers on the package. Check out http://www.gromacs.org.
37 * Implements force provider for density fitting
39 * \author Christian Blau <blau@kth.se>
40 * \ingroup module_applied_forces
44 #include "densityfittingforceprovider.h"
48 #include "gromacs/compat/optional.h"
49 #include "gromacs/gmxlib/network.h"
50 #include "gromacs/math/densityfit.h"
51 #include "gromacs/math/densityfittingforce.h"
52 #include "gromacs/math/exponentialmovingaverage.h"
53 #include "gromacs/math/gausstransform.h"
54 #include "gromacs/mdtypes/commrec.h"
55 #include "gromacs/mdtypes/enerdata.h"
56 #include "gromacs/mdtypes/forceoutput.h"
57 #include "gromacs/mdtypes/iforceprovider.h"
58 #include "gromacs/pbcutil/pbc.h"
60 #include "densityfittingamplitudelookup.h"
61 #include "densityfittingparameters.h"
69 /*! \internal \brief Generate the spread kernel from Gaussian parameters.
71 * \param[in] sigma the width of the Gaussian to be spread
72 * \param[in] nSigma the range of the Gaussian in multiples of sigma
73 * \param[in] scaleToLattice the coordinate transformation into the spreading lattice
74 * \returns A Gauss-transform kernel shape
76 GaussianSpreadKernelParameters::Shape
77 makeSpreadKernel(real sigma
, real nSigma
, const ScaleCoordinates
&scaleToLattice
)
79 RVec sigmaInLatticeCoordinates
{
82 scaleToLattice( { &sigmaInLatticeCoordinates
, &sigmaInLatticeCoordinates
+ 1 });
85 sigmaInLatticeCoordinates
[XX
], sigmaInLatticeCoordinates
[YY
], sigmaInLatticeCoordinates
[ZZ
]
92 /********************************************************************
93 * DensityFittingForceProvider::Impl
96 class DensityFittingForceProvider::Impl
99 //! \copydoc DensityFittingForceProvider(const DensityFittingParameters ¶meters)
100 Impl(const DensityFittingParameters
¶meters
,
101 basic_mdspan
<const float, dynamicExtents3D
> referenceDensity
,
102 const TranslateAndScale
&transformationToDensityLattice
,
103 const LocalAtomSet
&localAtomSet
,
105 double simulationTimeStep
,
106 const DensityFittingForceProviderState
&state
);
108 void calculateForces(const ForceProviderInput
&forceProviderInput
, ForceProviderOutput
*forceProviderOutput
);
110 DensityFittingForceProviderState
state();
113 const DensityFittingParameters
¶meters_
;
114 DensityFittingForceProviderState state_
;
115 LocalAtomSet localAtomSet_
;
117 GaussianSpreadKernelParameters::Shape spreadKernel_
;
118 GaussTransform3D gaussTransform_
;
119 DensitySimilarityMeasure measure_
;
120 DensityFittingForce densityFittingForce_
;
121 //! the local atom coordinates transformed into the grid coordinate system
122 std::vector
<RVec
> transformedCoordinates_
;
123 std::vector
<RVec
> forces_
;
124 DensityFittingAmplitudeLookup amplitudeLookup_
;
125 TranslateAndScale transformationToDensityLattice_
;
126 RVec referenceDensityCenter_
;
129 //! Optionally scale the force according to a moving average of the similarity
130 compat::optional
<ExponentialMovingAverage
> expAverageSimilarity_
;
133 DensityFittingForceProvider::Impl::~Impl() = default;
135 DensityFittingForceProvider::Impl::Impl(const DensityFittingParameters
¶meters
,
136 basic_mdspan
<const float, dynamicExtents3D
> referenceDensity
,
137 const TranslateAndScale
&transformationToDensityLattice
,
138 const LocalAtomSet
&localAtomSet
,
140 double simulationTimeStep
,
141 const DensityFittingForceProviderState
&state
) :
142 parameters_(parameters
),
144 localAtomSet_(localAtomSet
),
145 spreadKernel_(makeSpreadKernel(parameters_
.gaussianTransformSpreadingWidth_
,
146 parameters_
.gaussianTransformSpreadingRangeInMultiplesOfWidth_
,
147 transformationToDensityLattice
.scaleOperationOnly())),
148 gaussTransform_(referenceDensity
.extents(), spreadKernel_
),
149 measure_(parameters
.similarityMeasureMethod_
, referenceDensity
),
150 densityFittingForce_(spreadKernel_
),
151 transformedCoordinates_(localAtomSet_
.numAtomsLocal()),
152 amplitudeLookup_(parameters_
.amplitudeLookupMethod_
),
153 transformationToDensityLattice_(transformationToDensityLattice
),
155 expAverageSimilarity_(compat::nullopt
)
157 if (parameters_
.adaptiveForceScaling_
)
159 GMX_ASSERT(simulationTimeStep
> 0, "Simulation time step must be larger than zero for adaptive for scaling.");
160 expAverageSimilarity_
.emplace(ExponentialMovingAverage(
161 parameters_
.adaptiveForceScalingTimeConstant_
162 / (simulationTimeStep
* parameters_
.calculationIntervalInSteps_
),
163 state
.exponentialMovingAverageState_
));
165 referenceDensityCenter_
= {
166 real(referenceDensity
.extent(XX
))/2,
167 real(referenceDensity
.extent(YY
))/2,
168 real(referenceDensity
.extent(ZZ
))/2
170 transformationToDensityLattice_
.scaleOperationOnly().inverseIgnoringZeroScale(
171 { &referenceDensityCenter_
, &referenceDensityCenter_
+ 1 });
172 // correct the reference density center for a shift
173 // if the reference density does not have its origin at (0,0,0)
174 RVec
referenceDensityOriginShift(0, 0, 0);
175 transformationToDensityLattice_({ &referenceDensityOriginShift
, &referenceDensityOriginShift
+ 1 });
176 transformationToDensityLattice_
.scaleOperationOnly().inverseIgnoringZeroScale(
177 { &referenceDensityOriginShift
, &referenceDensityOriginShift
+ 1 });
178 referenceDensityCenter_
-= referenceDensityOriginShift
;
181 void DensityFittingForceProvider::Impl::calculateForces(const ForceProviderInput
&forceProviderInput
,
182 ForceProviderOutput
*forceProviderOutput
)
184 // do nothing but count number of steps when not in density fitting step
185 if (state_
.stepsSinceLastCalculation_
% parameters_
.calculationIntervalInSteps_
!= 0)
187 ++(state_
.stepsSinceLastCalculation_
);
191 state_
.stepsSinceLastCalculation_
= 0;
193 // do nothing if there are no density fitting atoms on this node
194 if (localAtomSet_
.numAtomsLocal() == 0)
198 transformedCoordinates_
.resize(localAtomSet_
.numAtomsLocal());
199 // pick and copy atom coordinates
200 std::transform(std::cbegin(localAtomSet_
.localIndex()),
201 std::cend(localAtomSet_
.localIndex()),
202 std::begin(transformedCoordinates_
),
203 [&forceProviderInput
](int index
) { return forceProviderInput
.x_
[index
]; });
205 // pick periodic image that is closest to the center of the reference density
208 set_pbc(&pbc
, pbcType_
, forceProviderInput
.box_
);
209 for (RVec
&x
: transformedCoordinates_
)
212 pbc_dx(&pbc
, x
, referenceDensityCenter_
, dx
);
213 x
= referenceDensityCenter_
+ dx
;
217 // transform local atom coordinates to density grid coordinates
218 transformationToDensityLattice_(transformedCoordinates_
);
220 // spread atoms on grid
221 gaussTransform_
.setZero();
223 std::vector
<real
> amplitudes
= amplitudeLookup_(forceProviderInput
.mdatoms_
, localAtomSet_
.localIndex());
225 if (parameters_
.normalizeDensities_
)
227 real sum
= std::accumulate(std::begin(amplitudes
), std::end(amplitudes
), 0.);
228 if (PAR(&forceProviderInput
.cr_
))
230 gmx_sum(1, &sum
, &forceProviderInput
.cr_
);
232 for (real
&litude
: amplitudes
)
238 auto amplitudeIterator
= amplitudes
.cbegin();
240 for (const auto &r
: transformedCoordinates_
)
242 gaussTransform_
.add({ r
, *amplitudeIterator
});
247 if (PAR(&forceProviderInput
.cr_
))
249 // \todo update to real once GaussTransform class returns real
250 gmx_sumf(gaussTransform_
.view().mapping().required_span_size(),
251 gaussTransform_
.view().data(), &forceProviderInput
.cr_
);
254 // calculate grid derivative
255 const DensitySimilarityMeasure::density
&densityDerivative
=
256 measure_
.gradient(gaussTransform_
.constView());
258 forces_
.resize(localAtomSet_
.numAtomsLocal());
260 std::begin(transformedCoordinates_
),
261 std::end(transformedCoordinates_
),
262 std::begin(amplitudes
),
264 [&densityDerivative
, this](const RVec r
, real amplitude
)
266 return densityFittingForce_
.evaluateForce({r
, amplitude
}, densityDerivative
);
270 transformationToDensityLattice_
.scaleOperationOnly().inverseIgnoringZeroScale(forces_
);
272 auto densityForceIterator
= forces_
.cbegin();
273 const real effectiveForceConstant
= state_
.adaptiveForceConstantScale_
*
274 parameters_
.calculationIntervalInSteps_
* parameters_
.forceConstant_
;
275 for (const auto localAtomIndex
: localAtomSet_
.localIndex())
277 forceProviderOutput
->forceWithVirial_
.force_
[localAtomIndex
]
278 += effectiveForceConstant
* *densityForceIterator
;
279 ++densityForceIterator
;
282 // calculate corresponding potential energy
283 const float similarity
= measure_
.similarity(gaussTransform_
.constView());
284 const real energy
= -similarity
* parameters_
.forceConstant_
* state_
.adaptiveForceConstantScale_
;
285 forceProviderOutput
->enerd_
.term
[F_DENSITYFITTING
] += energy
;
287 if (expAverageSimilarity_
.has_value())
289 expAverageSimilarity_
->updateWithDataPoint(similarity
);
291 if (expAverageSimilarity_
->increasing())
293 state_
.adaptiveForceConstantScale_
/= 1._real
+ expAverageSimilarity_
->inverseTimeConstant();
297 state_
.adaptiveForceConstantScale_
*= 1._real
+ expAverageSimilarity_
->inverseTimeConstant();
302 DensityFittingForceProviderState
303 DensityFittingForceProvider::Impl::state()
305 if (expAverageSimilarity_
.has_value())
307 state_
.exponentialMovingAverageState_
= expAverageSimilarity_
->state();
312 /********************************************************************
313 * DensityFittingForceProvider
316 DensityFittingForceProvider::~DensityFittingForceProvider() = default;
318 DensityFittingForceProvider::DensityFittingForceProvider(const DensityFittingParameters
¶meters
,
319 basic_mdspan
<const float, dynamicExtents3D
> referenceDensity
,
320 const TranslateAndScale
&transformationToDensityLattice
,
321 const LocalAtomSet
&localAtomSet
,
323 double simulationTimeStep
,
324 const DensityFittingForceProviderState
&state
)
325 : impl_(new Impl(parameters
, referenceDensity
, transformationToDensityLattice
, localAtomSet
, pbcType
, simulationTimeStep
, state
))
328 void DensityFittingForceProvider::calculateForces(const ForceProviderInput
&forceProviderInput
,
329 ForceProviderOutput
* forceProviderOutput
)
331 impl_
->calculateForces(forceProviderInput
, forceProviderOutput
);
334 DensityFittingForceProviderState
DensityFittingForceProvider::state()
336 return impl_
->state();