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 : #include <algorithm>
28 : #include <cmath>
29 :
30 : namespace PLMD {
31 : namespace isdb {
32 :
33 : //+PLUMEDOC ISDB_BIAS CALIBER
34 : /*
35 : Add a time-dependent, harmonic restraint on one or more variables.
36 :
37 : This allows implementing a maximum caliber restraint on one or more experimental time series by replica-averaged restrained simulations.
38 :
39 : The time resolved experiments are read from a text file and intermediate values are obtained by splines.
40 :
41 : ## Examples
42 :
43 : In the following example a restraint is applied on the time evolution of a saxs spectrum
44 :
45 : ```plumed
46 : MOLINFO STRUCTURE=first.pdb
47 :
48 : # Define saxs variable
49 : SAXS ...
50 : LABEL=saxs
51 : ATOMISTIC
52 : ATOMS=1-436
53 : QVALUE1=0.02 # Q-value at which calculate the scattering
54 : QVALUE2=0.0808
55 : QVALUE3=0.1264
56 : QVALUE4=0.1568
57 : QVALUE5=0.172
58 : QVALUE6=0.1872
59 : QVALUE7=0.2176
60 : QVALUE8=0.2328
61 : QVALUE9=0.248
62 : QVALUE10=0.2632
63 : QVALUE11=0.2936
64 : QVALUE12=0.3088
65 : QVALUE13=0.324
66 : QVALUE14=0.3544
67 : QVALUE15=0.4
68 : ... SAXS
69 :
70 :
71 : #define the caliber restraint
72 : CALIBER ...
73 : ARG=(saxs\.q_.*)
74 : FILE=expsaxs.dat
75 : KAPPA=10
76 : LABEL=cal0
77 : STRIDE=10
78 : REGRES_ZERO=200
79 : AVERAGING=200
80 : ... CALIBER
81 : ```
82 :
83 : 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.
84 : 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
85 : Because [SAXS](SAXS.md) 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.
86 : The bias is here applied every tenth step.
87 :
88 : */
89 : //+ENDPLUMEDOC
90 :
91 :
92 : class Caliber : public bias::Bias {
93 : public:
94 : explicit Caliber(const ActionOptions&);
95 : void calculate();
96 : static void registerKeywords( Keywords& keys );
97 : private:
98 : std::vector<double> time;
99 : std::vector< std::vector<double> > var;
100 : std::vector< std::vector<double> > dvar;
101 : double mult;
102 : double scale_;
103 : bool master;
104 : unsigned replica_;
105 : unsigned nrep_;
106 : // scale and offset regression
107 : bool doregres_zero_;
108 : int nregres_zero_;
109 : // force constant
110 : unsigned optsigmamean_stride_;
111 : std::vector<double> sigma_mean2_;
112 : std::vector< std::vector<double> > sigma_mean2_last_;
113 : std::vector<Value*> x0comp;
114 : std::vector<Value*> kcomp;
115 : std::vector<Value*> mcomp;
116 : Value* valueScale;
117 :
118 : void get_sigma_mean(const double fact, const std::vector<double> &mean);
119 : void replica_averaging(const double fact, std::vector<double> &mean);
120 : double getSpline(const unsigned iarg);
121 : void do_regression_zero(const std::vector<double> &mean);
122 : };
123 :
124 : PLUMED_REGISTER_ACTION(Caliber,"CALIBER")
125 :
126 20 : void Caliber::registerKeywords( Keywords& keys ) {
127 20 : Bias::registerKeywords(keys);
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 40 : keys.addOutputComponent("x0","default","scalar","the instantaneous value of the center of the potential");
136 40 : keys.addOutputComponent("mean","default","scalar","the current average value of the calculated observable");
137 40 : keys.addOutputComponent("kappa","default","scalar","the current force constant");
138 40 : keys.addOutputComponent("scale","REGRES_ZERO","scalar","the current scaling constant");
139 20 : keys.addDOI("10.1063/1.5030339");
140 20 : }
141 :
142 16 : Caliber::Caliber(const ActionOptions&ao):
143 : PLUMED_BIAS_INIT(ao),
144 16 : mult(0),
145 16 : scale_(1),
146 16 : doregres_zero_(false),
147 16 : nregres_zero_(0),
148 16 : optsigmamean_stride_(0) {
149 32 : parse("KAPPA",mult);
150 : std::string filename;
151 32 : parse("FILE",filename);
152 16 : if( filename.length()==0 ) {
153 0 : error("No external variable file was specified");
154 : }
155 16 : unsigned averaging=0;
156 16 : parse("AVERAGING", averaging);
157 16 : if(averaging>0) {
158 4 : optsigmamean_stride_ = averaging;
159 : }
160 16 : double tscale=1.0;
161 16 : parse("TSCALE", tscale);
162 16 : if(tscale<=0.) {
163 0 : error("The time scale factor must be greater than 0.");
164 : }
165 16 : parse("SCALE", scale_);
166 16 : if(scale_==0.) {
167 0 : error("The time scale factor cannot be 0.");
168 : }
169 : // regression with zero intercept
170 16 : parse("REGRES_ZERO", nregres_zero_);
171 16 : if(nregres_zero_>0) {
172 : // set flag
173 8 : doregres_zero_=true;
174 8 : log.printf(" doing regression with zero intercept with stride: %d\n", nregres_zero_);
175 : }
176 :
177 :
178 16 : bool noensemble = false;
179 16 : parseFlag("NOENSEMBLE", noensemble);
180 :
181 16 : checkRead();
182 :
183 : // set up replica stuff
184 16 : master = (comm.Get_rank()==0);
185 16 : nrep_ = getNumberOfReplicas();
186 16 : replica_ = getReplicaIndex();
187 16 : if(noensemble) {
188 0 : nrep_ = 1;
189 : }
190 :
191 : const unsigned narg = getNumberOfArguments();
192 16 : sigma_mean2_.resize(narg,1);
193 16 : sigma_mean2_last_.resize(narg);
194 36 : for(unsigned j=0; j<narg; j++) {
195 20 : sigma_mean2_last_[j].push_back(0.000001);
196 : }
197 :
198 16 : log.printf(" Time resolved data from file %s\n",filename.c_str());
199 16 : std::ifstream varfile(filename.c_str());
200 16 : if(varfile.fail()) {
201 0 : error("Cannot open "+filename);
202 : }
203 16 : var.resize(narg);
204 16 : dvar.resize(narg);
205 : double tempT;
206 8032 : while(varfile >> tempT) {
207 8016 : time.push_back(tempT/tscale);
208 18036 : for(unsigned i=0; i<narg; i++) {
209 : double tempVar;
210 10020 : if(!(varfile >> tempVar)) {
211 0 : error("the last line of "+filename+" does not contain a value for each argument");
212 : }
213 10020 : var[i].push_back(tempVar);
214 : }
215 : }
216 16 : if(!varfile.eof()) {
217 0 : error("cannot read "+filename+": only numbers are allowed");
218 : }
219 16 : varfile.close();
220 16 : if(time.size()<2) {
221 0 : error("the file "+filename+" must contain at least two times");
222 : }
223 :
224 16 : const double deltat = time[1] - time[0];
225 16 : if(!(deltat>0.)) {
226 0 : error("the times in "+filename+" must be increasing");
227 : }
228 : // getSpline assumes equally spaced times
229 8016 : for(unsigned j=1; j<time.size(); j++) {
230 8000 : if(std::abs(time[j]-time[j-1]-deltat)>1.e-6*deltat) {
231 0 : error("the times in "+filename+" must be equally spaced");
232 : }
233 : }
234 36 : for(unsigned i=0; i<narg; i++) {
235 10040 : for(unsigned j=0; j<var[i].size(); j++) {
236 10020 : if(j==0) {
237 20 : dvar[i].push_back((var[i][j+1] - var[i][j])/(deltat));
238 10000 : } else if(j==var[i].size()-1) {
239 20 : dvar[i].push_back((var[i][j] - var[i][j-1])/(deltat));
240 : } else {
241 9980 : dvar[i].push_back((var[i][j+1] - var[i][j-1])/(2.*deltat));
242 : }
243 : }
244 : }
245 :
246 36 : for(unsigned i=0; i<narg; i++) {
247 : std::string num;
248 20 : Tools::convert(i,num);
249 40 : addComponent("x0-"+num);
250 20 : componentIsNotPeriodic("x0-"+num);
251 40 : x0comp.push_back(getPntrToComponent("x0-"+num));
252 40 : addComponent("kappa-"+num);
253 20 : componentIsNotPeriodic("kappa-"+num);
254 40 : kcomp.push_back(getPntrToComponent("kappa-"+num));
255 40 : addComponent("mean-"+num);
256 20 : componentIsNotPeriodic("mean-"+num);
257 40 : mcomp.push_back(getPntrToComponent("mean-"+num));
258 : }
259 :
260 16 : if(doregres_zero_) {
261 16 : addComponent("scale");
262 8 : componentIsNotPeriodic("scale");
263 8 : valueScale=getPntrToComponent("scale");
264 : }
265 :
266 16 : log<<" Bibliography "<<plumed.cite("Capelli, Tiana, Camilloni, J Chem Phys, 148, 184114");
267 32 : }
268 :
269 2004 : void Caliber::get_sigma_mean(const double fact, const std::vector<double> &mean) {
270 2004 : const unsigned narg = getNumberOfArguments();
271 2004 : const double dnrep = static_cast<double>(nrep_);
272 :
273 2004 : if(sigma_mean2_last_[0].size()==optsigmamean_stride_)
274 5904 : for(unsigned i=0; i<narg; ++i) {
275 3936 : sigma_mean2_last_[i].erase(sigma_mean2_last_[i].begin());
276 : }
277 2004 : std::vector<double> sigma_mean2_now(narg,0);
278 2004 : if(master) {
279 6012 : for(unsigned i=0; i<narg; ++i) {
280 4008 : double tmp = getArgument(i)-mean[i];
281 4008 : sigma_mean2_now[i] = fact*tmp*tmp;
282 : }
283 2004 : if(nrep_>1) {
284 2004 : multi_sim_comm.Sum(&sigma_mean2_now[0], narg);
285 : }
286 : }
287 2004 : comm.Sum(&sigma_mean2_now[0], narg);
288 :
289 6012 : for(unsigned i=0; i<narg; ++i) {
290 4008 : sigma_mean2_last_[i].push_back(sigma_mean2_now[i]/dnrep);
291 4008 : sigma_mean2_[i] = *max_element(sigma_mean2_last_[i].begin(), sigma_mean2_last_[i].end());
292 : }
293 2004 : }
294 :
295 8016 : void Caliber::replica_averaging(const double fact, std::vector<double> &mean) {
296 8016 : const unsigned narg = getNumberOfArguments();
297 8016 : if(master) {
298 18036 : for(unsigned i=0; i<narg; ++i) {
299 10020 : mean[i] = fact*getArgument(i);
300 : }
301 8016 : if(nrep_>1) {
302 8016 : multi_sim_comm.Sum(&mean[0], narg);
303 : }
304 : }
305 8016 : comm.Sum(&mean[0], narg);
306 8016 : }
307 :
308 12832 : double Caliber::getSpline(const unsigned iarg) {
309 12832 : const double deltat = time[1] - time[0];
310 : // index of the last data point before the current time (times are equally spaced)
311 12832 : const double tpos = (getTime()-time[0])/deltat;
312 12832 : if(tpos<-1.e-9 || tpos>static_cast<double>(time.size()-1)+1.e-9) {
313 0 : error("the simulation time is outside the time range of the experimental data");
314 : }
315 12832 : const unsigned tindex = std::min(static_cast<unsigned>(std::max(tpos,0.)),
316 12832 : static_cast<unsigned>(time.size()-1));
317 :
318 : unsigned start, end;
319 : start=tindex;
320 12832 : if(tindex+1<var[iarg].size()) {
321 12832 : end=tindex+2;
322 : } else {
323 0 : end=var[iarg].size();
324 : }
325 :
326 : double value=0;
327 38496 : for(unsigned ipoint=start; ipoint<end; ++ipoint) {
328 25664 : double grid=var[iarg][ipoint];
329 25664 : double dder=dvar[iarg][ipoint];
330 : double yy=0.;
331 25664 : if(std::abs(grid)>0.0000001) {
332 25664 : yy=-dder/grid;
333 : }
334 :
335 : int x0=1;
336 25664 : if(ipoint==tindex) {
337 : x0=0;
338 : }
339 :
340 25664 : double X=std::abs((getTime()-time[tindex])/deltat-(double)x0);
341 25664 : double X2=X*X;
342 25664 : double X3=X2*X;
343 25664 : double C=(1.0-3.0*X2+2.0*X3) - (x0?-1.0:1.0)*yy*(X-2.0*X2+X3)*deltat;
344 :
345 25664 : value+=grid*C;
346 : }
347 12832 : return value;
348 : }
349 :
350 2408 : void Caliber::do_regression_zero(const std::vector<double> &mean) {
351 : // parameters[i] = scale_ * mean[i]: find scale_ with linear regression
352 : double num = 0.0;
353 : double den = 0.0;
354 5220 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
355 2812 : num += mean[i] * getSpline(i);
356 2812 : den += mean[i] * mean[i];
357 : }
358 2408 : if(den>0) {
359 2408 : scale_ = num / den;
360 : } else {
361 0 : scale_ = 1.0;
362 : }
363 2408 : }
364 :
365 8016 : void Caliber::calculate() {
366 8016 : const unsigned narg = getNumberOfArguments();
367 8016 : const double dnrep = static_cast<double>(nrep_);
368 8016 : const double fact = 1.0/dnrep;
369 :
370 8016 : std::vector<double> mean(narg,0);
371 8016 : std::vector<double> dmean_x(narg,fact);
372 8016 : replica_averaging(fact, mean);
373 8016 : if(optsigmamean_stride_>0) {
374 2004 : get_sigma_mean(fact, mean);
375 : }
376 :
377 : // in case of regression with zero intercept, calculate scale
378 8016 : if(doregres_zero_ && getStep()%nregres_zero_==0) {
379 2408 : do_regression_zero(mean);
380 : }
381 :
382 : double ene=0;
383 18036 : for(unsigned i=0; i<narg; ++i) {
384 10020 : const double x0 = getSpline(i);
385 10020 : const double kappa = mult*dnrep/sigma_mean2_[i];
386 10020 : const double cv=difference(i,x0,scale_*mean[i]);
387 : // ene = 0.5*kappa*cv^2 with cv = scale*mean-x0 and d(mean)/d(arg) = dmean_x
388 10020 : const double f=-kappa*cv*dmean_x[i]*scale_;
389 10020 : setOutputForce(i,f);
390 10020 : ene+=0.5*kappa*cv*cv;
391 10020 : x0comp[i]->set(x0);
392 10020 : kcomp[i]->set(kappa);
393 10020 : mcomp[i]->set(mean[i]);
394 : }
395 :
396 8016 : if(doregres_zero_) {
397 4008 : valueScale->set(scale_);
398 : }
399 :
400 8016 : setBias(ene);
401 8016 : }
402 :
403 : }
404 : }
405 :
406 :
|
