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