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_CUSTOM
26 : /*
27 : Use some given CVs as a set of expansion collective variables (ECVs).
28 :
29 : This can be useful e.g. for quickly testing new ECVs, but from a performance point of view it is probably better to implement a new ECV class.
30 :
31 : By default the ARGs are expeted to be energies, \f$\Delta U_i\f$, and are then multiplied by the inverse temperature \f$\beta\f$
32 : \f[
33 : \Delta u_i=\beta \Delta U_i\, .
34 : \f]
35 : Use the DIMENSIONLESS flag to avoid this multiplication.
36 :
37 : The flag ADD_P0 adds also the unbiased distribution to the target.
38 : It is possible to specify a BARRIER as in \ref ECV_UMBRELLAS_LINE, to avoid a too high initial bias.
39 :
40 : \par Examples
41 :
42 : \plumedfile
43 : ene: ENERGY
44 : t1: CUSTOM PERIODIC=NO ARG=ene FUNC=(300/500-1)*x
45 : t2: CUSTOM PERIODIC=NO ARG=ene FUNC=(300/1000-1)*x
46 : ecv: ECV_CUSTOM ARG=t1,t2 TEMP=300 ADD_P0
47 : opes: OPES_EXPANDED ARG=ecv.* PACE=500
48 : \endplumedfile
49 :
50 : It is equivalent to the following:
51 :
52 : \plumedfile
53 : ene: ENERGY
54 : ecv: ECV_MULTITHERMAL ARG=ene TEMP=300 TEMP_SET_ALL=300,500,1000
55 : opes: OPES_EXPANDED ARG=ecv.* PACE=500
56 : \endplumedfile
57 :
58 : */
59 : //+ENDPLUMEDOC
60 :
61 : class ECVcustom :
62 : public ExpansionCVs {
63 : private:
64 : unsigned P0_contribution_;
65 : double barrier_;
66 : double beta0_;
67 :
68 : std::vector< std::vector<double> > ECVs_;
69 : std::vector< std::vector<double> > derECVs_;
70 : void initECVs();
71 :
72 : public:
73 : explicit ECVcustom(const ActionOptions&);
74 : static void registerKeywords(Keywords& keys);
75 : void calculateECVs(const double *) override;
76 : const double * getPntrToECVs(unsigned) override;
77 : const double * getPntrToDerECVs(unsigned) override;
78 : std::vector< std::vector<unsigned> > getIndex_k() const override;
79 : std::vector<std::string> getLambdas() const override;
80 : void initECVs_observ(const std::vector<double>&,const unsigned,const unsigned) override;
81 : void initECVs_restart(const std::vector<std::string>&) override;
82 : };
83 :
84 : PLUMED_REGISTER_ACTION(ECVcustom,"ECV_CUSTOM")
85 :
86 4 : void ECVcustom::registerKeywords(Keywords& keys) {
87 4 : ExpansionCVs::registerKeywords(keys);
88 4 : keys.remove("ARG");
89 8 : keys.add("compulsory","ARG","the labels of the single ECVs. Delta U_i, in energy units");
90 8 : keys.addFlag("ADD_P0",false,"add the unbiased Boltzmann distribution to the target distribution, to make sure to sample it");
91 8 : keys.addFlag("DIMENSIONLESS",false,"consider ARG as dimensionless rather than an energy, thus do not multiply it by beta");
92 8 : keys.add("optional","BARRIER","a guess of the free energy barrier to be overcome (better to stay higher than lower)");
93 4 : }
94 :
95 2 : ECVcustom::ECVcustom(const ActionOptions&ao)
96 : : Action(ao)
97 : , ExpansionCVs(ao)
98 2 : , beta0_(1./kbt_) {
99 : //set beta0_
100 : bool dimensionless;
101 2 : parseFlag("DIMENSIONLESS",dimensionless);
102 2 : if(dimensionless) {
103 0 : beta0_=1;
104 : }
105 :
106 : //set P0_contribution_
107 2 : bool add_P0=false;
108 2 : parseFlag("ADD_P0",add_P0);
109 2 : if(add_P0) {
110 2 : P0_contribution_=1;
111 : } else {
112 0 : P0_contribution_=0;
113 : }
114 :
115 : //set barrier_
116 2 : barrier_=std::numeric_limits<double>::infinity();
117 2 : parse("BARRIER",barrier_);
118 :
119 2 : checkRead();
120 :
121 : //set ECVs stuff
122 2 : totNumECVs_=getNumberOfArguments()+P0_contribution_;
123 2 : ECVs_.resize(getNumberOfArguments(),std::vector<double>(totNumECVs_));
124 2 : derECVs_.resize(getNumberOfArguments(),std::vector<double>(totNumECVs_));
125 6 : for(unsigned j=0; j<getNumberOfArguments(); j++) {
126 4 : derECVs_[j][j+P0_contribution_]=beta0_; //always constant
127 : }
128 :
129 : //print some info
130 2 : if(dimensionless) {
131 0 : log.printf(" -- DIMENSIONLESS: the ARG is not multiplied by beta\n");
132 : }
133 2 : if(barrier_!=std::numeric_limits<double>::infinity()) {
134 0 : log.printf(" guess for free energy BARRIER = %g\n",barrier_);
135 0 : if(dimensionless) {
136 0 : log.printf(" also the BARRIER is considered to be DIMENSIONLESS\n");
137 : }
138 : }
139 2 : if(P0_contribution_==1) {
140 2 : log.printf(" -- ADD_P0: the target includes also the unbiased probability itself\n");
141 : }
142 2 : }
143 :
144 91 : void ECVcustom::calculateECVs(const double * cv) {
145 273 : for(unsigned j=0; j<getNumberOfArguments(); j++) {
146 182 : ECVs_[j][j+P0_contribution_]=beta0_*cv[j];
147 : }
148 : //derivative is constant
149 91 : }
150 :
151 4 : const double * ECVcustom::getPntrToECVs(unsigned j) {
152 4 : plumed_massert(isReady_,"cannot access ECVs before initialization");
153 4 : plumed_massert(j<getNumberOfArguments(),getName()+" has fewer CVs");
154 4 : return &ECVs_[j][0];
155 : }
156 :
157 4 : const double * ECVcustom::getPntrToDerECVs(unsigned j) {
158 4 : plumed_massert(isReady_,"cannot access ECVs before initialization");
159 4 : plumed_massert(j<getNumberOfArguments(),getName()+" has fewer CVs");
160 4 : return &derECVs_[j][0];
161 : }
162 :
163 2 : std::vector< std::vector<unsigned> > ECVcustom::getIndex_k() const {
164 2 : plumed_massert(isReady_ && totNumECVs_>0,"cannot access getIndex_k() of ECV before initialization");
165 2 : std::vector< std::vector<unsigned> > index_k(totNumECVs_,std::vector<unsigned>(getNumberOfArguments()));
166 8 : for(unsigned k=0; k<totNumECVs_; k++)
167 18 : for(unsigned j=0; j<getNumberOfArguments(); j++)
168 12 : if(k==j+P0_contribution_) {
169 4 : index_k[k][j]=k;
170 : }
171 2 : return index_k;
172 0 : }
173 :
174 2 : std::vector<std::string> ECVcustom::getLambdas() const {
175 2 : std::vector<std::string> lambdas(totNumECVs_);
176 2 : if(P0_contribution_==1) {
177 2 : std::ostringstream subs;
178 2 : subs<<"P0";
179 4 : for(unsigned j=1; j<getNumberOfArguments(); j++) {
180 2 : subs<<"_P0";
181 : }
182 2 : lambdas[0]=subs.str();
183 2 : }
184 6 : for(unsigned k=P0_contribution_; k<totNumECVs_; k++) {
185 4 : const unsigned kk=k-P0_contribution_;
186 4 : std::ostringstream subs;
187 : //the getLambdas method is const, so it complains if one tries to access a non-const pointer, hence the const_cast
188 4 : if(kk==0) {
189 : subs<<const_cast<ECVcustom *>(this)->getPntrToArgument(kk)->getName();
190 : } else {
191 2 : subs<<"NaN";
192 : }
193 8 : for(unsigned j=1; j<getNumberOfArguments(); j++) {
194 4 : if(kk==j) {
195 2 : subs<<"_"<<const_cast<ECVcustom *>(this)->getPntrToArgument(kk)->getName();
196 : } else {
197 2 : subs<<"_NaN";
198 : }
199 : }
200 4 : lambdas[k]=subs.str();
201 4 : }
202 2 : return lambdas;
203 0 : }
204 :
205 2 : void ECVcustom::initECVs() {
206 2 : plumed_massert(!isReady_,"initialization should not be called twice");
207 2 : isReady_=true;
208 2 : log.printf(" *%4u ECVs for %s\n",totNumECVs_,getName().c_str());
209 2 : }
210 :
211 1 : void ECVcustom::initECVs_observ(const std::vector<double>& all_obs_cvs,const unsigned ncv,const unsigned index_j) {
212 1 : initECVs();
213 1 : calculateECVs(&all_obs_cvs[index_j]);
214 3 : for(unsigned j=0; j<getNumberOfArguments(); j++) {
215 4 : ECVs_[j][j+P0_contribution_]=std::min(barrier_*beta0_,ECVs_[j][j+P0_contribution_]);
216 : }
217 1 : }
218 :
219 1 : void ECVcustom::initECVs_restart(const std::vector<std::string>& lambdas) {
220 : std::size_t pos=0;
221 2 : for(unsigned j=0; j<getNumberOfArguments()-1; j++) {
222 1 : pos=lambdas[0].find("_",pos+1); //checking only lambdas[0] is hopefully enough
223 : }
224 1 : plumed_massert(pos<lambdas[0].length(),"this should not happen, fewer '_' than expected in "+getName());
225 1 : pos=lambdas[0].find("_",pos+1);
226 1 : plumed_massert(pos>lambdas[0].length(),"this should not happen, more '_' than expected in "+getName());
227 :
228 1 : std::vector<std::string> myLambdas=getLambdas();
229 1 : plumed_massert(myLambdas.size()==lambdas.size(),"RESTART - mismatch in number of "+getName());
230 1 : plumed_massert(std::equal(myLambdas.begin(),myLambdas.end(),lambdas.begin()),"RESTART - mismatch in lambda values of "+getName());
231 :
232 1 : initECVs();
233 1 : }
234 :
235 : }
236 : }
|
