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 : using std::vector;
45 : using std::string;
46 : using std::begin;
47 : using std::end;
48 :
49 : /// This class can store the minimum, maximum and bins of a dimension(axis).
50 : class DRRAxis {
51 : public:
52 : /// Default empty constructor
53 64 : DRRAxis() {
54 64 : min = max = 0.0;
55 64 : nbins = 0;
56 64 : periodic = false;
57 64 : domainMax = domainMin = 0.0;
58 64 : binWidth = 0.0;
59 : }
60 : /// Constructor using maximum value, minimum value and the number of bins(No
61 : /// pbc)
62 : DRRAxis(double l, double h, size_t n)
63 : : min(l), max(h), nbins(n), periodic(false), domainMax(0), domainMin(0),
64 : binWidth((max - min) / double(nbins)) {}
65 : /// PBC-aware constructor
66 : DRRAxis(double l, double h, size_t n, bool pbc, double dMax, double dMin)
67 4 : : min(l), max(h), nbins(n), periodic(pbc), domainMax(dMax),
68 4 : domainMin(dMin), binWidth((max - min) / double(nbins)) {}
69 : /// Set values
70 : void set(double l, double h, size_t n, bool pbc = false, double dmin = 0,
71 : double dmax = 0) {
72 38 : min = l;
73 38 : max = h;
74 38 : nbins = n;
75 38 : periodic = pbc;
76 38 : domainMax = dmax;
77 38 : domainMin = dmin;
78 38 : binWidth = (max - min) / nbins;
79 : }
80 : /// Set PBC data
81 : void setPeriodicity(double dmin, double dmax) {
82 16 : domainMax = dmax;
83 16 : domainMin = dmin;
84 16 : periodic = true;
85 : }
86 : /// Getters
87 : double getMin() const {
88 41 : return this->min;
89 : }
90 : double getMax() const {
91 40 : return this->max;
92 : }
93 : double getWidth() const {
94 31 : return binWidth;
95 : }
96 : double getDomainMax() const {
97 : return this->domainMax;
98 : }
99 : double getDomainMin() const {
100 : return this->domainMin;
101 : }
102 : size_t getBins() const {
103 : return this->nbins;
104 : }
105 :
106 : /// Check periodicity
107 : bool isPeriodic() const {
108 60 : return this->periodic;
109 : }
110 : /// Check real periodicity, i.e. the maximum == the domain maximum
111 : bool isRealPeriodic() const;
112 :
113 : /// Check whether x is in this axis
114 : bool isInBoundary(double x) const;
115 : /// Get an array of middle points of each bins
116 : vector<double> getMiddlePoints();
117 :
118 : /// Combine two axes if they share the same bin width.
119 : static DRRAxis merge(const DRRAxis &d1, const DRRAxis &d2);
120 :
121 : friend class DRRForceGrid;
122 :
123 : protected:
124 : double min; // Minimum value of the axis
125 : double max; // Maximum value of the axis
126 : size_t nbins; // Number of bins
127 : bool periodic; // Periodicity
128 : double domainMax; // Maximum value of the CV domain
129 : double domainMin; // Minimum value of the CV domain
130 : friend class boost::serialization::access;
131 : /// Use boost serialization
132 : template <typename Archive>
133 52 : void save(Archive &ar, const unsigned int version) const {
134 52 : ar &min;
135 52 : ar &max;
136 52 : ar &nbins;
137 52 : ar &periodic;
138 52 : ar &domainMax;
139 52 : ar &domainMin;
140 52 : }
141 : /// Split save and load. The bin width is calculated after initialization.
142 : template <typename Archive>
143 22 : void load(Archive &ar, const unsigned int version) {
144 22 : ar &min;
145 22 : ar &max;
146 22 : ar &nbins;
147 22 : ar &periodic;
148 22 : ar &domainMax;
149 22 : ar &domainMin;
150 22 : binWidth = (max - min) / double(nbins);
151 22 : }
152 : template <typename Archive>
153 : void serialize(Archive &ar, const unsigned int version) {
154 : boost::serialization::split_member(ar, *this, version);
155 : }
156 :
157 : private:
158 : double binWidth; // bin width
159 : };
160 :
161 : /// A class for collecting instantaneous forces, calculating average forces and
162 : /// build CV histogram.
163 : class DRRForceGrid {
164 : public:
165 : /// Empty constructor
166 : DRRForceGrid();
167 : /// "Real" constructor
168 : /// The 2D table vector is mainly used for print grid points in grad and count
169 : /// file.
170 : /// So when use binary output we can set initializeTable to false to save
171 : /// memory.
172 : explicit DRRForceGrid(const vector<DRRAxis> &p_dimensions,
173 : const string &p_suffix,
174 : bool initializeTable = true);
175 : /// Check whether a point is in this grid
176 : bool isInBoundary(const vector<double> &pos) const;
177 : // /// Get internal indices of a point
178 : // vector<size_t> index(const vector<double> &pos) const;
179 : /// Get internal counts address of a point
180 : size_t sampleAddress(const vector<double> &pos) const;
181 : /// Store instantaneous forces of a point
182 : /// nsamples > 1 is useful for merging windows
183 : bool store(const vector<double> &pos, const vector<double> &f,
184 : unsigned long int nsamples = 1);
185 : /// Get accumulated forces of a point
186 : vector<double>
187 : getAccumulatedForces(const vector<double> &pos) const;
188 : /// Get counts of a point
189 : unsigned long int getCount(const vector<double> &pos,
190 : bool SkipCheck = false) const;
191 : /// Virtual function! get gradients of a point
192 : /// CZAR and naive(ABF) have different gradient formulae
193 : virtual vector<double> getGradient(const vector<double> &pos,
194 : bool SkipCheck = false) const;
195 : /// Calculate divergence of the mean force field (experimental)
196 : double getDivergence(const vector<double> &pos) const;
197 : /// Calculate dln(ρ)/dz, useful for CZAR
198 : /// This function may be moved to CZAR class in the future
199 : vector<double>
200 : getCountsLogDerivative(const vector<double> &pos) const;
201 : /// Write grad file
202 : // void writeGrad(string filename) const;
203 : /// Write 1D pmf file on one dimensional occasion
204 : void write1DPMF(string filename, const string &fmt="%.9f") const;
205 : /// Write count file
206 : // void writeCount(string filename) const;
207 : /// Write necessary output file in one function (.grad and .count)
208 : void writeAll(const string &filename, const string &fmt="%.9f", bool addition = false) const;
209 : /// Output divergence (.div) (experimental)
210 : void writeDivergence(const string &filename, const string &fmt="%.9f") const;
211 : /// merge windows
212 : static vector<DRRAxis> merge(const vector<DRRAxis> &dA,
213 : const vector<DRRAxis> &dB);
214 : /// Get suffix
215 : string getSuffix() const {
216 : return suffix;
217 : }
218 : /// Set unit for .grad output
219 : void setOutputUnit(double unit) {
220 4 : outputunit = unit;
221 : }
222 : /// Destructor
223 228 : virtual ~DRRForceGrid() {}
224 :
225 : protected:
226 : /// The output suffix appended before .grad(.czar.grad) and
227 : /// .count(.czar.count)
228 : string suffix;
229 : /// Number of dimensions
230 : size_t ndims;
231 : /// Store each axes
232 : vector<DRRAxis> dimensions;
233 : /// Size of samples
234 : size_t sampleSize;
235 : /// The header lines of .grad and .count files
236 : string headers;
237 : /// A table stores the middle points of all dimensions.
238 : /// For output in .grad and .count files
239 : vector<vector<double>> table;
240 : /// Store the average force of each bins
241 : vector<double> forces;
242 : /// Store counts of each bins
243 : vector<unsigned long int> samples;
244 : /// Only for 1D pmf output
245 : vector<double> endpoints;
246 : /// For faster indexing
247 : /// shifts[0] = 1, shifts[n+1] = shifts[n] * dimensions[n].nbins
248 : vector<size_t> shifts;
249 : /// For set different output units
250 : double outputunit;
251 :
252 : /// Miscellaneous helper functions
253 : static size_t index1D(const DRRAxis &c, double x);
254 : void fillTable(const vector<vector<double>> &in);
255 :
256 : /// Boost serialization functions
257 : friend class boost::serialization::access;
258 : template <class Archive>
259 38 : void save(Archive &ar, const unsigned int version) const {
260 : // Don't save all members.
261 38 : ar << suffix;
262 38 : ar << dimensions;
263 38 : ar << forces;
264 38 : ar << samples;
265 38 : }
266 16 : template <class Archive> void load(Archive &ar, const unsigned int version) {
267 16 : ar >> suffix;
268 16 : ar >> dimensions;
269 16 : ar >> forces;
270 16 : ar >> samples;
271 : // Restore other members.
272 16 : ndims = dimensions.size();
273 16 : sampleSize = samples.size();
274 16 : std::stringstream ss;
275 16 : ss << "# " << ndims << '\n';
276 16 : vector<vector<double>> mp(ndims);
277 16 : shifts.resize(ndims, 0);
278 16 : shifts[0] = 1;
279 38 : for (size_t i = 0; i < ndims; ++i) {
280 22 : mp[i] = dimensions[i].getMiddlePoints();
281 22 : if (i > 0) {
282 6 : shifts[i] = shifts[i - 1] * dimensions[i - 1].nbins;
283 : }
284 : ss.precision(std::numeric_limits<double>::max_digits10);
285 44 : ss << std::fixed << "# " << dimensions[i].min << ' '
286 44 : << dimensions[i].binWidth << ' ' << dimensions[i].nbins;
287 22 : if (dimensions[i].isPeriodic()) {
288 24 : ss << " 1" << '\n';
289 : } else {
290 20 : ss << " 0" << '\n';
291 : }
292 : }
293 16 : fillTable(mp);
294 16 : headers = ss.str();
295 16 : outputunit = 1.0;
296 : // For 1D pmf
297 16 : if (ndims == 1) {
298 10 : endpoints.resize(dimensions[0].nbins + 1, 0);
299 10 : double ep = dimensions[0].min;
300 10 : double stride = dimensions[0].binWidth;
301 1020 : for (auto it = begin(endpoints); it != end(endpoints); ++it) {
302 1010 : (*it) = ep;
303 1010 : ep += stride;
304 : }
305 : }
306 16 : }
307 : template <typename Archive>
308 : void serialize(Archive &ar, const unsigned int version) {
309 : boost::serialization::split_member(ar, *this, version);
310 : }
311 : };
312 :
313 : class ABF : public DRRForceGrid {
314 : public:
315 16 : ABF() {}
316 2 : ABF(const vector<DRRAxis> &p_dimensions, const string &p_suffix,
317 : double fullSamples = 500.0, double maxFactor = 1.0,
318 : bool initializeTable = true)
319 2 : : DRRForceGrid(p_dimensions, p_suffix, initializeTable),
320 2 : mFullSamples(fullSamples), mMaxFactors(p_dimensions.size(), maxFactor) {}
321 11 : ABF(const vector<DRRAxis> &p_dimensions, const string &p_suffix,
322 : double fullSamples, const vector<double>& maxFactors,
323 : bool initializeTable = true)
324 11 : : DRRForceGrid(p_dimensions, p_suffix, initializeTable),
325 11 : mFullSamples(fullSamples), mMaxFactors(maxFactors) {}
326 : // Provide a setter for ABF parametres (fullsamples, maxfactor)
327 : void setParameters(double fullSamples, const vector<double>& maxFactors) {
328 1 : mFullSamples = fullSamples;
329 1 : mMaxFactors = maxFactors;
330 1 : }
331 : // Store the "instantaneous" spring force of a point and get ABF bias forces.
332 : bool getbias(const vector<double> &pos,
333 : const vector<double> &f,
334 : vector<double> &fbias) const;
335 : void store_force(const vector<double> &pos,
336 : const vector<double> &f);
337 : static ABF mergewindow(const ABF &aWA, const ABF &aWB);
338 38 : ~ABF() {}
339 :
340 : private:
341 : // Parametres for calculate bias force
342 : double mFullSamples;
343 : vector<double> mMaxFactors;
344 : // Boost serialization
345 : friend class boost::serialization::access;
346 : template <typename Archive>
347 : void serialize(Archive &ar, const unsigned int version) {
348 27 : ar &boost::serialization::base_object<DRRForceGrid>(*this);
349 : }
350 : };
351 :
352 28 : class CZAR : public DRRForceGrid {
353 : public:
354 16 : CZAR() : kbt(0) {}
355 : CZAR(const vector<DRRAxis> &p_dimensions, const string &p_suffix,
356 : double p_kbt, bool initializeTable = true)
357 13 : : DRRForceGrid(p_dimensions, p_suffix, initializeTable), kbt(p_kbt) {}
358 : vector<double> getGradient(const vector<double> &pos,
359 : bool SkipCheck = false) const;
360 : double getkbt() const {
361 : return kbt;
362 : }
363 : void setkbt(double p_kbt) {
364 : kbt = p_kbt;
365 : }
366 : static CZAR mergewindow(const CZAR &cWA, const CZAR &cWB);
367 : void writeZCountZGrad(const string &filename, bool addition = false) const;
368 23 : ~CZAR() {}
369 :
370 : private:
371 : double kbt;
372 : friend class boost::serialization::access;
373 : template <typename Archive>
374 27 : void serialize(Archive &ar, const unsigned int version) {
375 27 : ar &boost::serialization::base_object<DRRForceGrid>(*this);
376 27 : ar &kbt;
377 27 : }
378 : };
379 : }
380 : }
381 :
382 : #endif
383 : #endif
|
