src/gromacs/applied_forces/densityfittingforceprovider.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   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.
   8  *
   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.
  13  *
  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.
  18  *
  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.
  23  *
  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.
  31  *
  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.
  34  */
  35 /*! \internal \file
  36  * \brief
  37  * Implements force provider for density fitting
  38  *
  39  * \author Christian Blau <blau@kth.se>
  40  * \ingroup module_applied_forces
  41  */
  42 #include "gmxpre.h"
  43
  44 #include "densityfittingforceprovider.h"
  45
  46 #include "gromacs/gmxlib/network.h"
  47 #include "gromacs/math/densityfit.h"
  48 #include "gromacs/math/densityfittingforce.h"
  49 #include "gromacs/math/gausstransform.h"
  50 #include "gromacs/mdtypes/commrec.h"
  51 #include "gromacs/mdtypes/enerdata.h"
  52 #include "gromacs/mdtypes/forceoutput.h"
  53 #include "gromacs/mdtypes/iforceprovider.h"
  54 #include "gromacs/pbcutil/pbc.h"
  55
  56 #include "densityfittingamplitudelookup.h"
  57 #include "densityfittingparameters.h"
  58
  59 namespace gmx
  60 {
  61
  62 namespace
  63 {
  64
  65 /*! \internal \brief Generate the spread kernal from Gaussian parameters.
  66  *
  67  * \param[in] sigma the width of the Gaussian to be spread
  68  * \param[in] nSigma the range of the Gaussian in multiples of sigma
  69  * \param[in] scaleToLattice the coordinate transformation into the spreading lattice
  70  * \returns A Gauss-transform kernel shape
  71  */
  72 GaussianSpreadKernelParameters::Shape
  73 makeSpreadKernel(real sigma, real nSigma, const ScaleCoordinates &scaleToLattice)
  74 {
  75     RVec sigmaInLatticeCoordinates {
  76         sigma, sigma, sigma
  77     };
  78     scaleToLattice( { &sigmaInLatticeCoordinates, &sigmaInLatticeCoordinates + 1 });
  79     return {
  80                DVec {
  81                    sigmaInLatticeCoordinates[XX], sigmaInLatticeCoordinates[YY], sigmaInLatticeCoordinates[ZZ]
  82                }, nSigma
  83     };
  84 }
  85
  86 }   // namespace
  87
  88 /********************************************************************
  89  * DensityFittingForceProvider::Impl
  90  */
  91
  92 class DensityFittingForceProvider::Impl
  93 {
  94     public:
  95         //! \copydoc DensityFittingForceProvider(const DensityFittingParameters &parameters)
  96         Impl(const DensityFittingParameters &parameters,
  97              basic_mdspan<const float, dynamicExtents3D> referenceDensity,
  98              const TranslateAndScale &transformationToDensityLattice,
  99              const LocalAtomSet &localAtomSet,
 100              int pbcType);
 101         ~Impl();
 102         void calculateForces(const ForceProviderInput &forceProviderInput, ForceProviderOutput *forceProviderOutput);
 103     private:
 104         const DensityFittingParameters       &parameters_;
 105         LocalAtomSet                          localAtomSet_;
 106
 107         GaussianSpreadKernelParameters::Shape spreadKernel_;
 108         GaussTransform3D                      gaussTransform_;
 109         DensitySimilarityMeasure              measure_;
 110         DensityFittingForce                   densityFittingForce_;
 111         //! the local atom coordinates transformed into the grid coordinate system
 112         std::vector<RVec>                     transformedCoordinates_;
 113         std::vector<RVec>                     forces_;
 114         DensityFittingAmplitudeLookup         amplitudeLookup_;
 115         TranslateAndScale                     transformationToDensityLattice_;
 116         RVec                                  referenceDensityCenter_;
 117         int                                   pbcType_;
 118 };
 119
 120 DensityFittingForceProvider::Impl::~Impl() = default;
 121
 122 DensityFittingForceProvider::Impl::Impl(const DensityFittingParameters &parameters,
 123                                         basic_mdspan<const float, dynamicExtents3D> referenceDensity,
 124                                         const TranslateAndScale &transformationToDensityLattice,
 125                                         const LocalAtomSet &localAtomSet,
 126                                         int pbcType) :
 127     parameters_(parameters),
 128     localAtomSet_(localAtomSet),
 129     spreadKernel_(makeSpreadKernel(parameters_.gaussianTransformSpreadingWidth_,
 130                                    parameters_.gaussianTransformSpreadingRangeInMultiplesOfWidth_,
 131                                    transformationToDensityLattice.scaleOperationOnly())),
 132     gaussTransform_(referenceDensity.extents(), spreadKernel_),
 133     measure_(parameters.similarityMeasureMethod_, referenceDensity),
 134     densityFittingForce_(spreadKernel_),
 135     transformedCoordinates_(localAtomSet_.numAtomsLocal()),
 136     amplitudeLookup_(parameters_.amplitudeLookupMethod_),
 137     transformationToDensityLattice_(transformationToDensityLattice),
 138     pbcType_(pbcType)
 139 {
 140     referenceDensityCenter_  = {
 141         real(referenceDensity.extent(XX))/2,
 142         real(referenceDensity.extent(YY))/2,
 143         real(referenceDensity.extent(ZZ))/2
 144     };
 145     transformationToDensityLattice_.scaleOperationOnly().inverseIgnoringZeroScale(
 146             { &referenceDensityCenter_, &referenceDensityCenter_ + 1 });
 147 }
 148
 149 void DensityFittingForceProvider::Impl::calculateForces(const ForceProviderInput &forceProviderInput,
 150                                                         ForceProviderOutput      *forceProviderOutput)
 151 {
 152     // do nothing if there are no density fitting atoms on this node
 153     if (localAtomSet_.numAtomsLocal() == 0)
 154     {
 155         return;
 156     }
 157     transformedCoordinates_.resize(localAtomSet_.numAtomsLocal());
 158     // pick and copy atom coordinates
 159     std::transform(std::cbegin(localAtomSet_.localIndex()),
 160                    std::cend(localAtomSet_.localIndex()),
 161                    std::begin(transformedCoordinates_),
 162                    [&forceProviderInput](int index) { return forceProviderInput.x_[index]; });
 163
 164     // pick periodic image that is closest to the center of the reference density
 165     {
 166         t_pbc pbc;
 167         set_pbc(&pbc, pbcType_, forceProviderInput.box_);
 168         for (RVec &x : transformedCoordinates_)
 169         {
 170             rvec dx;
 171             pbc_dx(&pbc, x, referenceDensityCenter_, dx);
 172             x = referenceDensityCenter_ + dx;
 173         }
 174     }
 175
 176     // transform local atom coordinates to density grid coordinates
 177     transformationToDensityLattice_(transformedCoordinates_);
 178
 179     // spread atoms on grid
 180     gaussTransform_.setZero();
 181
 182     const std::vector<real> &amplitudes        = amplitudeLookup_(forceProviderInput.mdatoms_, localAtomSet_.localIndex());
 183     auto                     amplitudeIterator = amplitudes.cbegin();
 184     for (const auto &r : transformedCoordinates_)
 185     {
 186         gaussTransform_.add({ r, *amplitudeIterator });
 187         ++amplitudeIterator;
 188     }
 189
 190     // communicate grid
 191     if (PAR(&forceProviderInput.cr_))
 192     {
 193         // \todo update to real once GaussTransform class returns real
 194         gmx_sumf(gaussTransform_.view().mapping().required_span_size(),
 195                  gaussTransform_.view().data(), &forceProviderInput.cr_);
 196     }
 197
 198     // calculate grid derivative
 199     const DensitySimilarityMeasure::density &densityDerivative =
 200         measure_.gradient(gaussTransform_.constView());
 201     // calculate forces
 202     forces_.resize(localAtomSet_.numAtomsLocal());
 203     std::transform(
 204             std::begin(transformedCoordinates_),
 205             std::end(transformedCoordinates_),
 206             std::begin(amplitudes),
 207             std::begin(forces_),
 208             [&densityDerivative, this](const RVec r, real amplitude)
 209             {
 210                 return densityFittingForce_.evaluateForce({r, amplitude}, densityDerivative);
 211             }
 212             );
 213
 214     transformationToDensityLattice_.scaleOperationOnly().inverseIgnoringZeroScale(forces_);
 215
 216     auto densityForceIterator = forces_.cbegin();
 217     for (const auto localAtomIndex : localAtomSet_.localIndex())
 218     {
 219         forceProviderOutput->forceWithVirial_.force_[localAtomIndex] +=
 220             parameters_.forceConstant_ * *densityForceIterator;
 221         ++densityForceIterator;
 222     }
 223
 224     // calculate corresponding potential energy
 225     const float similarity  = measure_.similarity(gaussTransform_.constView());
 226     const real  energy      = -similarity * parameters_.forceConstant_;
 227     forceProviderOutput->enerd_.term[F_DENSITYFITTING] += energy;
 228 }
 229
 230 /********************************************************************
 231  * DensityFittingForceProvider
 232  */
 233
 234 DensityFittingForceProvider::~DensityFittingForceProvider() = default;
 235
 236 DensityFittingForceProvider::DensityFittingForceProvider(const DensityFittingParameters &parameters,
 237                                                          basic_mdspan<const float, dynamicExtents3D> referenceDensity,
 238                                                          const TranslateAndScale &transformationToDensityLattice,
 239                                                          const LocalAtomSet &localAtomSet,
 240                                                          int pbcType)
 241     : impl_(new Impl(parameters, referenceDensity, transformationToDensityLattice, localAtomSet, pbcType))
 242 {}
 243
 244 void DensityFittingForceProvider::calculateForces(const ForceProviderInput  &forceProviderInput,
 245                                                   ForceProviderOutput      * forceProviderOutput)
 246 {
 247     impl_->calculateForces(forceProviderInput, forceProviderOutput);
 248 }
 249
 250 } // namespace gmx