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 66 : 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 106 : for(unsigned i=0; i<paction->getNumberOfArguments(); ++i) {
42 84 : 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 0 : 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 : vector<double> delta;
102 : vector<double> cv;
103 778 : double decay=1./sigma;
104 : // this is done all the times from scratch. It is not an accumulator
105 : // here update the flexible bin according to the needs
106 778 : switch (type) {
107 : // This should be called every time
108 586 : case diffusion:
109 : // if you use this below then the decay is in time units
110 : //double decay=paction->getTimeStep()/sigma;
111 : // to be consistent with the rest of the program: everything is better to be in timesteps
112 : // THE AVERAGE VALUE
113 : // beware: the pbc
114 586 : delta.resize(ncv);
115 2344 : for(unsigned i=0; i<ncv; i++) cv.push_back(paction->getArgument(i));
116 586 : if(average.size()==0) { // initial time: just set the initial vector
117 2 : average.resize(ncv);
118 6 : for(unsigned i=0; i<ncv; i++) average[i]=cv[i];
119 : } else { // accumulate
120 1752 : for(unsigned i=0; i<ncv; i++) {
121 2336 : delta[i]=paction->difference(i,average[i],cv[i]);
122 1168 : average[i]+=decay*delta[i];
123 1168 : average[i]=paction->bringBackInPbc(i,average[i]); // equation 8 of "Metadynamics with adaptive Gaussians"
124 : }
125 : }
126 : // THE VARIANCE
127 586 : if(variance.size()==0) {
128 2 : variance.resize(dimension,0.); // nonredundant members dimension=ncv*(ncv+1)/2;
129 : } else {
130 : unsigned k=0;
131 1752 : for(unsigned i=0; i<ncv; i++) {
132 1752 : for(unsigned j=i; j<ncv; j++) { // upper diagonal loop
133 2336 : variance[k]+=decay*(delta[i]*delta[j]-variance[k]);
134 584 : k++;
135 : }
136 : }
137 : }
138 : break;
139 192 : case geometry:
140 : //this calculates in variance the \nabla CV_i \dot \nabla CV_j
141 192 : variance.resize(dimension);
142 : // now the signal for retrieving the gradients should be already given by checkNeedsGradients.
143 : // here just do the projections
144 : // note that the call checkNeedsGradients() in BiasMetaD takes care of switching on the call to gradients
145 192 : if (nowAddAHill) { // geometry is sync with hill deposition
146 : unsigned k=0;
147 415 : for(unsigned i=0; i<ncv; i++) {
148 664 : for(unsigned j=i; j<ncv; j++) {
149 : // eq 12 of "Metadynamics with adaptive Gaussians"
150 498 : variance[k]=sigma*sigma*(paction->getProjection(i,j));
151 249 : k++;
152 : }
153 : }
154 : }
155 : break;
156 0 : default:
157 0 : plumed_merror("This flexible bin method is not recognized");
158 : }
159 778 : }
160 :
161 0 : vector<double> FlexibleBin::getMatrix() const {
162 0 : return variance;
163 : }
164 :
165 : /// Update the flexible bin for PBMetaD like FlexBin
166 : /// in case of diffusion based: update at every step
167 : /// in case of gradient based: update only when you add the hill
168 0 : void FlexibleBin::update(bool nowAddAHill, unsigned iarg) {
169 : // this is done all the times from scratch. It is not an accumulator
170 : // here update the flexible bin according to the needs
171 : vector<double> cv;
172 : vector<double> delta;
173 : // if you use this below then the decay is in time units
174 : // to be consistent with the rest of the program: everything is better to be in timesteps
175 0 : double decay=1./sigma;
176 0 : switch (type) {
177 : // This should be called every time
178 0 : case diffusion:
179 : // THE AVERAGE VALUE
180 0 : delta.resize(1);
181 0 : cv.push_back(paction->getArgument(iarg));
182 0 : if(average.size()==0) { // initial time: just set the initial vector
183 0 : average.resize(1);
184 0 : average[0]=cv[0];
185 : } else { // accumulate
186 0 : delta[0]=paction->difference(iarg,average[0],cv[0]);
187 0 : average[0]+=decay*delta[0];
188 0 : average[0]=paction->bringBackInPbc(iarg,average[0]); // equation 8 of "Metadynamics with adaptive Gaussians"
189 : }
190 : // THE VARIANCE
191 0 : if(variance.size()==0) {
192 0 : variance.resize(1,0.); // nonredundant members dimension=ncv*(ncv+1)/2;
193 : } else {
194 0 : variance[0]+=decay*(delta[0]*delta[0]-variance[0]);
195 : }
196 : break;
197 0 : case geometry:
198 : //this calculates in variance the \nabla CV_i \dot \nabla CV_j
199 0 : variance.resize(1);
200 : // now the signal for retrieving the gradients should be already given by checkNeedsGradients.
201 : // here just do the projections
202 : // note that the call checkNeedsGradients() in BiasMetaD takes care of switching on the call to gradients
203 0 : if (nowAddAHill) { // geometry is sync with hill deposition
204 : // eq 12 of "Metadynamics with adaptive Gaussians"
205 0 : variance[0]=sigma*sigma*(paction->getProjection(iarg,iarg));
206 : }
207 : break;
208 0 : default:
209 0 : plumed_merror("This flexible bin is not recognized");
210 : }
211 0 : }
212 :
213 : ///
214 : /// Calculate the matrix of (dcv_i/dx)*(dcv_j/dx)^-1
215 : /// that is needed for the metrics in metadynamics
216 : ///
217 : ///
218 374 : vector<double> FlexibleBin::getInverseMatrix() const {
219 374 : unsigned ncv=paction->getNumberOfArguments();
220 : Matrix<double> matrix(ncv,ncv);
221 : unsigned i,j,k;
222 : k=0;
223 : //paction->log<<"------------ GET INVERSE MATRIX ---------------\n";
224 : // place the matrix in a complete matrix for compatibility
225 1288 : for (i=0; i<ncv; i++) {
226 1537 : for (j=i; j<ncv; j++) {
227 1620 : matrix(j,i)=matrix(i,j)=variance[k];
228 540 : k++;
229 : }
230 : }
231 : #define NEWFLEX
232 : #ifdef NEWFLEX
233 : // diagonalize to impose boundaries (only if boundaries are set)
234 : Matrix<double> eigenvecs(ncv,ncv);
235 374 : vector<double> eigenvals(ncv);
236 :
237 : //eigenvecs: first is eigenvec number, second is eigenvec component
238 374 : if(diagMat( matrix, eigenvals, eigenvecs )!=0) {plumed_merror("diagonalization in FlexibleBin failed! This matrix is weird\n");};
239 :
240 1288 : for (i=0; i<ncv; i++) { //loop on the dimension
241 914 : if( limitmax[i] ) {
242 : //limit every component that is larger
243 0 : for (j=0; j<ncv; j++) { //loop on components
244 0 : if(pow(eigenvals[j]*eigenvecs[j][i],2)>pow(sigmamax[i],2) ) {
245 0 : eigenvals[j]=sqrt(pow(sigmamax[i]/(eigenvecs[j][i]),2))*copysign(1.,eigenvals[j]);
246 : }
247 : }
248 : }
249 : }
250 1288 : for (i=0; i<ncv; i++) { //loop on the dimension
251 : // find the largest one: if it is smaller than min then rescale
252 914 : if( limitmin[i] ) {
253 : unsigned imax=0;
254 : double fmax=-1.e10;
255 0 : for (j=0; j<ncv; j++) { //loop on components
256 0 : double fact=pow(eigenvals[j]*eigenvecs[j][i],2);
257 0 : if(fact>fmax) {
258 : fmax=fact; imax=j;
259 : }
260 : }
261 0 : if(fmax<pow(sigmamin[i],2) ) {
262 0 : eigenvals[imax]=sqrt(pow(sigmamin[i]/(eigenvecs[imax][i]),2))*copysign(1.,eigenvals[imax]);
263 : }
264 : }
265 : }
266 :
267 : // normalize eigenvecs
268 : Matrix<double> newinvmatrix(ncv,ncv);
269 1288 : for (i=0; i<ncv; i++) {
270 1703 : for (j=0; j<ncv; j++) {
271 1246 : newinvmatrix[j][i]=eigenvecs[j][i]/eigenvals[j];
272 : }
273 : }
274 :
275 374 : vector<double> uppervec(ncv*(ncv+1)/2);
276 : k=0;
277 1288 : for (i=0; i<ncv; i++) {
278 1537 : for (j=i; j<ncv; j++) {
279 : double scal=0;
280 2118 : for(unsigned l=0; l<ncv; ++l) {
281 789 : scal+=eigenvecs[l][i]*newinvmatrix[l][j];
282 : }
283 1080 : uppervec[k]=scal; k++;
284 : }
285 : }
286 : #else
287 : // get the inverted matrix
288 : Matrix<double> invmatrix(ncv,ncv);
289 : Invert(matrix,invmatrix);
290 : vector<double> uppervec(ncv*(ncv+1)/2);
291 : // upper diagonal of the inverted matrix (that is symmetric)
292 : k=0;
293 : for (i=0; i<ncv; i++) {
294 : for (j=i; j<ncv; j++) {
295 : uppervec[k]=invmatrix(i,j);
296 : //paction->log<<"VV "<<i<<" "<<j<<" "<<uppervec[k]<<"\n";
297 : k++;
298 : }
299 : }
300 : #endif
301 374 : return uppervec;
302 : }
303 :
304 : ///
305 : /// Calculate the matrix of (dcv_i/dx)*(dcv_j/dx)^-1
306 : /// that is needed for the metrics in metadynamics
307 : /// for PBMetaD like FlexBin
308 : ///
309 0 : vector<double> FlexibleBin::getInverseMatrix(unsigned iarg) const {
310 : // diagonalize to impose boundaries (only if boundaries are set)
311 0 : vector<double> eigenvals(1, variance[0]);
312 0 : if( limitmax[0] ) {
313 0 : if(eigenvals[0]>sigmamax[0]) {
314 0 : eigenvals[0]=sigmamax[0];
315 : }
316 : }
317 : // find the largest one: if it is smaller than min then rescale
318 0 : if( limitmin[0] ) {
319 : double fmax=-1.e10;
320 0 : double fact=eigenvals[0];
321 0 : if(fact>fmax) {
322 : fmax=fact;
323 : }
324 0 : if(fmax<sigmamin[0]) {
325 0 : eigenvals[0]=sigmamin[0];
326 : }
327 : }
328 0 : vector<double> uppervec(1,1./eigenvals[0]);
329 :
330 0 : return uppervec;
331 : }
332 :
333 4839 : }
|
