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 : #include <stdexcept>
21 :
22 : namespace PLMD {
23 : namespace drr {
24 :
25 : using std::vector;
26 : using std::string;
27 : using std::begin;
28 : using std::end;
29 :
30 0 : bool DRRAxis::isInBoundary(double x) const {
31 0 : if (x < min || x > max) {
32 : return false;
33 : } else {
34 0 : return true;
35 : }
36 : }
37 :
38 5860 : bool DRRAxis::isRealPeriodic() const {
39 5860 : if (periodic == true) {
40 5808 : if (std::abs(domainMax - max) < binWidth &&
41 5796 : std::abs(domainMin - min) < binWidth) {
42 : return true;
43 : } else {
44 12 : return false;
45 : }
46 : } else {
47 : return false;
48 : }
49 : }
50 :
51 4 : DRRAxis DRRAxis::merge(const DRRAxis &d1, const DRRAxis &d2) {
52 4 : const double newmin = std::min(d1.min, d2.min);
53 4 : const double newmax = std::max(d1.max, d2.max);
54 4 : const double newWidth = d1.binWidth;
55 4 : const size_t newbins = size_t(std::nearbyint((newmax - newmin) / newWidth));
56 4 : const bool newpbc = d1.periodic;
57 4 : const double newdmin = std::min(d1.domainMin, d2.domainMin);
58 4 : const double newdmax = std::max(d1.domainMax, d2.domainMax);
59 : DRRAxis result(newmin, newmax, newbins, newpbc, newdmin, newdmax);
60 4 : return result;
61 : }
62 :
63 60 : vector<double> DRRAxis::getMiddlePoints() {
64 60 : vector<double> result(nbins, 0);
65 60 : const double width = binWidth;
66 60 : double temp = min - width / 2;
67 : std::generate(begin(result), end(result), [&temp, &width]() {
68 3382 : temp += width;
69 : return temp;
70 : });
71 60 : return result;
72 : }
73 :
74 2407170 : size_t DRRForceGrid::index1D(const DRRAxis &c, double x) {
75 2407170 : size_t idx = size_t(std::floor((x - c.min) / c.binWidth));
76 2407170 : idx = (idx == c.nbins) ? (c.nbins - 1) : idx;
77 2407170 : return idx;
78 : }
79 :
80 34 : void DRRForceGrid::fillTable(const vector<vector<double>> &in) {
81 34 : table.resize(ndims, vector<double>(sampleSize, 0));
82 82 : for (size_t i = 0; i < ndims; ++i) {
83 : size_t repeatAll = 1, repeatOne = 1;
84 62 : for (size_t j = i + 1; j < ndims; ++j) {
85 14 : repeatOne *= in[j].size();
86 : }
87 62 : for (size_t j = 0; j < i; ++j) {
88 14 : repeatAll *= in[j].size();
89 : }
90 : size_t in_i_sz = in[i].size();
91 3390 : for (size_t l = 0; l < in_i_sz; ++l) {
92 3342 : std::fill_n(begin(table[i]) + l * repeatOne, repeatOne, in[i][l]);
93 : }
94 784 : for (size_t k = 0; k < repeatAll - 1; ++k)
95 736 : std::copy_n(begin(table[i]), repeatOne * in_i_sz,
96 736 : begin(table[i]) + repeatOne * in_i_sz * (k + 1));
97 : }
98 34 : }
99 :
100 38 : DRRForceGrid::DRRForceGrid()
101 38 : : suffix(""), ndims(0), dimensions(0), sampleSize(0),
102 38 : headers(""), table(0), forces(0), samples(0), endpoints(0), shifts(0),
103 76 : outputunit(1.0) {}
104 :
105 26 : DRRForceGrid::DRRForceGrid(const vector<DRRAxis> &p_dimensions,
106 26 : const string &p_suffix, bool initializeTable)
107 26 : : suffix(p_suffix), ndims(p_dimensions.size()), dimensions(p_dimensions) {
108 26 : sampleSize = 1;
109 26 : vector<vector<double>> mp(ndims);
110 26 : std::stringstream ss;
111 26 : ss << "# " << ndims << '\n';
112 26 : shifts.resize(ndims, 0);
113 26 : shifts[0] = 1;
114 64 : for (size_t i = 0; i < ndims; ++i) {
115 38 : sampleSize = dimensions[i].nbins * sampleSize;
116 38 : mp[i] = dimensions[i].getMiddlePoints();
117 38 : if (i > 0) {
118 12 : shifts[i] = shifts[i - 1] * dimensions[i - 1].nbins;
119 : }
120 : ss.precision(std::numeric_limits<double>::max_digits10);
121 76 : ss << std::fixed << "# " << dimensions[i].min << ' '
122 76 : << dimensions[i].binWidth << ' ' << dimensions[i].nbins;
123 38 : if (dimensions[i].isPeriodic()) {
124 56 : ss << " 1" << '\n';
125 : } else {
126 20 : ss << " 0" << '\n';
127 : }
128 : }
129 26 : headers = ss.str();
130 26 : if (initializeTable) {
131 18 : fillTable(mp);
132 : }
133 26 : forces.resize(sampleSize * ndims, 0.0);
134 26 : samples.resize(sampleSize, 0);
135 26 : outputunit = 1.0;
136 : // For 1D pmf
137 26 : if (ndims == 1) {
138 14 : endpoints.resize(dimensions[0].nbins + 1, 0);
139 14 : double ep = dimensions[0].min;
140 14 : double stride = dimensions[0].binWidth;
141 882 : for (auto &i : endpoints) {
142 868 : i = ep;
143 868 : ep += stride;
144 : }
145 : }
146 26 : }
147 :
148 393846 : bool DRRForceGrid::isInBoundary(const vector<double> &pos) const {
149 : bool result = true;
150 1175042 : for (size_t i = 0; i < ndims; ++i) {
151 782805 : if (pos[i] < dimensions[i].min || pos[i] > dimensions[i].max) {
152 : return false;
153 : }
154 : }
155 : return result;
156 : }
157 :
158 1012383 : size_t DRRForceGrid::sampleAddress(const vector<double> &pos) const {
159 : size_t saddr = 0;
160 3029541 : for (size_t i = 0; i < ndims; ++i) {
161 2017158 : saddr += shifts[i] * index1D(dimensions[i], pos[i]);
162 : }
163 1012383 : return saddr;
164 : }
165 :
166 1723 : bool DRRForceGrid::store(const vector<double> &pos, const vector<double> &f,
167 : unsigned long int nsamples) {
168 1723 : if (isInBoundary(pos)) {
169 1714 : if (nsamples == 0) {
170 : return true;
171 : }
172 914 : const size_t baseaddr = sampleAddress(pos);
173 914 : samples[baseaddr] += nsamples;
174 914 : auto it_fa = begin(forces) + baseaddr * ndims;
175 : std::transform(begin(f), end(f), it_fa, it_fa, std::plus<double>());
176 914 : return true;
177 : } else {
178 : return false;
179 : }
180 : }
181 :
182 4 : vector<DRRAxis> DRRForceGrid::merge(const vector<DRRAxis> &dA,
183 : const vector<DRRAxis> &dB) {
184 4 : vector<DRRAxis> dR(dA.size());
185 4 : std::transform(begin(dA), end(dA), begin(dB), begin(dR), DRRAxis::merge);
186 4 : return dR;
187 : }
188 :
189 : vector<double>
190 1600 : DRRForceGrid::getAccumulatedForces(const vector<double> &pos) const {
191 1600 : vector<double> result(ndims, 0);
192 1600 : if (!isInBoundary(pos)) {
193 : return result;
194 : }
195 800 : const size_t force_addr = sampleAddress(pos) * ndims;
196 800 : std::copy(begin(forces) + force_addr, begin(forces) + force_addr + ndims,
197 : begin(result));
198 800 : return result;
199 : }
200 :
201 67612 : unsigned long int DRRForceGrid::getCount(const vector<double> &pos,
202 : bool SkipCheck) const {
203 67612 : if (!SkipCheck) {
204 1600 : if (!isInBoundary(pos)) {
205 : return 0;
206 : }
207 : }
208 66812 : return samples[sampleAddress(pos)];
209 : }
210 :
211 195576 : vector<double> DRRForceGrid::getGradient(const vector<double> &pos,
212 : bool SkipCheck) const {
213 195576 : vector<double> result(ndims, 0);
214 195576 : if (!SkipCheck) {
215 162000 : if (!isInBoundary(pos)) {
216 : return result;
217 : }
218 : }
219 195576 : const size_t baseaddr = sampleAddress(pos);
220 : const unsigned long int &count = samples[baseaddr];
221 195576 : if (count == 0) {
222 : return result;
223 : }
224 195411 : auto it_fa = begin(forces) + baseaddr * ndims;
225 : std::transform(it_fa, it_fa + ndims, begin(result),
226 : [&count](double fa) {
227 389802 : return (-1.0) * fa / count;
228 : });
229 : return result;
230 : }
231 :
232 64800 : double DRRForceGrid::getDivergence(const vector<double> &pos) const {
233 : double div = 0.0;
234 64800 : vector<double> grad_deriv(ndims, 0.0);
235 64800 : if (!isInBoundary(pos)) {
236 : return div;
237 : }
238 64800 : const size_t force_addr = sampleAddress(pos) * ndims;
239 64800 : vector<double> grad = getGradient(pos);
240 194400 : for (size_t i = 0; i < ndims; ++i) {
241 129600 : const double binWidth = dimensions[i].binWidth;
242 129600 : vector<double> first = pos;
243 129600 : first[i] = dimensions[i].min + binWidth * 0.5;
244 129600 : vector<double> last = pos;
245 129600 : last[i] = dimensions[i].max - binWidth * 0.5;
246 129600 : const size_t force_addr_first = sampleAddress(first) * ndims;
247 129600 : const size_t force_addr_last = sampleAddress(last) * ndims;
248 129600 : if (force_addr == force_addr_first) {
249 720 : if (dimensions[i].isRealPeriodic() == true) {
250 720 : vector<double> next = pos;
251 720 : next[i] += binWidth;
252 720 : const vector<double> grad_next = getGradient(next);
253 720 : const vector<double> grad_prev = getGradient(last);
254 720 : grad_deriv[i] = (grad_next[i] - grad_prev[i]) / (2.0 * binWidth);
255 : } else {
256 0 : vector<double> next = pos;
257 0 : next[i] += binWidth;
258 0 : vector<double> next_2 = next;
259 0 : next_2[i] += binWidth;
260 0 : const vector<double> grad_next = getGradient(next);
261 0 : const vector<double> grad_next_2 = getGradient(next_2);
262 0 : grad_deriv[i] =
263 0 : (grad_next_2[i] * -1.0 + grad_next[i] * 4.0 - grad[i] * 3.0) /
264 0 : (2.0 * binWidth);
265 : }
266 128880 : } else if (force_addr == force_addr_last) {
267 720 : if (dimensions[i].isRealPeriodic() == true) {
268 720 : vector<double> prev = pos;
269 720 : prev[i] -= binWidth;
270 720 : const vector<double> grad_next = getGradient(first);
271 720 : const vector<double> grad_prev = getGradient(prev);
272 720 : grad_deriv[i] = (grad_next[i] - grad_prev[i]) / (2.0 * binWidth);
273 : } else {
274 0 : vector<double> prev = pos;
275 0 : prev[i] -= binWidth;
276 0 : vector<double> prev_2 = prev;
277 0 : prev_2[i] -= binWidth;
278 0 : const vector<double> grad_prev = getGradient(prev);
279 0 : const vector<double> grad_prev_2 = getGradient(prev_2);
280 0 : grad_deriv[i] =
281 0 : (grad[i] * 3.0 - grad_prev[i] * 4.0 + grad_prev_2[i] * 1.0) /
282 0 : (2.0 * binWidth);
283 : }
284 : } else {
285 128160 : vector<double> prev = pos;
286 128160 : prev[i] -= binWidth;
287 128160 : vector<double> next = pos;
288 128160 : next[i] += binWidth;
289 128160 : const vector<double> grad_next = getGradient(next);
290 128160 : const vector<double> grad_prev = getGradient(prev);
291 128160 : grad_deriv[i] = (grad_next[i] - grad_prev[i]) / (2.0 * binWidth);
292 : }
293 : }
294 : div = std::accumulate(begin(grad_deriv), end(grad_deriv), 0.0);
295 : return div;
296 : }
297 :
298 : vector<double>
299 195576 : DRRForceGrid::getCountsLogDerivative(const vector<double> &pos) const {
300 195576 : const size_t addr = sampleAddress(pos);
301 195576 : const unsigned long int count_this = samples[addr];
302 195576 : vector<double> result(ndims, 0);
303 585588 : for (size_t i = 0; i < ndims; ++i) {
304 390012 : const double binWidth = dimensions[i].binWidth;
305 : const size_t addr_first =
306 390012 : addr - shifts[i] * index1D(dimensions[i], pos[i]) + 0;
307 390012 : const size_t addr_last = addr_first + shifts[i] * (dimensions[i].nbins - 1);
308 390012 : if (addr == addr_first) {
309 2210 : if (dimensions[i].isRealPeriodic() == true) {
310 2178 : const unsigned long int &count_next = samples[addr + shifts[i]];
311 : const unsigned long int &count_prev = samples[addr_last];
312 2178 : if (count_next != 0 && count_prev != 0)
313 2160 : result[i] =
314 2160 : (std::log(count_next) - std::log(count_prev)) / (2 * binWidth);
315 : } else {
316 32 : const unsigned long int &count_next = samples[addr + shifts[i]];
317 32 : const unsigned long int &count_next2 = samples[addr + shifts[i] * 2];
318 32 : if (count_next != 0 && count_this != 0 && count_next2 != 0)
319 6 : result[i] =
320 6 : (std::log(count_next2) * (-1.0) + std::log(count_next) * 4.0 -
321 6 : std::log(count_this) * 3.0) /
322 6 : (2.0 * binWidth);
323 : }
324 387802 : } else if (addr == addr_last) {
325 2210 : if (dimensions[i].isRealPeriodic() == true) {
326 2178 : const unsigned long int &count_prev = samples[addr - shifts[i]];
327 : const unsigned long int &count_next = samples[addr_first];
328 2178 : if (count_next != 0 && count_prev != 0)
329 2160 : result[i] =
330 2160 : (std::log(count_next) - std::log(count_prev)) / (2 * binWidth);
331 : } else {
332 32 : const unsigned long int &count_prev = samples[addr - shifts[i]];
333 32 : const unsigned long int &count_prev2 = samples[addr - shifts[i] * 2];
334 32 : if (count_prev != 0 && count_this != 0 && count_prev2 != 0)
335 6 : result[i] = (std::log(count_this) * 3.0 - std::log(count_prev) * 4.0 +
336 6 : std::log(count_prev2)) /
337 6 : (2.0 * binWidth);
338 : }
339 : } else {
340 385592 : const unsigned long int &count_prev = samples[addr - shifts[i]];
341 385592 : const unsigned long int &count_next = samples[addr + shifts[i]];
342 385592 : if (count_next != 0 && count_prev != 0)
343 385432 : result[i] =
344 385432 : (std::log(count_next) - std::log(count_prev)) / (2 * binWidth);
345 : }
346 : }
347 195576 : return result;
348 : }
349 :
350 26 : void DRRForceGrid::write1DPMF(string filename, const string &fmt) const {
351 26 : filename += suffix + ".pmf";
352 52 : std::string fmtv=fmt+" "+fmt+"\n";
353 : FILE *ppmf;
354 26 : ppmf = fopen(filename.c_str(), "w");
355 26 : if (!ppmf) {
356 0 : throw std::runtime_error("cannot open " + filename + " for writing");
357 : }
358 26 : const double w = dimensions[0].binWidth;
359 : double pmf = 0;
360 26 : fprintf(ppmf, fmtv.c_str(), endpoints[0], pmf);
361 1166 : for (size_t i = 0; i < dimensions[0].nbins; ++i) {
362 1140 : vector<double> pos(1, 0);
363 1140 : pos[0] = table[0][i];
364 1140 : const vector<double> f = getGradient(pos, true);
365 1140 : pmf += f[0] * w / outputunit;
366 1140 : fprintf(ppmf, fmtv.c_str(), endpoints[i + 1], pmf);
367 : }
368 26 : fclose(ppmf);
369 26 : }
370 :
371 36 : void DRRForceGrid::writeAll(const string &filename, const string &fmt, bool addition) const {
372 36 : const std::string countname = filename + suffix + ".count";
373 36 : const std::string gradname = filename + suffix + ".grad";
374 36 : std::string fmtv=" "+fmt;
375 36 : vector<double> pos(ndims, 0);
376 : FILE *pGrad, *pCount;
377 36 : if (addition) {
378 8 : pGrad = fopen(gradname.c_str(), "a");
379 8 : pCount = fopen(countname.c_str(), "a");
380 : } else {
381 28 : pGrad = fopen(gradname.c_str(), "w");
382 28 : pCount = fopen(countname.c_str(), "w");
383 : }
384 36 : if (!pGrad || !pCount) {
385 0 : if (pGrad) {
386 0 : fclose(pGrad);
387 : }
388 0 : if (pCount) {
389 0 : fclose(pCount);
390 : }
391 0 : throw std::runtime_error("cannot open " + (!pGrad ? gradname : countname) +
392 0 : " for writing");
393 : }
394 :
395 : char *buffer1, *buffer2;
396 36 : buffer1 = (char *)malloc((sizeof(double)) * sampleSize * ndims);
397 36 : buffer2 = (char *)malloc((sizeof(double)) * sampleSize * ndims);
398 36 : setvbuf(pGrad, buffer1, _IOFBF, (sizeof(double)) * sampleSize * ndims);
399 36 : setvbuf(pCount, buffer2, _IOFBF, (sizeof(double)) * sampleSize * ndims);
400 36 : fwrite(headers.c_str(), sizeof(char), strlen(headers.c_str()), pGrad);
401 36 : fwrite(headers.c_str(), sizeof(char), strlen(headers.c_str()), pCount);
402 66048 : for (size_t i = 0; i < sampleSize; ++i) {
403 196896 : for (size_t j = 0; j < ndims; ++j) {
404 130884 : pos[j] = table[j][i];
405 130884 : fprintf(pGrad, fmtv.c_str(), table[j][i]);
406 130884 : fprintf(pCount, fmtv.c_str(), table[j][i]);
407 : }
408 66012 : fprintf(pCount, " %lu\n", getCount(pos, true));
409 66012 : vector<double> f = getGradient(pos, true);
410 196896 : for (size_t j = 0; j < ndims; ++j) {
411 130884 : fprintf(pGrad, fmtv.c_str(), (f[j] / outputunit));
412 : }
413 : fprintf(pGrad, "\n");
414 : }
415 36 : fclose(pGrad);
416 36 : fclose(pCount);
417 36 : free(buffer1);
418 36 : free(buffer2);
419 36 : if (ndims == 1) {
420 52 : write1DPMF(filename, fmt);
421 : }
422 36 : }
423 :
424 2 : void DRRForceGrid::writeDivergence(const string &filename, const string &fmt) const {
425 2 : const string divname = filename + suffix + ".div";
426 2 : std::string fmtv=" "+fmt;
427 2 : vector<double> pos(ndims, 0);
428 : FILE *pDiv;
429 2 : pDiv = fopen(divname.c_str(), "w");
430 2 : if (!pDiv) {
431 0 : throw std::runtime_error("cannot open " + divname + " for writing");
432 : }
433 2 : fwrite(headers.c_str(), sizeof(char), strlen(headers.c_str()), pDiv);
434 64802 : for (size_t i = 0; i < sampleSize; ++i) {
435 194400 : for (size_t j = 0; j < ndims; ++j) {
436 129600 : pos[j] = table[j][i];
437 129600 : fprintf(pDiv, fmtv.c_str(), table[j][i]);
438 : }
439 64800 : const double divergence = getDivergence(pos);
440 64800 : fprintf(pDiv, fmtv.c_str(), (divergence / outputunit));
441 : fprintf(pDiv, "\n");
442 : }
443 2 : fclose(pDiv);
444 2 : }
445 :
446 123 : bool ABF::store_getbias(const vector<double> &pos, const vector<double> &f,
447 : vector<double> &fbias) {
448 123 : if (!isInBoundary(pos)) {
449 : std::fill(begin(fbias), end(fbias), 0.0);
450 : return false;
451 : }
452 123 : const size_t baseaddr = sampleAddress(pos);
453 : unsigned long int &count = samples[baseaddr];
454 123 : ++count;
455 123 : double factor = 2 * (static_cast<double>(count)) / mFullSamples - 1;
456 123 : auto it_fa = begin(forces) + baseaddr * ndims;
457 : auto it_fb = begin(fbias);
458 : auto it_f = begin(f);
459 : auto it_maxFactor = begin(mMaxFactors);
460 : do {
461 : // Clamp to [0,maxFactors]
462 207 : factor = factor < 0 ? 0 : factor > (*it_maxFactor) ? (*it_maxFactor) : factor;
463 207 : (*it_fa) += (*it_f); // Accumulate instantaneous force
464 207 : (*it_fb) = factor * (*it_fa) * (-1.0) /
465 207 : static_cast<double>(count); // Calculate bias force
466 : ++it_fa;
467 : ++it_fb;
468 : ++it_f;
469 : ++it_maxFactor;
470 207 : } while (it_f != end(f));
471 :
472 : return true;
473 : }
474 :
475 2 : ABF ABF::mergewindow(const ABF &aWA, const ABF &aWB) {
476 2 : const vector<DRRAxis> dR = merge(aWA.dimensions, aWB.dimensions);
477 2 : const string suffix = ".abf";
478 2 : ABF result(dR, suffix);
479 2 : const size_t nrows = result.sampleSize;
480 2 : const size_t ncols = result.ndims;
481 2 : vector<double> pos(ncols, 0);
482 402 : for (size_t i = 0; i < nrows; ++i) {
483 800 : for (size_t j = 0; j < ncols; ++j) {
484 400 : pos[j] = result.table[j][i];
485 : }
486 400 : const unsigned long int countA = aWA.getCount(pos);
487 400 : const unsigned long int countB = aWB.getCount(pos);
488 400 : const vector<double> aForceA = aWA.getAccumulatedForces(pos);
489 400 : const vector<double> aForceB = aWB.getAccumulatedForces(pos);
490 400 : result.store(pos, aForceA, countA);
491 400 : result.store(pos, aForceB, countB);
492 : }
493 2 : return result;
494 0 : }
495 :
496 195576 : vector<double> CZAR::getGradient(const vector<double> &pos,
497 : bool SkipCheck) const {
498 195576 : vector<double> result(ndims, 0);
499 195576 : if (!SkipCheck) {
500 162000 : if (!isInBoundary(pos)) {
501 : return result;
502 : }
503 : }
504 195576 : if (kbt <= std::numeric_limits<double>::epsilon()) {
505 : std::cerr << "ERROR! The kbt shouldn't be zero when use CZAR estimator. "
506 0 : << '\n';
507 0 : std::abort();
508 : }
509 195576 : const size_t baseaddr = sampleAddress(pos);
510 195576 : const vector<double> log_deriv(getCountsLogDerivative(pos));
511 : const unsigned long int &count = samples[baseaddr];
512 195576 : if (count == 0) {
513 : return result;
514 : }
515 195421 : auto it_fa = begin(forces) + baseaddr * ndims;
516 195421 : std::transform(it_fa, it_fa + ndims, begin(log_deriv), begin(result),
517 389822 : [&count, this](double fa, double ld) {
518 389822 : return fa * (-1.0) / count - kbt * ld;
519 : });
520 195421 : return result;
521 : }
522 :
523 2 : CZAR CZAR::mergewindow(const CZAR &cWA, const CZAR &cWB) {
524 2 : const vector<DRRAxis> dR = merge(cWA.dimensions, cWB.dimensions);
525 2 : const double newkbt = cWA.kbt;
526 2 : const string suffix = ".czar";
527 : CZAR result(dR, suffix, newkbt);
528 2 : const size_t nrows = result.sampleSize;
529 2 : const size_t ncols = result.ndims;
530 2 : vector<double> pos(ncols, 0);
531 402 : for (size_t i = 0; i < nrows; ++i) {
532 800 : for (size_t j = 0; j < ncols; ++j) {
533 400 : pos[j] = result.table[j][i];
534 : }
535 400 : const unsigned long int countA = cWA.getCount(pos);
536 400 : const unsigned long int countB = cWB.getCount(pos);
537 400 : const vector<double> aForceA = cWA.getAccumulatedForces(pos);
538 400 : const vector<double> aForceB = cWB.getAccumulatedForces(pos);
539 400 : result.store(pos, aForceA, countA);
540 400 : result.store(pos, aForceB, countB);
541 : }
542 2 : return result;
543 : }
544 :
545 18 : void CZAR:: writeZCountZGrad(const string &filename, bool addition) const {
546 18 : const string countname = filename + ".zcount";
547 18 : const string gradname = filename + ".zgrad";
548 18 : vector<double> pos(ndims, 0);
549 : FILE *pCount;
550 : FILE *pGrad;
551 18 : if (addition) {
552 4 : pCount = fopen(countname.c_str(), "a");
553 4 : pGrad = fopen(gradname.c_str(), "a");
554 : } else {
555 14 : pCount = fopen(countname.c_str(), "w");
556 14 : pGrad = fopen(gradname.c_str(), "w");
557 : }
558 18 : if (!pCount || !pGrad) {
559 0 : if (pCount) {
560 0 : fclose(pCount);
561 : }
562 0 : if (pGrad) {
563 0 : fclose(pGrad);
564 : }
565 0 : throw std::runtime_error("cannot open " + (!pCount ? countname : gradname) +
566 0 : " for writing");
567 : }
568 18 : fwrite(headers.c_str(), sizeof(char), strlen(headers.c_str()), pCount);
569 18 : fwrite(headers.c_str(), sizeof(char), strlen(headers.c_str()), pGrad);
570 33024 : for (size_t i = 0; i < sampleSize; ++i) {
571 98448 : for (size_t j = 0; j < ndims; ++j) {
572 65442 : pos[j] = table[j][i];
573 65442 : fprintf(pCount, " %.9f", table[j][i]);
574 65442 : fprintf(pGrad, " %.9f", table[j][i]);
575 : }
576 33006 : const size_t baseaddr = sampleAddress(pos);
577 : const auto& current_sample = samples[baseaddr];
578 33006 : fprintf(pCount, " %lu\n", current_sample);
579 33006 : if (current_sample == 0) {
580 195 : for (size_t j = 0; j < ndims; ++j) {
581 : fprintf(pGrad, " %.9f", 0.0);
582 : }
583 : } else {
584 98253 : for (size_t j = 0; j < ndims; ++j) {
585 65332 : const double grad = -1.0 * forces[baseaddr * ndims + j] / current_sample;
586 65332 : fprintf(pGrad, " %.9f", grad / outputunit);
587 : }
588 : }
589 : fprintf(pGrad, "\n");
590 : }
591 18 : fclose(pCount);
592 18 : fclose(pGrad);
593 18 : }
594 :
595 : }
596 : }
597 :
598 : #endif
|
