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