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.use("ARG");
100 10 : keys.add("compulsory","FILE","the name of the file containing the umbrellas");
101 10 : keys.add("optional","BARRIER","a guess of the free energy barrier to be overcome (better to stay higher than lower)");
102 10 : keys.addFlag("ADD_P0",false,"add the unbiased Boltzmann distribution to the target distribution, to make sure to sample it");
103 10 : keys.addFlag("LOWER_HALF_ONLY",false,"use only the lower half of each umbrella potentials");
104 5 : }
105 :
106 3 : ECVumbrellasFile::ECVumbrellasFile(const ActionOptions&ao):
107 : Action(ao),
108 3 : ExpansionCVs(ao) {
109 : //get number of CVs
110 : const unsigned ncv=getNumberOfArguments();
111 3 : centers_.resize(ncv);
112 3 : sigmas_.resize(ncv);
113 :
114 : //set P0_contribution_
115 3 : bool add_P0=false;
116 3 : parseFlag("ADD_P0",add_P0);
117 3 : if(add_P0) {
118 2 : P0_contribution_=1;
119 : } else {
120 1 : P0_contribution_=0;
121 : }
122 :
123 : //set barrier_
124 3 : barrier_=std::numeric_limits<double>::infinity();
125 3 : parse("BARRIER",barrier_);
126 6 : parseFlag("LOWER_HALF_ONLY",lower_only_);
127 :
128 : //set umbrellas
129 : std::string umbrellasFileName;
130 3 : parse("FILE",umbrellasFileName);
131 3 : IFile ifile;
132 3 : ifile.link(*this);
133 3 : if(ifile.FileExist(umbrellasFileName)) {
134 3 : log.printf(" reading from FILE '%s'\n",umbrellasFileName.c_str());
135 3 : ifile.open(umbrellasFileName);
136 3 : ifile.allowIgnoredFields();
137 : double time; //first field is ignored
138 1332 : while(ifile.scanField("time",time)) {
139 1989 : for(unsigned j=0; j<ncv; j++) {
140 : double centers_j;
141 1326 : ifile.scanField(getPntrToArgument(j)->getName(),centers_j);
142 1326 : centers_[j].push_back(centers_j); //this might be slow
143 : }
144 1989 : for(unsigned j=0; j<ncv; j++) {
145 : double sigmas_j;
146 2652 : ifile.scanField("sigma_"+getPntrToArgument(j)->getName(),sigmas_j);
147 1326 : sigmas_[j].push_back(sigmas_j);
148 : }
149 663 : ifile.scanField();
150 : }
151 : } else {
152 0 : plumed_merror("Umbrellas FILE '"+umbrellasFileName+"' not found");
153 : }
154 :
155 3 : checkRead();
156 :
157 : //extra consistency checks
158 3 : const unsigned sizeUmbrellas=centers_[0].size();
159 9 : for(unsigned j=0; j<ncv; j++) {
160 6 : plumed_massert(centers_[j].size()==sizeUmbrellas,"mismatch in the number of centers read from file");
161 6 : plumed_massert(sigmas_[j].size()==sizeUmbrellas,"mismatch in the number of sigmas read from file");
162 : }
163 :
164 : //set ECVs stuff
165 3 : totNumECVs_=sizeUmbrellas+P0_contribution_;
166 3 : ECVs_.resize(ncv,std::vector<double>(totNumECVs_));
167 3 : derECVs_.resize(ncv,std::vector<double>(totNumECVs_));
168 :
169 : //printing some info
170 3 : log.printf(" total number of umbrellas = %u\n",sizeUmbrellas);
171 3 : if(barrier_!=std::numeric_limits<double>::infinity()) {
172 1 : log.printf(" guess for free energy BARRIER = %g\n",barrier_);
173 : }
174 3 : if(P0_contribution_==1) {
175 2 : log.printf(" -- ADD_P0: the target includes also the unbiased probability itself\n");
176 : }
177 3 : if(lower_only_) {
178 1 : log.printf(" -- LOWER_HALF_ONLY: the ECVs are set to zero for values of the CV above the respective center\n");
179 : }
180 6 : }
181 :
182 131 : void ECVumbrellasFile::calculateECVs(const double * cv) {
183 131 : if(lower_only_) {
184 153 : for(unsigned j=0; j<getNumberOfArguments(); j++) {
185 22644 : for(unsigned k=P0_contribution_; k<totNumECVs_; k++) { //if ADD_P0, the first ECVs=0
186 22542 : const unsigned kk=k-P0_contribution_;
187 22542 : const double dist_jk=difference(j,centers_[j][kk],cv[j])/sigmas_[j][kk]; //PBC might be present
188 22542 : if(dist_jk>=0) {
189 11483 : ECVs_[j][k]=0;
190 11483 : derECVs_[j][k]=0;
191 : } else {
192 11059 : ECVs_[j][k]=0.5*std::pow(dist_jk,2);
193 11059 : derECVs_[j][k]=dist_jk/sigmas_[j][kk];
194 : }
195 : }
196 : }
197 : } else {
198 240 : for(unsigned j=0; j<getNumberOfArguments(); j++) {
199 35520 : for(unsigned k=P0_contribution_; k<totNumECVs_; k++) { //if ADD_P0, the first ECVs=0
200 35360 : const unsigned kk=k-P0_contribution_;
201 35360 : const double dist_jk=difference(j,centers_[j][kk],cv[j])/sigmas_[j][kk]; //PBC might be present
202 35360 : ECVs_[j][k]=0.5*std::pow(dist_jk,2);
203 35360 : derECVs_[j][k]=dist_jk/sigmas_[j][kk];
204 : }
205 : }
206 : }
207 131 : }
208 :
209 6 : const double * ECVumbrellasFile::getPntrToECVs(unsigned j) {
210 6 : plumed_massert(isReady_,"cannot access ECVs before initialization");
211 6 : plumed_massert(j<getNumberOfArguments(),getName()+" has fewer CVs");
212 6 : return &ECVs_[j][0];
213 : }
214 :
215 6 : const double * ECVumbrellasFile::getPntrToDerECVs(unsigned j) {
216 6 : plumed_massert(isReady_,"cannot access ECVs before initialization");
217 6 : plumed_massert(j<getNumberOfArguments(),getName()+" has fewer CVs");
218 6 : return &derECVs_[j][0];
219 : }
220 :
221 3 : std::vector<std::string> ECVumbrellasFile::getLambdas() const {
222 : //notice that sigmas are not considered!
223 3 : std::vector<std::string> lambdas(totNumECVs_);
224 3 : if(P0_contribution_==1) {
225 2 : std::ostringstream subs;
226 2 : subs<<"P0";
227 4 : for(unsigned j=1; j<getNumberOfArguments(); j++) {
228 2 : subs<<"_P0";
229 : }
230 2 : lambdas[0]=subs.str();
231 2 : }
232 666 : for(unsigned k=P0_contribution_; k<totNumECVs_; k++) {
233 663 : const unsigned kk=k-P0_contribution_;
234 663 : std::ostringstream subs;
235 663 : subs<<centers_[0][kk];
236 1326 : for(unsigned j=1; j<getNumberOfArguments(); j++) {
237 663 : subs<<"_"<<centers_[j][kk];
238 : }
239 663 : lambdas[k]=subs.str();
240 663 : }
241 3 : return lambdas;
242 0 : }
243 :
244 3 : void ECVumbrellasFile::initECVs() {
245 3 : plumed_massert(!isReady_,"initialization should not be called twice");
246 3 : isReady_=true;
247 3 : log.printf(" *%4u windows for %s\n",totNumECVs_,getName().c_str());
248 3 : }
249 :
250 2 : void ECVumbrellasFile::initECVs_observ(const std::vector<double>& all_obs_cvs,const unsigned ncv,const unsigned index_j) {
251 : //this non-linear exansion never uses automatic initialization
252 2 : initECVs();
253 2 : calculateECVs(&all_obs_cvs[index_j]); //use only first obs point
254 6 : for(unsigned j=0; j<getNumberOfArguments(); j++)
255 888 : for(unsigned k=P0_contribution_; k<totNumECVs_; k++) {
256 1680 : ECVs_[j][k]=std::min(barrier_/kbt_,ECVs_[j][k]);
257 : }
258 2 : }
259 :
260 1 : void ECVumbrellasFile::initECVs_restart(const std::vector<std::string>& lambdas) {
261 : std::size_t pos=0;
262 2 : for(unsigned j=0; j<getNumberOfArguments()-1; j++) {
263 1 : pos=lambdas[0].find("_",pos+1); //checking only lambdas[0] is hopefully enough
264 : }
265 1 : plumed_massert(pos<lambdas[0].length(),"this should not happen, fewer '_' than expected in "+getName());
266 1 : pos=lambdas[0].find("_",pos+1);
267 1 : plumed_massert(pos>lambdas[0].length(),"this should not happen, more '_' than expected in "+getName());
268 :
269 1 : std::vector<std::string> myLambdas=getLambdas();
270 1 : plumed_massert(myLambdas.size()==lambdas.size(),"RESTART - mismatch in number of "+getName());
271 1 : plumed_massert(std::equal(myLambdas.begin(),myLambdas.end(),lambdas.begin()),"RESTART - mismatch in lambda values of "+getName());
272 :
273 1 : initECVs();
274 1 : }
275 :
276 : }
277 : }
|
