eabf2D.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 EABF2D_CPP
00014 #define EABF2D_CPP
00015 
00016 #include <cmath>
00017 #include <fstream>
00018 #include <string>
00019 #include <memory>
00020 #include <algorithm>
00021 #include <iomanip>
00022 #include "fstream_namd.h"
00023 
00024 #include "eabf2D.h"
00025 #include "eabffunc.h"
00026 
00027 // constructor
00028 eABF2D::eABF2D(const double _lowerboundary,
00029         const double _upperboundary,
00030         const double _width,
00031         const double _krestr,
00032         const double _lowerboundary2,
00033         const double _upperboundary2,
00034         const double _width2,
00035         const double _krestr2,
00036         const std::string& _outputfile,
00037         const int _outputfreq,
00038         const bool _restart,
00039         const std::string& _inputfile,
00040         const bool _outputgrad,
00041         const int _gradfreq,
00042         const double _temperature) :
00043         eABF(_outputfile, _outputfreq, _restart, _inputfile, _outputgrad, _gradfreq, _temperature)
00044 {
00045         lowerboundary.push_back(_lowerboundary);
00046         lowerboundary.push_back(_lowerboundary2);
00047         upperboundary.push_back(_upperboundary);
00048         upperboundary.push_back(_upperboundary2);
00049         width.push_back(_width);
00050         width.push_back(_width2);
00051         krestr.push_back(_krestr);
00052         krestr.push_back(_krestr2);
00053         this->initialize();
00054         this->restart = _restart;
00055     if(restart == true)
00056         this->readfile();
00057 }
00058 
00059 bool eABF2D::initialize()
00060 {
00061         for (int i = 0; i < 2; i++)
00062         {
00063                 min.push_back(eabffunc::doubletoint(lowerboundary[i] / width[i]) - 10);
00064                 max.push_back(eabffunc::doubletoint(upperboundary[i] / width[i]) + 10);
00065                 bins.push_back(max[i] - min[i]);
00066         }
00067 
00068     // initialize the distribution of "countall"
00069     // be very careful about it
00070         countall = new int*[bins[0] * bins[0]];
00071         for (int i = 0; i < bins[0] * bins[0]; i++)
00072         {
00073                 countall[i] = new int[bins[1] * bins[1]];
00074                 for (int j = 0; j < bins[1] * bins[1]; j++)
00075                         countall[i][j] = 0;
00076     }
00077 
00078     // initialize other variables
00079         for (int i = 0; i < bins[0]; i++)
00080         {
00081                 sumx1.push_back(std::vector<double>(bins[1], 0));
00082                 sumx21.push_back(std::vector<double>(bins[1], 0));
00083                 sumx2.push_back(std::vector<double>(bins[1], 0));
00084                 sumx22.push_back(std::vector<double>(bins[1], 0));
00085                 county.push_back(std::vector<int>(bins[1], 0));
00086         }
00087     line = 0;
00088 
00089     return true;
00090 }
00091 
00092 bool eABF2D::update(const std::string& colvarstring)
00093 {
00094     // this is the string obtained from colvar each step
00095     // the format of this string is the same as those in colvar.traj file
00096 
00097         std::vector<std::string> data;
00098         eabffunc::split(colvarstring, data);
00099 
00100     int step;
00101 
00102     double abf_force;
00103     double eabf_force;
00104     double abf_force2;
00105     double eabf_force2;
00106 
00107         step = eabffunc::chartoint(data[0]);
00108         abf_force = eabffunc::chartodouble(data[col[0]]);
00109         eabf_force = eabffunc::chartodouble(data[col[1]]);
00110         abf_force2 = eabffunc::chartodouble(data[col[2]]);
00111         eabf_force2 = eabffunc::chartodouble(data[col[3]]);
00112 
00113 
00114         if (abf_force > 150 && eabf_force < -150)
00115                 eabf_force += 360;
00116         else if (abf_force < -150 && eabf_force > 150)
00117                 eabf_force -= 360;
00118 
00119         if (abf_force2 > 150 && eabf_force2 < -150)
00120                 eabf_force2 += 360;
00121         else if (abf_force2 < -150 && eabf_force2 > 150)
00122         eabf_force2 -= 360;
00123 
00124         int binx = int(floor(abf_force / width[0])) - min[0];
00125         int biny = int(floor(eabf_force / width[0])) - min[0];
00126         int binx2 = int(floor(abf_force2 / width[1])) - min[1];
00127         int biny2 = int(floor(eabf_force2 / width[1])) - min[1];
00128 
00129         if (step%outputfreq == 0)
00130     {
00131         writefile();
00132     }
00133 
00134         if (outputgrad == true && (step%gradfreq) == 0)
00135     {
00136         calpmf();
00137     }
00138 
00139 
00140         if (binx < 0 || binx >= bins[0] || biny < 0 || biny >= bins[0] || binx2 < 0 || binx2 >= bins[1] || biny2 < 0 || biny2 >= bins[1])
00141                 return false;
00142     else
00143     {
00144         // distribution
00145                 countall[binx*bins[0] + biny][binx2*bins[1] + biny2]++;
00146 
00147         // if this is the first time add "countall", then the line in outputfile will increase
00148                 if (countall[binx*bins[0] + biny][binx2*bins[1] + biny2] == 1)
00149                         line++;
00150         county[biny][biny2]++;
00151     }
00152 
00153     sumx1[biny][biny2] += abf_force;
00154     sumx21[biny][biny2] += abf_force * abf_force;
00155     sumx2[biny][biny2] += abf_force2;
00156         sumx22[biny][biny2] += abf_force2 * abf_force2;
00157 
00158 
00159     // TODO: add convergence judgement
00160     return true;
00161 }
00162 
00163 
00164 // WARNGING: do not test whether this function works !!!
00165 bool eABF2D::readfile()
00166 {
00167     std::ifstream file(inputfile.c_str(), std::ios::in);
00168 
00169     // output file format, the first line is "this->line"
00170         int temp_line, temp_bins, temp_bins2;
00171         double temp_lowerboundary, temp_width, temp_lowerboundary2, temp_width2;
00172         file >> temp_line >> temp_lowerboundary >> temp_width >> temp_bins >> krestr[0] >> temp_lowerboundary2 >> temp_width2 >> temp_bins2 >> krestr[1] >> temperature;
00173 
00174 
00175         temp_bins = temp_bins + 20;
00176         temp_bins2 = temp_bins2 + 20;
00177 
00178         int x = 0, y = 0, m = 0, n = 0;
00179         int pos0 = 0, pos1 = 0, temp_countall = 0;
00180         for (int i = 0; i < temp_line; i++)
00181     {
00182                 file >> x >> y >> m >> n;
00183                 //if (x > temp_bins - 10 || x < 10 || y > temp_bins2 - 10 || y < 10)
00184                 //      continue;
00185                 file >> temp_countall;
00186                 ;
00187                 pos0 = convertscale(temp_lowerboundary, x, 1) * bins[0] + convertscale(temp_lowerboundary, y, 1);
00188                 pos1 = convertscale(temp_lowerboundary2, m, 2) * bins[1] + convertscale(temp_lowerboundary2, n, 2);
00189                 // find the correct place to increase countall
00190                 if (countall[pos0][pos1] == 0)
00191                         line++;
00192                 countall[pos0][pos1] += temp_countall;
00193     }
00194 
00195         double temp_sumx1, temp_sumx21, temp_sumx2, temp_sumx22;
00196         int temp_county;
00197         for (int i = 0; i < temp_bins; i++)
00198     {
00199                 for (int j = 0; j < temp_bins2; j++)
00200         {
00201                         file >> x >> y >> temp_sumx1 >> temp_sumx21 >> temp_sumx2 >> temp_sumx22 >> temp_county;
00202                         //if (x > temp_bins - 10 || x < 10 || y > temp_bins2 - 10 || y < 10)
00203                         //      continue;
00204                         ;
00205                         sumx1[convertscale(temp_lowerboundary, x, 1)][convertscale(temp_lowerboundary2, y, 2)] += temp_sumx1;
00206                         sumx21[convertscale(temp_lowerboundary, x, 1)][convertscale(temp_lowerboundary2, y, 2)] += temp_sumx21;
00207                         sumx2[convertscale(temp_lowerboundary, x, 1)][convertscale(temp_lowerboundary2, y, 2)] += temp_sumx2;
00208                         sumx22[convertscale(temp_lowerboundary, x, 1)][convertscale(temp_lowerboundary2, y, 2)] += temp_sumx22;
00209                         county[convertscale(temp_lowerboundary, x, 1)][convertscale(temp_lowerboundary2, y, 2)] += temp_county;
00210         }
00211     }
00212 
00213     file.close();
00214     return true;
00215 }
00216 
00217 bool eABF2D::writefile() const
00218 {
00219     ofstream_namd file(outputfile.c_str(),".old");
00220     file<<std::setprecision(15);
00221 
00222         file << line << " " << lowerboundary[0] << " " << width[0] << " " << bins[0] - 20 << " " << krestr[0] << " ";
00223         file << lowerboundary[1] << " " << width[1] << " " << bins[1] - 20 << " " << krestr[1] << " " << temperature << std::endl;
00224 
00225         for (int i = 0; i < bins[0]; i++)
00226                 for (int j = 0; j < bins[0]; j++)
00227                         for (int m = 0; m < bins[1]; m++)
00228                                 for (int n = 0; n < bins[1]; n++)
00229                                         if (countall[i*bins[0] + j][m*bins[1] + n] != 0)
00230                                                 // the number of line should be "this->line", I guess :)
00231                                                 file << i << " " << j << " " << m << " " << n << " " << countall[i*bins[0] + j][m*bins[1] + n] << std::endl;
00232 
00233         for (int i = 0; i < bins[0]; i++)
00234                 for (int j = 0; j < bins[1]; j++)
00235                         file << i << " " << j << " " << sumx1[i][j] << " " << sumx21[i][j] << " " << sumx2[i][j] << " " << sumx22[i][j] << " " << county[i][j] << " " << std::endl;
00236 
00237     file.close();
00238     return true;
00239 }
00240 
00241 
00242 //bool eABF2D::writecount(ofstream_namd& os) const
00243 //{
00244 //      double x1,x2;
00245 //    writehead(os);
00246 //    for(int binx=0; binx<bins; binx++)
00247 //    {
00248 //       x1 = (binx + _min + 0.5) * width;
00249 //        for(int binx2=0; binx2<bins2; binx2++)
00250 //       {
00251 //            x2 = (binx2 + _min2 + 0.5) * width2;
00252 //                      if(x1>=lowerboundary && x1<upperboundary && x2>=lowerboundary2 && x2<upperboundary2)
00253 //                              os<<x1<<" "<<x2<<" "<<county[binx][binx2]<<std::endl;
00254 //       }
00255 //      }
00256 //    return true;
00257 //}
00258 
00259 bool eABF2D::writehead(ofstream_namd& os) const
00260 {
00261     os<<"# 2"<<std::endl;
00262         // the "real" bin numbers is "bins - 20"
00263         os << "# " << lowerboundary[0] << " " << width[0] << " " << (bins[0] - 20) << " " << 0 << std::endl;
00264         os << "# " << lowerboundary[1] << " " << width[1] << " " << (bins[1] - 20) << " " << 0 << std::endl;
00265     os<<std::endl;
00266     return true;
00267 }
00268 
00269 bool eABF2D::calpmf() const
00270 {
00271         // double integration
00272         // this is just a conversion from Chipot's script
00273         int norm;
00274         std::vector<std::vector<double> > x_av(bins[0]);
00275         std::vector<std::vector<double> > x_av2(bins[0]);
00276         std::vector<std::vector<double> > sigma21(bins[0]);
00277         std::vector<std::vector<double> > sigma22(bins[0]);
00278         for (int i = 0; i < bins[0]; i++)
00279         {
00280                 x_av[i] = (std::vector<double>(bins[1], 0));
00281                 x_av2[i] = (std::vector<double>(bins[1], 0));
00282                 sigma21[i] = (std::vector<double>(bins[1], 0));
00283                 sigma22[i] = (std::vector<double>(bins[1], 0));
00284         }
00285         for (int biny = 0; biny < bins[0]; biny++)
00286         {
00287                 for (int biny2 = 0; biny2 < bins[1]; biny2++)
00288                 {
00289                         norm = county[biny][biny2] > 0 ? county[biny][biny2] : 1;
00290                         x_av[biny][biny2] = sumx1[biny][biny2] / norm;
00291                         x_av2[biny][biny2] = sumx2[biny][biny2] / norm;
00292                         sigma21[biny][biny2] = sumx21[biny][biny2] / norm - x_av[biny][biny2] * x_av[biny][biny2];
00293                         sigma22[biny][biny2] = sumx22[biny][biny2] / norm - x_av2[biny][biny2] * x_av2[biny][biny2];
00294                 }
00295 
00296         }
00297 
00298         static std::string gridfilename = outputfile + ".grad";
00299         static std::string histfilename = outputfile + ".hist.grad";
00300         static std::string countfilename = outputfile + ".count";
00301         //static std::string pmffilename = outputfile + ".pmf";
00302 
00303         ofstream_namd ofile(gridfilename.c_str(),".old");
00304         ofstream_namd ofile_hist(histfilename.c_str(), std::ios::app);
00305         ofstream_namd ofile_count(countfilename.c_str(),".old");
00306 
00307         writehead(ofile);
00308         writehead(ofile_hist);
00309         writehead(ofile_count);
00310         //writecount(ofile_count);
00311         //ofstream_namd ofile_pmf(pmffilename.c_str(),".old");
00312 
00313         //ofile_pmf<<lowerboundary<<" 0.0"<<std::endl;
00314 
00315         //double pmf = 0, pmf_x = 0;
00316         double x1, x2, y1, y2, av1, av2, diff_av1, diff_av2, grad1, grad2;
00317         for (int binx = 0; binx < bins[0]; binx++)
00318         {
00319                 x1 = (binx + min[0] + 0.5) * width[0];
00320                 for (int binx2 = 0; binx2 < bins[1]; binx2++)
00321                 {
00322                         x2 = (binx2 + min[1] + 0.5) * width[1];
00323                         norm = 0;
00324                         av1 = 0;
00325                         av2 = 0;
00326                         diff_av1 = 0;
00327                         diff_av2 = 0;
00328 
00329                         for (int biny = 0; biny < bins[0]; biny++)
00330                         {
00331                                 y1 = (biny + min[0] + 0.5) * width[0];
00332                                 for (int biny2 = 0; biny2 < bins[1]; biny2++)
00333                                 {
00334                                         y2 = (biny2 + min[1] + 0.5) * width[1];
00335                                         norm += countall[binx * bins[0] + biny][binx2 * bins[1] + biny2];
00336 
00337                                         if (sigma21[biny][biny2]>0.00001 || sigma21[biny][biny2] < -0.00001)
00338                                         {
00339                                                 av1 += countall[binx * bins[0] + biny][binx2 * bins[1] + biny2] * (x1 - x_av[biny][biny2]) / sigma21[biny][biny2];
00340                                         }
00341                                         else if (countall[binx * bins[0] + biny][binx2 * bins[1] + biny2] > 1)
00342                                         {
00343                                                 // TODOC:\Users\ChinaFu\eabf\eabf2d.cpp|165|error: passing 'const eABF2D' as 'this' argument of 'int eABF2D::position(int, int, int, int)' discards qualifiers [-fpermissive]|, print error
00344                                         }
00345                                         diff_av1 += countall[binx * bins[0] + biny][binx2 * bins[1] + biny2] * (x1 - y1);
00346 
00347                                         if (sigma22[biny][biny2]>0.00001 || sigma22[biny][biny2] < -0.00001)
00348                                         {
00349                                                 av2 += countall[binx * bins[0] + biny][binx2 * bins[1] + biny2] * (x2 - x_av2[biny][biny2]) / sigma22[biny][biny2];
00350                                         }
00351                                         else if (countall[binx * bins[0] + biny][binx2 * bins[1] + biny2] > 1)
00352                                         {
00353                                                 // TODO, print error
00354                                         }
00355 
00356                                         diff_av2 += countall[binx * bins[0] + biny][binx2 * bins[1] + biny2] * (x2 - y2);
00357                                 }
00358                         }
00359 
00360                         diff_av1 /= (norm > 0 ? norm : 1);
00361                         diff_av2 /= (norm > 0 ? norm : 1);
00362 
00363                         av1 = eabffunc::BOLTZMANN * temperature * av1 / (norm > 0 ? norm : 1);
00364                         av2 = eabffunc::BOLTZMANN * temperature * av2 / (norm > 0 ? norm : 1);
00365 
00366                         grad1 = av1 - krestr[0] * diff_av1;
00367                         grad2 = av2 - krestr[1] * diff_av2;
00368 
00369                         if (x1 >= lowerboundary[0] && x1 < upperboundary[0] && x2 >= lowerboundary[1] && x2 < upperboundary[1])
00370                         {
00371                                 ofile << x1 << " " << x2 << " " << grad1 << " " << grad2 << std::endl;
00372                                 ofile_hist << x1 << " " << x2 << " " << grad1 << " " << grad2 << std::endl;
00373                                 ofile_count << x1 << " " << x2 << " " << norm << std::endl;
00374                                 //pmf += grad * width;
00375                                 //pmf_x = x + 0.5 * width;
00376                                 //ofile_pmf<<pmf_x<<" "<<pmf<<std::endl;
00377                                 // TODO, rescale the pmf file to prevent negetive free energy
00378                         }
00379                 }
00380                 if (x1 > lowerboundary[0] && x1 < upperboundary[0])
00381                 {
00382                         ofile << std::endl;
00383                         ofile_count << std::endl;
00384                 }
00385         }
00386         ofile.close();
00387         ofile_hist.close();
00388         ofile_count.close();
00389         //ofile_pmf.close();
00390         return true;
00391 }
00392 
00393 #endif // EABF2D

Generated on Tue Sep 19 01:17:12 2017 for NAMD by  doxygen 1.4.7