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 : #ifdef __PLUMED_HAS_BOOST_SERIALIZATION
19 : #include "DRR.h"
20 :
21 : namespace PLMD {
22 : namespace drr {
23 :
24 0 : bool DRRAxis::isInBoundary(double x) const {
25 0 : if (x < min || x > max)
26 : return false;
27 : else
28 0 : return true;
29 : }
30 :
31 736 : bool DRRAxis::isRealPeriodic() const {
32 736 : if (periodic == true) {
33 2160 : if(std::abs(domainMax - max) < binWidth && std::abs(domainMin - min) < binWidth) {
34 : return true;
35 : }
36 : else {
37 : return false;
38 : }
39 : }
40 : else {
41 : return false;
42 : }
43 : }
44 :
45 2 : DRRAxis DRRAxis::merge(const DRRAxis &d1, const DRRAxis &d2) {
46 6 : const double newmin = std::min(d1.getMin(), d2.getMin());
47 6 : const double newmax = std::max(d1.getMax(), d2.getMax());
48 : const double newWidth = d1.getWidth();
49 2 : const size_t newbins = size_t(std::nearbyint((newmax - newmin) / newWidth));
50 : const bool newpbc = d1.isPeriodic();
51 6 : const double newdmin = std::min(d1.getDomainMin(), d2.getDomainMin());
52 6 : const double newdmax = std::max(d1.getDomainMax(), d2.getDomainMax());
53 : DRRAxis result(newmin, newmax, newbins, newpbc, newdmin, newdmax);
54 2 : return result;
55 : }
56 :
57 24 : std::vector<double> DRRAxis::getMiddlePoints() {
58 24 : std::vector<double> result(nbins, 0);
59 24 : const double width = getWidth();
60 24 : double temp = min - width / 2;
61 48 : std::generate(std::begin(result), std::end(result), [&]() {
62 1772 : temp += width;
63 : return temp;
64 : });
65 24 : return result;
66 : }
67 :
68 393293 : size_t DRRForceGrid::index1D(const DRRAxis &c, double x) {
69 : #ifdef DEBUG_DRR
70 : if (x < c.min || x > c.max) {
71 : std::cerr << "This is a bug!" << '\n';
72 : std::cerr << "CV should be larger than minimal value or smaller than the "
73 : "maximum value of dimension."
74 : << '\n';
75 : std::cerr << "min = " << c.min << '\n';
76 : std::cerr << "max = " << c.max << '\n';
77 : std::cerr << "x = " << x << std::endl;
78 : std::abort();
79 : }
80 : #endif
81 393293 : size_t idx = size_t(std::floor((x - c.min) / c.binWidth));
82 393293 : idx = (idx == c.nbins) ? (c.nbins - 1) : idx;
83 393293 : return idx;
84 : }
85 :
86 16 : void DRRForceGrid::fillTable(const std::vector<std::vector<double>> &in) {
87 32 : table.resize(ndims, std::vector<double>(sampleSize, 0));
88 36 : for (size_t i = 0; i < ndims; ++i) {
89 : size_t repeatAll = 1, repeatOne = 1;
90 24 : for (size_t j = i + 1; j < ndims; ++j)
91 4 : repeatOne *= in[j].size();
92 28 : for (size_t j = 0; j < i; ++j)
93 4 : repeatAll *= in[j].size();
94 : size_t in_i_sz = in[i].size();
95 3540 : for (size_t l = 0; l < in_i_sz; ++l)
96 1760 : std::fill_n(std::begin(table[i]) + l * repeatOne, repeatOne, in[i][l]);
97 362 : for (size_t k = 0; k < repeatAll - 1; ++k)
98 : std::copy_n(std::begin(table[i]), repeatOne * in_i_sz,
99 362 : std::begin(table[i]) + repeatOne * in_i_sz * (k + 1));
100 : }
101 16 : }
102 :
103 16 : DRRForceGrid::DRRForceGrid()
104 : : suffix(""), ndims(0), dimensions(0), sampleSize(0), forceSize(0),
105 16 : headers(""), table(0), forces(0), samples(0), endpoints(0), shifts(0), outputunit(1.0) {}
106 :
107 8 : DRRForceGrid::DRRForceGrid(const std::vector<DRRAxis> &p_dimensions,
108 8 : const std::string &p_suffix, bool initializeTable)
109 56 : : suffix(p_suffix), ndims(p_dimensions.size()), dimensions(p_dimensions) {
110 8 : sampleSize = 1;
111 16 : std::vector<std::vector<double>> mp(ndims);
112 16 : std::stringstream ss;
113 8 : ss << "# " << ndims << '\n';
114 8 : shifts.resize(ndims, 0);
115 28 : for (size_t i = 0; i < ndims; ++i) {
116 10 : sampleSize = dimensions[i].nbins * sampleSize;
117 20 : mp[i] = dimensions[i].getMiddlePoints();
118 10 : shifts[i] = std::accumulate(
119 : std::begin(dimensions), std::begin(dimensions) + i, size_t(1),
120 2 : [](size_t k, const DRRAxis &d) { return k * d.getBins(); });
121 : ss.precision(std::numeric_limits<double>::max_digits10);
122 10 : ss << std::fixed << "# " << dimensions[i].min << ' '
123 10 : << dimensions[i].getWidth() << ' ' << dimensions[i].nbins;
124 10 : if (dimensions[i].isPeriodic())
125 : ss << " 1" << '\n';
126 : else
127 : ss << " 0" << '\n';
128 : }
129 16 : headers = ss.str();
130 8 : if (initializeTable)
131 6 : fillTable(mp);
132 8 : forceSize = sampleSize * ndims;
133 8 : forces.resize(forceSize, 0.0);
134 8 : samples.resize(sampleSize, 0);
135 8 : outputunit = 1.0;
136 : // For 1D pmf
137 8 : if (ndims == 1) {
138 12 : endpoints.resize(dimensions[0].nbins + 1, 0);
139 6 : double ep = dimensions[0].min;
140 : double stride = dimensions[0].getWidth();
141 440 : for (auto &i : endpoints) {
142 434 : i = ep;
143 434 : ep += stride;
144 : }
145 : }
146 8 : }
147 :
148 2463 : bool DRRForceGrid::isInBoundary(const std::vector<double> &pos) const {
149 : bool result = true;
150 5885 : for (size_t i = 0; i < ndims; ++i) {
151 5022 : if (pos[i] < dimensions[i].min || pos[i] > dimensions[i].max)
152 : return false;
153 : }
154 : return result;
155 : }
156 :
157 165817 : size_t DRRForceGrid::sampleAddress(const std::vector<double> &pos) const {
158 : size_t saddr = 0;
159 821547 : for (size_t i = 0; i < ndims; ++i) {
160 983595 : saddr += shifts[i] * index1D(dimensions[i], pos[i]);
161 : }
162 165817 : return saddr;
163 : }
164 :
165 829 : bool DRRForceGrid::store(const std::vector<double> &pos,
166 : const std::vector<double> &f,
167 : unsigned long int nsamples) {
168 829 : if (isInBoundary(pos)) {
169 829 : if (nsamples == 0)
170 : return true;
171 429 : const size_t baseaddr = sampleAddress(pos) * ndims;
172 858 : samples[baseaddr / ndims] += nsamples;
173 429 : auto it_fa = std::begin(forces) + baseaddr;
174 429 : std::transform(std::begin(f), std::end(f), it_fa, it_fa,
175 : std::plus<double>());
176 429 : return true;
177 : } else {
178 : return false;
179 : }
180 : }
181 :
182 2 : std::vector<DRRAxis> DRRForceGrid::merge(const std::vector<DRRAxis> &dA,
183 : const std::vector<DRRAxis> &dB) {
184 2 : std::vector<DRRAxis> dR(dA.size());
185 2 : std::transform(std::begin(dA), std::end(dA), std::begin(dB), std::begin(dR),
186 : DRRAxis::merge);
187 2 : return dR;
188 : }
189 :
190 : std::vector<double>
191 800 : DRRForceGrid::getAccumulatedForces(const std::vector<double> &pos) const {
192 800 : std::vector<double> result(ndims, 0);
193 800 : if (!isInBoundary(pos))
194 : return result;
195 400 : const size_t baseaddr = sampleAddress(pos) * ndims;
196 : std::copy(std::begin(forces) + baseaddr,
197 : std::begin(forces) + baseaddr + ndims, std::begin(result));
198 400 : return result;
199 : }
200 :
201 66242 : unsigned long int DRRForceGrid::getCount(const std::vector<double> &pos,
202 : bool SkipCheck) const {
203 66242 : if (!SkipCheck) {
204 800 : if (!isInBoundary(pos)) {
205 : return 0;
206 : }
207 : }
208 131684 : return samples[sampleAddress(pos)];
209 : }
210 :
211 33056 : std::vector<double> DRRForceGrid::getGradient(const std::vector<double> &pos,
212 : bool SkipCheck) const {
213 33056 : std::vector<double> result(ndims, 0);
214 33056 : if (!SkipCheck) {
215 0 : if (!isInBoundary(pos)) {
216 : return result;
217 : }
218 : }
219 33056 : const size_t baseaddr = sampleAddress(pos) * ndims;
220 66112 : if (samples[baseaddr / ndims] == 0)
221 : return result;
222 : auto it_fa = std::begin(forces) + baseaddr;
223 33005 : std::transform(it_fa, it_fa + ndims, std::begin(result), [&](double fa) {
224 130804 : return (-1.0) * fa / samples[baseaddr / ndims];
225 98407 : });
226 33005 : return result;
227 : }
228 :
229 : std::vector<double>
230 33028 : DRRForceGrid::getCountsLogDerivative(const std::vector<double> &pos) const {
231 33028 : const size_t addr = sampleAddress(pos);
232 33028 : const unsigned long int count_this = samples[addr];
233 33028 : std::vector<double> result(ndims, 0);
234 163884 : for (size_t i = 0; i < ndims; ++i) {
235 : const double binWidth = dimensions[i].getWidth();
236 : const size_t addr_first =
237 130856 : addr - shifts[i] * index1D(dimensions[i], pos[i]) + 0;
238 130856 : const size_t addr_last = addr_first + shifts[i] * (dimensions[i].nbins - 1);
239 65428 : if (addr == addr_first) {
240 368 : if (dimensions[i].isRealPeriodic() == true) {
241 360 : const unsigned long int &count_next = samples[addr + shifts[i]];
242 : const unsigned long int &count_prev = samples[addr_last];
243 360 : if (count_next != 0 && count_prev != 0)
244 360 : result[i] =
245 1080 : (std::log(count_next) - std::log(count_prev)) / (2 * binWidth);
246 : } else {
247 8 : const unsigned long int &count_next = samples[addr + shifts[i]];
248 8 : const unsigned long int &count_next2 = samples[addr + shifts[i] * 2];
249 8 : if (count_next != 0 && count_this != 0 && count_next2 != 0)
250 4 : result[i] =
251 12 : (std::log(count_next2) * (-1.0) + std::log(count_next) * 4.0 -
252 8 : std::log(count_this) * 3.0) /
253 4 : (2.0 * binWidth);
254 : }
255 65060 : } else if (addr == addr_last) {
256 368 : if (dimensions[i].isRealPeriodic() == true) {
257 360 : const unsigned long int &count_prev = samples[addr - shifts[i]];
258 : const unsigned long int &count_next = samples[addr_first];
259 360 : if (count_next != 0 && count_prev != 0)
260 360 : result[i] =
261 1080 : (std::log(count_next) - std::log(count_prev)) / (2 * binWidth);
262 : } else {
263 8 : const unsigned long int &count_prev = samples[addr - shifts[i]];
264 8 : const unsigned long int &count_prev2 = samples[addr - shifts[i] * 2];
265 8 : if (count_prev != 0 && count_this != 0 && count_prev2 != 0)
266 16 : result[i] = (std::log(count_this) * 3.0 - std::log(count_prev) * 4.0 +
267 8 : std::log(count_prev2)) /
268 4 : (2.0 * binWidth);
269 : }
270 : } else {
271 64692 : const unsigned long int &count_prev = samples[addr - shifts[i]];
272 64692 : const unsigned long int &count_next = samples[addr + shifts[i]];
273 64692 : if (count_next != 0 && count_prev != 0)
274 64664 : result[i] =
275 193992 : (std::log(count_next) - std::log(count_prev)) / (2 * binWidth);
276 : }
277 : }
278 33028 : return result;
279 : }
280 :
281 : // Write the gradients to a .grad file.
282 : // void DRRForceGrid::writeGrad(std::string filename) const {
283 : // std::stringstream ssg;
284 : // std::fstream fsgrad;
285 : // filename += suffix + ".grad";
286 : // ssg << headers;
287 : // ssg << std::left << std::fixed << std::setprecision(OUTPUTPRECISION);
288 : // for (size_t i = 0; i < sampleSize; ++i) {
289 : // std::vector<double> pos(ndims, 0);
290 : // for (size_t j = 0; j < ndims; ++j) {
291 : // pos[j] = table[j][i];
292 : // ssg << ' ' << table[j][i];
293 : // }
294 : // const std::vector<double> f = getGradient(pos, true);
295 : // for (const auto &i : f)
296 : // ssg << ' ' << i;
297 : // ssg << '\n';
298 : // }
299 : // fsgrad.open(filename.c_str(), std::ios_base::out);
300 : // fsgrad.write(ssg.str().c_str(), ssg.str().length());
301 : // fsgrad.close();
302 : // }
303 :
304 10 : void DRRForceGrid::write1DPMF(std::string filename) const {
305 20 : filename += suffix + ".pmf";
306 : FILE *ppmf;
307 10 : ppmf = fopen(filename.c_str(), "w");
308 : const double w = dimensions[0].getWidth();
309 : double pmf = 0;
310 10 : fprintf(ppmf, "%.9f %.9f\n", endpoints[0], pmf);
311 1294 : for (size_t i = 0; i < dimensions[0].nbins; ++i) {
312 642 : std::vector<double> pos(1, 0);
313 642 : pos[0] = table[0][i];
314 642 : const std::vector<double> f = getGradient(pos, true);
315 642 : pmf += f[0] * w;
316 1284 : fprintf(ppmf, "%.9f %.9f\n", endpoints[i + 1], pmf);
317 : }
318 10 : fclose(ppmf);
319 10 : }
320 :
321 : // Write the gradients to a .count file.
322 : // void DRRForceGrid::writeCount(std::string filename) const {
323 : // std::stringstream ssc;
324 : // std::fstream fscount;
325 : // filename += suffix + ".count";
326 : // ssc << headers;
327 : // ssc << std::left << std::fixed << std::setprecision(OUTPUTPRECISION);
328 : // std::vector<double> pos(ndims, 0);
329 : // for (size_t i = 0; i < sampleSize; ++i) {
330 : // for (size_t j = 0; j < ndims; ++j) {
331 : // pos[j] = table[j][i];
332 : // ssc << ' ' << table[j][i];
333 : // }
334 : // ssc << ' ' << getCount(pos, true) << '\n';
335 : // }
336 : // fscount.open(filename.c_str(), std::ios_base::out);
337 : // fscount.write(ssc.str().c_str(), ssc.str().length());
338 : // fscount.close();
339 : // }
340 :
341 12 : void DRRForceGrid::writeAll(const std::string &filename) const {
342 24 : std::string countname = filename + suffix + ".count";
343 24 : std::string gradname = filename + suffix + ".grad";
344 12 : std::vector<double> pos(ndims, 0);
345 : FILE *pGrad, *pCount;
346 12 : pGrad = fopen(gradname.c_str(), "w");
347 12 : pCount = fopen(countname.c_str(), "w");
348 : char *buffer1, *buffer2;
349 12 : buffer1 = (char*)malloc((sizeof(double))*sampleSize*ndims);
350 12 : buffer2 = (char*)malloc((sizeof(double))*sampleSize*ndims);
351 12 : setvbuf(pGrad, buffer1, _IOFBF, (sizeof(double))*sampleSize*ndims);
352 12 : setvbuf(pCount, buffer2, _IOFBF, (sizeof(double))*sampleSize*ndims);
353 12 : fwrite(headers.c_str(), sizeof(char), strlen(headers.c_str()), pGrad);
354 12 : fwrite(headers.c_str(), sizeof(char), strlen(headers.c_str()), pCount);
355 130896 : for (size_t i = 0; i < sampleSize; ++i) {
356 325926 : for (size_t j = 0; j < ndims; ++j) {
357 130242 : pos[j] = table[j][i];
358 130242 : fprintf(pGrad, " %.9f", table[j][i]);
359 130242 : fprintf(pCount, " %.9f", table[j][i]);
360 : }
361 65442 : fprintf(pCount, " %lu\n", getCount(pos, true));
362 65442 : std::vector<double> f = getGradient(pos, true);
363 325926 : for (size_t j = 0; j < ndims; ++j) {
364 130242 : fprintf(pGrad, " %.9f", (f[j] / outputunit));
365 : }
366 : fprintf(pGrad, "\n");
367 : }
368 12 : fclose(pGrad);
369 12 : fclose(pCount);
370 12 : free(buffer1);
371 12 : free(buffer2);
372 12 : if (ndims == 1) {
373 20 : write1DPMF(filename);
374 : }
375 12 : }
376 :
377 34 : bool ABF::store_getbias(const std::vector<double> &pos,
378 : const std::vector<double> &f,
379 : std::vector<double> &fbias, double fullsamples) {
380 34 : if (!isInBoundary(pos)) {
381 0 : std::fill(std::begin(fbias), std::end(fbias), 0.0);
382 0 : return false;
383 : }
384 34 : const size_t baseaddr = sampleAddress(pos);
385 : unsigned long int &count = samples[baseaddr];
386 34 : ++count;
387 34 : double factor = 2 * (static_cast<double>(count)) / fullsamples - 1;
388 34 : factor = factor < 0 ? 0 : factor > 1 ? 1 : factor; // Clamp to [0,1]
389 34 : auto it_fa = std::begin(forces) + baseaddr * ndims;
390 : auto it_fb = std::begin(fbias);
391 : auto it_f = std::begin(f);
392 : do {
393 58 : (*it_fa) += (*it_f); // Accumulate instantaneous force
394 58 : (*it_fb) =
395 58 : factor * (*it_fa) * (-1.0) / static_cast<double>(count); // Calculate bias force
396 : ++it_fa;
397 : ++it_fb;
398 : ++it_f;
399 58 : } while (it_f != std::end(f));
400 :
401 : return true;
402 : }
403 :
404 1 : ABF ABF::mergewindow(const ABF &aWA, const ABF &aWB) {
405 : const std::vector<DRRAxis> dA = aWA.getDimensions();
406 : const std::vector<DRRAxis> dB = aWB.getDimensions();
407 1 : const std::vector<DRRAxis> dR = merge(dA, dB);
408 1 : const std::string suffix = ".abf";
409 : ABF result(dR, suffix);
410 1 : const std::vector<std::vector<double>> table = result.getTable();
411 : const size_t nrows = result.getSampleSize();
412 : const size_t ncols = result.getNumberOfDimension();
413 1 : std::vector<double> pos(ncols, 0);
414 401 : for (size_t i = 0; i < nrows; ++i) {
415 600 : for (size_t j = 0; j < ncols; ++j) {
416 200 : pos[j] = table[j][i];
417 : }
418 200 : const unsigned long int countA = aWA.getCount(pos, false);
419 200 : const unsigned long int countB = aWB.getCount(pos, false);
420 200 : const std::vector<double> aForceA = aWA.getAccumulatedForces(pos);
421 200 : const std::vector<double> aForceB = aWB.getAccumulatedForces(pos);
422 200 : result.store(pos, aForceA, countA);
423 200 : result.store(pos, aForceB, countB);
424 : }
425 1 : return result;
426 : }
427 :
428 33028 : std::vector<double> CZAR::getGradient(const std::vector<double> &pos,
429 : bool SkipCheck) const {
430 33028 : std::vector<double> result(ndims, 0);
431 33028 : if (!SkipCheck) {
432 0 : if (!isInBoundary(pos)) {
433 : return result;
434 : }
435 : }
436 33028 : if (kbt <= std::numeric_limits<double>::epsilon()) {
437 : std::cerr << "ERROR! The kbt shouldn't be zero when use CZAR estimator. "
438 : << '\n';
439 0 : std::abort();
440 : }
441 33028 : const size_t baseaddr = sampleAddress(pos) * ndims;
442 33028 : const std::vector<double> log_deriv(getCountsLogDerivative(pos));
443 66056 : if (samples[baseaddr / ndims] == 0)
444 : return result;
445 : auto it_fa = std::begin(forces) + baseaddr;
446 33002 : std::transform(it_fa, it_fa + ndims, std::begin(log_deriv),
447 : std::begin(result), [&](double fa, double ld) {
448 130800 : return fa * (-1.0) / samples[baseaddr / ndims] - kbt * ld;
449 98402 : });
450 33002 : return result;
451 : }
452 :
453 1 : CZAR CZAR::mergewindow(const CZAR &cWA, const CZAR &cWB) {
454 : const std::vector<DRRAxis> dA = cWA.getDimensions();
455 : const std::vector<DRRAxis> dB = cWB.getDimensions();
456 1 : const std::vector<DRRAxis> dR = merge(dA, dB);
457 : const double newkbt = cWA.getkbt();
458 1 : const std::string suffix = ".czar";
459 : CZAR result(dR, suffix, newkbt);
460 1 : const std::vector<std::vector<double>> table = result.getTable();
461 : const size_t nrows = result.getSampleSize();
462 : const size_t ncols = result.getNumberOfDimension();
463 1 : std::vector<double> pos(ncols, 0);
464 401 : for (size_t i = 0; i < nrows; ++i) {
465 600 : for (size_t j = 0; j < ncols; ++j) {
466 200 : pos[j] = table[j][i];
467 : }
468 200 : const unsigned long int countA = cWA.getCount(pos);
469 200 : const unsigned long int countB = cWB.getCount(pos);
470 200 : const std::vector<double> aForceA = cWA.getAccumulatedForces(pos);
471 200 : const std::vector<double> aForceB = cWB.getAccumulatedForces(pos);
472 200 : result.store(pos, aForceA, countA);
473 200 : result.store(pos, aForceB, countB);
474 : }
475 1 : return result;
476 : }
477 :
478 : }
479 4839 : }
480 :
481 : #endif
|
