Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2018-2023 The plumed team
3 : (see the PEOPLE file at the root of the distribution for a list of names)
4 :
5 : See http://www.plumed.org for more information.
6 :
7 : This file is part of plumed, version 2.
8 :
9 : plumed is free software: you can redistribute it and/or modify
10 : it under the terms of the GNU Lesser General Public License as published by
11 : the Free Software Foundation, either version 3 of the License, or
12 : (at your option) any later version.
13 :
14 : plumed is distributed in the hope that it will be useful,
15 : but WITHOUT ANY WARRANTY; without even the implied warranty of
16 : MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 : GNU Lesser General Public License for more details.
18 :
19 : You should have received a copy of the GNU Lesser General Public License
20 : along with plumed. If not, see <http://www.gnu.org/licenses/>.
21 : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
22 : #include "bias/Bias.h"
23 : #include "core/ActionRegister.h"
24 : #include "core/PlumedMain.h"
25 : #include "tools/Communicator.h"
26 : #include <fstream>
27 :
28 : namespace PLMD {
29 : namespace isdb {
30 :
31 : //+PLUMEDOC ISDB_BIAS CALIBER
32 : /*
33 : Add a time-dependent, harmonic restraint on one or more variables.
34 :
35 : This allows implementing a maximum caliber restraint on one or more experimental time series by replica-averaged restrained simulations.
36 : See \cite Capelli:2018jt .
37 :
38 : The time resolved experiments are read from a text file and intermediate values are obtained by splines.
39 :
40 : \par Examples
41 :
42 : In the following example a restraint is applied on the time evolution of a saxs spectrum
43 :
44 : \plumedfile
45 : MOLINFO STRUCTURE=first.pdb
46 :
47 : # Define saxs variable
48 : SAXS ...
49 : LABEL=saxs
50 : ATOMISTIC
51 : ATOMS=1-436
52 : QVALUE1=0.02 # Q-value at which calculate the scattering
53 : QVALUE2=0.0808
54 : QVALUE3=0.1264
55 : QVALUE4=0.1568
56 : QVALUE5=0.172
57 : QVALUE6=0.1872
58 : QVALUE7=0.2176
59 : QVALUE8=0.2328
60 : QVALUE9=0.248
61 : QVALUE10=0.2632
62 : QVALUE11=0.2936
63 : QVALUE12=0.3088
64 : QVALUE13=0.324
65 : QVALUE14=0.3544
66 : QVALUE15=0.4
67 : ... SAXS
68 :
69 :
70 : #define the caliber restraint
71 : CALIBER ...
72 : ARG=(saxs\.q_.*)
73 : FILE=expsaxs.dat
74 : KAPPA=10
75 : LABEL=cal0
76 : STRIDE=10
77 : REGRES_ZERO=200
78 : AVERAGING=200
79 : ... CALIBER
80 : \endplumedfile
81 :
82 : In particular the file expsaxs.dat contains the time traces for the 15 intensities at the selected scattering lengths, organized as time, q_1, etc.
83 : The strength of the bias is automatically evaluated from the standard error of the mean over AVERAGING steps and multiplied by KAPPA. This is useful when working with multiple experimental data
84 : Because \ref SAXS is usually defined in a manner that is irrespective of a scaling factor the scaling is evaluated from a linear fit every REGRES_ZERO step. Alternatively it can be given as a fixed constant as SCALE.
85 : The bias is here applied every tenth step.
86 :
87 : */
88 : //+ENDPLUMEDOC
89 :
90 :
91 : class Caliber : public bias::Bias {
92 : public:
93 : explicit Caliber(const ActionOptions&);
94 : void calculate();
95 : static void registerKeywords( Keywords& keys );
96 : private:
97 : std::vector<double> time;
98 : std::vector< std::vector<double> > var;
99 : std::vector< std::vector<double> > dvar;
100 : double mult;
101 : double scale_;
102 : bool master;
103 : unsigned replica_;
104 : unsigned nrep_;
105 : // scale and offset regression
106 : bool doregres_zero_;
107 : int nregres_zero_;
108 : // force constant
109 : unsigned optsigmamean_stride_;
110 : std::vector<double> sigma_mean2_;
111 : std::vector< std::vector<double> > sigma_mean2_last_;
112 : std::vector<Value*> x0comp;
113 : std::vector<Value*> kcomp;
114 : std::vector<Value*> mcomp;
115 : Value* valueScale;
116 :
117 : void get_sigma_mean(const double fact, const std::vector<double> &mean);
118 : void replica_averaging(const double fact, std::vector<double> &mean);
119 : double getSpline(const unsigned iarg);
120 : void do_regression_zero(const std::vector<double> &mean);
121 : };
122 :
123 : PLUMED_REGISTER_ACTION(Caliber,"CALIBER")
124 :
125 6 : void Caliber::registerKeywords( Keywords& keys ) {
126 6 : Bias::registerKeywords(keys);
127 12 : keys.addFlag("NOENSEMBLE",false,"don't perform any replica-averaging");
128 12 : keys.add("compulsory","FILE","the name of the file containing the time-resolved values");
129 12 : keys.add("compulsory","KAPPA","a force constant, this can be use to scale a constant estimated on-the-fly using AVERAGING");
130 12 : keys.add("optional","AVERAGING", "Stride for calculation of the optimum kappa, if 0 only KAPPA is used.");
131 12 : keys.add("compulsory","TSCALE","1.0","Apply a time scaling on the experimental time scale");
132 12 : keys.add("compulsory","SCALE","1.0","Apply a constant scaling on the data provided as arguments");
133 12 : keys.add("optional","REGRES_ZERO","stride for regression with zero offset");
134 12 : keys.addOutputComponent("x0","default","scalar","the instantaneous value of the center of the potential");
135 12 : keys.addOutputComponent("mean","default","scalar","the current average value of the calculated observable");
136 12 : keys.addOutputComponent("kappa","default","scalar","the current force constant");
137 12 : keys.addOutputComponent("scale","REGRES_ZERO","scalar","the current scaling constant");
138 6 : }
139 :
140 4 : Caliber::Caliber(const ActionOptions&ao):
141 : PLUMED_BIAS_INIT(ao),
142 4 : mult(0),
143 4 : scale_(1),
144 4 : doregres_zero_(false),
145 4 : nregres_zero_(0),
146 4 : optsigmamean_stride_(0)
147 : {
148 8 : parse("KAPPA",mult);
149 : std::string filename;
150 8 : parse("FILE",filename);
151 4 : if( filename.length()==0 ) error("No external variable file was specified");
152 4 : unsigned averaging=0;
153 4 : parse("AVERAGING", averaging);
154 4 : if(averaging>0) optsigmamean_stride_ = averaging;
155 4 : double tscale=1.0;
156 4 : parse("TSCALE", tscale);
157 4 : if(tscale<=0.) error("The time scale factor must be greater than 0.");
158 4 : parse("SCALE", scale_);
159 4 : if(scale_==0.) error("The time scale factor cannot be 0.");
160 : // regression with zero intercept
161 4 : parse("REGRES_ZERO", nregres_zero_);
162 4 : if(nregres_zero_>0) {
163 : // set flag
164 0 : doregres_zero_=true;
165 0 : log.printf(" doing regression with zero intercept with stride: %d\n", nregres_zero_);
166 : }
167 :
168 :
169 4 : bool noensemble = false;
170 4 : parseFlag("NOENSEMBLE", noensemble);
171 :
172 4 : checkRead();
173 :
174 : // set up replica stuff
175 4 : master = (comm.Get_rank()==0);
176 4 : if(master) {
177 4 : nrep_ = multi_sim_comm.Get_size();
178 4 : replica_ = multi_sim_comm.Get_rank();
179 4 : if(noensemble) nrep_ = 1;
180 : } else {
181 0 : nrep_ = 0;
182 0 : replica_ = 0;
183 : }
184 4 : comm.Sum(&nrep_,1);
185 4 : comm.Sum(&replica_,1);
186 :
187 : const unsigned narg = getNumberOfArguments();
188 4 : sigma_mean2_.resize(narg,1);
189 4 : sigma_mean2_last_.resize(narg);
190 8 : for(unsigned j=0; j<narg; j++) sigma_mean2_last_[j].push_back(0.000001);
191 :
192 4 : log.printf(" Time resolved data from file %s\n",filename.c_str());
193 4 : std::ifstream varfile(filename.c_str());
194 4 : if(varfile.fail()) error("Cannot open "+filename);
195 4 : var.resize(narg);
196 4 : dvar.resize(narg);
197 2012 : while (!varfile.eof()) {
198 : double tempT, tempVar;
199 : varfile >> tempT;
200 2008 : time.push_back(tempT/tscale);
201 4016 : for(unsigned i=0; i<narg; i++) {
202 : varfile >> tempVar;
203 2008 : var[i].push_back(tempVar);
204 : }
205 : }
206 4 : varfile.close();
207 :
208 4 : const double deltat = time[1] - time[0];
209 8 : for(unsigned i=0; i<narg; i++) {
210 2012 : for(unsigned j=0; j<var[i].size(); j++) {
211 2008 : if(j==0) dvar[i].push_back((var[i][j+1] - var[i][j])/(deltat));
212 2004 : else if(j==var[i].size()-1) dvar[i].push_back((var[i][j] - var[i][j-1])/(deltat));
213 2000 : else dvar[i].push_back((var[i][j+1] - var[i][j-1])/(2.*deltat));
214 : }
215 : }
216 :
217 8 : for(unsigned i=0; i<narg; i++) {
218 4 : std::string num; Tools::convert(i,num);
219 12 : addComponent("x0-"+num); componentIsNotPeriodic("x0-"+num); x0comp.push_back(getPntrToComponent("x0-"+num));
220 12 : addComponent("kappa-"+num); componentIsNotPeriodic("kappa-"+num); kcomp.push_back(getPntrToComponent("kappa-"+num));
221 12 : addComponent("mean-"+num); componentIsNotPeriodic("mean-"+num); mcomp.push_back(getPntrToComponent("mean-"+num));
222 : }
223 :
224 4 : if(doregres_zero_) {
225 0 : addComponent("scale");
226 0 : componentIsNotPeriodic("scale");
227 0 : valueScale=getPntrToComponent("scale");
228 : }
229 :
230 8 : log<<" Bibliography "<<plumed.cite("Capelli, Tiana, Camilloni, J Chem Phys, 148, 184114");
231 8 : }
232 :
233 0 : void Caliber::get_sigma_mean(const double fact, const std::vector<double> &mean)
234 : {
235 0 : const unsigned narg = getNumberOfArguments();
236 0 : const double dnrep = static_cast<double>(nrep_);
237 :
238 0 : if(sigma_mean2_last_[0].size()==optsigmamean_stride_) for(unsigned i=0; i<narg; ++i) sigma_mean2_last_[i].erase(sigma_mean2_last_[i].begin());
239 0 : std::vector<double> sigma_mean2_now(narg,0);
240 0 : if(master) {
241 0 : for(unsigned i=0; i<narg; ++i) {
242 0 : double tmp = getArgument(i)-mean[i];
243 0 : sigma_mean2_now[i] = fact*tmp*tmp;
244 : }
245 0 : if(nrep_>1) multi_sim_comm.Sum(&sigma_mean2_now[0], narg);
246 : }
247 0 : comm.Sum(&sigma_mean2_now[0], narg);
248 :
249 0 : for(unsigned i=0; i<narg; ++i) {
250 0 : sigma_mean2_last_[i].push_back(sigma_mean2_now[i]/dnrep);
251 0 : sigma_mean2_[i] = *max_element(sigma_mean2_last_[i].begin(), sigma_mean2_last_[i].end());
252 : }
253 0 : }
254 :
255 2004 : void Caliber::replica_averaging(const double fact, std::vector<double> &mean)
256 : {
257 2004 : const unsigned narg = getNumberOfArguments();
258 2004 : if(master) {
259 4008 : for(unsigned i=0; i<narg; ++i) mean[i] = fact*getArgument(i);
260 2004 : if(nrep_>1) multi_sim_comm.Sum(&mean[0], narg);
261 : }
262 2004 : comm.Sum(&mean[0], narg);
263 2004 : }
264 :
265 2004 : double Caliber::getSpline(const unsigned iarg)
266 : {
267 2004 : const double deltat = time[1] - time[0];
268 2004 : const int tindex = static_cast<int>(getTime()/deltat);
269 :
270 : unsigned start, end;
271 2004 : start=tindex;
272 2004 : if(tindex+1<var[iarg].size()) end=tindex+2;
273 0 : else end=var[iarg].size();
274 :
275 : double value=0;
276 6012 : for(unsigned ipoint=start; ipoint<end; ++ipoint) {
277 4008 : double grid=var[iarg][ipoint];
278 4008 : double dder=dvar[iarg][ipoint];
279 : double yy=0.;
280 4008 : if(std::abs(grid)>0.0000001) yy=-dder/grid;
281 :
282 : int x0=1;
283 4008 : if(ipoint==tindex) x0=0;
284 :
285 4008 : double X=std::abs((getTime()-time[tindex])/deltat-(double)x0);
286 4008 : double X2=X*X;
287 4008 : double X3=X2*X;
288 4008 : double C=(1.0-3.0*X2+2.0*X3) - (x0?-1.0:1.0)*yy*(X-2.0*X2+X3)*deltat;
289 :
290 4008 : value+=grid*C;
291 : }
292 2004 : return value;
293 : }
294 :
295 0 : void Caliber::do_regression_zero(const std::vector<double> &mean)
296 : {
297 : // parameters[i] = scale_ * mean[i]: find scale_ with linear regression
298 : double num = 0.0;
299 : double den = 0.0;
300 0 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
301 0 : num += mean[i] * getSpline(i);
302 0 : den += mean[i] * mean[i];
303 : }
304 0 : if(den>0) {
305 0 : scale_ = num / den;
306 : } else {
307 0 : scale_ = 1.0;
308 : }
309 0 : }
310 :
311 2004 : void Caliber::calculate()
312 : {
313 2004 : const unsigned narg = getNumberOfArguments();
314 2004 : const double dnrep = static_cast<double>(nrep_);
315 2004 : const double fact = 1.0/dnrep;
316 :
317 2004 : std::vector<double> mean(narg,0);
318 2004 : std::vector<double> dmean_x(narg,fact);
319 2004 : replica_averaging(fact, mean);
320 2004 : if(optsigmamean_stride_>0) get_sigma_mean(fact, mean);
321 :
322 : // in case of regression with zero intercept, calculate scale
323 2004 : if(doregres_zero_ && getStep()%nregres_zero_==0) do_regression_zero(mean);
324 :
325 : double ene=0;
326 4008 : for(unsigned i=0; i<narg; ++i) {
327 2004 : const double x0 = getSpline(i);
328 2004 : const double kappa = mult*dnrep/sigma_mean2_[i];
329 2004 : const double cv=difference(i,x0,scale_*mean[i]);
330 2004 : const double f=-kappa*cv*dmean_x[i]/scale_;
331 2004 : setOutputForce(i,f);
332 2004 : ene+=0.5*kappa*cv*cv;
333 2004 : x0comp[i]->set(x0);
334 2004 : kcomp[i]->set(kappa);
335 2004 : mcomp[i]->set(mean[i]);
336 : }
337 :
338 2004 : if(doregres_zero_) valueScale->set(scale_);
339 :
340 2004 : setBias(ene);
341 2004 : }
342 :
343 : }
344 : }
345 :
346 :
|
