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 10 : keys.use("ARG");
103 20 : keys.add("compulsory","CV_MIN","the minimum of the CV range to be explored");
104 20 : keys.add("compulsory","CV_MAX","the maximum of the CV range to be explored");
105 20 : keys.add("compulsory","SIGMA","sigma of the umbrella Gaussians");
106 20 : keys.add("compulsory","SPACING","1","the distance between umbrellas, in units of SIGMA");
107 20 : keys.add("optional","BARRIER","a guess of the free energy barrier to be overcome (better to stay higher than lower)");
108 20 : keys.addFlag("ADD_P0",false,"add the unbiased Boltzmann distribution to the target distribution, to make sure to sample it");
109 20 : keys.addFlag("LOWER_HALF_ONLY",false,"use only the lower half of each umbrella potentials");
110 10 : }
111 :
112 8 : ECVumbrellasLine::ECVumbrellasLine(const ActionOptions&ao):
113 : Action(ao),
114 8 : ExpansionCVs(ao)
115 : {
116 : //set P0_contribution_
117 8 : bool add_P0=false;
118 8 : parseFlag("ADD_P0",add_P0);
119 8 : if(add_P0)
120 2 : P0_contribution_=1;
121 : else
122 6 : P0_contribution_=0;
123 :
124 : //set barrier_
125 8 : barrier_=std::numeric_limits<double>::infinity();
126 8 : parse("BARRIER",barrier_);
127 8 : parseFlag("LOWER_HALF_ONLY",lower_only_);
128 :
129 : //set umbrellas
130 16 : parse("SIGMA",sigma_);
131 : std::vector<double> cv_min;
132 : std::vector<double> cv_max;
133 8 : parseVector("CV_MIN",cv_min);
134 16 : parseVector("CV_MAX",cv_max);
135 8 : plumed_massert(cv_min.size()==getNumberOfArguments(),"wrong number of CV_MINs");
136 8 : plumed_massert(cv_max.size()==getNumberOfArguments(),"wrong number of CV_MAXs");
137 : double spacing;
138 8 : parse("SPACING",spacing);
139 : double length=0;
140 18 : for(unsigned j=0; j<getNumberOfArguments(); j++)
141 10 : length+=std::pow(cv_max[j]-cv_min[j],2);
142 8 : length=std::sqrt(length);
143 8 : unsigned sizeUmbrellas=1+std::round(length/(sigma_*spacing));
144 8 : centers_.resize(getNumberOfArguments()); //centers_[cv][umbrellas]
145 : unsigned full_period=0;
146 18 : for(unsigned j=0; j<getNumberOfArguments(); j++)
147 : {
148 10 : centers_[j].resize(sizeUmbrellas);
149 10 : if(sizeUmbrellas>1)
150 140 : for(unsigned k=0; k<sizeUmbrellas; k++)
151 130 : centers_[j][k]=cv_min[j]+k*(cv_max[j]-cv_min[j])/(sizeUmbrellas-1);
152 : else
153 0 : centers_[j][0]=(cv_min[j]+cv_max[j])/2.;
154 10 : if(getPntrToArgument(j)->isPeriodic())
155 : {
156 : double min,max;
157 : std::string min_str,max_str;
158 10 : getPntrToArgument(j)->getDomain(min,max);
159 10 : getPntrToArgument(j)->getDomain(min_str,max_str);
160 10 : plumed_massert(cv_min[j]>=min,"ARG "+std::to_string(j)+": CV_MIN cannot be smaller than the periodic bound "+min_str);
161 10 : plumed_massert(cv_max[j]<=max,"ARG "+std::to_string(j)+": CV_MAX cannot be greater than the periodic bound "+max_str);
162 10 : if(cv_min[j]==min && cv_max[j]==max)
163 6 : full_period++;
164 : }
165 : }
166 8 : if(full_period==getNumberOfArguments() && sizeUmbrellas>1) //first and last are the same point
167 : {
168 6 : sizeUmbrellas--;
169 12 : for(unsigned j=0; j<getNumberOfArguments(); j++)
170 6 : centers_[j].pop_back();
171 : }
172 :
173 8 : checkRead();
174 :
175 : //set ECVs stuff
176 8 : totNumECVs_=sizeUmbrellas+P0_contribution_;
177 8 : ECVs_.resize(getNumberOfArguments(),std::vector<double>(totNumECVs_));
178 8 : derECVs_.resize(getNumberOfArguments(),std::vector<double>(totNumECVs_));
179 :
180 : //printing some info
181 8 : log.printf(" total number of umbrellas = %u\n",sizeUmbrellas);
182 8 : log.printf(" with SIGMA = %g\n",sigma_);
183 8 : log.printf(" and SPACING = %g\n",spacing);
184 8 : if(barrier_!=std::numeric_limits<double>::infinity())
185 2 : log.printf(" guess for free energy BARRIER = %g\n",barrier_);
186 8 : if(P0_contribution_==1)
187 2 : log.printf(" -- ADD_P0: the target includes also the unbiased probability itself\n");
188 8 : if(lower_only_)
189 1 : log.printf(" -- LOWER_HALF_ONLY: the ECVs are set to zero for values of the CV above the respective center\n");
190 8 : }
191 :
192 433 : void ECVumbrellasLine::calculateECVs(const double * cv)
193 : {
194 433 : if(lower_only_)
195 : {
196 153 : for(unsigned j=0; j<getNumberOfArguments(); j++)
197 : {
198 2040 : for(unsigned k=P0_contribution_; k<totNumECVs_; k++) //if ADD_P0, the first ECVs=0
199 : {
200 1938 : const unsigned kk=k-P0_contribution_;
201 1938 : const double dist_jk=difference(j,centers_[j][kk],cv[j])/sigma_; //PBC might be present
202 1938 : if(dist_jk>=0)
203 : {
204 933 : ECVs_[j][k]=0;
205 933 : derECVs_[j][k]=0;
206 : }
207 : else
208 : {
209 1005 : ECVs_[j][k]=0.5*std::pow(dist_jk,2);
210 1005 : derECVs_[j][k]=dist_jk/sigma_;
211 : }
212 : }
213 : }
214 : }
215 : else
216 : {
217 804 : for(unsigned j=0; j<getNumberOfArguments(); j++)
218 : {
219 4678 : for(unsigned k=P0_contribution_; k<totNumECVs_; k++) //if ADD_P0, the first ECVs=0
220 : {
221 4256 : const unsigned kk=k-P0_contribution_;
222 4256 : const double dist_jk=difference(j,centers_[j][kk],cv[j])/sigma_; //PBC might be present
223 4256 : ECVs_[j][k]=0.5*std::pow(dist_jk,2);
224 4256 : derECVs_[j][k]=dist_jk/sigma_;
225 : }
226 : }
227 : }
228 433 : }
229 :
230 10 : const double * ECVumbrellasLine::getPntrToECVs(unsigned j)
231 : {
232 10 : plumed_massert(isReady_,"cannot access ECVs before initialization");
233 10 : plumed_massert(j<getNumberOfArguments(),getName()+" has fewer CVs");
234 10 : return &ECVs_[j][0];
235 : }
236 :
237 10 : const double * ECVumbrellasLine::getPntrToDerECVs(unsigned j)
238 : {
239 10 : plumed_massert(isReady_,"cannot access ECVs before initialization");
240 10 : plumed_massert(j<getNumberOfArguments(),getName()+" has fewer CVs");
241 10 : return &derECVs_[j][0];
242 : }
243 :
244 8 : std::vector<std::string> ECVumbrellasLine::getLambdas() const
245 : {
246 8 : std::vector<std::string> lambdas(totNumECVs_);
247 8 : if(P0_contribution_==1)
248 : {
249 2 : std::ostringstream subs;
250 2 : subs<<"P0";
251 4 : for(unsigned j=1; j<getNumberOfArguments(); j++)
252 2 : subs<<"_P0";
253 2 : lambdas[0]=subs.str();
254 2 : }
255 94 : for(unsigned k=P0_contribution_; k<totNumECVs_; k++)
256 : {
257 86 : const unsigned kk=k-P0_contribution_;
258 86 : std::ostringstream subs;
259 86 : subs<<centers_[0][kk];
260 124 : for(unsigned j=1; j<getNumberOfArguments(); j++)
261 38 : subs<<"_"<<centers_[j][kk];
262 86 : lambdas[k]=subs.str();
263 86 : }
264 8 : return lambdas;
265 0 : }
266 :
267 8 : void ECVumbrellasLine::initECVs()
268 : {
269 8 : plumed_massert(!isReady_,"initialization should not be called twice");
270 8 : isReady_=true;
271 8 : log.printf(" *%4u windows for %s\n",totNumECVs_,getName().c_str());
272 8 : }
273 :
274 5 : void ECVumbrellasLine::initECVs_observ(const std::vector<double>& all_obs_cvs,const unsigned ncv,const unsigned index_j)
275 : {
276 : //this non-linear exansion never uses automatic initialization
277 5 : initECVs();
278 5 : calculateECVs(&all_obs_cvs[index_j]); //use only first obs point
279 11 : for(unsigned j=0; j<getNumberOfArguments(); j++)
280 76 : for(unsigned k=P0_contribution_; k<totNumECVs_; k++)
281 123 : ECVs_[j][k]=std::min(barrier_/kbt_,ECVs_[j][k]);
282 5 : }
283 :
284 3 : void ECVumbrellasLine::initECVs_restart(const std::vector<std::string>& lambdas)
285 : {
286 : std::size_t pos=0;
287 4 : for(unsigned j=0; j<getNumberOfArguments()-1; j++)
288 1 : pos=lambdas[0].find("_",pos+1); //checking only lambdas[0] is hopefully enough
289 3 : plumed_massert(pos<lambdas[0].length(),"this should not happen, fewer '_' than expected in "+getName());
290 3 : pos=lambdas[0].find("_",pos+1);
291 3 : plumed_massert(pos>lambdas[0].length(),"this should not happen, more '_' than expected in "+getName());
292 :
293 3 : std::vector<std::string> myLambdas=getLambdas();
294 3 : plumed_massert(myLambdas.size()==lambdas.size(),"RESTART - mismatch in number of "+getName());
295 3 : plumed_massert(std::equal(myLambdas.begin(),myLambdas.end(),lambdas.begin()),"RESTART - mismatch in lambda values of "+getName());
296 :
297 3 : initECVs();
298 3 : }
299 :
300 : }
301 : }
|
