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