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