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