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 2407170 : size_t DRRForceGrid::index1D(const DRRAxis &c, double x) {
74 2407170 : size_t idx = size_t(std::floor((x - c.min) / c.binWidth));
75 2407170 : idx = (idx == c.nbins) ? (c.nbins - 1) : idx;
76 2407170 : 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 393846 : bool DRRForceGrid::isInBoundary(const vector<double> &pos) const {
148 : bool result = true;
149 1175042 : for (size_t i = 0; i < ndims; ++i) {
150 782805 : if (pos[i] < dimensions[i].min || pos[i] > dimensions[i].max) {
151 : return false;
152 : }
153 : }
154 : return result;
155 : }
156 :
157 1012383 : size_t DRRForceGrid::sampleAddress(const vector<double> &pos) const {
158 : size_t saddr = 0;
159 3029541 : for (size_t i = 0; i < ndims; ++i) {
160 2017158 : saddr += shifts[i] * index1D(dimensions[i], pos[i]);
161 : }
162 1012383 : 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 {
350 52 : filename += suffix + ".pmf";
351 : FILE *ppmf;
352 26 : ppmf = fopen(filename.c_str(), "w");
353 26 : const double w = dimensions[0].binWidth;
354 : double pmf = 0;
355 26 : fprintf(ppmf, "%.9f %.9f\n", endpoints[0], pmf);
356 1166 : for (size_t i = 0; i < dimensions[0].nbins; ++i) {
357 1140 : vector<double> pos(1, 0);
358 1140 : pos[0] = table[0][i];
359 1140 : const vector<double> f = getGradient(pos, true);
360 1140 : pmf += f[0] * w / outputunit;
361 1140 : fprintf(ppmf, "%.9f %.9f\n", endpoints[i + 1], pmf);
362 : }
363 26 : fclose(ppmf);
364 26 : }
365 :
366 36 : void DRRForceGrid::writeAll(const string &filename, bool addition) const {
367 36 : const string countname = filename + suffix + ".count";
368 36 : const string gradname = filename + suffix + ".grad";
369 36 : vector<double> pos(ndims, 0);
370 : FILE *pGrad, *pCount;
371 36 : if (addition) {
372 8 : pGrad = fopen(gradname.c_str(), "a");
373 8 : pCount = fopen(countname.c_str(), "a");
374 : } else {
375 28 : pGrad = fopen(gradname.c_str(), "w");
376 28 : pCount = fopen(countname.c_str(), "w");
377 : }
378 :
379 : char *buffer1, *buffer2;
380 36 : buffer1 = (char *)malloc((sizeof(double)) * sampleSize * ndims);
381 36 : buffer2 = (char *)malloc((sizeof(double)) * sampleSize * ndims);
382 36 : setvbuf(pGrad, buffer1, _IOFBF, (sizeof(double)) * sampleSize * ndims);
383 36 : setvbuf(pCount, buffer2, _IOFBF, (sizeof(double)) * sampleSize * ndims);
384 36 : fwrite(headers.c_str(), sizeof(char), strlen(headers.c_str()), pGrad);
385 36 : fwrite(headers.c_str(), sizeof(char), strlen(headers.c_str()), pCount);
386 66048 : for (size_t i = 0; i < sampleSize; ++i) {
387 196896 : for (size_t j = 0; j < ndims; ++j) {
388 130884 : pos[j] = table[j][i];
389 130884 : fprintf(pGrad, " %.9f", table[j][i]);
390 130884 : fprintf(pCount, " %.9f", table[j][i]);
391 : }
392 66012 : fprintf(pCount, " %lu\n", getCount(pos, true));
393 66012 : vector<double> f = getGradient(pos, true);
394 196896 : for (size_t j = 0; j < ndims; ++j) {
395 130884 : fprintf(pGrad, " %.9f", (f[j] / outputunit));
396 : }
397 : fprintf(pGrad, "\n");
398 : }
399 36 : fclose(pGrad);
400 36 : fclose(pCount);
401 36 : free(buffer1);
402 36 : free(buffer2);
403 36 : if (ndims == 1) {
404 52 : write1DPMF(filename);
405 : }
406 36 : }
407 :
408 2 : void DRRForceGrid::writeDivergence(const string &filename) const {
409 2 : const string divname = filename + suffix + ".div";
410 2 : vector<double> pos(ndims, 0);
411 : FILE *pDiv;
412 2 : pDiv = fopen(divname.c_str(), "w");
413 2 : fwrite(headers.c_str(), sizeof(char), strlen(headers.c_str()), pDiv);
414 64802 : for (size_t i = 0; i < sampleSize; ++i) {
415 194400 : for (size_t j = 0; j < ndims; ++j) {
416 129600 : pos[j] = table[j][i];
417 129600 : fprintf(pDiv, " %.9f", table[j][i]);
418 : }
419 64800 : const double divergence = getDivergence(pos);
420 64800 : fprintf(pDiv, " %.9f", (divergence / outputunit));
421 : fprintf(pDiv, "\n");
422 : }
423 2 : fclose(pDiv);
424 2 : }
425 :
426 123 : bool ABF::store_getbias(const vector<double> &pos, const vector<double> &f,
427 : vector<double> &fbias) {
428 123 : if (!isInBoundary(pos)) {
429 : std::fill(begin(fbias), end(fbias), 0.0);
430 : return false;
431 : }
432 123 : const size_t baseaddr = sampleAddress(pos);
433 : unsigned long int &count = samples[baseaddr];
434 123 : ++count;
435 123 : double factor = 2 * (static_cast<double>(count)) / mFullSamples - 1;
436 123 : auto it_fa = begin(forces) + baseaddr * ndims;
437 : auto it_fb = begin(fbias);
438 : auto it_f = begin(f);
439 : auto it_maxFactor = begin(mMaxFactors);
440 : do {
441 : // Clamp to [0,maxFactors]
442 207 : factor = factor < 0 ? 0 : factor > (*it_maxFactor) ? (*it_maxFactor) : factor;
443 207 : (*it_fa) += (*it_f); // Accumulate instantaneous force
444 207 : (*it_fb) = factor * (*it_fa) * (-1.0) /
445 207 : static_cast<double>(count); // Calculate bias force
446 : ++it_fa;
447 : ++it_fb;
448 : ++it_f;
449 : ++it_maxFactor;
450 207 : } while (it_f != end(f));
451 :
452 : return true;
453 : }
454 :
455 2 : ABF ABF::mergewindow(const ABF &aWA, const ABF &aWB) {
456 2 : const vector<DRRAxis> dR = merge(aWA.dimensions, aWB.dimensions);
457 2 : const string suffix = ".abf";
458 2 : ABF result(dR, suffix);
459 2 : const size_t nrows = result.sampleSize;
460 2 : const size_t ncols = result.ndims;
461 2 : vector<double> pos(ncols, 0);
462 402 : for (size_t i = 0; i < nrows; ++i) {
463 800 : for (size_t j = 0; j < ncols; ++j) {
464 400 : pos[j] = result.table[j][i];
465 : }
466 400 : const unsigned long int countA = aWA.getCount(pos);
467 400 : const unsigned long int countB = aWB.getCount(pos);
468 400 : const vector<double> aForceA = aWA.getAccumulatedForces(pos);
469 400 : const vector<double> aForceB = aWB.getAccumulatedForces(pos);
470 400 : result.store(pos, aForceA, countA);
471 400 : result.store(pos, aForceB, countB);
472 : }
473 2 : return result;
474 0 : }
475 :
476 195576 : vector<double> CZAR::getGradient(const vector<double> &pos,
477 : bool SkipCheck) const {
478 195576 : vector<double> result(ndims, 0);
479 195576 : if (!SkipCheck) {
480 162000 : if (!isInBoundary(pos)) {
481 : return result;
482 : }
483 : }
484 195576 : if (kbt <= std::numeric_limits<double>::epsilon()) {
485 : std::cerr << "ERROR! The kbt shouldn't be zero when use CZAR estimator. "
486 0 : << '\n';
487 0 : std::abort();
488 : }
489 195576 : const size_t baseaddr = sampleAddress(pos);
490 195576 : const vector<double> log_deriv(getCountsLogDerivative(pos));
491 : const unsigned long int &count = samples[baseaddr];
492 195576 : if (count == 0) {
493 : return result;
494 : }
495 195421 : auto it_fa = begin(forces) + baseaddr * ndims;
496 195421 : std::transform(it_fa, it_fa + ndims, begin(log_deriv), begin(result),
497 389822 : [&count, this](double fa, double ld) {
498 389822 : return fa * (-1.0) / count - kbt * ld;
499 : });
500 195421 : return result;
501 : }
502 :
503 2 : CZAR CZAR::mergewindow(const CZAR &cWA, const CZAR &cWB) {
504 2 : const vector<DRRAxis> dR = merge(cWA.dimensions, cWB.dimensions);
505 2 : const double newkbt = cWA.kbt;
506 2 : const string suffix = ".czar";
507 : CZAR result(dR, suffix, newkbt);
508 2 : const size_t nrows = result.sampleSize;
509 2 : const size_t ncols = result.ndims;
510 2 : vector<double> pos(ncols, 0);
511 402 : for (size_t i = 0; i < nrows; ++i) {
512 800 : for (size_t j = 0; j < ncols; ++j) {
513 400 : pos[j] = result.table[j][i];
514 : }
515 400 : const unsigned long int countA = cWA.getCount(pos);
516 400 : const unsigned long int countB = cWB.getCount(pos);
517 400 : const vector<double> aForceA = cWA.getAccumulatedForces(pos);
518 400 : const vector<double> aForceB = cWB.getAccumulatedForces(pos);
519 400 : result.store(pos, aForceA, countA);
520 400 : result.store(pos, aForceB, countB);
521 : }
522 2 : return result;
523 : }
524 :
525 18 : void CZAR:: writeZCountZGrad(const string &filename, bool addition) const {
526 18 : const string countname = filename + ".zcount";
527 18 : const string gradname = filename + ".zgrad";
528 18 : vector<double> pos(ndims, 0);
529 : FILE *pCount;
530 : FILE *pGrad;
531 18 : if (addition) {
532 4 : pCount = fopen(countname.c_str(), "a");
533 4 : pGrad = fopen(gradname.c_str(), "a");
534 : } else {
535 14 : pCount = fopen(countname.c_str(), "w");
536 14 : pGrad = fopen(gradname.c_str(), "w");
537 : }
538 18 : fwrite(headers.c_str(), sizeof(char), strlen(headers.c_str()), pCount);
539 18 : fwrite(headers.c_str(), sizeof(char), strlen(headers.c_str()), pGrad);
540 33024 : for (size_t i = 0; i < sampleSize; ++i) {
541 98448 : for (size_t j = 0; j < ndims; ++j) {
542 65442 : pos[j] = table[j][i];
543 65442 : fprintf(pCount, " %.9f", table[j][i]);
544 65442 : fprintf(pGrad, " %.9f", table[j][i]);
545 : }
546 33006 : const size_t baseaddr = sampleAddress(pos);
547 : const auto& current_sample = samples[baseaddr];
548 33006 : fprintf(pCount, " %lu\n", current_sample);
549 33006 : if (current_sample == 0) {
550 195 : for (size_t j = 0; j < ndims; ++j) {
551 : fprintf(pGrad, " %.9f", 0.0);
552 : }
553 : } else {
554 98253 : for (size_t j = 0; j < ndims; ++j) {
555 65332 : const double grad = -1.0 * forces[baseaddr * ndims + j] / current_sample;
556 65332 : fprintf(pGrad, " %.9f", grad / outputunit);
557 : }
558 : }
559 : fprintf(pGrad, "\n");
560 : }
561 18 : fclose(pCount);
562 18 : fclose(pGrad);
563 18 : }
564 :
565 : }
566 : }
567 :
568 : #endif
|
