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 :
21 : #include "tools/OpenMP.h"
22 : #include "core/PlumedMain.h"
23 : #include "core/Atoms.h"
24 :
25 : namespace PLMD {
26 : namespace opes {
27 :
28 43 : void ExpansionCVs::registerKeywords(Keywords& keys)
29 : {
30 43 : Action::registerKeywords(keys);
31 43 : ActionWithValue::registerKeywords(keys);
32 43 : ActionWithArguments::registerKeywords(keys);
33 43 : ActionWithValue::useCustomisableComponents(keys);
34 86 : keys.add("compulsory","TEMP","-1","temperature. If not specified tries to get it from MD engine");
35 43 : }
36 :
37 37 : ExpansionCVs::ExpansionCVs(const ActionOptions&ao)
38 : : Action(ao)
39 : , ActionWithValue(ao)
40 : , ActionWithArguments(ao)
41 37 : , isReady_(false)
42 37 : , totNumECVs_(0)
43 : {
44 : //set kbt_
45 37 : const double kB=plumed.getAtoms().getKBoltzmann();
46 37 : kbt_=plumed.getAtoms().getKbT();
47 37 : double temp=-1;
48 37 : parse("TEMP",temp);
49 37 : if(temp>0)
50 : {
51 37 : if(kbt_>0 && std::abs(kbt_-kB*temp)>1e-4)
52 0 : log.printf(" +++ WARNING +++ using TEMP=%g while MD engine uses %g\n",temp,kbt_/kB);
53 37 : kbt_=kB*temp;
54 : }
55 37 : plumed_massert(kbt_>0,"your MD engine does not pass the temperature to plumed, you must specify it using TEMP");
56 37 : log.printf(" temperature = %g, beta = %g\n",kbt_/kB,1./kbt_);
57 :
58 : //set components
59 37 : plumed_massert( getNumberOfArguments()!=0, "you must specify the underlying CV");
60 90 : for(unsigned j=0; j<getNumberOfArguments(); j++)
61 : {
62 53 : std::string name_j=getPntrToArgument(j)->getName();
63 53 : ActionWithValue::addComponentWithDerivatives(name_j);
64 53 : getPntrToComponent(j)->resizeDerivatives(1);
65 53 : if(getPntrToArgument(j)->isPeriodic()) //it should not be necessary, but why not
66 : {
67 : std::string min,max;
68 17 : getPntrToArgument(j)->getDomain(min,max);
69 17 : getPntrToComponent(j)->setDomain(min,max);
70 : }
71 : else
72 36 : getPntrToComponent(j)->setNotPeriodic();
73 : }
74 37 : plumed_massert((int)getNumberOfArguments()==getNumberOfComponents(),"Expansion CVs have same number of arguments and components");
75 37 : }
76 :
77 1847 : void ExpansionCVs::calculate()
78 : {
79 1847 : std::vector<double> args(getNumberOfArguments());
80 4470 : for(unsigned j=0; j<getNumberOfArguments(); j++)
81 : {
82 2623 : args[j]=getArgument(j);
83 2623 : getPntrToComponent(j)->set(args[j]); //components are equal to arguments
84 2623 : getPntrToComponent(j)->addDerivative(0,1.); //the derivative of the identity is 1
85 : }
86 1847 : if(isReady_)
87 1417 : calculateECVs(&args[0]);
88 1847 : }
89 :
90 1847 : void ExpansionCVs::apply()
91 : {
92 4470 : for(unsigned j=0; j<getNumberOfArguments(); j++)
93 : {
94 2623 : std::vector<double> force_j(1);
95 2623 : if(getPntrToComponent(j)->applyForce(force_j)) //a bias is applied?
96 2623 : getPntrToArgument(j)->addForce(force_j[0]); //just tell it to the CV!
97 : }
98 1847 : }
99 :
100 26 : std::vector< std::vector<unsigned> > ExpansionCVs::getIndex_k() const
101 : {
102 26 : plumed_massert(isReady_ && totNumECVs_>0,"cannot access getIndex_k() of ECV before initialization");
103 26 : std::vector< std::vector<unsigned> > index_k(totNumECVs_,std::vector<unsigned>(getNumberOfArguments()));
104 869 : for(unsigned k=0; k<totNumECVs_; k++)
105 2391 : for(unsigned j=0; j<getNumberOfArguments(); j++)
106 1548 : index_k[k][j]=k; //each CV gives rise to the same number of ECVs
107 26 : return index_k;
108 0 : }
109 :
110 : //following methods are meant to be used only in case of linear expansions
111 24 : std::vector<double> ExpansionCVs::getSteps(double lambda_min,double lambda_max,const unsigned lambda_steps,const std::string& msg,const bool geom_spacing, const double shift) const
112 : {
113 24 : plumed_massert(!(lambda_min==lambda_max && lambda_steps>1),"cannot have multiple "+msg+"_STEPS if "+msg+"_MIN=="+msg+"_MAX");
114 24 : std::vector<double> lambda(lambda_steps);
115 24 : if(lambda_steps==1)
116 : {
117 0 : lambda[0]=(lambda_min+lambda_max)/2.;
118 0 : log.printf(" +++ WARNING +++ using one single %s as target = %g\n",msg.c_str(),lambda[0]);
119 : }
120 : else
121 : {
122 24 : if(geom_spacing) //geometric spacing
123 : { //this way lambda[k]/lambda[k+1] is constant
124 14 : lambda_min+=shift;
125 14 : lambda_max+=shift;
126 14 : plumed_massert(lambda_min>0,"cannot use GEOM_SPACING when %s_MIN is not greater than zero");
127 14 : plumed_massert(lambda_max>0,"cannot use GEOM_SPACING when %s_MAX is not greater than zero");
128 14 : const double log_lambda_min=std::log(lambda_min);
129 14 : const double log_lambda_max=std::log(lambda_max);
130 196 : for(unsigned k=0; k<lambda.size(); k++)
131 182 : lambda[k]=std::exp(log_lambda_min+k*(log_lambda_max-log_lambda_min)/(lambda_steps-1))-shift;
132 : }
133 : else //linear spacing
134 108 : for(unsigned k=0; k<lambda.size(); k++)
135 98 : lambda[k]=lambda_min+k*(lambda_max-lambda_min)/(lambda_steps-1);
136 : }
137 24 : return lambda;
138 : }
139 :
140 6 : unsigned ExpansionCVs::estimateNumSteps(const double left_side,const double right_side,const std::vector<double>& obs,const std::string& msg) const
141 : { //for linear expansions only, it uses effective sample size (Neff) to estimate the grid spacing
142 6 : if(left_side==0 && right_side==0)
143 : {
144 0 : log.printf(" +++ WARNING +++ %s_MIN and %s_MAX are equal to %s, using single step\n",msg.c_str(),msg.c_str(),msg.c_str());
145 0 : return 1;
146 : }
147 9 : auto get_neff_HWHM=[](const double side,const std::vector<double>& obs) //HWHM = half width at half maximum. neff is in general not symmetric
148 : {
149 : //func: Neff/N-0.5 is a function between -0.5 and 0.5
150 109 : auto func=[](const double delta,const std::vector<double>& obs)
151 : {
152 : double sum_w=0;
153 : double sum_w2=0;
154 : //we could avoid recomputing safe_shift every time, but here speed is not a concern
155 109 : const double safe_shift=delta<0?*std::max_element(obs.begin(),obs.end()):*std::min_element(obs.begin(),obs.end());
156 899 : for(unsigned t=0; t<obs.size(); t++)
157 : {
158 790 : const double w=std::exp(-delta*(obs[t]-safe_shift)); //robust to overflow
159 790 : sum_w+=w;
160 790 : sum_w2+=w*w;
161 : }
162 109 : return sum_w*sum_w/sum_w2/obs.size()-0.5;
163 : };
164 : //here we find the root of func using the regula falsi (false position) method
165 : //but any method would be OK, not much precision is needed. src/tools/RootFindingBase.h looked complicated
166 : const double tolerance=1e-4; //seems to be a good default
167 : double a=0; //default is right side case
168 : double func_a=0.5;
169 : double b=side;
170 9 : double func_b=func(side,obs);
171 9 : if(func_b>=0)
172 : return 0.0; //no zero is present!
173 9 : if(b<0) //left side case
174 : {
175 : std::swap(a,b);
176 : std::swap(func_a,func_b);
177 : }
178 : double c=a;
179 : double func_c=func_a;
180 109 : while(std::abs(func_c)>tolerance)
181 : {
182 100 : if(func_a*func_c>0)
183 : {
184 : a=c;
185 : func_a=func_c;
186 : }
187 : else
188 : {
189 : b=c;
190 : func_b=func_c;
191 : }
192 100 : c=(a*func_b-b*func_a)/(func_b-func_a);
193 100 : func_c=func(c,obs); //func is evaluated only here, it might be expensive
194 : }
195 9 : return std::abs(c);
196 : };
197 :
198 : //estimation
199 : double left_HWHM=0;
200 6 : if(left_side!=0)
201 4 : left_HWHM=get_neff_HWHM(left_side,obs);
202 : double right_HWHM=0;
203 6 : if(right_side!=0)
204 5 : right_HWHM=get_neff_HWHM(right_side,obs);
205 6 : if(left_HWHM==0)
206 : {
207 2 : right_HWHM*=2;
208 2 : if(left_side==0)
209 2 : log.printf(" --- %s_MIN is equal to %s\n",msg.c_str(),msg.c_str());
210 : else
211 0 : log.printf(" +++ WARNING +++ %s_MIN is very close to %s\n",msg.c_str(),msg.c_str());
212 : }
213 6 : if(right_HWHM==0)
214 : {
215 1 : left_HWHM*=2;
216 1 : if(right_side==0)
217 1 : log.printf(" --- %s_MAX is equal to %s\n",msg.c_str(),msg.c_str());
218 : else
219 0 : log.printf(" +++ WARNING +++ %s_MAX is very close to %s\n",msg.c_str(),msg.c_str());
220 : }
221 6 : const double grid_spacing=left_HWHM+right_HWHM;
222 6 : log.printf(" estimated %s spacing = %g\n",msg.c_str(),grid_spacing);
223 6 : unsigned steps=std::ceil(std::abs(right_side-left_side)/grid_spacing);
224 6 : if(steps<2 || grid_spacing==0)
225 : {
226 0 : log.printf(" +++ WARNING +++ %s range is very narrow, using %s_MIN and %s_MAX as only steps\n",msg.c_str(),msg.c_str(),msg.c_str());
227 : steps=2;
228 : }
229 : return steps;
230 : }
231 :
232 : }
233 : }
|
