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