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 : using std::vector;
25 : using std::string;
26 : using std::begin;
27 : using std::end;
28 :
29 0 : bool DRRAxis::isInBoundary(double x) const {
30 0 : if (x < min || x > max)
31 : return false;
32 : else
33 0 : return true;
34 : }
35 :
36 5820 : bool DRRAxis::isRealPeriodic() const {
37 5820 : if (periodic == true) {
38 17388 : if (std::abs(domainMax - max) < binWidth &&
39 5796 : std::abs(domainMin - min) < binWidth) {
40 : return true;
41 : } else {
42 : return false;
43 : }
44 : } else {
45 : return false;
46 : }
47 : }
48 :
49 2 : DRRAxis DRRAxis::merge(const DRRAxis &d1, const DRRAxis &d2) {
50 4 : const double newmin = std::min(d1.min, d2.min);
51 4 : const double newmax = std::max(d1.max, d2.max);
52 2 : const double newWidth = d1.binWidth;
53 2 : const size_t newbins = size_t(std::nearbyint((newmax - newmin) / newWidth));
54 2 : const bool newpbc = d1.periodic;
55 4 : const double newdmin = std::min(d1.domainMin, d2.domainMin);
56 4 : const double newdmax = std::max(d1.domainMax, d2.domainMax);
57 : DRRAxis result(newmin, newmax, newbins, newpbc, newdmin, newdmax);
58 2 : return result;
59 : }
60 :
61 48 : vector<double> DRRAxis::getMiddlePoints() {
62 48 : vector<double> result(nbins, 0);
63 48 : const double width = binWidth;
64 48 : double temp = min - width / 2;
65 2652 : std::generate(begin(result), end(result), [&temp, &width]() {
66 5112 : temp += width;
67 : return temp;
68 : });
69 48 : return result;
70 : }
71 :
72 2337960 : size_t DRRForceGrid::index1D(const DRRAxis &c, double x) {
73 2337960 : size_t idx = size_t(std::floor((x - c.min) / c.binWidth));
74 2337960 : idx = (idx == c.nbins) ? (c.nbins - 1) : idx;
75 2337960 : return idx;
76 : }
77 :
78 24 : void DRRForceGrid::fillTable(const vector<vector<double>> &in) {
79 48 : table.resize(ndims, vector<double>(sampleSize, 0));
80 60 : for (size_t i = 0; i < ndims; ++i) {
81 : size_t repeatAll = 1, repeatOne = 1;
82 48 : for (size_t j = i + 1; j < ndims; ++j)
83 12 : repeatOne *= in[j].size();
84 60 : for (size_t j = 0; j < i; ++j)
85 12 : repeatAll *= in[j].size();
86 : size_t in_i_sz = in[i].size();
87 5068 : for (size_t l = 0; l < in_i_sz; ++l)
88 2516 : std::fill_n(begin(table[i]) + l * repeatOne, repeatOne, in[i][l]);
89 732 : for (size_t k = 0; k < repeatAll - 1; ++k)
90 : std::copy_n(begin(table[i]), repeatOne * in_i_sz,
91 732 : begin(table[i]) + repeatOne * in_i_sz * (k + 1));
92 : }
93 24 : }
94 :
95 30 : DRRForceGrid::DRRForceGrid()
96 : : suffix(""), ndims(0), dimensions(0), sampleSize(0),
97 : headers(""), table(0), forces(0), samples(0), endpoints(0), shifts(0),
98 30 : outputunit(1.0) {}
99 :
100 20 : DRRForceGrid::DRRForceGrid(const vector<DRRAxis> &p_dimensions,
101 20 : const string &p_suffix, bool initializeTable)
102 140 : : suffix(p_suffix), ndims(p_dimensions.size()), dimensions(p_dimensions) {
103 20 : sampleSize = 1;
104 40 : vector<vector<double>> mp(ndims);
105 40 : std::stringstream ss;
106 20 : ss << "# " << ndims << '\n';
107 20 : shifts.resize(ndims, 0);
108 20 : shifts[0] = 1;
109 80 : for (size_t i = 0; i < ndims; ++i) {
110 30 : sampleSize = dimensions[i].nbins * sampleSize;
111 60 : mp[i] = dimensions[i].getMiddlePoints();
112 30 : if (i > 0) {
113 40 : shifts[i] = shifts[i - 1] * dimensions[i - 1].nbins;
114 : }
115 : ss.precision(std::numeric_limits<double>::max_digits10);
116 30 : ss << std::fixed << "# " << dimensions[i].min << ' '
117 60 : << dimensions[i].binWidth << ' ' << dimensions[i].nbins;
118 30 : if (dimensions[i].isPeriodic())
119 : ss << " 1" << '\n';
120 : else
121 : ss << " 0" << '\n';
122 : }
123 40 : headers = ss.str();
124 20 : if (initializeTable)
125 12 : fillTable(mp);
126 20 : forces.resize(sampleSize * ndims, 0.0);
127 20 : samples.resize(sampleSize, 0);
128 20 : outputunit = 1.0;
129 : // For 1D pmf
130 20 : if (ndims == 1) {
131 20 : endpoints.resize(dimensions[0].nbins + 1, 0);
132 10 : double ep = dimensions[0].min;
133 10 : double stride = dimensions[0].binWidth;
134 460 : for (auto &i : endpoints) {
135 450 : i = ep;
136 450 : ep += stride;
137 : }
138 : }
139 20 : }
140 :
141 391412 : bool DRRForceGrid::isInBoundary(const vector<double> &pos) const {
142 : bool result = true;
143 1950524 : for (size_t i = 0; i < ndims; ++i) {
144 1560712 : if (pos[i] < dimensions[i].min || pos[i] > dimensions[i].max)
145 : return false;
146 : }
147 : return result;
148 : }
149 :
150 976171 : size_t DRRForceGrid::sampleAddress(const vector<double> &pos) const {
151 : size_t saddr = 0;
152 4873071 : for (size_t i = 0; i < ndims; ++i) {
153 5845350 : saddr += shifts[i] * index1D(dimensions[i], pos[i]);
154 : }
155 976171 : return saddr;
156 : }
157 :
158 906 : bool DRRForceGrid::store(const vector<double> &pos, const vector<double> &f,
159 : unsigned long int nsamples) {
160 906 : if (isInBoundary(pos)) {
161 906 : if (nsamples == 0)
162 : return true;
163 506 : const size_t baseaddr = sampleAddress(pos) * ndims;
164 1012 : samples[baseaddr / ndims] += nsamples;
165 506 : auto it_fa = begin(forces) + baseaddr;
166 506 : std::transform(begin(f), end(f), it_fa, it_fa, std::plus<double>());
167 506 : return true;
168 : } else {
169 : return false;
170 : }
171 : }
172 :
173 2 : vector<DRRAxis> DRRForceGrid::merge(const vector<DRRAxis> &dA,
174 : const vector<DRRAxis> &dB) {
175 2 : vector<DRRAxis> dR(dA.size());
176 2 : std::transform(begin(dA), end(dA), begin(dB), begin(dR), DRRAxis::merge);
177 2 : return dR;
178 : }
179 :
180 : vector<double>
181 800 : DRRForceGrid::getAccumulatedForces(const vector<double> &pos) const {
182 800 : vector<double> result(ndims, 0);
183 800 : if (!isInBoundary(pos))
184 : return result;
185 400 : const size_t baseaddr = sampleAddress(pos) * ndims;
186 : std::copy(begin(forces) + baseaddr, begin(forces) + baseaddr + ndims,
187 : begin(result));
188 400 : return result;
189 : }
190 :
191 66310 : unsigned long int DRRForceGrid::getCount(const vector<double> &pos,
192 : bool SkipCheck) const {
193 66310 : if (!SkipCheck) {
194 800 : if (!isInBoundary(pos)) {
195 : return 0;
196 : }
197 : }
198 131820 : return samples[sampleAddress(pos)];
199 : }
200 :
201 195083 : vector<double> DRRForceGrid::getGradient(const vector<double> &pos,
202 : bool SkipCheck) const {
203 195083 : vector<double> result(ndims, 0);
204 195083 : if (!SkipCheck) {
205 162000 : if (!isInBoundary(pos)) {
206 : return result;
207 : }
208 : }
209 195083 : const size_t baseaddr = sampleAddress(pos);
210 : const unsigned long int &count = samples[baseaddr];
211 195083 : if (count == 0)
212 : return result;
213 194998 : auto it_fa = begin(forces) + baseaddr * ndims;
214 194998 : std::transform(it_fa, it_fa + ndims, begin(result),
215 584386 : [&count](double fa) { return (-1.0) * fa / count; });
216 194998 : return result;
217 : }
218 :
219 64800 : double DRRForceGrid::getDivergence(const vector<double> &pos) const {
220 : double div = 0.0;
221 64800 : vector<double> grad_deriv(ndims, 0.0);
222 64800 : if (!isInBoundary(pos)) {
223 : return div;
224 : }
225 64800 : const size_t force_addr = sampleAddress(pos) * ndims;
226 64800 : vector<double> grad = getGradient(pos);
227 324000 : for (size_t i = 0; i < ndims; ++i) {
228 129600 : const double binWidth = dimensions[i].binWidth;
229 129600 : vector<double> first = pos;
230 259200 : first[i] = dimensions[i].min + binWidth * 0.5;
231 129600 : vector<double> last = pos;
232 259200 : last[i] = dimensions[i].max - binWidth * 0.5;
233 129600 : const size_t force_addr_first = sampleAddress(first) * ndims;
234 129600 : const size_t force_addr_last = sampleAddress(last) * ndims;
235 129600 : if (force_addr == force_addr_first) {
236 720 : if (dimensions[i].isRealPeriodic() == true) {
237 720 : vector<double> next = pos;
238 720 : next[i] += binWidth;
239 720 : const vector<double> grad_next = getGradient(next);
240 720 : const vector<double> grad_prev = getGradient(last);
241 2160 : grad_deriv[i] = (grad_next[i] - grad_prev[i]) / (2.0 * binWidth);
242 : } else {
243 0 : vector<double> next = pos;
244 0 : next[i] += binWidth;
245 0 : vector<double> next_2 = next;
246 0 : next_2[i] += binWidth;
247 0 : const vector<double> grad_next = getGradient(next);
248 0 : const vector<double> grad_next_2 = getGradient(next_2);
249 0 : grad_deriv[i] =
250 0 : (grad_next_2[i] * -1.0 + grad_next[i] * 4.0 - grad[i] * 3.0) /
251 0 : (2.0 * binWidth);
252 : }
253 128880 : } else if (force_addr == force_addr_last) {
254 720 : if (dimensions[i].isRealPeriodic() == true) {
255 720 : vector<double> prev = pos;
256 720 : prev[i] -= binWidth;
257 720 : const vector<double> grad_next = getGradient(first);
258 720 : const vector<double> grad_prev = getGradient(prev);
259 2160 : grad_deriv[i] = (grad_next[i] - grad_prev[i]) / (2.0 * binWidth);
260 : } else {
261 0 : vector<double> prev = pos;
262 0 : prev[i] -= binWidth;
263 0 : vector<double> prev_2 = prev;
264 0 : prev_2[i] -= binWidth;
265 0 : const vector<double> grad_prev = getGradient(prev);
266 0 : const vector<double> grad_prev_2 = getGradient(prev_2);
267 0 : grad_deriv[i] =
268 0 : (grad[i] * 3.0 - grad_prev[i] * 4.0 + grad_prev_2[i] * 1.0) /
269 0 : (2.0 * binWidth);
270 : }
271 : } else {
272 128160 : vector<double> prev = pos;
273 128160 : prev[i] -= binWidth;
274 128160 : vector<double> next = pos;
275 128160 : next[i] += binWidth;
276 128160 : const vector<double> grad_next = getGradient(next);
277 128160 : const vector<double> grad_prev = getGradient(prev);
278 384480 : grad_deriv[i] = (grad_next[i] - grad_prev[i]) / (2.0 * binWidth);
279 : }
280 : }
281 : div = std::accumulate(begin(grad_deriv), end(grad_deriv), 0.0);
282 : return div;
283 : }
284 :
285 : vector<double>
286 195083 : DRRForceGrid::getCountsLogDerivative(const vector<double> &pos) const {
287 195083 : const size_t addr = sampleAddress(pos);
288 195083 : const unsigned long int count_this = samples[addr];
289 195083 : vector<double> result(ndims, 0);
290 974103 : for (size_t i = 0; i < ndims; ++i) {
291 389510 : const double binWidth = dimensions[i].binWidth;
292 : const size_t addr_first =
293 779020 : addr - shifts[i] * index1D(dimensions[i], pos[i]) + 0;
294 779020 : const size_t addr_last = addr_first + shifts[i] * (dimensions[i].nbins - 1);
295 389510 : if (addr == addr_first) {
296 2190 : if (dimensions[i].isRealPeriodic() == true) {
297 2178 : const unsigned long int &count_next = samples[addr + shifts[i]];
298 : const unsigned long int &count_prev = samples[addr_last];
299 2178 : if (count_next != 0 && count_prev != 0)
300 2160 : result[i] =
301 6480 : (std::log(count_next) - std::log(count_prev)) / (2 * binWidth);
302 : } else {
303 12 : const unsigned long int &count_next = samples[addr + shifts[i]];
304 12 : const unsigned long int &count_next2 = samples[addr + shifts[i] * 2];
305 12 : if (count_next != 0 && count_this != 0 && count_next2 != 0)
306 4 : result[i] =
307 12 : (std::log(count_next2) * (-1.0) + std::log(count_next) * 4.0 -
308 8 : std::log(count_this) * 3.0) /
309 4 : (2.0 * binWidth);
310 : }
311 387320 : } else if (addr == addr_last) {
312 2190 : if (dimensions[i].isRealPeriodic() == true) {
313 2178 : const unsigned long int &count_prev = samples[addr - shifts[i]];
314 : const unsigned long int &count_next = samples[addr_first];
315 2178 : if (count_next != 0 && count_prev != 0)
316 2160 : result[i] =
317 6480 : (std::log(count_next) - std::log(count_prev)) / (2 * binWidth);
318 : } else {
319 12 : const unsigned long int &count_prev = samples[addr - shifts[i]];
320 12 : const unsigned long int &count_prev2 = samples[addr - shifts[i] * 2];
321 12 : if (count_prev != 0 && count_this != 0 && count_prev2 != 0)
322 16 : result[i] = (std::log(count_this) * 3.0 - std::log(count_prev) * 4.0 +
323 8 : std::log(count_prev2)) /
324 4 : (2.0 * binWidth);
325 : }
326 : } else {
327 385130 : const unsigned long int &count_prev = samples[addr - shifts[i]];
328 385130 : const unsigned long int &count_next = samples[addr + shifts[i]];
329 385130 : if (count_next != 0 && count_prev != 0)
330 385036 : result[i] =
331 1155108 : (std::log(count_next) - std::log(count_prev)) / (2 * binWidth);
332 : }
333 : }
334 195083 : return result;
335 : }
336 :
337 12 : void DRRForceGrid::write1DPMF(string filename) const {
338 24 : filename += suffix + ".pmf";
339 : FILE *ppmf;
340 12 : ppmf = fopen(filename.c_str(), "w");
341 12 : const double w = dimensions[0].binWidth;
342 : double pmf = 0;
343 12 : fprintf(ppmf, "%.9f %.9f\n", endpoints[0], pmf);
344 1324 : for (size_t i = 0; i < dimensions[0].nbins; ++i) {
345 656 : vector<double> pos(1, 0);
346 656 : pos[0] = table[0][i];
347 656 : const vector<double> f = getGradient(pos, true);
348 656 : pmf += f[0] * w / outputunit;
349 1312 : fprintf(ppmf, "%.9f %.9f\n", endpoints[i + 1], pmf);
350 : }
351 12 : fclose(ppmf);
352 12 : }
353 :
354 20 : void DRRForceGrid::writeAll(const string &filename) const {
355 40 : string countname = filename + suffix + ".count";
356 40 : string gradname = filename + suffix + ".grad";
357 20 : vector<double> pos(ndims, 0);
358 : FILE *pGrad, *pCount;
359 20 : pGrad = fopen(gradname.c_str(), "w");
360 20 : pCount = fopen(countname.c_str(), "w");
361 : char *buffer1, *buffer2;
362 20 : buffer1 = (char *)malloc((sizeof(double)) * sampleSize * ndims);
363 20 : buffer2 = (char *)malloc((sizeof(double)) * sampleSize * ndims);
364 20 : setvbuf(pGrad, buffer1, _IOFBF, (sizeof(double)) * sampleSize * ndims);
365 20 : setvbuf(pCount, buffer2, _IOFBF, (sizeof(double)) * sampleSize * ndims);
366 20 : fwrite(headers.c_str(), sizeof(char), strlen(headers.c_str()), pGrad);
367 20 : fwrite(headers.c_str(), sizeof(char), strlen(headers.c_str()), pCount);
368 131040 : for (size_t i = 0; i < sampleSize; ++i) {
369 326238 : for (size_t j = 0; j < ndims; ++j) {
370 130364 : pos[j] = table[j][i];
371 130364 : fprintf(pGrad, " %.9f", table[j][i]);
372 130364 : fprintf(pCount, " %.9f", table[j][i]);
373 : }
374 65510 : fprintf(pCount, " %lu\n", getCount(pos, true));
375 65510 : vector<double> f = getGradient(pos, true);
376 326238 : for (size_t j = 0; j < ndims; ++j) {
377 130364 : fprintf(pGrad, " %.9f", (f[j] / outputunit));
378 : }
379 : fprintf(pGrad, "\n");
380 : }
381 20 : fclose(pGrad);
382 20 : fclose(pCount);
383 20 : free(buffer1);
384 20 : free(buffer2);
385 20 : if (ndims == 1) {
386 24 : write1DPMF(filename);
387 : }
388 20 : }
389 :
390 2 : void DRRForceGrid::writeDivergence(const string &filename) const {
391 4 : string divname = filename + suffix + ".div";
392 2 : vector<double> pos(ndims, 0);
393 : FILE *pDiv;
394 2 : pDiv = fopen(divname.c_str(), "w");
395 2 : fwrite(headers.c_str(), sizeof(char), strlen(headers.c_str()), pDiv);
396 129602 : for (size_t i = 0; i < sampleSize; ++i) {
397 324000 : for (size_t j = 0; j < ndims; ++j) {
398 129600 : pos[j] = table[j][i];
399 129600 : fprintf(pDiv, " %.9f", table[j][i]);
400 : }
401 64800 : const double divergence = getDivergence(pos);
402 64800 : fprintf(pDiv, " %.9f", (divergence / outputunit));
403 : fprintf(pDiv, "\n");
404 : }
405 2 : fclose(pDiv);
406 2 : }
407 :
408 106 : bool ABF::store_getbias(const vector<double> &pos, const vector<double> &f,
409 : vector<double> &fbias) {
410 106 : if (!isInBoundary(pos)) {
411 0 : std::fill(begin(fbias), end(fbias), 0.0);
412 0 : return false;
413 : }
414 106 : const size_t baseaddr = sampleAddress(pos);
415 : unsigned long int &count = samples[baseaddr];
416 106 : ++count;
417 106 : double factor = 2 * (static_cast<double>(count)) / mFullSamples - 1;
418 : // Clamp to [0,maxFactor]
419 106 : factor = factor < 0 ? 0 : factor > mMaxFactor ? mMaxFactor : factor;
420 106 : auto it_fa = begin(forces) + baseaddr * ndims;
421 : auto it_fb = begin(fbias);
422 : auto it_f = begin(f);
423 : do {
424 178 : (*it_fa) += (*it_f); // Accumulate instantaneous force
425 356 : (*it_fb) = factor * (*it_fa) * (-1.0) /
426 178 : static_cast<double>(count); // Calculate bias force
427 : ++it_fa;
428 : ++it_fb;
429 : ++it_f;
430 178 : } while (it_f != end(f));
431 :
432 : return true;
433 : }
434 :
435 1 : ABF ABF::mergewindow(const ABF &aWA, const ABF &aWB) {
436 1 : const vector<DRRAxis> dR = merge(aWA.dimensions, aWB.dimensions);
437 1 : const string suffix = ".abf";
438 : ABF result(dR, suffix);
439 1 : const size_t nrows = result.sampleSize;
440 1 : const size_t ncols = result.ndims;
441 1 : vector<double> pos(ncols, 0);
442 401 : for (size_t i = 0; i < nrows; ++i) {
443 600 : for (size_t j = 0; j < ncols; ++j) {
444 200 : pos[j] = result.table[j][i];
445 : }
446 200 : const unsigned long int countA = aWA.getCount(pos);
447 200 : const unsigned long int countB = aWB.getCount(pos);
448 200 : const vector<double> aForceA = aWA.getAccumulatedForces(pos);
449 200 : const vector<double> aForceB = aWB.getAccumulatedForces(pos);
450 200 : result.store(pos, aForceA, countA);
451 200 : result.store(pos, aForceB, countB);
452 : }
453 1 : return result;
454 : }
455 :
456 195083 : vector<double> CZAR::getGradient(const vector<double> &pos,
457 : bool SkipCheck) const {
458 195083 : vector<double> result(ndims, 0);
459 195083 : if (!SkipCheck) {
460 162000 : if (!isInBoundary(pos)) {
461 : return result;
462 : }
463 : }
464 195083 : if (kbt <= std::numeric_limits<double>::epsilon()) {
465 : std::cerr << "ERROR! The kbt shouldn't be zero when use CZAR estimator. "
466 : << '\n';
467 0 : std::abort();
468 : }
469 195083 : const size_t baseaddr = sampleAddress(pos);
470 195083 : const vector<double> log_deriv(getCountsLogDerivative(pos));
471 : const unsigned long int &count = samples[baseaddr];
472 195083 : if (count == 0)
473 : return result;
474 195008 : auto it_fa = begin(forces) + baseaddr * ndims;
475 195008 : std::transform(it_fa, it_fa + ndims, begin(log_deriv), begin(result),
476 778816 : [&count, this](double fa, double ld) {
477 1168224 : return fa * (-1.0) / count - kbt * ld;
478 584416 : });
479 195008 : return result;
480 : }
481 :
482 1 : CZAR CZAR::mergewindow(const CZAR &cWA, const CZAR &cWB) {
483 1 : const vector<DRRAxis> dR = merge(cWA.dimensions, cWB.dimensions);
484 1 : const double newkbt = cWA.kbt;
485 1 : const string suffix = ".czar";
486 : CZAR result(dR, suffix, newkbt);
487 1 : const size_t nrows = result.sampleSize;
488 1 : const size_t ncols = result.ndims;
489 1 : vector<double> pos(ncols, 0);
490 401 : for (size_t i = 0; i < nrows; ++i) {
491 600 : for (size_t j = 0; j < ncols; ++j) {
492 200 : pos[j] = result.table[j][i];
493 : }
494 200 : const unsigned long int countA = cWA.getCount(pos);
495 200 : const unsigned long int countB = cWB.getCount(pos);
496 200 : const vector<double> aForceA = cWA.getAccumulatedForces(pos);
497 200 : const vector<double> aForceB = cWB.getAccumulatedForces(pos);
498 200 : result.store(pos, aForceA, countA);
499 200 : result.store(pos, aForceB, countB);
500 : }
501 1 : return result;
502 : }
503 : }
504 5517 : }
505 :
506 : #endif
|
