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 : #ifndef __PLUMED_bias_MetaDHills_h
23 : #define __PLUMED_bias_MetaDHills_h
24 :
25 : #include "MetaDKernel.h"
26 : #include "tools/Communicator.h"
27 : #include <deque>
28 : #include <functional>
29 : #include <limits>
30 : #include <memory>
31 : #include <string>
32 : #include <vector>
33 :
34 : namespace PLMD {
35 :
36 : class ActionWithValue;
37 : class IFile;
38 : class OFile;
39 : class Value;
40 :
41 : namespace bias {
42 :
43 : class MetaDBiasGrid;
44 : class MetaDDiffusionRestart;
45 :
46 : /// shortest representation of x that is converted back to the same double
47 : std::string exactString(double x);
48 :
49 : /// broadcast a string from process 0: the length is sent first, as on the other processes the string
50 : /// can have a different length (Communicator::Bcast of a string uses the local length)
51 : void bcastString(Communicator& c, std::string& s);
52 :
53 : /// remove the last line of a hills file if it is incomplete (left by a run that was interrupted while
54 : /// writing it); returns the message for the log (empty if nothing was done)
55 : std::string removeIncompleteLastLine(const std::string& path);
56 :
57 : /// The hills files of METAD and PBMETAD for a set of arguments: writing of the hills, and reading of the files
58 : /// of a restart or of the other walkers.
59 : ///
60 : /// With several MPI processes the files are opened and read by the first process only, which also decides if
61 : /// they exist, and the hills are broadcast to the others: the file of another walker can grow while it is read,
62 : /// and if each process read it they could read a different number of hills, and their biases would differ.
63 : /// The file operations do not synchronize the processes (the files are linked to a communicator of the process
64 : /// alone). Everything written on the log is written by all the processes.
65 : class MetaDHillsFiles {
66 : public:
67 : /// called for each hill that is read: time, multivariate, height, center, sigma (the widths, or the upper
68 : /// triangle of the inverse covariance matrix for multivariate hills)
69 : using Callback=std::function<void(double,bool,double,const std::vector<double>&,const std::vector<double>&)>;
70 : /// The bias to which the hills that are read are added (see readHills)
71 1728 : struct Bias {
72 : /// add a hill to the bias
73 : std::function<void(const MetaDGaussian&)> add;
74 : /// the grid of the bias, if any (it completes the line of the log when it is extended)
75 : MetaDBiasGrid* grid=nullptr;
76 : /// ADAPTIVE=DIFF: when it is pending, the hills of the own file are kept here for the statistics of the fluctuations
77 : MetaDDiffusionRestart* diffrestart=nullptr;
78 : /// the height of a hill is the one in the file times heightscale (well-tempered metadynamics)
79 : double heightscale=1.0;
80 : /// times closer than this are the same time (times are written with a finite precision): half a PACE
81 : double tolerance=0.0;
82 : };
83 : private:
84 : ActionWithValue& action_;
85 : Communicator& comm_;
86 : /// communicator of the process alone, linked to the input files
87 : Communicator selfcomm_;
88 : std::vector<Value*> args_;
89 : std::vector<std::unique_ptr<IFile>> ifiles_;
90 : std::vector<std::string> ifilesnames_;
91 : /// number of hills read from each file
92 : std::vector<unsigned long long> nread_;
93 : Bias bias_;
94 : bool noStretchWarningDone_=false;
95 : /// a hill with a Gaussian (not stretched) kernel was read
96 : bool gaussianKernelsFound_=false;
97 : void noStretchWarning();
98 : /// an error in a hills file
99 : [[noreturn]] void inputError(const std::string& msg) const;
100 : /// the numbers of hills of the other walkers in a grid file (on one process)
101 : static std::vector<long long> readGridComment(const std::string& path, unsigned nwalkers);
102 : bool scanOneHill(IFile* ifile, std::vector<Value>& v, std::vector<double>& center, std::vector<double>& sigma, double& height,
103 : bool& multivariate, double& time);
104 : public:
105 : MetaDHillsFiles(ActionWithValue& action, Communicator& comm, const std::vector<Value*>& args);
106 : /// the name of the hills file of a walker: name.walker with WALKERS_N, in the directory dir if given
107 : /// (with WALKERS_N or WALKERS_MPI)
108 : static std::string getWalkerFileName(const std::string& name, const std::string& dir, int nwalkers, int walker, bool mpi);
109 : ~MetaDHillsFiles();
110 : /// add an input file, which is opened if it exists; returns its index. With nosuffix, the name is used without
111 : /// the suffix of the replica (with multiple walkers, as for the files that are written)
112 : unsigned addInput(const std::string& name, bool nosuffix=false);
113 : const std::string& getInputName(unsigned i) const {
114 6715 : return ifilesnames_[i];
115 : }
116 : /// true if the file is open (the same on all the processes)
117 : bool isOpen(unsigned i);
118 : /// path of an open file (on the first process only, empty on the others)
119 : std::string getPath(unsigned i) const;
120 : /// open the file if it exists now; returns true if it has been opened
121 : bool openIfExists(unsigned i);
122 : /// read the hills that are in the file (an open file) up to now and call f for each one; returns their number.
123 : /// An error while reading (e.g. a wrong field) is thrown by all the processes.
124 : unsigned read(unsigned i, const Callback& f);
125 : /// With WALKERS_N, the grid written by a walker contains the hills it has read from the files of the other walkers:
126 : /// their numbers are written in a comment line after the grid, so that when restarting from the grid only the
127 : /// hills that are not in it are added. The comment line (own is the index of this walker)
128 : std::string getGridComment(unsigned own) const;
129 : /// the numbers of hills of the other walkers in the grid file path (written by getGridComment), -1 for this walker;
130 : /// empty if the file does not contain them
131 : /// (read by the first process of comm and broadcast)
132 : static std::vector<long long> readGridComment(Communicator& comm, const std::string& path, unsigned nwalkers);
133 : void setBias(const Bias& bias) {
134 432 : bias_=bias;
135 432 : }
136 : const Bias& getBias() const {
137 : return bias_;
138 : }
139 : /// Read the hills of file i that are in the file up to now (an open file) and add them to the bias, writing their
140 : /// number on the log (the caller has started the line). Returns their number.
141 : /// lasttime: if given, the time of the last hill is returned, and the hills are kept for ADAPTIVE=DIFF (own file);
142 : /// addtobias: false to only scan the file; the hills up to time skipupto, and the first skipcount, are not added
143 : /// (they are already in the grid read from file)
144 : unsigned readHills(unsigned i, double* lasttime=nullptr, bool addtobias=true,
145 : double skipupto=-std::numeric_limits<double>::infinity(), unsigned long long skipcount=0);
146 : /// prepare an open file for the next read (the hills read up to now are not read again)
147 : void reset(unsigned i);
148 : void close(unsigned i);
149 : /// write a hill; height is the one written in the file
150 : void write(OFile& file, double time, const MetaDGaussian& hill, double height, double biasf, bool clock) const;
151 : };
152 :
153 : /// ADAPTIVE=DIFF: when restarting, the statistics of the fluctuations of the variables (running average and
154 : /// covariance) are initialized from the last hills of this walker in its hills file. The hills are kept here while the
155 : /// file is read, until the first update, when the current time is known.
156 : class MetaDDiffusionRestart {
157 : std::deque<double> time_;
158 : std::deque<std::vector<double>> center_;
159 : /// as stored in the hills: upper triangle of the inverse covariance matrix
160 : std::deque<std::vector<double>> sigma_;
161 : std::string file_;
162 : public:
163 : /// the statistics are still to be initialized
164 : bool pending=false;
165 : /// maximum number of hills to keep (the older ones have a negligible weight)
166 : unsigned maxhills=0;
167 18 : void clear() {
168 18 : time_.clear();
169 18 : center_.clear();
170 18 : sigma_.clear();
171 18 : }
172 : /// a hill read from file (hills that are not multivariate cannot be used: the previous ones are discarded)
173 : void add(const std::string& file, double time, bool multivariate, const std::vector<double>& center, const std::vector<double>& sigma);
174 : const std::string& getFile() const {
175 : return file_;
176 : }
177 : /// The statistics from the hills of this walker. The hills with times within tol form blocks (one deposition):
178 : /// with blocksize 0 all the hills are used, otherwise the hill at position offset of the blocks with blocksize hills.
179 : /// Only the hills up to now, and the last sequence with increasing times, are used. The covariance is that of the
180 : /// last hill; the average is the average of the centers with the exponential weights of the running average
181 : /// (decay per update, updatetime is the time between updates). difference(i,a,b) and bringBackInPbc(i,x) are those
182 : /// of variable i. Returns the number of hills used (0: no statistics), and the time of the last one.
183 : unsigned compute(double tol, double now, unsigned blocksize, unsigned offset, double decay, double updatetime,
184 : const std::function<double(unsigned,double,double)>& difference,
185 : const std::function<double(unsigned,double)>& bringBackInPbc,
186 : std::vector<double>& average, std::vector<double>& variance, double& lasttime) const;
187 : };
188 :
189 : /// The steps when a hill is added: every stride steps, except the first step of a run. At step 0 no hill is
190 : /// added (also when the bias of a previous simulation is read with RESTART). When continuing a simulation, the
191 : /// hill is not added only if the hills file already contains it, i.e. if the previous run already performed this
192 : /// step (some MD codes compute again the last step of the previous run when restarting, others do not).
193 350 : class MetaDHillSchedule {
194 : bool runStarted_=false;
195 : long long int runStartStep_=0;
196 : /// the time of the last hill in the hills file read when restarting
197 : bool lastHillTimeKnown_=false;
198 : double lastHillTime_=0.0;
199 : public:
200 : void setLastHillTime(double time) {
201 62 : lastHillTimeKnown_=true;
202 62 : lastHillTime_=time;
203 62 : }
204 : /// to be called at every calculate (the first call starts the run)
205 : void start(long long int step) {
206 52496 : if(!runStarted_) {
207 330 : runStarted_=true;
208 330 : runStartStep_=step;
209 : }
210 : }
211 : /// true before the first calculate and at the first step of the run
212 : bool isRunFirstStep(long long int step) const {
213 14689 : return !runStarted_ || step==runStartStep_;
214 : }
215 : /// tolerance: times are written with a finite precision (FMT), so that a hill within tolerance (half a PACE)
216 : /// is considered to be the hill of this step; a later hill means that the previous run went beyond this step
217 53156 : bool isHillStep(long long int step, long long int stride, double time, double tolerance) const {
218 53156 : if(step%stride!=0) {
219 : return false;
220 : }
221 : if(isRunFirstStep(step)) {
222 348 : if(step==0 || !lastHillTimeKnown_) {
223 : return false;
224 : }
225 39 : return lastHillTime_<time-tolerance;
226 : }
227 : return true;
228 : }
229 : };
230 :
231 : }
232 : }
233 :
234 : #endif
|
