Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2011-2026 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 "MetaDHills.h"
23 : #include "MetaDGrid.h"
24 : #include "core/ActionWithValue.h"
25 : #include "core/Value.h"
26 : #include "tools/Exception.h"
27 : #include "tools/File.h"
28 : #include "tools/Log.h"
29 : #include "tools/Matrix.h"
30 : #include "tools/Tools.h"
31 : #include <cstdio>
32 : #include <algorithm>
33 : #include <cmath>
34 : #include <cstdlib>
35 : #include <ctime>
36 : #include <exception>
37 : #include <filesystem>
38 : #include <fstream>
39 : #include <system_error>
40 :
41 : namespace PLMD {
42 : namespace bias {
43 :
44 211 : std::string exactString(double x) {
45 : char buffer[64];
46 211 : for(int precision=15; precision<17; precision++) {
47 : std::snprintf(buffer,sizeof(buffer),"%.*g",precision,x);
48 211 : if(std::strtod(buffer,nullptr)==x) {
49 211 : return buffer;
50 : }
51 : }
52 : std::snprintf(buffer,sizeof(buffer),"%.17g",x);
53 0 : return buffer;
54 : }
55 :
56 363 : void bcastString(Communicator& c, std::string& s) {
57 363 : unsigned n=s.size();
58 363 : c.Bcast(n,0);
59 363 : s.resize(n);
60 363 : if(n>0) {
61 173 : c.Bcast(s,0);
62 : }
63 363 : }
64 :
65 63 : std::string removeIncompleteLastLine(const std::string& path) {
66 : // compressed files cannot be fixed in place
67 90 : if(path.length()>3 && path.substr(path.length()-3)==".gz") {
68 3 : return "";
69 : }
70 60 : std::ifstream f(path,std::ios::binary);
71 60 : if(!f) {
72 0 : return "";
73 : }
74 60 : f.seekg(0,std::ios::end);
75 60 : const std::streamoff size=f.tellg();
76 60 : if(size<=0) {
77 0 : return "";
78 : }
79 : char last;
80 60 : f.seekg(size-1);
81 60 : f.get(last);
82 60 : if(last=='\n') {
83 56 : return "";
84 : }
85 : // look backward for the end of the last complete line
86 : std::streamoff pos=size-1;
87 255 : while(pos>0) {
88 : char c;
89 255 : f.seekg(pos-1);
90 255 : f.get(c);
91 255 : if(c=='\n') {
92 : break;
93 : }
94 : pos--;
95 : }
96 4 : std::string tail(size-pos,' ');
97 4 : f.seekg(pos);
98 4 : f.read(&tail[0],size-pos);
99 4 : f.close();
100 4 : std::string message=" WARNING: the last line of "+path+" is incomplete (the previous run was probably interrupted while writing it) and is removed:\n";
101 4 : message+=" "+tail+"\n";
102 : std::error_code ec;
103 8 : std::filesystem::resize_file(path,pos,ec);
104 4 : if(ec) {
105 0 : message+=" WARNING: the incomplete last line of "+path+" could not be removed ("+ec.message()+")\n";
106 : }
107 4 : return message;
108 64 : }
109 :
110 :
111 432 : MetaDHillsFiles::MetaDHillsFiles(ActionWithValue& action, Communicator& comm, const std::vector<Value*>& args):
112 432 : action_(action),
113 432 : comm_(comm),
114 432 : args_(args) {
115 432 : }
116 :
117 432 : MetaDHillsFiles::~MetaDHillsFiles() {
118 864 : }
119 :
120 468 : std::string MetaDHillsFiles::getWalkerFileName(const std::string& name, const std::string& dir, int nwalkers, int walker, bool mpi) {
121 468 : const std::string suffix=(nwalkers>1 ? "."+std::to_string(walker) : std::string());
122 468 : if(dir!="" && (nwalkers>1 || mpi)) {
123 334 : return dir+"/"+name+suffix;
124 : }
125 301 : return name+suffix;
126 : }
127 :
128 468 : unsigned MetaDHillsFiles::addInput(const std::string& name, bool nosuffix) {
129 468 : ifiles_.emplace_back(Tools::make_unique<IFile>());
130 468 : ifilesnames_.push_back(name);
131 468 : nread_.push_back(0);
132 : IFile* ifile=ifiles_.back().get();
133 468 : ifile->link(action_);
134 468 : if(nosuffix) {
135 448 : ifile->enforceSuffix("");
136 : }
137 : // the file operations must not synchronize the processes, as with comm they would do
138 468 : ifile->link(selfcomm_);
139 468 : const unsigned i=ifiles_.size()-1;
140 468 : openIfExists(i);
141 468 : return i;
142 : }
143 :
144 6692 : bool MetaDHillsFiles::isOpen(unsigned i) {
145 6692 : int isopen=ifiles_[i]->isOpen();
146 6692 : if(comm_.Get_size()>1) {
147 402 : comm_.Bcast(isopen,0);
148 : }
149 6692 : return isopen;
150 : }
151 :
152 6359 : std::string MetaDHillsFiles::getPath(unsigned i) const {
153 6359 : return ifiles_[i]->isOpen() ? ifiles_[i]->getPath() : std::string();
154 : }
155 :
156 501 : bool MetaDHillsFiles::openIfExists(unsigned i) {
157 501 : int exists=0;
158 501 : if(comm_.Get_rank()==0) {
159 395 : exists=ifiles_[i]->FileExist(ifilesnames_[i]);
160 : }
161 501 : if(comm_.Get_size()>1) {
162 210 : comm_.Bcast(exists,0);
163 : }
164 501 : if(exists && comm_.Get_rank()==0) {
165 110 : ifiles_[i]->open(ifilesnames_[i]);
166 : }
167 501 : return exists;
168 : }
169 :
170 6283 : unsigned MetaDHillsFiles::readHills(unsigned i, double* lasttime, bool addtobias, double skipupto, unsigned long long skipcount) {
171 6283 : Log& log=action_.log;
172 : // the callers have started a line of the log
173 6283 : if(bias_.grid) {
174 : bias_.grid->setMidline(true);
175 : }
176 6283 : const std::string path=getPath(i);
177 6283 : unsigned long long nskipped=0;
178 : unsigned nhills=0;
179 : try {
180 6283 : nhills=read(i,[&](double time,bool multivariate,double height,const std::vector<double>& center,const std::vector<double>& sigma) {
181 6308 : if(lasttime) {
182 3110 : *lasttime=time;
183 : }
184 6308 : if(lasttime && bias_.diffrestart && bias_.diffrestart->pending) {
185 : // the own hills file is being read: keep the last hills for the statistics of ADAPTIVE=DIFF
186 1068 : bias_.diffrestart->add(path,time,multivariate,center,sigma);
187 : }
188 6308 : if(!addtobias) {
189 : return;
190 : }
191 5142 : if(nskipped<skipcount || time<=skipupto+bias_.tolerance) {
192 313 : nskipped++;
193 313 : return;
194 : }
195 4829 : bias_.add(MetaDGaussian(multivariate,height*bias_.heightscale,center,sigma));
196 : });
197 2 : } catch(...) {
198 2 : if(bias_.grid) {
199 : bias_.grid->setMidline(false);
200 : }
201 2 : throw;
202 2 : }
203 6281 : if(bias_.grid) {
204 : bias_.grid->setMidline(false);
205 : }
206 6281 : if(!addtobias) {
207 18 : log.printf(" %u Gaussians scanned (the bias is read from the grid)\n",nhills);
208 6263 : } else if(nskipped>0) {
209 8 : log.printf(" %u Gaussians read, %llu of them already in the grid\n",nhills,nskipped);
210 : } else {
211 6255 : log.printf(" %u Gaussians read\n",nhills);
212 : }
213 6281 : return nhills;
214 : }
215 :
216 6315 : void MetaDHillsFiles::reset(unsigned i) {
217 6315 : ifiles_[i]->reset(false);
218 6315 : }
219 :
220 88 : void MetaDHillsFiles::close(unsigned i) {
221 88 : ifiles_[i]->close();
222 88 : }
223 :
224 0 : void MetaDHillsFiles::inputError(const std::string& msg) const {
225 : // as Action::error, but without writing on the log: the files are read by the first process only, and
226 : // the log requires the same output on all the processes (read throws the error on all of them)
227 0 : plumed_merror("ERROR in input to action "+action_.getName()+" with label "+action_.getLabel()+" : "+msg);
228 : }
229 :
230 3 : void MetaDHillsFiles::noStretchWarning() {
231 3 : if(!noStretchWarningDone_) {
232 3 : action_.log<<"\nWARNING: you are using a HILLS file with Gaussian kernels, PLUMED 2.8 uses stretched Gaussians by default\n";
233 : }
234 3 : noStretchWarningDone_=true;
235 3 : }
236 :
237 6283 : unsigned MetaDHillsFiles::read(unsigned i, const Callback& f) {
238 6283 : IFile* ifile=ifiles_[i].get();
239 6283 : const unsigned ncv=args_.size();
240 6283 : std::vector<double> center(ncv);
241 6285 : std::vector<double> sigma(ncv);
242 : double height;
243 : unsigned nhills=0;
244 6283 : bool multivariate=false;
245 :
246 6283 : std::vector<Value> tmpvalues;
247 18632 : for(unsigned j=0; j<ncv; ++j) {
248 24698 : tmpvalues.push_back( Value( &action_, args_[j]->getName(), false ) );
249 : }
250 :
251 : // The hills are read and broadcast in chunks; each one is stored as
252 : // time, multivariate, height, center (ncv), sigma (ncv or ncv*(ncv+1)/2).
253 6283 : const unsigned nsigmamax=ncv*(ncv+1)/2;
254 6283 : const std::size_t recordsize=3+ncv+nsigmamax;
255 : const unsigned chunksize=4096;
256 6283 : const bool parallel=(comm_.Get_size()>1);
257 6283 : const bool reader=(comm_.Get_rank()==0);
258 6283 : std::vector<double> buffer;
259 : bool more=true;
260 12564 : while(more) {
261 6283 : unsigned nread=0;
262 6283 : if(reader) {
263 6180 : buffer.clear();
264 6180 : int failed=0;
265 : std::string message;
266 : std::exception_ptr exception;
267 : try {
268 : double time;
269 12077 : while(nread<chunksize && (more=scanOneHill(ifile,tmpvalues,center,sigma,height,multivariate,time))) {
270 5897 : buffer.push_back(time);
271 10555 : buffer.push_back(multivariate ? 1.0 : 0.0);
272 5897 : buffer.push_back(height);
273 5897 : buffer.insert(buffer.end(),center.begin(),center.end());
274 5897 : buffer.insert(buffer.end(),sigma.begin(),sigma.end());
275 5897 : buffer.resize(buffer.size()+nsigmamax-sigma.size(),0.0);
276 5897 : nread++;
277 : }
278 1 : } catch(const std::exception& e) {
279 1 : failed=1;
280 1 : message=e.what();
281 1 : exception=std::current_exception();
282 1 : }
283 6180 : if(parallel) {
284 : // the other processes must know if reading failed, so that all of them stop
285 103 : comm_.Bcast(failed,0);
286 103 : if(failed) {
287 1 : bcastString(comm_,message);
288 : // all the processes throw the same exception: its message can be written on the log,
289 : // which requires the same output on all the processes
290 2 : throw ExceptionError(message);
291 : }
292 : }
293 6179 : if(failed) {
294 1 : std::rethrow_exception(exception);
295 : }
296 : } else {
297 103 : int failed=0;
298 103 : comm_.Bcast(failed,0);
299 103 : if(failed) {
300 : std::string message;
301 1 : bcastString(comm_,message);
302 1 : throw ExceptionError(message);
303 : }
304 : }
305 : // the warning about Gaussian kernels is written on the log after reading, by all the processes
306 : // (the log requires the same output on all the processes)
307 6281 : int gaussiankernels=gaussianKernelsFound_;
308 6281 : if(parallel) {
309 204 : int imore=more;
310 204 : comm_.Bcast(nread,0);
311 204 : comm_.Bcast(imore,0);
312 204 : comm_.Bcast(gaussiankernels,0);
313 204 : more=imore;
314 204 : buffer.resize(nread*recordsize);
315 204 : if(nread>0) {
316 166 : comm_.Bcast(buffer,0);
317 : }
318 : }
319 6281 : if(gaussiankernels) {
320 3 : noStretchWarning();
321 : }
322 12589 : for(unsigned h=0; h<nread; h++) {
323 6308 : const double* record=&buffer[h*recordsize];
324 6308 : const double time=record[0];
325 6308 : multivariate=(record[1]!=0.0);
326 6308 : height=record[2];
327 6308 : center.assign(record+3,record+3+ncv);
328 6308 : sigma.assign(record+3+ncv,record+3+ncv+(multivariate ? nsigmamax : ncv));
329 6308 : nhills++;
330 6308 : f(time,multivariate,height,center,sigma);
331 : }
332 : }
333 6281 : nread_[i]+=nhills;
334 6281 : return nhills;
335 6283 : }
336 :
337 : namespace {
338 : const std::string gridCommentHeader="# WALKERS_N: hills read from the hills files of the walkers:";
339 : }
340 :
341 0 : std::string MetaDHillsFiles::getGridComment(unsigned own) const {
342 0 : std::string line=gridCommentHeader;
343 0 : for(unsigned i=0; i<nread_.size(); i++) {
344 0 : line+=" "+(i==own ? std::string("-") : std::to_string(nread_[i]));
345 : }
346 0 : return line+"\n";
347 : }
348 :
349 8 : std::vector<long long> MetaDHillsFiles::readGridComment(Communicator& comm, const std::string& path, unsigned nwalkers) {
350 8 : std::vector<long long> counts;
351 8 : if(comm.Get_rank()==0) {
352 16 : counts=readGridComment(path,nwalkers);
353 : }
354 8 : if(comm.Get_size()>1) {
355 : // (as doubles: they are exact up to 2^53)
356 0 : std::vector<double> buffer(counts.begin(),counts.end());
357 0 : unsigned n=buffer.size();
358 0 : comm.Bcast(n,0);
359 0 : buffer.resize(n);
360 0 : if(n>0) {
361 0 : comm.Bcast(buffer,0);
362 : }
363 0 : counts.assign(buffer.begin(),buffer.end());
364 : }
365 8 : return counts;
366 0 : }
367 :
368 8 : std::vector<long long> MetaDHillsFiles::readGridComment(const std::string& path, unsigned nwalkers) {
369 8 : std::vector<long long> counts;
370 : // compressed files are not read (the hills of the other walkers are then selected by time)
371 16 : if(path.length()>3 && path.substr(path.length()-3)==".gz") {
372 : return counts;
373 : }
374 8 : std::ifstream f(path);
375 : std::string line;
376 : // the comment follows the grid (with STORE_GRIDS, the first grid is the one that is read)
377 31422 : while(std::getline(f,line)) {
378 31420 : if(line.compare(0,gridCommentHeader.length(),gridCommentHeader)!=0) {
379 31414 : continue;
380 : }
381 12 : std::vector<std::string> words=Tools::getWords(line.substr(gridCommentHeader.length()));
382 6 : if(words.size()!=nwalkers) {
383 : return counts;
384 : }
385 18 : for(const auto & w : words) {
386 12 : long long n=-1;
387 12 : if(w!="-" && !Tools::convertNoexcept(w,n)) {
388 0 : counts.clear();
389 0 : return counts;
390 : }
391 12 : counts.push_back(n);
392 : }
393 : return counts;
394 6 : }
395 : return counts;
396 8 : }
397 :
398 : /// read the next hill: tmpvalues are used to read the CVs (with their periodicity); returns false at the end of the file
399 12077 : bool MetaDHillsFiles::scanOneHill(IFile* ifile, std::vector<Value>& tmpvalues, std::vector<double>& center, std::vector<double>& sigma,
400 : double& height, bool& multivariate, double& time) {
401 : double dummy;
402 12077 : multivariate=false;
403 24153 : if(ifile->scanField("time",time)) {
404 5897 : unsigned ncv=tmpvalues.size();
405 15447 : for(unsigned i=0; i<ncv; ++i) {
406 9550 : ifile->scanField( &tmpvalues[i] );
407 9550 : if( tmpvalues[i].isPeriodic() && ! args_[i]->isPeriodic() ) {
408 0 : inputError("in hills file periodicity for variable " + tmpvalues[i].getName() + " does not match periodicity in input");
409 9550 : } else if( tmpvalues[i].isPeriodic() ) {
410 : std::string imin, imax;
411 3263 : tmpvalues[i].getDomain( imin, imax );
412 : std::string rmin, rmax;
413 3263 : args_[i]->getDomain( rmin, rmax );
414 3263 : if( imin!=rmin || imax!=rmax ) {
415 0 : inputError("in hills file periodicity for variable " + tmpvalues[i].getName() + " does not match periodicity in input");
416 : }
417 : }
418 9550 : center[i]=tmpvalues[i].get();
419 : }
420 : // scan for kerneltype
421 5897 : std::string ktype="stretched-gaussian";
422 11794 : if( ifile->FieldExist("kerneltype") ) {
423 11776 : ifile->scanField("kerneltype",ktype);
424 : }
425 5897 : if( ktype=="gaussian" ) {
426 : // the warning is written by read
427 12 : gaussianKernelsFound_=true;
428 5885 : } else if( ktype!="stretched-gaussian") {
429 0 : inputError("non Gaussian kernels are not supported in MetaD");
430 : }
431 : // scan for multivariate label
432 : std::string sss;
433 5897 : ifile->scanField("multivariate",sss);
434 5897 : if(sss=="true") {
435 1239 : multivariate=true;
436 4658 : } else if(sss=="false") {
437 4658 : multivariate=false;
438 : } else {
439 0 : plumed_merror("cannot parse multivariate = "+ sss);
440 : }
441 5897 : if(multivariate) {
442 1239 : sigma.resize(ncv*(ncv+1)/2);
443 1239 : Matrix<double> upper(ncv,ncv);
444 1239 : Matrix<double> lower(ncv,ncv);
445 3126 : for(unsigned i=0; i<ncv; i++) {
446 4422 : for(unsigned j=0; j<ncv-i; j++) {
447 7605 : ifile->scanField("sigma_"+args_[j+i]->getName()+"_"+args_[j]->getName(),lower(j+i,j));
448 2535 : upper(j,j+i)=lower(j+i,j);
449 : }
450 : }
451 1239 : Matrix<double> mymult(ncv,ncv);
452 1239 : Matrix<double> invmatrix(ncv,ncv);
453 1239 : mult(lower,upper,mymult);
454 : // now invert and get the sigmas
455 1239 : Invert(mymult,invmatrix);
456 : // put the sigmas in the usual order: upper diagonal (this time in normal form and not in band form)
457 : unsigned k=0;
458 3126 : for(unsigned i=0; i<ncv; i++) {
459 4422 : for(unsigned j=i; j<ncv; j++) {
460 2535 : sigma[k]=invmatrix(i,j);
461 2535 : k++;
462 : }
463 : }
464 : } else {
465 12321 : for(unsigned i=0; i<ncv; ++i) {
466 15326 : ifile->scanField("sigma_"+args_[i]->getName(),sigma[i]);
467 : }
468 : }
469 :
470 5897 : ifile->scanField("height",height);
471 5897 : ifile->scanField("biasf",dummy);
472 11794 : if(ifile->FieldExist("clock")) {
473 7736 : ifile->scanField("clock",dummy);
474 : }
475 11794 : if(ifile->FieldExist("lower_int")) {
476 0 : ifile->scanField("lower_int",dummy);
477 : }
478 11794 : if(ifile->FieldExist("upper_int")) {
479 0 : ifile->scanField("upper_int",dummy);
480 : }
481 5897 : ifile->scanField();
482 : return true;
483 : } else {
484 : return false;
485 : }
486 : }
487 :
488 19377 : void MetaDHillsFiles::write(OFile& file, double time, const MetaDGaussian& hill, double height, double biasf, bool clock) const {
489 19377 : const unsigned ncv=args_.size();
490 19377 : file.printField("time",time);
491 46853 : for(unsigned i=0; i<ncv; ++i) {
492 27476 : file.printField(args_[i],hill.center[i]);
493 : }
494 38754 : file.printField("kerneltype","stretched-gaussian");
495 19377 : if(hill.multivariate) {
496 7310 : file.printField("multivariate","true");
497 3655 : Matrix<double> mymatrix(ncv,ncv);
498 : unsigned k=0;
499 8969 : for(unsigned i=0; i<ncv; i++) {
500 12287 : for(unsigned j=i; j<ncv; j++) {
501 : // recompose the full inverse matrix
502 6973 : mymatrix(i,j)=mymatrix(j,i)=hill.sigma[k];
503 6973 : k++;
504 : }
505 : }
506 : // invert it
507 3655 : Matrix<double> invmatrix(ncv,ncv);
508 3655 : Invert(mymatrix,invmatrix);
509 : // enforce symmetry
510 8969 : for(unsigned i=0; i<ncv; i++) {
511 12287 : for(unsigned j=i; j<ncv; j++) {
512 6973 : invmatrix(i,j)=invmatrix(j,i);
513 : }
514 : }
515 :
516 : // do cholesky so to have a "sigma like" number
517 3655 : Matrix<double> lower(ncv,ncv);
518 3655 : cholesky(invmatrix,lower);
519 : // loop in band form
520 8969 : for(unsigned i=0; i<ncv; i++) {
521 12287 : for(unsigned j=0; j<ncv-i; j++) {
522 20919 : file.printField("sigma_"+args_[j+i]->getName()+"_"+args_[j]->getName(),lower(j+i,j));
523 : }
524 : }
525 : } else {
526 31444 : file.printField("multivariate","false");
527 37884 : for(unsigned i=0; i<ncv; ++i) {
528 44324 : file.printField("sigma_"+args_[i]->getName(),hill.sigma[i]);
529 : }
530 : }
531 38754 : file.printField("height",height).printField("biasf",biasf);
532 19377 : if(clock) {
533 4442 : file.printField("clock",int(std::time(0)));
534 : }
535 19377 : file.printField();
536 19377 : }
537 :
538 1068 : void MetaDDiffusionRestart::add(const std::string& file, double time, bool multivariate, const std::vector<double>& center,
539 : const std::vector<double>& sigma) {
540 1068 : if(!multivariate) {
541 : // e.g. a previous run without ADAPTIVE=DIFF: these hills cannot be used
542 0 : clear();
543 : } else {
544 1068 : time_.push_back(time);
545 1068 : center_.push_back(center);
546 1068 : sigma_.push_back(sigma);
547 1068 : if(time_.size()>maxhills) {
548 0 : time_.pop_front();
549 0 : center_.pop_front();
550 0 : sigma_.pop_front();
551 : }
552 : }
553 1068 : file_=file;
554 1068 : }
555 :
556 18 : unsigned MetaDDiffusionRestart::compute(double tol, double now, unsigned blocksize, unsigned offset, double decay, double updatetime,
557 : const std::function<double(unsigned,double,double)>& difference,
558 : const std::function<double(unsigned,double)>& bringBackInPbc,
559 : std::vector<double>& average, std::vector<double>& variance, double& lasttime) const {
560 : // the hills of this walker
561 18 : std::vector<unsigned> own;
562 968 : for(unsigned k=0; k<time_.size();) {
563 : unsigned l=k;
564 2294 : while(l<time_.size() && std::abs(time_[l]-time_[k])<tol) {
565 1344 : l++;
566 : }
567 950 : if(blocksize==0) {
568 836 : for(unsigned j=k; j<l; j++) {
569 418 : own.push_back(j);
570 : }
571 532 : } else if(l-k==blocksize) {
572 532 : own.push_back(k+offset);
573 : }
574 : k=l;
575 : }
576 : // only the hills deposited up to now (when restarting from a grid, or from a checkpoint older than the
577 : // end of the hills file, the file can go further)
578 : {
579 18 : std::vector<unsigned> before;
580 968 : for(const auto j : own) {
581 950 : if(time_[j]<=now+tol) {
582 836 : before.push_back(j);
583 : }
584 : }
585 : own.swap(before);
586 : }
587 : // only the last sequence of increasing times: the file can contain several runs (e.g. a run restarted
588 : // from an older checkpoint, or a new simulation appended to an old file)
589 18 : if(!own.empty()) {
590 15 : unsigned first=own.size()-1;
591 836 : while(first>0 && time_[own[first-1]]<time_[own[first]]+tol) {
592 : first--;
593 : }
594 : own.erase(own.begin(),own.begin()+first);
595 : }
596 18 : if(own.empty()) {
597 : return 0;
598 : }
599 : // covariance: the one of the last hill, i.e. the running estimate of the previous run at that time
600 15 : const unsigned last=own.back();
601 15 : const unsigned ncv=center_[last].size();
602 15 : Matrix<double> mymatrix(ncv,ncv);
603 : unsigned k=0;
604 39 : for(unsigned i=0; i<ncv; i++) {
605 57 : for(unsigned j=i; j<ncv; j++) {
606 33 : mymatrix(i,j)=mymatrix(j,i)=sigma_[last][k];
607 33 : k++;
608 : }
609 : }
610 15 : Matrix<double> mycov(ncv,ncv);
611 15 : Invert(mymatrix,mycov);
612 15 : variance.resize(ncv*(ncv+1)/2);
613 : k=0;
614 39 : for(unsigned i=0; i<ncv; i++) {
615 57 : for(unsigned j=i; j<ncv; j++) {
616 33 : variance[k]=mycov(i,j);
617 33 : k++;
618 : }
619 : }
620 : // running average: the center of a hill is the instantaneous value of the CVs, which differs from
621 : // the running average by about one standard deviation, so the centers of the last hills are averaged
622 : // with the same exponential weights used by the running average
623 : const std::vector<double>& ref=center_[last];
624 15 : average.assign(ncv,0.0);
625 : double wsum=0.0;
626 851 : for(const auto j : own) {
627 : // age in number of updates of the running average
628 836 : const double age=std::max(0.0,(time_[last]-time_[j])/updatetime);
629 836 : const double w=std::pow(1.0-decay,age);
630 2151 : for(unsigned i=0; i<ncv; i++) {
631 1315 : average[i]+=w*difference(i,ref[i],center_[j][i]);
632 : }
633 836 : wsum+=w;
634 : }
635 39 : for(unsigned i=0; i<ncv; i++) {
636 24 : average[i]=bringBackInPbc(i,ref[i]+average[i]/wsum);
637 : }
638 15 : lasttime=time_[last];
639 15 : return own.size();
640 : }
641 :
642 : }
643 : }
|
