Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2020-2021 of Michele Invernizzi.
3 :
4 : This file is part of the OPES plumed module.
5 :
6 : The OPES plumed module is free software: you can redistribute it and/or modify
7 : it under the terms of the GNU Lesser General Public License as published by
8 : the Free Software Foundation, either version 3 of the License, or
9 : (at your option) any later version.
10 :
11 : The OPES plumed module is distributed in the hope that it will be useful,
12 : but WITHOUT ANY WARRANTY; without even the implied warranty of
13 : MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 : GNU Lesser General Public License for more details.
15 :
16 : You should have received a copy of the GNU Lesser General Public License
17 : along with plumed. If not, see <http://www.gnu.org/licenses/>.
18 : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
19 : #include "ExpansionCVs.h"
20 : #include "core/ActionRegister.h"
21 : #include "tools/File.h"
22 :
23 : namespace PLMD {
24 : namespace opes {
25 :
26 : //+PLUMEDOC OPES_EXPANSION_CV ECV_UMBRELLAS_FILE
27 : /*
28 : Target a multiumbrella ensemble, by combining systems each with a parabolic bias potential at a different location.
29 :
30 : Any set of collective variables \f$\mathbf{s}\f$ can be used as ARG.
31 : The positions \f$\mathbf{s}_i\f$ and dimension \f$\mathbf{\sigma}_i\f$ of the umbrellas are read from file.
32 : \f[
33 : \Delta u_{\mathbf{s}_i}(\mathbf{s})=\sum_j^{\text{dim}}\frac{([s]_j-[s_i]_j)^2}{2[\sigma_i]_j^2}\, .
34 : \f]
35 : Notice that \f$\mathbf{\sigma}_i\f$ is diagonal, thus only one SIGMA per CV has to be specified for each umbrella.
36 : You can choose the umbrellas manually, or place them on a grid, or along a path, similar to \ref PATH.
37 : They must cover all the CV space that one wishes to sample.
38 :
39 : The first column of the umbrellas file is always ignored and must be called "time".
40 : You can also use as input file a STATE file from an earlier \ref OPES_METAD run (or an \ref OPES_MEAD_EXPLORE run, if you combine it with other ECVs).
41 :
42 : Similarly to \ref ECV_UMBRELLAS_LINE, you should set the flag ADD_P0 if you think your umbrellas might not properly cover all the CV region relevant for the unbiased distribution.
43 : You can also use BARRIER to set the maximum barrier height to be explored, and avoid huge biases at the beginning of your simulation.
44 : See also Appendix B of Ref.\cite Invernizzi2020unified for more details on these last two options.
45 :
46 : \par Examples
47 :
48 : \plumedfile
49 : cv1: DISTANCE ATOMS=1,2
50 : cv2: DISTANCE ATOMS=3,4
51 : cv3: DISTANCE ATOMS=4,1
52 : ecv: ECV_UMBRELLAS_FILE ARG=cv1,cv2,cv3 FILE=Umbrellas.data ADD_P0 BARRIER=70
53 : opes: OPES_EXPANDED ARG=ecv.* PACE=500
54 : PRINT FILE=COLVAR STRIDE=500 ARG=cv1,cv2,cv3,opes.bias
55 : \endplumedfile
56 :
57 : The umbrellas file might look like this:
58 : \auxfile{Umbrellas.data}
59 : #! FIELDS time cv1 cv2 cv3 sigma_cv1 sigma_cv2 sigma_cv3
60 : 1 1.17958 2.93697 1.06109 0.19707 0.28275 0.32427
61 : 2 2.04023 2.69714 1.84770 0.22307 0.25933 0.31783
62 : 3 1.99693 1.10299 1.13351 0.19517 0.26260 0.37427
63 : 4 1.15954 1.37447 2.25975 0.20096 0.27168 0.33353
64 : 5 1.10126 2.45936 2.40260 0.19747 0.24215 0.35523
65 : \endauxfile
66 :
67 : */
68 : //+ENDPLUMEDOC
69 :
70 : class ECVumbrellasFile :
71 : public ExpansionCVs {
72 : private:
73 : double barrier_;
74 : unsigned P0_contribution_;
75 : bool lower_only_;
76 :
77 : std::vector< std::vector<double> > centers_;
78 : std::vector< std::vector<double> > sigmas_;
79 :
80 : std::vector< std::vector<double> > ECVs_;
81 : std::vector< std::vector<double> > derECVs_;
82 : void initECVs();
83 :
84 : public:
85 : explicit ECVumbrellasFile(const ActionOptions&);
86 : static void registerKeywords(Keywords& keys);
87 : void calculateECVs(const double *) override;
88 : const double * getPntrToECVs(unsigned) override;
89 : const double * getPntrToDerECVs(unsigned) override;
90 : std::vector<std::string> getLambdas() const override;
91 : void initECVs_observ(const std::vector<double>&,const unsigned,const unsigned) override;
92 : void initECVs_restart(const std::vector<std::string>&) override;
93 : };
94 :
95 : PLUMED_REGISTER_ACTION(ECVumbrellasFile,"ECV_UMBRELLAS_FILE")
96 :
97 5 : void ECVumbrellasFile::registerKeywords(Keywords& keys) {
98 5 : ExpansionCVs::registerKeywords(keys);
99 5 : keys.add("compulsory","FILE","the name of the file containing the umbrellas");
100 5 : keys.add("optional","BARRIER","a guess of the free energy barrier to be overcome (better to stay higher than lower)");
101 5 : keys.addFlag("ADD_P0",false,"add the unbiased Boltzmann distribution to the target distribution, to make sure to sample it");
102 5 : keys.addFlag("LOWER_HALF_ONLY",false,"use only the lower half of each umbrella potentials");
103 5 : }
104 :
105 3 : ECVumbrellasFile::ECVumbrellasFile(const ActionOptions&ao):
106 : Action(ao),
107 3 : ExpansionCVs(ao) {
108 : //get number of CVs
109 : const unsigned ncv=getNumberOfArguments();
110 : centers_.resize(ncv);
111 3 : sigmas_.resize(ncv);
112 :
113 : //set P0_contribution_
114 3 : bool add_P0=false;
115 3 : parseFlag("ADD_P0",add_P0);
116 3 : if(add_P0) {
117 2 : P0_contribution_=1;
118 : } else {
119 1 : P0_contribution_=0;
120 : }
121 :
122 : //set barrier_
123 3 : barrier_=std::numeric_limits<double>::infinity();
124 3 : parse("BARRIER",barrier_);
125 6 : parseFlag("LOWER_HALF_ONLY",lower_only_);
126 :
127 : //set umbrellas
128 : std::string umbrellasFileName;
129 3 : parse("FILE",umbrellasFileName);
130 3 : IFile ifile;
131 3 : ifile.link(*this);
132 3 : if(ifile.FileExist(umbrellasFileName)) {
133 3 : log.printf(" reading from FILE '%s'\n",umbrellasFileName.c_str());
134 3 : ifile.open(umbrellasFileName);
135 3 : ifile.allowIgnoredFields();
136 : double time; //first field is ignored
137 1332 : while(ifile.scanField("time",time)) {
138 1989 : for(unsigned j=0; j<ncv; j++) {
139 : double centers_j;
140 1326 : ifile.scanField(getPntrToArgument(j)->getName(),centers_j);
141 1326 : centers_[j].push_back(centers_j); //this might be slow
142 : }
143 1989 : for(unsigned j=0; j<ncv; j++) {
144 : double sigmas_j;
145 2652 : ifile.scanField("sigma_"+getPntrToArgument(j)->getName(),sigmas_j);
146 1326 : sigmas_[j].push_back(sigmas_j);
147 : }
148 663 : ifile.scanField();
149 : }
150 : } else {
151 0 : plumed_merror("Umbrellas FILE '"+umbrellasFileName+"' not found");
152 : }
153 :
154 3 : checkRead();
155 :
156 : //extra consistency checks
157 3 : const unsigned sizeUmbrellas=centers_[0].size();
158 9 : for(unsigned j=0; j<ncv; j++) {
159 6 : plumed_massert(centers_[j].size()==sizeUmbrellas,"mismatch in the number of centers read from file");
160 6 : plumed_massert(sigmas_[j].size()==sizeUmbrellas,"mismatch in the number of sigmas read from file");
161 : }
162 :
163 : //set ECVs stuff
164 3 : totNumECVs_=sizeUmbrellas+P0_contribution_;
165 3 : ECVs_.resize(ncv,std::vector<double>(totNumECVs_));
166 3 : derECVs_.resize(ncv,std::vector<double>(totNumECVs_));
167 :
168 : //printing some info
169 3 : log.printf(" total number of umbrellas = %u\n",sizeUmbrellas);
170 3 : if(barrier_!=std::numeric_limits<double>::infinity()) {
171 1 : log.printf(" guess for free energy BARRIER = %g\n",barrier_);
172 : }
173 3 : if(P0_contribution_==1) {
174 2 : log.printf(" -- ADD_P0: the target includes also the unbiased probability itself\n");
175 : }
176 3 : if(lower_only_) {
177 1 : log.printf(" -- LOWER_HALF_ONLY: the ECVs are set to zero for values of the CV above the respective center\n");
178 : }
179 6 : }
180 :
181 131 : void ECVumbrellasFile::calculateECVs(const double * cv) {
182 131 : if(lower_only_) {
183 153 : for(unsigned j=0; j<getNumberOfArguments(); j++) {
184 22644 : for(unsigned k=P0_contribution_; k<totNumECVs_; k++) { //if ADD_P0, the first ECVs=0
185 22542 : const unsigned kk=k-P0_contribution_;
186 22542 : const double dist_jk=difference(j,centers_[j][kk],cv[j])/sigmas_[j][kk]; //PBC might be present
187 22542 : if(dist_jk>=0) {
188 11483 : ECVs_[j][k]=0;
189 11483 : derECVs_[j][k]=0;
190 : } else {
191 11059 : ECVs_[j][k]=0.5*std::pow(dist_jk,2);
192 11059 : derECVs_[j][k]=dist_jk/sigmas_[j][kk];
193 : }
194 : }
195 : }
196 : } else {
197 240 : for(unsigned j=0; j<getNumberOfArguments(); j++) {
198 35520 : for(unsigned k=P0_contribution_; k<totNumECVs_; k++) { //if ADD_P0, the first ECVs=0
199 35360 : const unsigned kk=k-P0_contribution_;
200 35360 : const double dist_jk=difference(j,centers_[j][kk],cv[j])/sigmas_[j][kk]; //PBC might be present
201 35360 : ECVs_[j][k]=0.5*std::pow(dist_jk,2);
202 35360 : derECVs_[j][k]=dist_jk/sigmas_[j][kk];
203 : }
204 : }
205 : }
206 131 : }
207 :
208 6 : const double * ECVumbrellasFile::getPntrToECVs(unsigned j) {
209 6 : plumed_massert(isReady_,"cannot access ECVs before initialization");
210 6 : plumed_massert(j<getNumberOfArguments(),getName()+" has fewer CVs");
211 6 : return &ECVs_[j][0];
212 : }
213 :
214 6 : const double * ECVumbrellasFile::getPntrToDerECVs(unsigned j) {
215 6 : plumed_massert(isReady_,"cannot access ECVs before initialization");
216 6 : plumed_massert(j<getNumberOfArguments(),getName()+" has fewer CVs");
217 6 : return &derECVs_[j][0];
218 : }
219 :
220 3 : std::vector<std::string> ECVumbrellasFile::getLambdas() const {
221 : //notice that sigmas are not considered!
222 3 : std::vector<std::string> lambdas(totNumECVs_);
223 3 : if(P0_contribution_==1) {
224 2 : std::ostringstream subs;
225 2 : subs<<"P0";
226 4 : for(unsigned j=1; j<getNumberOfArguments(); j++) {
227 2 : subs<<"_P0";
228 : }
229 2 : lambdas[0]=subs.str();
230 2 : }
231 666 : for(unsigned k=P0_contribution_; k<totNumECVs_; k++) {
232 663 : const unsigned kk=k-P0_contribution_;
233 663 : std::ostringstream subs;
234 663 : subs<<centers_[0][kk];
235 1326 : for(unsigned j=1; j<getNumberOfArguments(); j++) {
236 663 : subs<<"_"<<centers_[j][kk];
237 : }
238 663 : lambdas[k]=subs.str();
239 663 : }
240 3 : return lambdas;
241 0 : }
242 :
243 3 : void ECVumbrellasFile::initECVs() {
244 3 : plumed_massert(!isReady_,"initialization should not be called twice");
245 3 : isReady_=true;
246 3 : log.printf(" *%4u windows for %s\n",totNumECVs_,getName().c_str());
247 3 : }
248 :
249 2 : void ECVumbrellasFile::initECVs_observ(const std::vector<double>& all_obs_cvs,const unsigned ncv,const unsigned index_j) {
250 : //this non-linear exansion never uses automatic initialization
251 2 : initECVs();
252 2 : calculateECVs(&all_obs_cvs[index_j]); //use only first obs point
253 6 : for(unsigned j=0; j<getNumberOfArguments(); j++)
254 888 : for(unsigned k=P0_contribution_; k<totNumECVs_; k++) {
255 1680 : ECVs_[j][k]=std::min(barrier_/kbt_,ECVs_[j][k]);
256 : }
257 2 : }
258 :
259 1 : void ECVumbrellasFile::initECVs_restart(const std::vector<std::string>& lambdas) {
260 : std::size_t pos=0;
261 2 : for(unsigned j=0; j<getNumberOfArguments()-1; j++) {
262 1 : pos=lambdas[0].find("_",pos+1); //checking only lambdas[0] is hopefully enough
263 : }
264 1 : plumed_massert(pos<lambdas[0].length(),"this should not happen, fewer '_' than expected in "+getName());
265 1 : pos=lambdas[0].find("_",pos+1);
266 1 : plumed_massert(pos>lambdas[0].length(),"this should not happen, more '_' than expected in "+getName());
267 :
268 1 : std::vector<std::string> myLambdas=getLambdas();
269 1 : plumed_massert(myLambdas.size()==lambdas.size(),"RESTART - mismatch in number of "+getName());
270 1 : plumed_massert(std::equal(myLambdas.begin(),myLambdas.end(),lambdas.begin()),"RESTART - mismatch in lambda values of "+getName());
271 :
272 1 : initECVs();
273 1 : }
274 :
275 : }
276 : }
|
