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