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 :
22 : namespace PLMD {
23 : namespace opes {
24 :
25 : //+PLUMEDOC OPES_EXPANSION_CV ECV_UMBRELLAS_LINE
26 : /*
27 : Target a multiumbrella ensemble, by combining systems each with a parabolic bias potential at a different location.
28 :
29 : Any set of collective variables \f$\mathbf{s}\f$ can be used as ARG.
30 : \f[
31 : \Delta u_{\mathbf{s}_i}(\mathbf{s})=\sum_j^{\text{dim}}\frac{([s]_j-[s_i]_j)^2}{2\sigma^2}\, .
32 : \f]
33 : The Gaussian umbrellas are placed along a line, from CV_MIN to CV_MAX.
34 : The umbrellas are placed at a distance SIGMA*SPACING from each other, if you need more flexibility use \ref ECV_UMBRELLAS_FILE.
35 : The unbiased fluctuations in the basin usually are a reasonable guess for the value of SIGMA.
36 : Typically, a SPACING greater than 1 can lead to faster convergence, by reducing the total number of umbrellas.
37 : The umbrellas can be multidimensional, but the CVs dimensions should be rescaled since a single SIGMA must be used.
38 :
39 : The keyword BARRIER can be helpful to avoid breaking your system due to a too strong initial bias.
40 : If you think the placed umbrellas will not cover the whole unbiased probability distribution you should add it explicitly to the target, with the flag ADD_P0, for more robust convergence.
41 : See also Appendix B of Ref.\cite Invernizzi2020unified for more details on these last two options.
42 :
43 : The flag LOWER_HALF_ONLY modifies the ECVs so that they are set to zero when \f$\mathbf{s}>\mathbf{s}_i\f$, as in \ref LOWER_WALLS.
44 : This can be useful e.g. when the CV used is the \ref ENERGY and one wants to sample a broad range of high energy values, similar to \ref ECV_MULTITHERMAL but with a flat energy distribution as target.
45 : By pushing only from below one can avoid too extreme forces that could crash the simulation.
46 :
47 : \par Examples
48 :
49 : \plumedfile
50 : cv: DISTANCE ATOMS=1,2
51 : ecv: ECV_UMBRELLAS_LINE ARG=cv CV_MIN=1.2 CV_MAX=4.3 SIGMA=0.5 SPACING=1.5
52 : opes: OPES_EXPANDED ARG=ecv.* PACE=500
53 : \endplumedfile
54 :
55 : It is also possible to combine different ECV_UMBRELLAS_LINE to build a grid of CV values that will be sampled.
56 : For example the following code will sample a whole 2D region of cv1 and cv2.
57 :
58 : \plumedfile
59 : cv1: DISTANCE ATOMS=1,2
60 : ecv2: ECV_UMBRELLAS_LINE ARG=cv1 CV_MIN=1.2 CV_MAX=4.3 SIGMA=0.5
61 :
62 : cv2: DISTANCE ATOMS=3,4
63 : ecv1: ECV_UMBRELLAS_LINE ARG=cv2 CV_MIN=13.8 CV_MAX=21.4 SIGMA=4.3
64 :
65 : opes: OPES_EXPANDED ARG=ecv1.*,ecv2.* PACE=500
66 : \endplumedfile
67 :
68 : */
69 : //+ENDPLUMEDOC
70 :
71 : class ECVumbrellasLine :
72 : public ExpansionCVs
73 : {
74 : private:
75 : double barrier_;
76 : unsigned P0_contribution_;
77 : bool lower_only_;
78 :
79 : std::vector< std::vector<double> > centers_;
80 : double sigma_;
81 :
82 : std::vector< std::vector<double> > ECVs_;
83 : std::vector< std::vector<double> > derECVs_;
84 : void initECVs();
85 :
86 : public:
87 : explicit ECVumbrellasLine(const ActionOptions&);
88 : static void registerKeywords(Keywords& keys);
89 : void calculateECVs(const double *) override;
90 : const double * getPntrToECVs(unsigned) override;
91 : const double * getPntrToDerECVs(unsigned) override;
92 : std::vector<std::string> getLambdas() const override;
93 : void initECVs_observ(const std::vector<double>&,const unsigned,const unsigned) override;
94 : void initECVs_restart(const std::vector<std::string>&) override;
95 : };
96 :
97 : PLUMED_REGISTER_ACTION(ECVumbrellasLine,"ECV_UMBRELLAS_LINE")
98 :
99 10 : void ECVumbrellasLine::registerKeywords(Keywords& keys)
100 : {
101 10 : ExpansionCVs::registerKeywords(keys);
102 20 : keys.add("compulsory","CV_MIN","the minimum of the CV range to be explored");
103 20 : keys.add("compulsory","CV_MAX","the maximum of the CV range to be explored");
104 20 : keys.add("compulsory","SIGMA","sigma of the umbrella Gaussians");
105 20 : keys.add("compulsory","SPACING","1","the distance between umbrellas, in units of SIGMA");
106 20 : keys.add("optional","BARRIER","a guess of the free energy barrier to be overcome (better to stay higher than lower)");
107 20 : keys.addFlag("ADD_P0",false,"add the unbiased Boltzmann distribution to the target distribution, to make sure to sample it");
108 20 : keys.addFlag("LOWER_HALF_ONLY",false,"use only the lower half of each umbrella potentials");
109 10 : }
110 :
111 8 : ECVumbrellasLine::ECVumbrellasLine(const ActionOptions&ao):
112 : Action(ao),
113 8 : ExpansionCVs(ao)
114 : {
115 : //set P0_contribution_
116 8 : bool add_P0=false;
117 8 : parseFlag("ADD_P0",add_P0);
118 8 : if(add_P0)
119 2 : P0_contribution_=1;
120 : else
121 6 : P0_contribution_=0;
122 :
123 : //set barrier_
124 8 : barrier_=std::numeric_limits<double>::infinity();
125 8 : parse("BARRIER",barrier_);
126 8 : parseFlag("LOWER_HALF_ONLY",lower_only_);
127 :
128 : //set umbrellas
129 16 : parse("SIGMA",sigma_);
130 : std::vector<double> cv_min;
131 : std::vector<double> cv_max;
132 8 : parseVector("CV_MIN",cv_min);
133 16 : parseVector("CV_MAX",cv_max);
134 8 : plumed_massert(cv_min.size()==getNumberOfArguments(),"wrong number of CV_MINs");
135 8 : plumed_massert(cv_max.size()==getNumberOfArguments(),"wrong number of CV_MAXs");
136 : double spacing;
137 8 : parse("SPACING",spacing);
138 : double length=0;
139 18 : for(unsigned j=0; j<getNumberOfArguments(); j++)
140 10 : length+=std::pow(cv_max[j]-cv_min[j],2);
141 8 : length=std::sqrt(length);
142 8 : unsigned sizeUmbrellas=1+std::round(length/(sigma_*spacing));
143 8 : centers_.resize(getNumberOfArguments()); //centers_[cv][umbrellas]
144 : unsigned full_period=0;
145 18 : for(unsigned j=0; j<getNumberOfArguments(); j++)
146 : {
147 10 : centers_[j].resize(sizeUmbrellas);
148 10 : if(sizeUmbrellas>1)
149 140 : for(unsigned k=0; k<sizeUmbrellas; k++)
150 130 : centers_[j][k]=cv_min[j]+k*(cv_max[j]-cv_min[j])/(sizeUmbrellas-1);
151 : else
152 0 : centers_[j][0]=(cv_min[j]+cv_max[j])/2.;
153 10 : if(getPntrToArgument(j)->isPeriodic())
154 : {
155 : double min,max;
156 : std::string min_str,max_str;
157 10 : getPntrToArgument(j)->getDomain(min,max);
158 10 : getPntrToArgument(j)->getDomain(min_str,max_str);
159 10 : plumed_massert(cv_min[j]>=min,"ARG "+std::to_string(j)+": CV_MIN cannot be smaller than the periodic bound "+min_str);
160 10 : plumed_massert(cv_max[j]<=max,"ARG "+std::to_string(j)+": CV_MAX cannot be greater than the periodic bound "+max_str);
161 10 : if(cv_min[j]==min && cv_max[j]==max)
162 6 : full_period++;
163 : }
164 : }
165 8 : if(full_period==getNumberOfArguments() && sizeUmbrellas>1) //first and last are the same point
166 : {
167 6 : sizeUmbrellas--;
168 12 : for(unsigned j=0; j<getNumberOfArguments(); j++)
169 6 : centers_[j].pop_back();
170 : }
171 :
172 8 : checkRead();
173 :
174 : //set ECVs stuff
175 8 : totNumECVs_=sizeUmbrellas+P0_contribution_;
176 8 : ECVs_.resize(getNumberOfArguments(),std::vector<double>(totNumECVs_));
177 8 : derECVs_.resize(getNumberOfArguments(),std::vector<double>(totNumECVs_));
178 :
179 : //printing some info
180 8 : log.printf(" total number of umbrellas = %u\n",sizeUmbrellas);
181 8 : log.printf(" with SIGMA = %g\n",sigma_);
182 8 : log.printf(" and SPACING = %g\n",spacing);
183 8 : if(barrier_!=std::numeric_limits<double>::infinity())
184 2 : log.printf(" guess for free energy BARRIER = %g\n",barrier_);
185 8 : if(P0_contribution_==1)
186 2 : log.printf(" -- ADD_P0: the target includes also the unbiased probability itself\n");
187 8 : if(lower_only_)
188 1 : log.printf(" -- LOWER_HALF_ONLY: the ECVs are set to zero for values of the CV above the respective center\n");
189 8 : }
190 :
191 433 : void ECVumbrellasLine::calculateECVs(const double * cv)
192 : {
193 433 : if(lower_only_)
194 : {
195 153 : for(unsigned j=0; j<getNumberOfArguments(); j++)
196 : {
197 2040 : for(unsigned k=P0_contribution_; k<totNumECVs_; k++) //if ADD_P0, the first ECVs=0
198 : {
199 1938 : const unsigned kk=k-P0_contribution_;
200 1938 : const double dist_jk=difference(j,centers_[j][kk],cv[j])/sigma_; //PBC might be present
201 1938 : if(dist_jk>=0)
202 : {
203 933 : ECVs_[j][k]=0;
204 933 : derECVs_[j][k]=0;
205 : }
206 : else
207 : {
208 1005 : ECVs_[j][k]=0.5*std::pow(dist_jk,2);
209 1005 : derECVs_[j][k]=dist_jk/sigma_;
210 : }
211 : }
212 : }
213 : }
214 : else
215 : {
216 804 : for(unsigned j=0; j<getNumberOfArguments(); j++)
217 : {
218 4678 : for(unsigned k=P0_contribution_; k<totNumECVs_; k++) //if ADD_P0, the first ECVs=0
219 : {
220 4256 : const unsigned kk=k-P0_contribution_;
221 4256 : const double dist_jk=difference(j,centers_[j][kk],cv[j])/sigma_; //PBC might be present
222 4256 : ECVs_[j][k]=0.5*std::pow(dist_jk,2);
223 4256 : derECVs_[j][k]=dist_jk/sigma_;
224 : }
225 : }
226 : }
227 433 : }
228 :
229 10 : const double * ECVumbrellasLine::getPntrToECVs(unsigned j)
230 : {
231 10 : plumed_massert(isReady_,"cannot access ECVs before initialization");
232 10 : plumed_massert(j<getNumberOfArguments(),getName()+" has fewer CVs");
233 10 : return &ECVs_[j][0];
234 : }
235 :
236 10 : const double * ECVumbrellasLine::getPntrToDerECVs(unsigned j)
237 : {
238 10 : plumed_massert(isReady_,"cannot access ECVs before initialization");
239 10 : plumed_massert(j<getNumberOfArguments(),getName()+" has fewer CVs");
240 10 : return &derECVs_[j][0];
241 : }
242 :
243 8 : std::vector<std::string> ECVumbrellasLine::getLambdas() const
244 : {
245 8 : std::vector<std::string> lambdas(totNumECVs_);
246 8 : if(P0_contribution_==1)
247 : {
248 2 : std::ostringstream subs;
249 2 : subs<<"P0";
250 4 : for(unsigned j=1; j<getNumberOfArguments(); j++)
251 2 : subs<<"_P0";
252 2 : lambdas[0]=subs.str();
253 2 : }
254 94 : for(unsigned k=P0_contribution_; k<totNumECVs_; k++)
255 : {
256 86 : const unsigned kk=k-P0_contribution_;
257 86 : std::ostringstream subs;
258 86 : subs<<centers_[0][kk];
259 124 : for(unsigned j=1; j<getNumberOfArguments(); j++)
260 38 : subs<<"_"<<centers_[j][kk];
261 86 : lambdas[k]=subs.str();
262 86 : }
263 8 : return lambdas;
264 0 : }
265 :
266 8 : void ECVumbrellasLine::initECVs()
267 : {
268 8 : plumed_massert(!isReady_,"initialization should not be called twice");
269 8 : isReady_=true;
270 8 : log.printf(" *%4u windows for %s\n",totNumECVs_,getName().c_str());
271 8 : }
272 :
273 5 : void ECVumbrellasLine::initECVs_observ(const std::vector<double>& all_obs_cvs,const unsigned ncv,const unsigned index_j)
274 : {
275 : //this non-linear exansion never uses automatic initialization
276 5 : initECVs();
277 5 : calculateECVs(&all_obs_cvs[index_j]); //use only first obs point
278 11 : for(unsigned j=0; j<getNumberOfArguments(); j++)
279 76 : for(unsigned k=P0_contribution_; k<totNumECVs_; k++)
280 123 : ECVs_[j][k]=std::min(barrier_/kbt_,ECVs_[j][k]);
281 5 : }
282 :
283 3 : void ECVumbrellasLine::initECVs_restart(const std::vector<std::string>& lambdas)
284 : {
285 : std::size_t pos=0;
286 4 : for(unsigned j=0; j<getNumberOfArguments()-1; j++)
287 1 : pos=lambdas[0].find("_",pos+1); //checking only lambdas[0] is hopefully enough
288 3 : plumed_massert(pos<lambdas[0].length(),"this should not happen, fewer '_' than expected in "+getName());
289 3 : pos=lambdas[0].find("_",pos+1);
290 3 : plumed_massert(pos>lambdas[0].length(),"this should not happen, more '_' than expected in "+getName());
291 :
292 3 : std::vector<std::string> myLambdas=getLambdas();
293 3 : plumed_massert(myLambdas.size()==lambdas.size(),"RESTART - mismatch in number of "+getName());
294 3 : plumed_massert(std::equal(myLambdas.begin(),myLambdas.end(),lambdas.begin()),"RESTART - mismatch in lambda values of "+getName());
295 :
296 3 : initECVs();
297 3 : }
298 :
299 : }
300 : }
|
