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 52 : DRRAxis() {
54 52 : min = max = 0.0;
55 52 : nbins = 0;
56 52 : periodic = false;
57 52 : domainMax = domainMin = 0.0;
58 52 : 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 2 : : min(l), max(h), nbins(n), periodic(pbc), domainMax(dMax),
68 2 : 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 32 : min = l;
73 32 : max = h;
74 32 : nbins = n;
75 32 : periodic = pbc;
76 32 : domainMax = dmax;
77 32 : domainMin = dmin;
78 32 : binWidth = (max - min) / nbins;
79 : }
80 : /// Set PBC data
81 : void setPeriodicity(double dmin, double dmax) {
82 28 : domainMax = dmax;
83 28 : domainMin = dmin;
84 28 : periodic = true;
85 : }
86 : /// Getters
87 14 : double getMin() const { return this->min; }
88 14 : double getMax() const { return this->max; }
89 28 : 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 48 : 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 36 : void save(Archive &ar, const unsigned int version) const {
120 36 : ar &min;
121 36 : ar &max;
122 36 : ar &nbins;
123 36 : ar &periodic;
124 36 : ar &domainMax;
125 36 : ar &domainMin;
126 36 : }
127 : /// Split save and load. The bin width is calculated after initialization.
128 : template <typename Archive>
129 18 : void load(Archive &ar, const unsigned int version) {
130 18 : ar &min;
131 18 : ar &max;
132 18 : ar &nbins;
133 18 : ar &periodic;
134 18 : ar &domainMax;
135 18 : ar &domainMin;
136 18 : binWidth = (max - min) / double(nbins);
137 18 : }
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 72 : 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;
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;
195 : /// Output divergence (.div) (experimental)
196 : void writeDivergence(const string &filename) 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 5 : void setOutputUnit(double unit) { outputunit = unit; }
204 : /// Destructor
205 168 : 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 24 : void save(Archive &ar, const unsigned int version) const {
242 : // Don't save all members.
243 24 : ar << suffix;
244 24 : ar << dimensions;
245 24 : ar << forces;
246 24 : ar << samples;
247 24 : }
248 12 : template <class Archive> void load(Archive &ar, const unsigned int version) {
249 12 : ar >> suffix;
250 12 : ar >> dimensions;
251 12 : ar >> forces;
252 12 : ar >> samples;
253 : // Restore other members.
254 12 : ndims = dimensions.size();
255 12 : sampleSize = samples.size();
256 24 : std::stringstream ss;
257 12 : ss << "# " << ndims << '\n';
258 24 : vector<vector<double>> mp(ndims);
259 12 : shifts.resize(ndims, 0);
260 12 : shifts[0] = 1;
261 48 : for (size_t i = 0; i < ndims; ++i) {
262 36 : mp[i] = dimensions[i].getMiddlePoints();
263 18 : if (i > 0) {
264 24 : shifts[i] = shifts[i - 1] * dimensions[i - 1].nbins;
265 : }
266 : ss.precision(std::numeric_limits<double>::max_digits10);
267 18 : ss << std::fixed << "# " << dimensions[i].min << ' '
268 36 : << dimensions[i].binWidth << ' ' << dimensions[i].nbins;
269 18 : if (dimensions[i].isPeriodic())
270 : ss << " 1" << '\n';
271 : else
272 : ss << " 0" << '\n';
273 : }
274 12 : fillTable(mp);
275 24 : headers = ss.str();
276 12 : outputunit = 1.0;
277 : // For 1D pmf
278 12 : if (ndims == 1) {
279 12 : endpoints.resize(dimensions[0].nbins + 1, 0);
280 6 : double ep = dimensions[0].min;
281 6 : double stride = dimensions[0].binWidth;
282 612 : for (auto it = begin(endpoints); it != end(endpoints); ++it) {
283 606 : (*it) = ep;
284 606 : ep += stride;
285 : }
286 : }
287 12 : }
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 12 : class ABF : public DRRForceGrid {
295 : public:
296 15 : ABF() {}
297 : ABF(const vector<DRRAxis> &p_dimensions, const string &p_suffix,
298 : double fullSamples = 500.0, double maxFactor = 1.0,
299 : bool initializeTable = true)
300 9 : : DRRForceGrid(p_dimensions, p_suffix, initializeTable),
301 10 : mFullSamples(fullSamples), mMaxFactor(maxFactor) {}
302 : // Provide a setter for ABF parametres (fullsamples, maxfactor)
303 : void setParameters(double fullSamples, double maxFactor) {
304 1 : mFullSamples = fullSamples;
305 1 : mMaxFactor = maxFactor;
306 : }
307 : // Store the "instantaneous" spring force of a point and get ABF bias forces.
308 : bool store_getbias(const vector<double> &pos,
309 : const vector<double> &f,
310 : vector<double> &fbias);
311 : static ABF mergewindow(const ABF &aWA, const ABF &aWB);
312 28 : ~ABF() {}
313 :
314 : private:
315 : // Parametres for calculate bias force
316 : double mFullSamples;
317 : double mMaxFactor;
318 : // Boost serialization
319 : friend class boost::serialization::access;
320 : template <typename Archive>
321 18 : void serialize(Archive &ar, const unsigned int version) {
322 18 : ar &boost::serialization::base_object<DRRForceGrid>(*this);
323 18 : }
324 : };
325 :
326 12 : class CZAR : public DRRForceGrid {
327 : public:
328 15 : CZAR() : kbt(0) {}
329 : CZAR(const vector<DRRAxis> &p_dimensions, const string &p_suffix,
330 : double p_kbt, bool initializeTable = true)
331 10 : : DRRForceGrid(p_dimensions, p_suffix, initializeTable), kbt(p_kbt) {}
332 : vector<double> getGradient(const vector<double> &pos,
333 : bool SkipCheck = false) const;
334 : double getkbt() const { return kbt; }
335 : void setkbt(double p_kbt) { kbt = p_kbt; }
336 : static CZAR mergewindow(const CZAR &cWA, const CZAR &cWB);
337 28 : ~CZAR() {}
338 :
339 : private:
340 : double kbt;
341 : friend class boost::serialization::access;
342 : template <typename Archive>
343 18 : void serialize(Archive &ar, const unsigned int version) {
344 18 : ar &boost::serialization::base_object<DRRForceGrid>(*this);
345 18 : ar &kbt;
346 18 : }
347 : };
348 : }
349 : }
350 :
351 : #endif
352 : #endif
|
