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 $\mathbf{s}$ can be used as ARG.
31 : The positions $\mathbf{s}_i$ and dimension $\mathbf{\sigma}_i$ of the umbrellas are read from file.
32 :
33 : $$
34 : \Delta u_{\mathbf{s}_i}(\mathbf{s})=\sum_j^{\text{dim}}\frac{([s]_j-[s_i]_j)^2}{2[\sigma_i]_j^2}\, .
35 : $$
36 :
37 : Notice that $\mathbf{\sigma}_i$ is diagonal, thus only one SIGMA per CV has to be specified for each umbrella.
38 : You can choose the umbrellas manually, or place them on a grid, or along a path, similar to [PATH](PATH.md).
39 : They must cover all the CV space that one wishes to sample.
40 :
41 : The first column of the umbrellas file is always ignored and must be called "time".
42 : You can also use as input file a STATE file from an earlier [OPES_METAD](OPES_METAD.md) run (or an [OPES_METAD_EXPLORE](OPES_METAD_EXPLORE.md) run, if you combine it with other ECVs).
43 :
44 : Similarly to [ECV_UMBRELLAS_LINE](ECV_UMBRELLAS_LINE.md), 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.
45 : You can also use BARRIER to set the maximum barrier height to be explored, and avoid huge biases at the beginning of your simulation.
46 : See also Appendix B of the paper cited below for more details on these last two options.
47 :
48 : ## Examples
49 :
50 : ```plumed
51 : #SETTINGS INPUTFILES=extras/Umbrellas.data
52 :
53 : cv1: DISTANCE ATOMS=1,2
54 : cv2: DISTANCE ATOMS=3,4
55 : cv3: DISTANCE ATOMS=4,1
56 : ecv: ECV_UMBRELLAS_FILE ...
57 : ARG=cv1,cv2,cv3
58 : FILE=extras/Umbrellas.data
59 : ADD_P0 BARRIER=70
60 : ...
61 : opes: OPES_EXPANDED ARG=ecv.* PACE=500
62 : PRINT FILE=COLVAR STRIDE=500 ARG=cv1,cv2,cv3,opes.bias
63 : ```
64 :
65 : */
66 : //+ENDPLUMEDOC
67 :
68 : class ECVumbrellasFile :
69 : public ExpansionCVs {
70 : private:
71 : double barrier_;
72 : unsigned P0_contribution_;
73 : bool lower_only_;
74 :
75 : std::vector< std::vector<double> > centers_;
76 : std::vector< std::vector<double> > sigmas_;
77 :
78 : std::vector< std::vector<double> > ECVs_;
79 : std::vector< std::vector<double> > derECVs_;
80 : void initECVs();
81 :
82 : public:
83 : explicit ECVumbrellasFile(const ActionOptions&);
84 : static void registerKeywords(Keywords& keys);
85 : void calculateECVs(const double *) override;
86 : const double * getPntrToECVs(unsigned) override;
87 : const double * getPntrToDerECVs(unsigned) override;
88 : std::vector<std::string> getLambdas() const override;
89 : void initECVs_observ(const std::vector<double>&,const unsigned,const unsigned) override;
90 : void initECVs_restart(const std::vector<std::string>&) override;
91 : };
92 :
93 : PLUMED_REGISTER_ACTION(ECVumbrellasFile,"ECV_UMBRELLAS_FILE")
94 :
95 5 : void ECVumbrellasFile::registerKeywords(Keywords& keys) {
96 5 : ExpansionCVs::registerKeywords(keys);
97 5 : keys.add("compulsory","FILE","the name of the file containing the umbrellas");
98 5 : keys.add("optional","BARRIER","a guess of the free energy barrier to be overcome (better to stay higher than lower)");
99 5 : keys.addFlag("ADD_P0",false,"add the unbiased Boltzmann distribution to the target distribution, to make sure to sample it");
100 5 : keys.addFlag("LOWER_HALF_ONLY",false,"use only the lower half of each umbrella potentials");
101 5 : keys.addDOI("10.1103/PhysRevX.10.041034");
102 5 : }
103 :
104 3 : ECVumbrellasFile::ECVumbrellasFile(const ActionOptions&ao):
105 : Action(ao),
106 3 : ExpansionCVs(ao) {
107 : //get number of CVs
108 : const unsigned ncv=getNumberOfArguments();
109 3 : centers_.resize(ncv);
110 3 : sigmas_.resize(ncv);
111 :
112 : //set P0_contribution_
113 3 : bool add_P0=false;
114 3 : parseFlag("ADD_P0",add_P0);
115 3 : if(add_P0) {
116 2 : P0_contribution_=1;
117 : } else {
118 1 : P0_contribution_=0;
119 : }
120 :
121 : //set barrier_
122 3 : barrier_=std::numeric_limits<double>::infinity();
123 3 : parse("BARRIER",barrier_);
124 6 : parseFlag("LOWER_HALF_ONLY",lower_only_);
125 :
126 : //set umbrellas
127 : std::string umbrellasFileName;
128 3 : parse("FILE",umbrellasFileName);
129 3 : IFile ifile;
130 3 : ifile.link(*this);
131 3 : if(ifile.FileExist(umbrellasFileName)) {
132 3 : log.printf(" reading from FILE '%s'\n",umbrellasFileName.c_str());
133 3 : ifile.open(umbrellasFileName);
134 3 : ifile.allowIgnoredFields();
135 : double time; //first field is ignored
136 1332 : while(ifile.scanField("time",time)) {
137 1989 : for(unsigned j=0; j<ncv; j++) {
138 : double centers_j;
139 1326 : ifile.scanField(getPntrToArgument(j)->getName(),centers_j);
140 1326 : centers_[j].push_back(centers_j); //this might be slow
141 : }
142 1989 : for(unsigned j=0; j<ncv; j++) {
143 : double sigmas_j;
144 2652 : ifile.scanField("sigma_"+getPntrToArgument(j)->getName(),sigmas_j);
145 1326 : sigmas_[j].push_back(sigmas_j);
146 : }
147 663 : ifile.scanField();
148 : }
149 : } else {
150 0 : plumed_merror("Umbrellas FILE '"+umbrellasFileName+"' not found");
151 : }
152 :
153 3 : checkRead();
154 :
155 : //extra consistency checks
156 3 : const unsigned sizeUmbrellas=centers_[0].size();
157 9 : for(unsigned j=0; j<ncv; j++) {
158 6 : plumed_massert(centers_[j].size()==sizeUmbrellas,"mismatch in the number of centers read from file");
159 6 : plumed_massert(sigmas_[j].size()==sizeUmbrellas,"mismatch in the number of sigmas read from file");
160 : }
161 :
162 : //set ECVs stuff
163 3 : totNumECVs_=sizeUmbrellas+P0_contribution_;
164 3 : ECVs_.resize(ncv,std::vector<double>(totNumECVs_));
165 3 : derECVs_.resize(ncv,std::vector<double>(totNumECVs_));
166 :
167 : //printing some info
168 3 : log.printf(" total number of umbrellas = %u\n",sizeUmbrellas);
169 3 : if(barrier_!=std::numeric_limits<double>::infinity()) {
170 1 : log.printf(" guess for free energy BARRIER = %g\n",barrier_);
171 : }
172 3 : if(P0_contribution_==1) {
173 2 : log.printf(" -- ADD_P0: the target includes also the unbiased probability itself\n");
174 : }
175 3 : if(lower_only_) {
176 1 : log.printf(" -- LOWER_HALF_ONLY: the ECVs are set to zero for values of the CV above the respective center\n");
177 : }
178 6 : }
179 :
180 131 : void ECVumbrellasFile::calculateECVs(const double * cv) {
181 131 : if(lower_only_) {
182 153 : for(unsigned j=0; j<getNumberOfArguments(); j++) {
183 22644 : for(unsigned k=P0_contribution_; k<totNumECVs_; k++) { //if ADD_P0, the first ECVs=0
184 22542 : const unsigned kk=k-P0_contribution_;
185 22542 : const double dist_jk=difference(j,centers_[j][kk],cv[j])/sigmas_[j][kk]; //PBC might be present
186 22542 : if(dist_jk>=0) {
187 11483 : ECVs_[j][k]=0;
188 11483 : derECVs_[j][k]=0;
189 : } else {
190 11059 : ECVs_[j][k]=0.5*std::pow(dist_jk,2);
191 11059 : derECVs_[j][k]=dist_jk/sigmas_[j][kk];
192 : }
193 : }
194 : }
195 : } else {
196 240 : for(unsigned j=0; j<getNumberOfArguments(); j++) {
197 35520 : for(unsigned k=P0_contribution_; k<totNumECVs_; k++) { //if ADD_P0, the first ECVs=0
198 35360 : const unsigned kk=k-P0_contribution_;
199 35360 : const double dist_jk=difference(j,centers_[j][kk],cv[j])/sigmas_[j][kk]; //PBC might be present
200 35360 : ECVs_[j][k]=0.5*std::pow(dist_jk,2);
201 35360 : derECVs_[j][k]=dist_jk/sigmas_[j][kk];
202 : }
203 : }
204 : }
205 131 : }
206 :
207 6 : const double * ECVumbrellasFile::getPntrToECVs(unsigned j) {
208 6 : plumed_massert(isReady_,"cannot access ECVs before initialization");
209 6 : plumed_massert(j<getNumberOfArguments(),getName()+" has fewer CVs");
210 6 : return &ECVs_[j][0];
211 : }
212 :
213 6 : const double * ECVumbrellasFile::getPntrToDerECVs(unsigned j) {
214 6 : plumed_massert(isReady_,"cannot access ECVs before initialization");
215 6 : plumed_massert(j<getNumberOfArguments(),getName()+" has fewer CVs");
216 6 : return &derECVs_[j][0];
217 : }
218 :
219 3 : std::vector<std::string> ECVumbrellasFile::getLambdas() const {
220 : //notice that sigmas are not considered!
221 3 : std::vector<std::string> lambdas(totNumECVs_);
222 3 : if(P0_contribution_==1) {
223 2 : std::ostringstream subs;
224 2 : subs<<"P0";
225 4 : for(unsigned j=1; j<getNumberOfArguments(); j++) {
226 2 : subs<<"_P0";
227 : }
228 2 : lambdas[0]=subs.str();
229 2 : }
230 666 : for(unsigned k=P0_contribution_; k<totNumECVs_; k++) {
231 663 : const unsigned kk=k-P0_contribution_;
232 663 : std::ostringstream subs;
233 663 : subs<<centers_[0][kk];
234 1326 : for(unsigned j=1; j<getNumberOfArguments(); j++) {
235 663 : subs<<"_"<<centers_[j][kk];
236 : }
237 663 : lambdas[k]=subs.str();
238 663 : }
239 3 : return lambdas;
240 0 : }
241 :
242 3 : void ECVumbrellasFile::initECVs() {
243 3 : plumed_massert(!isReady_,"initialization should not be called twice");
244 3 : isReady_=true;
245 3 : log.printf(" *%4u windows for %s\n",totNumECVs_,getName().c_str());
246 3 : }
247 :
248 2 : void ECVumbrellasFile::initECVs_observ(const std::vector<double>& all_obs_cvs,const unsigned ncv,const unsigned index_j) {
249 : //this non-linear exansion never uses automatic initialization
250 2 : initECVs();
251 2 : calculateECVs(&all_obs_cvs[index_j]); //use only first obs point
252 6 : for(unsigned j=0; j<getNumberOfArguments(); j++)
253 888 : for(unsigned k=P0_contribution_; k<totNumECVs_; k++) {
254 1680 : ECVs_[j][k]=std::min(barrier_/kbt_,ECVs_[j][k]);
255 : }
256 2 : }
257 :
258 1 : void ECVumbrellasFile::initECVs_restart(const std::vector<std::string>& lambdas) {
259 : std::size_t pos=0;
260 2 : for(unsigned j=0; j<getNumberOfArguments()-1; j++) {
261 1 : pos=lambdas[0].find("_",pos+1); //checking only lambdas[0] is hopefully enough
262 : }
263 1 : plumed_massert(pos<lambdas[0].length(),"this should not happen, fewer '_' than expected in "+getName());
264 1 : pos=lambdas[0].find("_",pos+1);
265 1 : plumed_massert(pos>lambdas[0].length(),"this should not happen, more '_' than expected in "+getName());
266 :
267 1 : std::vector<std::string> myLambdas=getLambdas();
268 1 : plumed_massert(myLambdas.size()==lambdas.size(),"RESTART - mismatch in number of "+getName());
269 1 : plumed_massert(std::equal(myLambdas.begin(),myLambdas.end(),lambdas.begin()),"RESTART - mismatch in lambda values of "+getName());
270 :
271 1 : initECVs();
272 1 : }
273 :
274 : }
275 : }
|
