2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2019,2020, 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 density fitting forces.
39 * \author Christian Blau <blau@kth.se>
40 * \ingroup module_math
44 #include "densityfittingforce.h"
46 #include "gromacs/math/functions.h"
47 #include "gromacs/math/multidimarray.h"
52 /********************************************************************
53 * DensityFittingForce::Impl
57 * Private implementation class for DensityFittingForce.
59 class DensityFittingForce::Impl
62 /*! \brief Construct densityfitting force implementation from
63 * spread of function used to generate the density and spread range.
64 * \param[in] kernelShapeParameters determine the shape of the spreading kernel
66 explicit Impl(const GaussianSpreadKernelParameters::Shape
& kernelShapeParameters
);
67 //! \copydoc DensityFittingForce::evaluateForce(const DensitySpreadKernelParameters::PositionAndAmplitude & localParameters, basic_mdspan<const float, dynamicExtents3D> densityDerivative)
68 RVec
evaluateForce(const GaussianSpreadKernelParameters::PositionAndAmplitude
& localParameters
,
69 basic_mdspan
<const float, dynamicExtents3D
> densityDerivative
);
70 //! The width of the Gaussian in lattice spacing units
72 //! The spread range in lattice points
73 IVec latticeSpreadRange_
;
74 //! The three one-dimensional Gaussians that are used in the force calculation
75 std::array
<GaussianOn1DLattice
, DIM
> gauss1d_
;
76 //! The outer product of a Gaussian along the z and y dimension
77 OuterProductEvaluator outerProductZY_
;
80 DensityFittingForce::Impl::Impl(const GaussianSpreadKernelParameters::Shape
& kernelShapeParameters
) :
81 sigma_
{ kernelShapeParameters
.sigma_
},
82 latticeSpreadRange_
{ kernelShapeParameters
.latticeSpreadRange()[XX
],
83 kernelShapeParameters
.latticeSpreadRange()[YY
],
84 kernelShapeParameters
.latticeSpreadRange()[ZZ
] },
85 gauss1d_({ GaussianOn1DLattice(latticeSpreadRange_
[XX
], sigma_
[XX
]),
86 GaussianOn1DLattice(latticeSpreadRange_
[YY
], sigma_
[YY
]),
87 GaussianOn1DLattice(latticeSpreadRange_
[ZZ
], sigma_
[ZZ
]) })
91 RVec
DensityFittingForce::Impl::evaluateForce(const GaussianSpreadKernelParameters::PositionAndAmplitude
& localParameters
,
92 basic_mdspan
<const float, dynamicExtents3D
> densityDerivative
)
94 const IVec
closestLatticePoint(roundToInt(localParameters
.coordinate_
[XX
]),
95 roundToInt(localParameters
.coordinate_
[YY
]),
96 roundToInt(localParameters
.coordinate_
[ZZ
]));
97 const auto spreadRange
= spreadRangeWithinLattice(
98 closestLatticePoint
, densityDerivative
.extents(), latticeSpreadRange_
);
100 // do nothing if the added Gaussian will never reach the lattice
101 if (spreadRange
.empty())
106 for (int dimension
= XX
; dimension
<= ZZ
; ++dimension
)
108 // multiply with amplitude so that Gauss3D = (amplitude * Gauss_x) * Gauss_y * Gauss_z
109 const float gauss1DAmplitude
= dimension
> XX
? 1.0 : localParameters
.amplitude_
;
110 gauss1d_
[dimension
].spread(gauss1DAmplitude
, localParameters
.coordinate_
[dimension
]
111 - closestLatticePoint
[dimension
]);
114 const auto spreadZY
= outerProductZY_(gauss1d_
[ZZ
].view(), gauss1d_
[YY
].view());
115 const auto spreadX
= gauss1d_
[XX
].view();
116 const IVec
spreadGridOffset(latticeSpreadRange_
[XX
] - closestLatticePoint
[XX
],
117 latticeSpreadRange_
[YY
] - closestLatticePoint
[YY
],
118 latticeSpreadRange_
[ZZ
] - closestLatticePoint
[ZZ
]);
120 const DVec differenceVectorScale
= { 1. / (square(sigma_
[XX
])), 1. / (square(sigma_
[YY
])),
121 1. / (square(sigma_
[ZZ
])) };
122 const DVec differenceVectorOffset
= scaleByVector(
123 spreadRange
.begin().toDVec() - localParameters
.coordinate_
.toDVec(), differenceVectorScale
);
125 DVec differenceVector
= differenceVectorOffset
;
127 DVec force
= { 0., 0., 0. };
129 for (int zLatticeIndex
= spreadRange
.begin()[ZZ
]; zLatticeIndex
< spreadRange
.end()[ZZ
];
130 ++zLatticeIndex
, differenceVector
[ZZ
] += differenceVectorScale
[ZZ
])
132 auto zSliceOfDerivative
= densityDerivative
[zLatticeIndex
];
134 differenceVector
[YY
] = differenceVectorOffset
[YY
];
135 for (int yLatticeIndex
= spreadRange
.begin()[YY
]; yLatticeIndex
< spreadRange
.end()[YY
];
136 ++yLatticeIndex
, differenceVector
[YY
] += differenceVectorScale
[YY
])
138 auto ySliceOfDerivative
= zSliceOfDerivative
[yLatticeIndex
];
139 const auto zyPrefactor
= spreadZY(zLatticeIndex
+ spreadGridOffset
[ZZ
],
140 yLatticeIndex
+ spreadGridOffset
[YY
]);
142 differenceVector
[XX
] = differenceVectorOffset
[XX
];
143 for (int xLatticeIndex
= spreadRange
.begin()[XX
]; xLatticeIndex
< spreadRange
.end()[XX
];
144 ++xLatticeIndex
, differenceVector
[XX
] += differenceVectorScale
[XX
])
146 const double preFactor
= zyPrefactor
* spreadX
[xLatticeIndex
+ spreadGridOffset
[XX
]]
147 * ySliceOfDerivative
[xLatticeIndex
];
148 force
+= preFactor
* differenceVector
;
152 return localParameters
.amplitude_
* force
.toRVec();
155 /********************************************************************
156 * DensityFittingForce
159 DensityFittingForce::DensityFittingForce(const GaussianSpreadKernelParameters::Shape
& kernelShapeParameters
) :
160 impl_(new Impl(kernelShapeParameters
))
164 RVec
DensityFittingForce::evaluateForce(const GaussianSpreadKernelParameters::PositionAndAmplitude
& localParameters
,
165 basic_mdspan
<const float, dynamicExtents3D
> densityDerivative
)
167 return impl_
->evaluateForce(localParameters
, densityDerivative
);
170 DensityFittingForce::~DensityFittingForce() {}
172 DensityFittingForce::DensityFittingForce(const DensityFittingForce
& other
) :
173 impl_(new Impl(*other
.impl_
))
177 DensityFittingForce
& DensityFittingForce::operator=(const DensityFittingForce
& other
)
179 *impl_
= *other
.impl_
;
183 DensityFittingForce::DensityFittingForce(DensityFittingForce
&&) noexcept
= default;
185 DensityFittingForce
& DensityFittingForce::operator=(DensityFittingForce
&&) noexcept
= default;