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