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