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             : #ifndef __PLUMED_drr_DRR_h
      19             : #define __PLUMED_drr_DRR_h
      20             : // Build requirement: boost, c++11 compatible compiler.
      21             : #ifdef __PLUMED_HAS_BOOST_SERIALIZATION
      22             : 
      23             : #include <algorithm>
      24             : #include <cmath>
      25             : #include <cstddef>
      26             : #include <cstdlib>
      27             : #include <fstream>
      28             : #include <iomanip>
      29             : #include <iostream>
      30             : #include <iterator>
      31             : #include <limits>
      32             : #include <numeric>
      33             : #include <sstream>
      34             : 
      35             : // boost headers for serialization
      36             : #include <boost/archive/binary_iarchive.hpp>
      37             : #include <boost/archive/binary_oarchive.hpp>
      38             : #include <boost/serialization/string.hpp>
      39             : #include <boost/serialization/vector.hpp>
      40             : 
      41             : namespace PLMD {
      42             : namespace drr {
      43             : 
      44             : using std::vector;
      45             : using std::string;
      46             : using std::begin;
      47             : using std::end;
      48             : 
      49             : /// This class can store the minimum, maximum and bins of a dimension(axis).
      50             : class DRRAxis {
      51             : public:
      52             :   /// Default empty constructor
      53          64 :   DRRAxis() {
      54          64 :     min = max = 0.0;
      55          64 :     nbins = 0;
      56          64 :     periodic = false;
      57          64 :     domainMax = domainMin = 0.0;
      58          64 :     binWidth = 0.0;
      59             :   }
      60             :   /// Constructor using maximum value, minimum value and the number of bins(No
      61             :   /// pbc)
      62             :   DRRAxis(double l, double h, size_t n)
      63             :     : min(l), max(h), nbins(n), periodic(false), domainMax(0), domainMin(0),
      64             :       binWidth((max - min) / double(nbins)) {}
      65             :   /// PBC-aware constructor
      66             :   DRRAxis(double l, double h, size_t n, bool pbc, double dMax, double dMin)
      67           4 :     : min(l), max(h), nbins(n), periodic(pbc), domainMax(dMax),
      68           4 :       domainMin(dMin), binWidth((max - min) / double(nbins)) {}
      69             :   /// Set values
      70             :   void set(double l, double h, size_t n, bool pbc = false, double dmin = 0,
      71             :            double dmax = 0) {
      72          38 :     min = l;
      73          38 :     max = h;
      74          38 :     nbins = n;
      75          38 :     periodic = pbc;
      76          38 :     domainMax = dmax;
      77          38 :     domainMin = dmin;
      78          38 :     binWidth = (max - min) / nbins;
      79             :   }
      80             :   /// Set PBC data
      81             :   void setPeriodicity(double dmin, double dmax) {
      82          16 :     domainMax = dmax;
      83          16 :     domainMin = dmin;
      84          16 :     periodic = true;
      85             :   }
      86             :   /// Getters
      87          41 :   double getMin() const { return this->min; }
      88          40 :   double getMax() const { return this->max; }
      89          31 :   double getWidth() const { return binWidth; }
      90             :   double getDomainMax() const { return this->domainMax; }
      91             :   double getDomainMin() const { return this->domainMin; }
      92             :   size_t getBins() const { return this->nbins; }
      93             : 
      94             :   /// Check periodicity
      95          60 :   bool isPeriodic() const { return this->periodic; }
      96             :   /// Check real periodicity, i.e. the maximum == the domain maximum
      97             :   bool isRealPeriodic() const;
      98             : 
      99             :   /// Check whether x is in this axis
     100             :   bool isInBoundary(double x) const;
     101             :   /// Get an array of middle points of each bins
     102             :   vector<double> getMiddlePoints();
     103             : 
     104             :   /// Combine two axes if they share the same bin width.
     105             :   static DRRAxis merge(const DRRAxis &d1, const DRRAxis &d2);
     106             : 
     107             :   friend class DRRForceGrid;
     108             : 
     109             : protected:
     110             :   double min;       // Minimum value of the axis
     111             :   double max;       // Maximum value of the axis
     112             :   size_t nbins;     // Number of bins
     113             :   bool periodic;    // Periodicity
     114             :   double domainMax; // Maximum value of the CV domain
     115             :   double domainMin; // Minimum value of the CV domain
     116             :   friend class boost::serialization::access;
     117             :   /// Use boost serialization
     118             :   template <typename Archive>
     119          52 :   void save(Archive &ar, const unsigned int version) const {
     120          52 :     ar &min;
     121          52 :     ar &max;
     122          52 :     ar &nbins;
     123          52 :     ar &periodic;
     124          52 :     ar &domainMax;
     125          52 :     ar &domainMin;
     126          52 :   }
     127             :   /// Split save and load. The bin width is calculated after initialization.
     128             :   template <typename Archive>
     129          22 :   void load(Archive &ar, const unsigned int version) {
     130          22 :     ar &min;
     131          22 :     ar &max;
     132          22 :     ar &nbins;
     133          22 :     ar &periodic;
     134          22 :     ar &domainMax;
     135          22 :     ar &domainMin;
     136          22 :     binWidth = (max - min) / double(nbins);
     137          22 :   }
     138             :   template <typename Archive>
     139             :   void serialize(Archive &ar, const unsigned int version) {
     140             :     boost::serialization::split_member(ar, *this, version);
     141             :   }
     142             : 
     143             : private:
     144             :   double binWidth; // bin width
     145             : };
     146             : 
     147             : /// A class for collecting instantaneous forces, calculating average forces and
     148             : /// build CV histogram.
     149             : class DRRForceGrid {
     150             : public:
     151             :   /// Empty constructor
     152             :   DRRForceGrid();
     153             :   /// "Real" constructor
     154             :   /// The 2D table vector is mainly used for print grid points in grad and count
     155             :   /// file.
     156             :   /// So when use binary output we can set initializeTable to false to save
     157             :   /// memory.
     158             :   explicit DRRForceGrid(const vector<DRRAxis> &p_dimensions,
     159             :                         const string &p_suffix,
     160             :                         bool initializeTable = true);
     161             :   /// Check whether a point is in this grid
     162             :   bool isInBoundary(const vector<double> &pos) const;
     163             :   //  /// Get internal indices of a point
     164             :   //  vector<size_t> index(const vector<double> &pos) const;
     165             :   /// Get internal counts address of a point
     166             :   size_t sampleAddress(const vector<double> &pos) const;
     167             :   /// Store instantaneous forces of a point
     168             :   /// nsamples > 1 is useful for merging windows
     169             :   bool store(const vector<double> &pos, const vector<double> &f,
     170             :              unsigned long int nsamples = 1);
     171             :   /// Get accumulated forces of a point
     172             :   vector<double>
     173             :   getAccumulatedForces(const vector<double> &pos) const;
     174             :   /// Get counts of a point
     175             :   unsigned long int getCount(const vector<double> &pos,
     176             :                              bool SkipCheck = false) const;
     177             :   /// Virtual function! get gradients of a point
     178             :   /// CZAR and naive(ABF) have different gradient formulae
     179             :   virtual vector<double> getGradient(const vector<double> &pos,
     180             :                                      bool SkipCheck = false) const;
     181             :   /// Calculate divergence of the mean force field (experimental)
     182             :   double getDivergence(const vector<double> &pos) const;
     183             :   /// Calculate dln(ρ)/dz, useful for CZAR
     184             :   /// This function may be moved to CZAR class in the future
     185             :   vector<double>
     186             :   getCountsLogDerivative(const vector<double> &pos) const;
     187             :   /// Write grad file
     188             : //   void writeGrad(string filename) const;
     189             :   /// Write 1D pmf file on one dimensional occasion
     190             :   void write1DPMF(string filename, const string &fmt="%.9f") const;
     191             :   /// Write count file
     192             : //   void writeCount(string filename) const;
     193             :   /// Write necessary output file in one function (.grad and .count)
     194             :   void writeAll(const string &filename, const string &fmt="%.9f", bool addition = false) const;
     195             :   /// Output divergence (.div) (experimental)
     196             :   void writeDivergence(const string &filename, const string &fmt="%.9f") const;
     197             :   /// merge windows
     198             :   static vector<DRRAxis> merge(const vector<DRRAxis> &dA,
     199             :                                const vector<DRRAxis> &dB);
     200             :   /// Get suffix
     201             :   string getSuffix() const { return suffix; }
     202             :   /// Set unit for .grad output
     203           4 :   void setOutputUnit(double unit) { outputunit = unit; }
     204             :   /// Destructor
     205         228 :   virtual ~DRRForceGrid() {}
     206             : 
     207             : protected:
     208             :   /// The output suffix appended before .grad(.czar.grad) and
     209             :   /// .count(.czar.count)
     210             :   string suffix;
     211             :   /// Number of dimensions
     212             :   size_t ndims;
     213             :   /// Store each axes
     214             :   vector<DRRAxis> dimensions;
     215             :   /// Size of samples
     216             :   size_t sampleSize;
     217             :   /// The header lines of .grad and .count files
     218             :   string headers;
     219             :   /// A table stores the middle points of all dimensions.
     220             :   /// For output in .grad and .count files
     221             :   vector<vector<double>> table;
     222             :   /// Store the average force of each bins
     223             :   vector<double> forces;
     224             :   /// Store counts of each bins
     225             :   vector<unsigned long int> samples;
     226             :   /// Only for 1D pmf output
     227             :   vector<double> endpoints;
     228             :   /// For faster indexing
     229             :   /// shifts[0] = 1, shifts[n+1] = shifts[n] * dimensions[n].nbins
     230             :   vector<size_t> shifts;
     231             :   /// For set different output units
     232             :   double outputunit;
     233             : 
     234             :   /// Miscellaneous helper functions
     235             :   static size_t index1D(const DRRAxis &c, double x);
     236             :   void fillTable(const vector<vector<double>> &in);
     237             : 
     238             :   /// Boost serialization functions
     239             :   friend class boost::serialization::access;
     240             :   template <class Archive>
     241          38 :   void save(Archive &ar, const unsigned int version) const {
     242             :     // Don't save all members.
     243          38 :     ar << suffix;
     244          38 :     ar << dimensions;
     245          38 :     ar << forces;
     246          38 :     ar << samples;
     247          38 :   }
     248          16 :   template <class Archive> void load(Archive &ar, const unsigned int version) {
     249          16 :     ar >> suffix;
     250          16 :     ar >> dimensions;
     251          16 :     ar >> forces;
     252          16 :     ar >> samples;
     253             :     // Restore other members.
     254          16 :     ndims = dimensions.size();
     255          16 :     sampleSize = samples.size();
     256          16 :     std::stringstream ss;
     257          16 :     ss << "# " << ndims << '\n';
     258          16 :     vector<vector<double>> mp(ndims);
     259          16 :     shifts.resize(ndims, 0);
     260          16 :     shifts[0] = 1;
     261          38 :     for (size_t i = 0; i < ndims; ++i) {
     262          22 :       mp[i] = dimensions[i].getMiddlePoints();
     263          22 :       if (i > 0) {
     264           6 :         shifts[i] = shifts[i - 1] * dimensions[i - 1].nbins;
     265             :       }
     266             :       ss.precision(std::numeric_limits<double>::max_digits10);
     267          44 :       ss << std::fixed << "# " << dimensions[i].min << ' '
     268          44 :          << dimensions[i].binWidth << ' ' << dimensions[i].nbins;
     269          22 :       if (dimensions[i].isPeriodic())
     270          24 :         ss << " 1" << '\n';
     271             :       else
     272          20 :         ss << " 0" << '\n';
     273             :     }
     274          16 :     fillTable(mp);
     275          16 :     headers = ss.str();
     276          16 :     outputunit = 1.0;
     277             :     // For 1D pmf
     278          16 :     if (ndims == 1) {
     279          10 :       endpoints.resize(dimensions[0].nbins + 1, 0);
     280          10 :       double ep = dimensions[0].min;
     281          10 :       double stride = dimensions[0].binWidth;
     282        1020 :       for (auto it = begin(endpoints); it != end(endpoints); ++it) {
     283        1010 :         (*it) = ep;
     284        1010 :         ep += stride;
     285             :       }
     286             :     }
     287          16 :   }
     288             :   template <typename Archive>
     289             :   void serialize(Archive &ar, const unsigned int version) {
     290             :     boost::serialization::split_member(ar, *this, version);
     291             :   }
     292             : };
     293             : 
     294             : class ABF : public DRRForceGrid {
     295             : public:
     296          16 :   ABF() {}
     297           2 :   ABF(const vector<DRRAxis> &p_dimensions, const string &p_suffix,
     298             :       double fullSamples = 500.0, double maxFactor = 1.0,
     299             :       bool initializeTable = true)
     300           2 :     : DRRForceGrid(p_dimensions, p_suffix, initializeTable),
     301           2 :       mFullSamples(fullSamples), mMaxFactors(p_dimensions.size(), maxFactor) {}
     302          11 :   ABF(const vector<DRRAxis> &p_dimensions, const string &p_suffix,
     303             :       double fullSamples, const vector<double>& maxFactors,
     304             :       bool initializeTable = true)
     305          11 :     : DRRForceGrid(p_dimensions, p_suffix, initializeTable),
     306          11 :       mFullSamples(fullSamples), mMaxFactors(maxFactors) {}
     307             :   // Provide a setter for ABF parametres (fullsamples, maxfactor)
     308             :   void setParameters(double fullSamples, const vector<double>& maxFactors) {
     309           1 :     mFullSamples = fullSamples;
     310           1 :     mMaxFactors = maxFactors;
     311           1 :   }
     312             :   // Store the "instantaneous" spring force of a point and get ABF bias forces.
     313             :   bool getbias(const vector<double> &pos,
     314             :                const vector<double> &f,
     315             :                vector<double> &fbias) const;
     316             :   void store_force(const vector<double> &pos,
     317             :                    const vector<double> &f);
     318             :   static ABF mergewindow(const ABF &aWA, const ABF &aWB);
     319          38 :   ~ABF() {}
     320             : 
     321             : private:
     322             :   // Parametres for calculate bias force
     323             :   double mFullSamples;
     324             :   vector<double> mMaxFactors;
     325             :   // Boost serialization
     326             :   friend class boost::serialization::access;
     327             :   template <typename Archive>
     328             :   void serialize(Archive &ar, const unsigned int version) {
     329          27 :     ar &boost::serialization::base_object<DRRForceGrid>(*this);
     330             :   }
     331             : };
     332             : 
     333          28 : class CZAR : public DRRForceGrid {
     334             : public:
     335          16 :   CZAR() : kbt(0) {}
     336             :   CZAR(const vector<DRRAxis> &p_dimensions, const string &p_suffix,
     337             :        double p_kbt, bool initializeTable = true)
     338          13 :     : DRRForceGrid(p_dimensions, p_suffix, initializeTable), kbt(p_kbt) {}
     339             :   vector<double> getGradient(const vector<double> &pos,
     340             :                              bool SkipCheck = false) const;
     341             :   double getkbt() const { return kbt; }
     342             :   void setkbt(double p_kbt) { kbt = p_kbt; }
     343             :   static CZAR mergewindow(const CZAR &cWA, const CZAR &cWB);
     344             :   void writeZCountZGrad(const string &filename, bool addition = false) const;
     345          23 :   ~CZAR() {}
     346             : 
     347             : private:
     348             :   double kbt;
     349             :   friend class boost::serialization::access;
     350             :   template <typename Archive>
     351          27 :   void serialize(Archive &ar, const unsigned int version) {
     352          27 :     ar &boost::serialization::base_object<DRRForceGrid>(*this);
     353          27 :     ar &kbt;
     354          27 :   }
     355             : };
     356             : }
     357             : }
     358             : 
     359             : #endif
     360             : #endif