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 10 : void Caliber::registerKeywords( Keywords& keys ) {
126 10 : Bias::registerKeywords(keys);
127 10 : keys.use("ARG");
128 20 : keys.addFlag("NOENSEMBLE",false,"don't perform any replica-averaging");
129 20 : keys.add("compulsory","FILE","the name of the file containing the time-resolved values");
130 20 : keys.add("compulsory","KAPPA","a force constant, this can be use to scale a constant estimated on-the-fly using AVERAGING");
131 20 : keys.add("optional","AVERAGING", "Stride for calculation of the optimum kappa, if 0 only KAPPA is used.");
132 20 : keys.add("compulsory","TSCALE","1.0","Apply a time scaling on the experimental time scale");
133 20 : keys.add("compulsory","SCALE","1.0","Apply a constant scaling on the data provided as arguments");
134 20 : keys.add("optional","REGRES_ZERO","stride for regression with zero offset");
135 20 : keys.addOutputComponent("x0","default","the instantaneous value of the center of the potential");
136 20 : keys.addOutputComponent("mean","default","the current average value of the calculated observable");
137 20 : keys.addOutputComponent("kappa","default","the current force constant");
138 20 : keys.addOutputComponent("scale","REGRES_ZERO","the current scaling constant");
139 10 : }
140 :
141 8 : Caliber::Caliber(const ActionOptions&ao):
142 : PLUMED_BIAS_INIT(ao),
143 8 : mult(0),
144 8 : scale_(1),
145 8 : doregres_zero_(false),
146 8 : nregres_zero_(0),
147 8 : optsigmamean_stride_(0) {
148 16 : parse("KAPPA",mult);
149 : std::string filename;
150 16 : parse("FILE",filename);
151 8 : if( filename.length()==0 ) {
152 0 : error("No external variable file was specified");
153 : }
154 8 : unsigned averaging=0;
155 8 : parse("AVERAGING", averaging);
156 8 : if(averaging>0) {
157 0 : optsigmamean_stride_ = averaging;
158 : }
159 8 : double tscale=1.0;
160 8 : parse("TSCALE", tscale);
161 8 : if(tscale<=0.) {
162 0 : error("The time scale factor must be greater than 0.");
163 : }
164 8 : parse("SCALE", scale_);
165 8 : if(scale_==0.) {
166 0 : error("The time scale factor cannot be 0.");
167 : }
168 : // regression with zero intercept
169 8 : parse("REGRES_ZERO", nregres_zero_);
170 8 : if(nregres_zero_>0) {
171 : // set flag
172 4 : doregres_zero_=true;
173 4 : log.printf(" doing regression with zero intercept with stride: %d\n", nregres_zero_);
174 : }
175 :
176 :
177 8 : bool noensemble = false;
178 8 : parseFlag("NOENSEMBLE", noensemble);
179 :
180 8 : checkRead();
181 :
182 : // set up replica stuff
183 8 : master = (comm.Get_rank()==0);
184 8 : if(master) {
185 8 : nrep_ = multi_sim_comm.Get_size();
186 8 : replica_ = multi_sim_comm.Get_rank();
187 8 : if(noensemble) {
188 0 : nrep_ = 1;
189 : }
190 : } else {
191 0 : nrep_ = 0;
192 0 : replica_ = 0;
193 : }
194 8 : comm.Sum(&nrep_,1);
195 8 : comm.Sum(&replica_,1);
196 :
197 : const unsigned narg = getNumberOfArguments();
198 8 : sigma_mean2_.resize(narg,1);
199 8 : sigma_mean2_last_.resize(narg);
200 16 : for(unsigned j=0; j<narg; j++) {
201 8 : sigma_mean2_last_[j].push_back(0.000001);
202 : }
203 :
204 8 : log.printf(" Time resolved data from file %s\n",filename.c_str());
205 8 : std::ifstream varfile(filename.c_str());
206 8 : if(varfile.fail()) {
207 0 : error("Cannot open "+filename);
208 : }
209 8 : var.resize(narg);
210 8 : dvar.resize(narg);
211 4024 : while (!varfile.eof()) {
212 : double tempT, tempVar;
213 : varfile >> tempT;
214 4016 : time.push_back(tempT/tscale);
215 8032 : for(unsigned i=0; i<narg; i++) {
216 : varfile >> tempVar;
217 4016 : var[i].push_back(tempVar);
218 : }
219 : }
220 8 : varfile.close();
221 :
222 8 : const double deltat = time[1] - time[0];
223 16 : for(unsigned i=0; i<narg; i++) {
224 4024 : for(unsigned j=0; j<var[i].size(); j++) {
225 4016 : if(j==0) {
226 8 : dvar[i].push_back((var[i][j+1] - var[i][j])/(deltat));
227 4008 : } else if(j==var[i].size()-1) {
228 8 : dvar[i].push_back((var[i][j] - var[i][j-1])/(deltat));
229 : } else {
230 4000 : dvar[i].push_back((var[i][j+1] - var[i][j-1])/(2.*deltat));
231 : }
232 : }
233 : }
234 :
235 16 : for(unsigned i=0; i<narg; i++) {
236 : std::string num;
237 8 : Tools::convert(i,num);
238 16 : addComponent("x0-"+num);
239 8 : componentIsNotPeriodic("x0-"+num);
240 8 : x0comp.push_back(getPntrToComponent("x0-"+num));
241 16 : addComponent("kappa-"+num);
242 8 : componentIsNotPeriodic("kappa-"+num);
243 8 : kcomp.push_back(getPntrToComponent("kappa-"+num));
244 16 : addComponent("mean-"+num);
245 8 : componentIsNotPeriodic("mean-"+num);
246 8 : mcomp.push_back(getPntrToComponent("mean-"+num));
247 : }
248 :
249 8 : if(doregres_zero_) {
250 8 : addComponent("scale");
251 4 : componentIsNotPeriodic("scale");
252 4 : valueScale=getPntrToComponent("scale");
253 : }
254 :
255 16 : log<<" Bibliography "<<plumed.cite("Capelli, Tiana, Camilloni, J Chem Phys, 148, 184114");
256 16 : }
257 :
258 0 : void Caliber::get_sigma_mean(const double fact, const std::vector<double> &mean) {
259 0 : const unsigned narg = getNumberOfArguments();
260 0 : const double dnrep = static_cast<double>(nrep_);
261 :
262 0 : if(sigma_mean2_last_[0].size()==optsigmamean_stride_)
263 0 : for(unsigned i=0; i<narg; ++i) {
264 0 : sigma_mean2_last_[i].erase(sigma_mean2_last_[i].begin());
265 : }
266 0 : std::vector<double> sigma_mean2_now(narg,0);
267 0 : if(master) {
268 0 : for(unsigned i=0; i<narg; ++i) {
269 0 : double tmp = getArgument(i)-mean[i];
270 0 : sigma_mean2_now[i] = fact*tmp*tmp;
271 : }
272 0 : if(nrep_>1) {
273 0 : multi_sim_comm.Sum(&sigma_mean2_now[0], narg);
274 : }
275 : }
276 0 : comm.Sum(&sigma_mean2_now[0], narg);
277 :
278 0 : for(unsigned i=0; i<narg; ++i) {
279 0 : sigma_mean2_last_[i].push_back(sigma_mean2_now[i]/dnrep);
280 0 : sigma_mean2_[i] = *max_element(sigma_mean2_last_[i].begin(), sigma_mean2_last_[i].end());
281 : }
282 0 : }
283 :
284 4008 : void Caliber::replica_averaging(const double fact, std::vector<double> &mean) {
285 4008 : const unsigned narg = getNumberOfArguments();
286 4008 : if(master) {
287 8016 : for(unsigned i=0; i<narg; ++i) {
288 4008 : mean[i] = fact*getArgument(i);
289 : }
290 4008 : if(nrep_>1) {
291 4008 : multi_sim_comm.Sum(&mean[0], narg);
292 : }
293 : }
294 4008 : comm.Sum(&mean[0], narg);
295 4008 : }
296 :
297 6012 : double Caliber::getSpline(const unsigned iarg) {
298 6012 : const double deltat = time[1] - time[0];
299 6012 : const int tindex = static_cast<int>(getTime()/deltat);
300 :
301 : unsigned start, end;
302 6012 : start=tindex;
303 6012 : if(tindex+1<var[iarg].size()) {
304 6012 : end=tindex+2;
305 : } else {
306 0 : end=var[iarg].size();
307 : }
308 :
309 : double value=0;
310 18036 : for(unsigned ipoint=start; ipoint<end; ++ipoint) {
311 12024 : double grid=var[iarg][ipoint];
312 12024 : double dder=dvar[iarg][ipoint];
313 : double yy=0.;
314 12024 : if(std::abs(grid)>0.0000001) {
315 12024 : yy=-dder/grid;
316 : }
317 :
318 : int x0=1;
319 12024 : if(ipoint==tindex) {
320 : x0=0;
321 : }
322 :
323 12024 : double X=std::abs((getTime()-time[tindex])/deltat-(double)x0);
324 12024 : double X2=X*X;
325 12024 : double X3=X2*X;
326 12024 : double C=(1.0-3.0*X2+2.0*X3) - (x0?-1.0:1.0)*yy*(X-2.0*X2+X3)*deltat;
327 :
328 12024 : value+=grid*C;
329 : }
330 6012 : return value;
331 : }
332 :
333 2004 : void Caliber::do_regression_zero(const std::vector<double> &mean) {
334 : // parameters[i] = scale_ * mean[i]: find scale_ with linear regression
335 : double num = 0.0;
336 : double den = 0.0;
337 4008 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
338 2004 : num += mean[i] * getSpline(i);
339 2004 : den += mean[i] * mean[i];
340 : }
341 2004 : if(den>0) {
342 2004 : scale_ = num / den;
343 : } else {
344 0 : scale_ = 1.0;
345 : }
346 2004 : }
347 :
348 4008 : void Caliber::calculate() {
349 4008 : const unsigned narg = getNumberOfArguments();
350 4008 : const double dnrep = static_cast<double>(nrep_);
351 4008 : const double fact = 1.0/dnrep;
352 :
353 4008 : std::vector<double> mean(narg,0);
354 4008 : std::vector<double> dmean_x(narg,fact);
355 4008 : replica_averaging(fact, mean);
356 4008 : if(optsigmamean_stride_>0) {
357 0 : get_sigma_mean(fact, mean);
358 : }
359 :
360 : // in case of regression with zero intercept, calculate scale
361 4008 : if(doregres_zero_ && getStep()%nregres_zero_==0) {
362 2004 : do_regression_zero(mean);
363 : }
364 :
365 : double ene=0;
366 8016 : for(unsigned i=0; i<narg; ++i) {
367 4008 : const double x0 = getSpline(i);
368 4008 : const double kappa = mult*dnrep/sigma_mean2_[i];
369 4008 : const double cv=difference(i,x0,scale_*mean[i]);
370 4008 : const double f=-kappa*cv*dmean_x[i]/scale_;
371 4008 : setOutputForce(i,f);
372 4008 : ene+=0.5*kappa*cv*cv;
373 4008 : x0comp[i]->set(x0);
374 4008 : kcomp[i]->set(kappa);
375 4008 : mcomp[i]->set(mean[i]);
376 : }
377 :
378 4008 : if(doregres_zero_) {
379 2004 : valueScale->set(scale_);
380 : }
381 :
382 4008 : setBias(ene);
383 4008 : }
384 :
385 : }
386 : }
387 :
388 :
|
