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