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