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