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