Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2012-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 "FlexibleBin.h"
23 : #include "ActionWithArguments.h"
24 : #include <cmath>
25 : #include <vector>
26 : #include "tools/Matrix.h"
27 :
28 : namespace PLMD {
29 :
30 :
31 49 : FlexibleBin::FlexibleBin(int mytype, ActionWithArguments *awargs, double const &d, std::vector<double> &smin, const std::vector<double> &smax):
32 49 : type(mytype),
33 49 : paction(awargs),
34 49 : sigma(d),
35 49 : sigmamin(smin),
36 49 : sigmamax(smax) {
37 49 : paction->log<<" Limits for sigmas using adaptive hills: \n";
38 141 : for(unsigned i=0; i<paction->getNumberOfArguments(); ++i) {
39 92 : paction->log<<" CV "<<paction->getPntrToArgument(i)->getName()<<":\n";
40 92 : if(sigmamin[i]>0.) {
41 18 : limitmin.push_back(true);
42 18 : paction->log<<" Min "<<sigmamin[i];
43 18 : sigmamin[i]*=sigmamin[i]; // this is because the matrix which is calculated is the sigmasquared
44 : } else {
45 74 : limitmin.push_back(false);
46 74 : paction->log<<" Min No ";
47 : }
48 92 : if(sigmamax[i]>0.) {
49 4 : limitmax.push_back(true);
50 4 : paction->log<<" Max "<<sigmamax[i];
51 4 : sigmamax[i]*=sigmamax[i];
52 : } else {
53 88 : limitmax.push_back(false);
54 88 : paction->log<<" Max No ";
55 : }
56 92 : paction->log<<" \n";
57 : }
58 :
59 49 : }
60 :
61 : /// Constructure for 1D FB for PBMETAD
62 38 : FlexibleBin::FlexibleBin(int mytype, ActionWithArguments *awargs, unsigned iarg,
63 38 : double const &d, std::vector<double> &smin, const std::vector<double> &smax):
64 38 : type(mytype),
65 38 : paction(awargs),
66 38 : sigma(d),
67 38 : sigmamin(smin),
68 38 : sigmamax(smax) {
69 38 : paction->log<<" Limits for sigmas using adaptive hills: \n";
70 38 : paction->log<<" CV "<<paction->getPntrToArgument(iarg)->getName()<<":\n";
71 38 : if(sigmamin[0]>0.) {
72 17 : limitmin.push_back(true);
73 17 : paction->log<<" Min "<<sigmamin[0];
74 17 : sigmamin[0]*=sigmamin[0];
75 : } else {
76 21 : limitmin.push_back(false);
77 21 : paction->log<<" Min No ";
78 : }
79 38 : if(sigmamax[0]>0.) {
80 3 : limitmax.push_back(true);
81 3 : paction->log<<" Max "<<sigmamax[0];
82 3 : sigmamax[0]*=sigmamax[0];
83 : } else {
84 35 : limitmax.push_back(false);
85 35 : paction->log<<" Max No ";
86 : }
87 38 : paction->log<<" \n";
88 38 : }
89 :
90 : /// Update the flexible bin
91 : /// in case of diffusion based: update at every step
92 : /// in case of gradient based: update only when you add the hill
93 8488 : void FlexibleBin::update(bool nowAddAHill) {
94 8488 : unsigned ncv=paction->getNumberOfArguments();
95 8488 : unsigned dimension=ncv*(ncv+1)/2;
96 8488 : std::vector<double> delta;
97 8488 : std::vector<double> cv;
98 8488 : double decay=1./sigma;
99 : // this is done all the times from scratch. It is not an accumulator
100 : // here update the flexible bin according to the needs
101 8488 : switch (type) {
102 : // This should be called every time
103 5566 : case diffusion:
104 : // if you use this below then the decay is in time units
105 : //double decay=paction->getTimeStep()/sigma;
106 : // to be consistent with the rest of the program: everything is better to be in timesteps
107 : // THE AVERAGE VALUE
108 : // beware: the pbc
109 5566 : delta.resize(ncv);
110 14912 : for(unsigned i=0; i<ncv; i++) {
111 9346 : cv.push_back(paction->getArgument(i));
112 : }
113 : bool restarted;
114 : restarted=false;
115 5566 : if(average.size()==0) { // initial time: just set the initial vector
116 19 : average.resize(ncv);
117 47 : for(unsigned i=0; i<ncv; i++) {
118 28 : average[i]=cv[i];
119 : }
120 19 : restarted=(variance.size()>0); // after restartAverage() the variance is kept as it is
121 : } else { // accumulate
122 14865 : for(unsigned i=0; i<ncv; i++) {
123 9318 : delta[i]=paction->difference(i,average[i],cv[i]);
124 9318 : average[i]+=decay*delta[i];
125 9318 : average[i]=paction->bringBackInPbc(i,average[i]); // equation 8 of "Metadynamics with adaptive Gaussians"
126 : }
127 : }
128 5566 : nsamples++;
129 : // THE VARIANCE
130 5566 : if(variance.size()==0) {
131 15 : variance.resize(dimension,0.); // nonredundant members dimension=ncv*(ncv+1)/2;
132 5551 : } else if(!restarted) {
133 : unsigned k=0;
134 14865 : for(unsigned i=0; i<ncv; i++) {
135 22407 : for(unsigned j=i; j<ncv; j++) { // upper diagonal loop
136 13089 : variance[k]+=decay*(delta[i]*delta[j]-variance[k]);
137 13089 : k++;
138 : }
139 : }
140 : }
141 : break;
142 2922 : case geometry:
143 : //this calculates in variance the \nabla CV_i \dot \nabla CV_j
144 2922 : variance.resize(dimension);
145 : // now the signal for retrieving the gradients should be already given by checkNeedsGradients.
146 : // here just do the projections
147 : // note that the call checkNeedsGradients() in BiasMetaD takes care of switching on the call to gradients
148 2922 : if (nowAddAHill) { // geometry is sync with hill deposition
149 : unsigned k=0;
150 2709 : for(unsigned i=0; i<ncv; i++) {
151 4515 : for(unsigned j=i; j<ncv; j++) {
152 : // eq 12 of "Metadynamics with adaptive Gaussians"
153 2709 : variance[k]=sigma*sigma*(paction->getProjection(i,j));
154 2709 : k++;
155 : }
156 : }
157 : }
158 : break;
159 0 : default:
160 0 : plumed_merror("This flexible bin method is not recognized");
161 : }
162 8488 : }
163 :
164 12 : void FlexibleBin::restartAverage() {
165 12 : plumed_massert(type==diffusion,"the average can be restarted only for diffusion-based adaptive hills");
166 12 : average.clear();
167 12 : }
168 :
169 15 : void FlexibleBin::setStatistics(const std::vector<double>& av, const std::vector<double>& var, unsigned ns) {
170 15 : plumed_massert(type==diffusion,"the statistics can be set only for diffusion-based adaptive hills");
171 15 : const unsigned ncv=paction->getNumberOfArguments();
172 : // (one dimensional, for one argument of PBMETAD, or for all the arguments)
173 15 : plumed_assert((av.size()==1 && var.size()==1) || (av.size()==ncv && var.size()==ncv*(ncv+1)/2));
174 15 : average=av;
175 15 : variance=var;
176 15 : nsamples=ns;
177 15 : }
178 :
179 18 : std::vector<double> FlexibleBin::getMatrix() const {
180 18 : return variance;
181 : }
182 :
183 : /// Update the flexible bin for PBMetaD like FlexBin
184 : /// in case of diffusion based: update at every step
185 : /// in case of gradient based: update only when you add the hill
186 8644 : void FlexibleBin::update(bool nowAddAHill, unsigned iarg) {
187 : // this is done all the times from scratch. It is not an accumulator
188 : // here update the flexible bin according to the needs
189 8644 : std::vector<double> cv;
190 8644 : std::vector<double> delta;
191 : // if you use this below then the decay is in time units
192 : // to be consistent with the rest of the program: everything is better to be in timesteps
193 8644 : double decay=1./sigma;
194 8644 : switch (type) {
195 : // This should be called every time
196 7006 : case diffusion:
197 : // THE AVERAGE VALUE
198 7006 : delta.resize(1);
199 7006 : cv.push_back(paction->getArgument(iarg));
200 : bool restarted;
201 : restarted=false;
202 7006 : if(average.size()==0) { // initial time: just set the initial vector
203 37 : average.resize(1);
204 37 : average[0]=cv[0];
205 37 : restarted=(variance.size()>0); // after restartAverage() the variance is kept as it is
206 : } else { // accumulate
207 6969 : delta[0]=paction->difference(iarg,average[0],cv[0]);
208 6969 : average[0]+=decay*delta[0];
209 6969 : average[0]=paction->bringBackInPbc(iarg,average[0]); // equation 8 of "Metadynamics with adaptive Gaussians"
210 : }
211 7006 : nsamples++;
212 : // THE VARIANCE
213 7006 : if(variance.size()==0) {
214 29 : variance.resize(1,0.); // nonredundant members dimension=ncv*(ncv+1)/2;
215 6977 : } else if(!restarted) {
216 6969 : variance[0]+=decay*(delta[0]*delta[0]-variance[0]);
217 : }
218 : break;
219 1638 : case geometry:
220 : //this calculates in variance the \nabla CV_i \dot \nabla CV_j
221 1638 : variance.resize(1);
222 : // now the signal for retrieving the gradients should be already given by checkNeedsGradients.
223 : // here just do the projections
224 : // note that the call checkNeedsGradients() in BiasMetaD takes care of switching on the call to gradients
225 1638 : if (nowAddAHill) { // geometry is sync with hill deposition
226 : // eq 12 of "Metadynamics with adaptive Gaussians"
227 327 : variance[0]=sigma*sigma*(paction->getProjection(iarg,iarg));
228 : }
229 : break;
230 0 : default:
231 0 : plumed_merror("This flexible bin is not recognized");
232 : }
233 8644 : }
234 :
235 : ///
236 : /// Calculate the matrix of (dcv_i/dx)*(dcv_j/dx)^-1
237 : /// that is needed for the metrics in metadynamics
238 : ///
239 : ///
240 1954 : std::vector<double> FlexibleBin::getInverseMatrix() const {
241 1954 : unsigned ncv=paction->getNumberOfArguments();
242 1954 : Matrix<double> matrix(ncv,ncv);
243 : unsigned i,j,k;
244 : k=0;
245 : // place the matrix in a complete matrix for compatibility
246 5455 : for (i=0; i<ncv; i++) {
247 8549 : for (j=i; j<ncv; j++) {
248 5048 : matrix(j,i)=matrix(i,j)=variance[k];
249 5048 : k++;
250 : }
251 : }
252 : #define NEWFLEX
253 : #ifdef NEWFLEX
254 : // diagonalize to impose boundaries (only if boundaries are set)
255 1954 : Matrix<double> eigenvecs(ncv,ncv);
256 1954 : std::vector<double> eigenvals(ncv);
257 :
258 : //eigenvecs: first is eigenvec number, second is eigenvec component
259 1954 : if(diagMat( matrix, eigenvals, eigenvecs )!=0) {
260 0 : plumed_merror("diagonalization in FlexibleBin failed! This matrix is weird\n");
261 : };
262 :
263 5455 : for (i=0; i<ncv; i++) { //loop on the dimension
264 3501 : if( limitmax[i] ) {
265 : //limit every component that is larger
266 1308 : for (j=0; j<ncv; j++) { //loop on components
267 872 : if(std::pow(eigenvals[j]*eigenvecs[j][i],2)>std::pow(sigmamax[i],2) ) {
268 552 : eigenvals[j]=std::sqrt(std::pow(sigmamax[i]/(eigenvecs[j][i]),2))*copysign(1.,eigenvals[j]);
269 : }
270 : }
271 : }
272 : }
273 5455 : for (i=0; i<ncv; i++) { //loop on the dimension
274 : // find the largest one: if it is smaller than min then rescale
275 3501 : if( limitmin[i] ) {
276 : unsigned imax=0;
277 : double fmax=-1.e10;
278 3186 : for (j=0; j<ncv; j++) { //loop on components
279 2124 : double fact=std::pow(eigenvals[j]*eigenvecs[j][i],2);
280 2124 : if(fact>fmax) {
281 : fmax=fact;
282 : imax=j;
283 : }
284 : }
285 1062 : if(fmax<std::pow(sigmamin[i],2) ) {
286 226 : eigenvals[imax]=std::sqrt(std::pow(sigmamin[i]/(eigenvecs[imax][i]),2))*copysign(1.,eigenvals[imax]);
287 : }
288 : }
289 : }
290 :
291 : // normalize eigenvecs
292 1954 : Matrix<double> newinvmatrix(ncv,ncv);
293 5455 : for (i=0; i<ncv; i++) {
294 10096 : for (j=0; j<ncv; j++) {
295 6595 : newinvmatrix[j][i]=eigenvecs[j][i]/eigenvals[j];
296 : }
297 : }
298 :
299 1954 : std::vector<double> uppervec(ncv*(ncv+1)/2);
300 : k=0;
301 5455 : for (i=0; i<ncv; i++) {
302 8549 : for (j=i; j<ncv; j++) {
303 : double scal=0;
304 14737 : for(unsigned l=0; l<ncv; ++l) {
305 9689 : scal+=eigenvecs[l][i]*newinvmatrix[l][j];
306 : }
307 5048 : uppervec[k]=scal;
308 5048 : k++;
309 : }
310 : }
311 : #else
312 : // get the inverted matrix
313 : Matrix<double> invmatrix(ncv,ncv);
314 : Invert(matrix,invmatrix);
315 : std::vector<double> uppervec(ncv*(ncv+1)/2);
316 : // upper diagonal of the inverted matrix (that is symmetric)
317 : k=0;
318 : for (i=0; i<ncv; i++) {
319 : for (j=i; j<ncv; j++) {
320 : uppervec[k]=invmatrix(i,j);
321 : k++;
322 : }
323 : }
324 : #endif
325 1954 : return uppervec;
326 1954 : }
327 :
328 : ///
329 : /// Calculate the matrix of (dcv_i/dx)*(dcv_j/dx)^-1
330 : /// that is needed for the metrics in metadynamics
331 : /// for PBMetaD like FlexBin
332 : ///
333 1517 : std::vector<double> FlexibleBin::getInverseMatrix(unsigned iarg) const {
334 : // diagonalize to impose boundaries (only if boundaries are set)
335 1517 : std::vector<double> eigenvals(1, variance[0]);
336 1517 : if( limitmax[0] ) {
337 327 : if(eigenvals[0]>sigmamax[0]) {
338 218 : eigenvals[0]=sigmamax[0];
339 : }
340 : }
341 : // find the largest one: if it is smaller than min then rescale
342 1517 : if( limitmin[0] ) {
343 : double fmax=-1.e10;
344 783 : double fact=eigenvals[0];
345 783 : if(fact>fmax) {
346 : fmax=fact;
347 : }
348 783 : if(fmax<sigmamin[0]) {
349 102 : eigenvals[0]=sigmamin[0];
350 : }
351 : }
352 1517 : std::vector<double> uppervec(1,1./eigenvals[0]);
353 :
354 1517 : return uppervec;
355 1517 : }
356 :
357 : }
|
