a3bcbc1fb23b5b123c9792fb787cf03e6b66356a
[alexxy/gromacs.git] / src / gromacs / applied_forces / densityfitting / densityfittingamplitudelookup.cpp
1 /*
2  * This file is part of the GROMACS molecular simulation package.
3  *
4  * Copyright (c) 2019,2020,2021, 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 amplitude lookup for density fitting
38  *
39  * \author Christian Blau <blau@kth.se>
40  * \ingroup module_applied_forces
41  */
42 #include "gmxpre.h"
43
44 #include "densityfittingamplitudelookup.h"
45
46 #include <algorithm>
47 #include <vector>
48
49 #include "gromacs/utility/arrayref.h"
50
51 namespace gmx
52 {
53 /********************************************************************
54  * DensityFittingAmplitudeLookup::Impl
55  */
56
57 class DensityFittingAmplitudeLookupImpl
58 {
59 public:
60     DensityFittingAmplitudeLookupImpl()                                         = default;
61     DensityFittingAmplitudeLookupImpl(const DensityFittingAmplitudeLookupImpl&) = default;
62     virtual ~DensityFittingAmplitudeLookupImpl()                                = default;
63
64     virtual const std::vector<real>& operator()(ArrayRef<const real> /*chargeA*/,
65                                                 ArrayRef<const real> /*massT*/,
66                                                 ArrayRef<const int> localIndex) = 0;
67     virtual std::unique_ptr<DensityFittingAmplitudeLookupImpl> clone()          = 0;
68 };
69
70 namespace
71 {
72 class UnitAmplitudes final : public DensityFittingAmplitudeLookupImpl
73 {
74 public:
75     UnitAmplitudes()                      = default;
76     UnitAmplitudes(const UnitAmplitudes&) = default;
77     ~UnitAmplitudes() override            = default;
78     std::unique_ptr<DensityFittingAmplitudeLookupImpl> clone() override;
79     const std::vector<real>&                           operator()(ArrayRef<const real> /*chargeA*/,
80                                         ArrayRef<const real> /*massT*/,
81                                         ArrayRef<const int> localIndex) override;
82
83 private:
84     std::vector<real> amplitude_;
85 };
86
87 std::unique_ptr<DensityFittingAmplitudeLookupImpl> UnitAmplitudes::clone()
88 {
89     return std::make_unique<UnitAmplitudes>(*this);
90 };
91
92 const std::vector<real>& UnitAmplitudes::operator()(ArrayRef<const real> /*chargeA*/,
93                                                     ArrayRef<const real> /*massT*/,
94                                                     ArrayRef<const int> localIndex)
95 {
96     if (amplitude_.size() != localIndex.size())
97     {
98         amplitude_.resize(localIndex.size(), 1.);
99     }
100
101     return amplitude_;
102 }
103
104 class ChargesAsAmplitudes final : public DensityFittingAmplitudeLookupImpl
105 {
106 public:
107     ChargesAsAmplitudes()                           = default;
108     ChargesAsAmplitudes(const ChargesAsAmplitudes&) = default;
109     ~ChargesAsAmplitudes() override                 = default;
110     std::unique_ptr<DensityFittingAmplitudeLookupImpl> clone() override;
111     const std::vector<real>&                           operator()(ArrayRef<const real> chargeA,
112                                         ArrayRef<const real> /*massT*/,
113                                         ArrayRef<const int> localIndex) override;
114
115 private:
116     std::vector<real> amplitude_;
117 };
118
119 std::unique_ptr<DensityFittingAmplitudeLookupImpl> ChargesAsAmplitudes::clone()
120 {
121     return std::make_unique<ChargesAsAmplitudes>(*this);
122 };
123
124 const std::vector<real>& ChargesAsAmplitudes::operator()(ArrayRef<const real> chargeA,
125                                                          ArrayRef<const real> /*massT*/,
126                                                          ArrayRef<const int> localIndex)
127 {
128     if (amplitude_.size() != localIndex.size())
129     {
130         amplitude_.resize(localIndex.size());
131     }
132
133     std::transform(std::begin(localIndex),
134                    std::end(localIndex),
135                    std::begin(amplitude_),
136                    [&chargeA](gmx::index index) { return chargeA[index]; });
137     return amplitude_;
138 }
139
140 class MassesAsAmplitudes final : public DensityFittingAmplitudeLookupImpl
141 {
142 public:
143     MassesAsAmplitudes()                          = default;
144     MassesAsAmplitudes(const MassesAsAmplitudes&) = default;
145     ~MassesAsAmplitudes() override                = default;
146     std::unique_ptr<DensityFittingAmplitudeLookupImpl> clone() override;
147     const std::vector<real>&                           operator()(ArrayRef<const real> /*chargeA*/,
148                                         ArrayRef<const real> massT,
149                                         ArrayRef<const int>  localIndex) override;
150
151 private:
152     std::vector<real> amplitude_;
153 };
154
155 std::unique_ptr<DensityFittingAmplitudeLookupImpl> MassesAsAmplitudes::clone()
156 {
157     return std::make_unique<MassesAsAmplitudes>(*this);
158 };
159
160 const std::vector<real>& MassesAsAmplitudes::operator()(ArrayRef<const real> /*chargeA*/,
161                                                         ArrayRef<const real> massT,
162                                                         ArrayRef<const int>  localIndex)
163 {
164     if (amplitude_.size() != localIndex.size())
165     {
166         amplitude_.resize(localIndex.size());
167     }
168
169     std::transform(std::begin(localIndex),
170                    std::end(localIndex),
171                    std::begin(amplitude_),
172                    [&massT](gmx::index index) { return massT[index]; });
173     return amplitude_;
174 }
175
176 } // namespace
177
178 /********************************************************************
179  * DensityFittingAmplitudeLookup
180  */
181
182
183 DensityFittingAmplitudeLookup::DensityFittingAmplitudeLookup(const DensityFittingAmplitudeMethod& method)
184 {
185     switch (method)
186     {
187         case DensityFittingAmplitudeMethod::Unity:
188             impl_ = std::make_unique<UnitAmplitudes>();
189             break;
190         case DensityFittingAmplitudeMethod::Mass:
191             impl_ = std::make_unique<MassesAsAmplitudes>();
192             break;
193         case DensityFittingAmplitudeMethod::Charge:
194             impl_ = std::make_unique<ChargesAsAmplitudes>();
195             break;
196         default: break;
197     }
198 }
199
200 const std::vector<real>& DensityFittingAmplitudeLookup::operator()(ArrayRef<const real> chargeA,
201                                                                    ArrayRef<const real> massT,
202                                                                    ArrayRef<const int>  localIndex)
203 {
204     return (*impl_)(chargeA, massT, localIndex);
205 }
206
207 DensityFittingAmplitudeLookup::~DensityFittingAmplitudeLookup() = default;
208
209 DensityFittingAmplitudeLookup::DensityFittingAmplitudeLookup(const DensityFittingAmplitudeLookup& other) :
210     impl_(other.impl_->clone())
211 {
212 }
213
214 DensityFittingAmplitudeLookup& DensityFittingAmplitudeLookup::operator=(const DensityFittingAmplitudeLookup& other)
215 {
216     impl_ = other.impl_->clone();
217     return *this;
218 }
219
220 DensityFittingAmplitudeLookup::DensityFittingAmplitudeLookup(DensityFittingAmplitudeLookup&&) noexcept = default;
221
222 DensityFittingAmplitudeLookup& DensityFittingAmplitudeLookup::
223                                operator=(DensityFittingAmplitudeLookup&&) noexcept = default;
224
225 } // namespace gmx