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