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