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