eabf1D.C

Go to the documentation of this file.
00001 //
00002 // The extended adaptive biasing force method has been contributed to NAMD by the following authors:
00003 //
00004 //    Haohao Fu and Christophe Chipot
00005 //    Laboratoire International Associ\'e
00006 //    Centre National de la Recherche Scientifique et University of Illinois at Urbana--Champaign
00007 //    Unit\'e Mixte de Recherche No. 7565, Universit\'e de Lorraine
00008 //    B.P. 70239, 54506 Vand\oe uvre-lès-Nancy cedex, France
00009 //
00010 // Copyright 2016, Centre National de la Recherche Scientifique
00011 //
00012 
00013 #ifndef EABF1D_CPP
00014 #define EABF1D_CPP
00015 
00016 #include <cmath>
00017 #include <fstream>
00018 #include <string>
00019 #include <memory>
00020 #include <iomanip>
00021 #include "fstream_namd.h"
00022 
00023 #include "eabf1D.h"
00024 #include "eabffunc.h"
00025 
00026 
00027 // constructor
00028 eABF1D::eABF1D(const double _lowerboundary, const double _upperboundary, const double _width,
00029         const double _krestr,
00030         const std::string& _outputfile,
00031         const int _outputfreq,
00032         const bool _restart, const std::string& _inputfile,
00033         const bool _outputgrad, const int _gradfreq, const double _temperature) :
00034         eABF(_outputfile, _outputfreq, _restart, _inputfile, _outputgrad, _gradfreq, _temperature)
00035 {
00036         lowerboundary.push_back(_lowerboundary);
00037         upperboundary.push_back(_upperboundary);
00038         width.push_back(_width);
00039         krestr.push_back(_krestr);
00040         this->initialize();
00041         if (restart == true)
00042                 this->readfile();
00043 }
00044 
00045 bool eABF1D::initialize()
00046 {
00047         // calculate the number of bins in the calculation
00048         // the bins is equal to "real bins + 20" to prevent boundary effect
00049     min.push_back(eabffunc::doubletoint(lowerboundary[0] / width[0]) - 10);
00050     max.push_back(eabffunc::doubletoint(upperboundary[0] / width[0]) + 10);
00051         bins.push_back(max[0] - min[0]);
00052 
00053     // initialize the distribution of "countall"
00054         for (int i = 0; i < bins[0]; i++)
00055     {
00056                 countall.push_back(std::vector<int>(bins[0], 0));
00057     }
00058 
00059     // initialize other variables
00060         sumx = std::vector<double>(bins[0], 0);
00061         sumx2 = std::vector<double>(bins[0], 0);
00062         county = std::vector<int>(bins[0], 0);
00063     line = 0;
00064     return true;
00065 }
00066 
00067 bool eABF1D::update(const std::string& colvarstring)
00068 {
00069     // this is the string obtained from colvar each step
00070     // the format of this string is the same as those in colvar.traj file
00071 
00072         std::vector<std::string> data;
00073         eabffunc::split(colvarstring, data);  // split the string to a vector<string>
00074 
00075     int step;
00076 
00077     double abf_force;       // the value of the two columns characterizing the value of CV and extended CV
00078     double eabf_force;
00079 
00080         step = eabffunc::chartoint(data[0]);    // read the variables fron the "colvarstring"
00081         abf_force = eabffunc::chartodouble(data[col[0]]);
00082         eabf_force = eabffunc::chartodouble(data[col[1]]);
00083 
00084 
00085         // for dihedral RC, it is possible that abf_force = 179 and eabf_force = -179, should correct it
00086     if(abf_force > 150 && eabf_force < -150)
00087         eabf_force += 360;
00088     else if(abf_force < -150 && eabf_force > 150)
00089         eabf_force -= 360;
00090 
00091         // normalize the index of the value of CV and extended CV
00092         int binx = int(floor(abf_force / width[0])) - min[0];
00093         int biny = int(floor(eabf_force / width[0])) - min[0];
00094 
00095     // update the 
00096         if (binx < 0 || binx >= bins[0] || biny < 0 || biny >= bins[0])
00097         return false;
00098     else
00099     {
00100                 countall[binx][biny]++;
00101         // if this is the first time add "countall", then the line in outputfile will increase
00102         if(countall[binx][biny]==1)
00103             line++;
00104         county[biny]++;
00105     }
00106 
00107     sumx[biny] += abf_force;
00108     sumx2[biny] += abf_force * abf_force;
00109 
00110         if (step%outputfreq == 0)
00111     {
00112         writefile();
00113     }
00114 
00115     if(outputgrad==1&&(step%gradfreq)==0)
00116     {
00117         calpmf();
00118     }
00119 
00120     return true;
00121 }
00122 
00123 bool eABF1D::readfile()
00124 {
00125     std::ifstream file(inputfile.c_str(), std::ios::in);
00126 
00127     // output file format: (also see writefile function)
00128     // (number of lines) lowerboundary width (number of bins) krestr temperature
00129         // (normalized CV) (normalized extended CV) numbers
00130         // ......  "only non-zero value is recorded here"
00131         // (normalized CV) sumx sumx2 county
00132         // ......
00133     int temp_line, temp_bins;
00134     double temp_lowerboundary, temp_width;
00135         file >> temp_line >> temp_lowerboundary >> temp_width >> temp_bins >> krestr[0] >> temperature;
00136 
00137     temp_bins = temp_bins + 20;
00138 
00139         int x = 0, y = 0, temp_countall;
00140         for (int i = 0; i < temp_line; i++)
00141     {
00142                 file >> x >> y;
00143                 file >> temp_countall;
00144                 //if (x > temp_bins - 10 || x < 10)
00145                 //      continue;
00146         ;
00147                 if (countall[convertscale(temp_lowerboundary, x)][convertscale(temp_lowerboundary, y)] == 0)
00148                         line++;
00149         countall[convertscale(temp_lowerboundary,x)][convertscale(temp_lowerboundary,y)] += temp_countall;
00150     }
00151 
00152         double temp_sumx, temp_sumx2;
00153         int temp_county;
00154         for (int i = 0; i < temp_bins; i++)
00155     {
00156                 file >> x >> temp_sumx >> temp_sumx2 >> temp_county;
00157                 //if (x > temp_bins - 10 || x < 10)
00158                 //      continue;
00159         ;
00160                 sumx[convertscale(temp_lowerboundary, x)] += temp_sumx;
00161                 sumx2[convertscale(temp_lowerboundary, x)] += temp_sumx2;
00162                 county[convertscale(temp_lowerboundary, x)] += temp_county;
00163     }
00164 
00165     file.close();
00166     return true;
00167 }
00168 
00169 bool eABF1D::writefile() const
00170 {
00171     ofstream_namd file(outputfile.c_str(),".old");
00172     file<<std::setprecision(15);
00173 
00174         file << line << " " << lowerboundary[0] << " " << width[0] << " " << bins[0] - 20 << " " << krestr[0] << " " << temperature << std::endl;
00175 
00176         for (int i = 0; i < bins[0]; i++)
00177                 for (int j = 0; j < bins[0]; j++)
00178                         if (countall[i][j] != 0)
00179                                 // the number of line should be "this->line", I guess :)
00180                                 file << i << " " << j << " " << countall[i][j] << std::endl;
00181 
00182         for (int i = 0; i < bins[0]; i++)
00183                 file << i << " " << sumx[i] << " " << sumx2[i] << " " << county[i] << " " << std::endl;
00184 
00185     file.close();
00186     return true;
00187 }
00188 
00189 bool eABF1D::writehead(ofstream_namd& os) const
00190 {
00191         os << "# 1" << std::endl;
00192         os << "# " << lowerboundary[0] << " " << width[0] << " " << bins[0] - 20 << " " << 0 << std::endl;
00193         os << std::endl;
00194     return true;
00195 }
00196 
00197 bool eABF1D::calpmf() const
00198 {
00199         // implementation of YangWei's estimator
00200         // mimic Jerome's script
00201         int norm;
00202         std::vector<double> x_av(bins[0]);
00203         std::vector<double> sigma2(bins[0]);
00204         for (int biny = 0; biny < bins[0]; biny++)
00205         {
00206                 norm = county[biny] > 0 ? county[biny] : 1;
00207                 x_av[biny] = sumx[biny] / norm;
00208                 sigma2[biny] = sumx2[biny] / norm - x_av[biny] * x_av[biny];
00209         }
00210 
00211         static std::string gridfilename = outputfile + ".grad";
00212         static std::string histfilename = outputfile + ".hist.grad";
00213         static std::string pmffilename = outputfile + ".pmf";
00214         static std::string countfilename = outputfile + ".count";
00215 
00216         ofstream_namd ofile(gridfilename.c_str(),".old");
00217         ofstream_namd ofile_hist(histfilename.c_str(), std::ios::app);
00218         ofstream_namd ofile_pmf(pmffilename.c_str(),".old");
00219         ofstream_namd ofile_count(countfilename.c_str(),".old");
00220 
00221         writehead(ofile);
00222         writehead(ofile_hist);
00223         writehead(ofile_count);
00224 
00225         ofile_pmf << lowerboundary[0] << " 0.0" << std::endl;
00226 
00227         double pmf = 0, pmf_x = 0;
00228         double x, y, av, diff_av, grad;
00229         for (int binx = 0; binx < bins[0]; binx++)
00230         {
00231                 x = (binx + min[0] + 0.5) * width[0];
00232                 norm = 0;
00233                 av = 0;
00234                 diff_av = 0;
00235 
00236                 for (int biny = 0; biny < bins[0]; biny++)
00237                 {
00238                         y = (biny + min[0] + 0.5) * width[0];
00239                         norm += countall[binx][biny];
00240 
00241                         if (sigma2[biny] > 0.00001 || sigma2[biny] < -0.00001)
00242                         {
00243                                 av += countall[binx][biny] * (x - x_av[biny]) / sigma2[biny];
00244                         }
00245                         else if (countall[binx][biny] > 1)
00246                         {
00247                                 // TODO, print error
00248                         }
00249                         diff_av += countall[binx][biny] * (x - y);
00250                 }
00251 
00252                 diff_av /= (norm > 0 ? norm : 1);
00253                 av = eabffunc::BOLTZMANN * temperature * av / (norm > 0 ? norm : 1);
00254                 grad = av - krestr[0] * diff_av;
00255 
00256 
00257                 if (x >= lowerboundary[0] && x < upperboundary[0])
00258                 {
00259                         ofile << x << " " << grad << std::endl;
00260                         ofile_hist << x << " " << grad << std::endl;
00261                         pmf += grad * width[0];
00262                         pmf_x = x + 0.5 * width[0];
00263                         ofile_pmf << pmf_x << " " << pmf << std::endl;
00264                         ofile_count << x << " " << norm << std::endl;
00265                         // TODO, rescale the pmf file to prevent negetive free energy
00266                 }
00267         }
00268         ofile.close();
00269         ofile_pmf.close();
00270         ofile_hist.close();
00271         ofile_count.close();
00272         return true;
00273 }
00274 
00275 #endif // EABF

Generated on Thu Nov 23 01:17:12 2017 for NAMD by  doxygen 1.4.7