Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2017 of Haochuan Chen (excluding colvar_UIestimator.h)
3 : Copyright (c) 2017 of Haohao Fu (colvar_UIestimator.h)
4 :
5 : This program is free software: you can redistribute it and/or modify
6 : it under the terms of the GNU Lesser General Public License as published
7 : by the Free Software Foundation, either version 3 of the License, or
8 : (at your option) any later version.
9 :
10 : This program is distributed in the hope that it will be useful,
11 : but WITHOUT ANY WARRANTY; without even the implied warranty of
12 : MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 : GNU Lesser General Public License for more details.
14 :
15 : You should have received a copy of the GNU Lesser General Public License
16 : along with this program. If not, see <http://www.gnu.org/licenses/>.
17 : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
18 : #ifndef __PLUMED_drr_DRR_h
19 : #define __PLUMED_drr_DRR_h
20 : // Build requirement: boost, c++11 compatible compiler.
21 : #ifdef __PLUMED_HAS_BOOST_SERIALIZATION
22 :
23 : #include <algorithm>
24 : #include <cmath>
25 : #include <cstddef>
26 : #include <cstdlib>
27 : #include <fstream>
28 : #include <iomanip>
29 : #include <iostream>
30 : #include <iterator>
31 : #include <limits>
32 : #include <numeric>
33 : #include <sstream>
34 :
35 : // boost headers for serialization
36 : #include <boost/archive/binary_iarchive.hpp>
37 : #include <boost/archive/binary_oarchive.hpp>
38 : #include <boost/serialization/string.hpp>
39 : #include <boost/serialization/vector.hpp>
40 :
41 : namespace PLMD {
42 : namespace drr {
43 :
44 : /// This class can store the minimum, maximum and bins of a dimension(axis).
45 : class DRRAxis {
46 : public:
47 : /// Default empty constructor
48 22 : DRRAxis() {
49 22 : min = max = 0.0;
50 22 : nbins = 0;
51 22 : periodic = false;
52 22 : domainMax = domainMin = 0.0;
53 22 : binWidth = 0.0;
54 22 : }
55 : /// Constructor using maximum value, minimum value and the number of bins(No
56 : /// pbc)
57 : DRRAxis(double l, double h, size_t n)
58 : : min(l), max(h), nbins(n), periodic(false), domainMax(0), domainMin(0),
59 : binWidth((max - min) / double(nbins)) {}
60 : /// PBC-aware constructor
61 2 : DRRAxis(double l, double h, size_t n, bool pbc, double dMax, double dMin)
62 : : min(l), max(h), nbins(n), periodic(pbc), domainMax(dMax),
63 2 : domainMin(dMin), binWidth((max - min) / double(nbins)) {}
64 : /// Set values
65 6 : void set(double l, double h, size_t n, bool pbc = false, double dmin = 0,
66 : double dmax = 0) {
67 6 : min = l;
68 6 : max = h;
69 6 : nbins = n;
70 6 : periodic = pbc;
71 6 : domainMax = dmax;
72 6 : domainMin = dmin;
73 6 : binWidth = (max - min) / nbins;
74 6 : }
75 : /// Set PBC data
76 4 : void setPeriodicity(double dmin, double dmax) {
77 4 : domainMax = dmax;
78 4 : domainMin = dmin;
79 4 : periodic = true;
80 4 : }
81 : /// Getters
82 10 : double getMin() const { return this->min; }
83 10 : double getMax() const { return this->max; }
84 65514 : double getWidth() const { return binWidth; }
85 4 : double getDomainMax() const { return this->domainMax; }
86 4 : double getDomainMin() const { return this->domainMin; }
87 6 : size_t getBins() const { return this->nbins; }
88 :
89 : /// Check periodicity
90 26 : bool isPeriodic() const { return this->periodic; }
91 : /// Check real periodicity, i.e. the maximum == the domain maximum
92 : bool isRealPeriodic() const;
93 :
94 : /// Check whether x is in this axis
95 : bool isInBoundary(double x) const;
96 : /// Get an array of middle points of each bins
97 : std::vector<double> getMiddlePoints();
98 :
99 : /// Combine two axes if they share the same bin width.
100 : static DRRAxis merge(const DRRAxis &d1, const DRRAxis &d2);
101 :
102 : friend class DRRForceGrid;
103 :
104 : protected:
105 : double min; // Minimum value of the axis
106 : double max; // Maximum value of the axis
107 : size_t nbins; // Number of bins
108 : bool periodic; // Periodicity
109 : double domainMax; // Maximum value of the CV domain
110 : double domainMin; // Minimum value of the CV domain
111 : friend class boost::serialization::access;
112 : /// Use boost serialization
113 : template <typename Archive>
114 16 : void save(Archive &ar, const unsigned int version) const {
115 16 : ar &min;
116 16 : ar &max;
117 16 : ar &nbins;
118 16 : ar &periodic;
119 16 : ar &domainMax;
120 16 : ar &domainMin;
121 16 : }
122 : /// Split save and load. The bin width is calculated after initialization.
123 : template <typename Archive>
124 14 : void load(Archive &ar, const unsigned int version) {
125 14 : ar &min;
126 14 : ar &max;
127 14 : ar &nbins;
128 14 : ar &periodic;
129 14 : ar &domainMax;
130 14 : ar &domainMin;
131 14 : binWidth = (max - min) / double(nbins);
132 14 : }
133 : template <typename Archive>
134 30 : void serialize(Archive &ar, const unsigned int version) {
135 30 : boost::serialization::split_member(ar, *this, version);
136 30 : }
137 :
138 : private:
139 : double binWidth; // bin width
140 : };
141 :
142 : /// A class for collecting instantaneous forces, calculating average forces and
143 : /// build CV histogram.
144 12 : class DRRForceGrid {
145 : public:
146 : /// Empty constructor
147 : DRRForceGrid();
148 : /// "Real" constructor
149 : /// The 2D table vector is mainly used for print grid points in grad and count
150 : /// file.
151 : /// So when use binary output we can set initializeTable to false to save
152 : /// memory.
153 : explicit DRRForceGrid(const std::vector<DRRAxis> &p_dimensions,
154 : const std::string &p_suffix,
155 : bool initializeTable = true);
156 : /// Check whether a point is in this grid
157 : bool isInBoundary(const std::vector<double> &pos) const;
158 : // /// Get internal indices of a point
159 : // std::vector<size_t> index(const std::vector<double> &pos) const;
160 : /// Get internal counts address of a point
161 : size_t sampleAddress(const std::vector<double> &pos) const;
162 : /// Store instantaneous forces of a point
163 : /// nsamples > 1 is useful for merging windows
164 : bool store(const std::vector<double> &pos, const std::vector<double> &f,
165 : unsigned long int nsamples = 1);
166 : /// Get accumulated forces of a point
167 : std::vector<double>
168 : getAccumulatedForces(const std::vector<double> &pos) const;
169 : /// Get counts of a point
170 : unsigned long int getCount(const std::vector<double> &pos,
171 : bool SkipCheck = false) const;
172 : /// Virtual function! get gradients of a point
173 : /// CZAR and naive(ABF) have different gradient formulae
174 : virtual std::vector<double> getGradient(const std::vector<double> &pos,
175 : bool SkipCheck = false) const;
176 : /// Calculate dln(ρ)/dz, useful for CZAR
177 : /// This function may be moved to CZAR class in the future
178 : std::vector<double>
179 : getCountsLogDerivative(const std::vector<double> &pos) const;
180 : /// Write grad file
181 : // void writeGrad(std::string filename) const;
182 : /// Write 1D pmf file on one dimensional occasion
183 : void write1DPMF(std::string filename) const;
184 : /// Write count file
185 : // void writeCount(std::string filename) const;
186 : /// Write necessary output file in one function
187 : void writeAll(const std::string &filename) const;
188 : /// Miscellaneous getter functions, useful for merging windows
189 4 : std::vector<DRRAxis> getDimensions() const { return this->dimensions; }
190 2 : size_t getNumberOfDimension() const { return ndims; }
191 2 : size_t getSampleSize() const { return sampleSize; }
192 2 : std::vector<std::vector<double>> getTable() const { return table; }
193 : /// merge windows
194 : static std::vector<DRRAxis> merge(const std::vector<DRRAxis> &dA,
195 : const std::vector<DRRAxis> &dB);
196 : /// Get suffix
197 : std::string getSuffix() const { return suffix; }
198 : /// Set unit for .grad output
199 8 : void setOutputUnit(double unit) { outputunit = unit; }
200 : /// Destructor
201 30 : virtual ~DRRForceGrid() {}
202 :
203 : protected:
204 : /// The output suffix appended before .grad(.czar.grad) and
205 : /// .count(.czar.count)
206 : std::string suffix;
207 : /// Number of dimensions
208 : size_t ndims;
209 : /// Store each axes
210 : std::vector<DRRAxis> dimensions;
211 : /// Size of samples
212 : size_t sampleSize;
213 : /// Size of forces
214 : size_t forceSize;
215 : /// The header lines of .grad and .count files
216 : std::string headers;
217 : /// A table stores the middle points of all dimensions.
218 : /// For output in .grad and .count files
219 : std::vector<std::vector<double>> table;
220 : /// Store the average force of each bins
221 : std::vector<double> forces;
222 : /// Store counts of each bins
223 : std::vector<unsigned long int> samples;
224 : /// Only for 1D pmf output
225 : std::vector<double> endpoints;
226 : /// For (possibly) faster indexing
227 : std::vector<size_t> shifts;
228 : /// Output precision
229 : /// The abf_intergrate program has precision requirement.
230 : /// I test 9 and it just works.
231 : static const size_t OUTPUTPRECISION = 9;
232 : /// For set different output units
233 : double outputunit;
234 :
235 : /// Miscellaneous helper functions
236 : static size_t index1D(const DRRAxis &c, double x);
237 : void fillTable(const std::vector<std::vector<double>> &in);
238 :
239 : /// Boost serialization functions
240 : friend class boost::serialization::access;
241 : template <class Archive>
242 12 : void save(Archive &ar, const unsigned int version) const {
243 : // Don't save all members.
244 12 : ar << suffix;
245 12 : ar << dimensions;
246 12 : ar << forces;
247 12 : ar << samples;
248 12 : }
249 10 : template <class Archive> void load(Archive &ar, const unsigned int version) {
250 10 : ar >> suffix;
251 10 : ar >> dimensions;
252 10 : ar >> forces;
253 10 : ar >> samples;
254 : // Restore other members.
255 10 : ndims = dimensions.size();
256 10 : sampleSize = samples.size();
257 10 : forceSize = forces.size();
258 10 : std::stringstream ss;
259 10 : ss << "# " << ndims << '\n';
260 20 : std::vector<std::vector<double>> mp(ndims);
261 10 : shifts.resize(ndims, 0);
262 24 : for (size_t i = 0; i < ndims; ++i) {
263 14 : mp[i] = dimensions[i].getMiddlePoints();
264 28 : shifts[i] = std::accumulate(
265 14 : std::begin(dimensions), std::begin(dimensions) + i, size_t(1),
266 46 : [](size_t k, const DRRAxis &d) { return k * d.getBins(); });
267 14 : ss.precision(std::numeric_limits<double>::max_digits10);
268 28 : ss << std::fixed << "# " << dimensions[i].min << ' '
269 28 : << dimensions[i].getWidth() << ' ' << dimensions[i].nbins;
270 14 : if (dimensions[i].isPeriodic())
271 8 : ss << " 1" << '\n';
272 : else
273 6 : ss << " 0" << '\n';
274 : }
275 10 : fillTable(mp);
276 10 : headers = ss.str();
277 10 : outputunit = 1.0;
278 : // For 1D pmf
279 10 : if (ndims == 1) {
280 6 : endpoints.resize(dimensions[0].nbins + 1, 0);
281 6 : double ep = dimensions[0].min;
282 6 : double stride = dimensions[0].getWidth();
283 612 : for (auto it = std::begin(endpoints); it != std::end(endpoints); ++it) {
284 606 : (*it) = ep;
285 606 : ep += stride;
286 : }
287 10 : }
288 10 : }
289 : template <typename Archive>
290 22 : void serialize(Archive &ar, const unsigned int version) {
291 22 : boost::serialization::split_member(ar, *this, version);
292 22 : }
293 : };
294 :
295 6 : class ABF : public DRRForceGrid {
296 : public:
297 8 : ABF() {}
298 4 : ABF(const std::vector<DRRAxis> &p_dimensions, const std::string &p_suffix,
299 : bool initializeTable = true)
300 4 : : DRRForceGrid(p_dimensions, p_suffix, initializeTable) {}
301 : // Store the "instantaneous" spring force of a point and get ABF bias forces.
302 : bool store_getbias(const std::vector<double> &pos,
303 : const std::vector<double> &f, std::vector<double> &fbias,
304 : double fullsamples);
305 : static ABF mergewindow(const ABF &aWA, const ABF &aWB);
306 15 : ~ABF() {}
307 :
308 : private:
309 : // Boost serialization
310 : friend class boost::serialization::access;
311 : template <typename Archive>
312 11 : void serialize(Archive &ar, const unsigned int version) {
313 11 : ar &boost::serialization::base_object<DRRForceGrid>(*this);
314 11 : }
315 : };
316 :
317 6 : class CZAR : public DRRForceGrid {
318 : public:
319 8 : CZAR() : kbt(0) {}
320 4 : CZAR(const std::vector<DRRAxis> &p_dimensions, const std::string &p_suffix,
321 : double p_kbt, bool initializeTable = true)
322 4 : : DRRForceGrid(p_dimensions, p_suffix, initializeTable), kbt(p_kbt) {}
323 : std::vector<double> getGradient(const std::vector<double> &pos,
324 : bool SkipCheck = false) const;
325 1 : double getkbt() const { return kbt; }
326 : void setkbt(double p_kbt) { kbt = p_kbt; }
327 : static CZAR mergewindow(const CZAR &cWA, const CZAR &cWB);
328 15 : ~CZAR() {}
329 :
330 : private:
331 : double kbt;
332 : friend class boost::serialization::access;
333 : template <typename Archive>
334 11 : void serialize(Archive &ar, const unsigned int version) {
335 11 : ar &boost::serialization::base_object<DRRForceGrid>(*this);
336 11 : ar &kbt;
337 11 : }
338 : };
339 : }
340 : }
341 :
342 : #endif
343 : #endif
|
