Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2012-2018 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 21 : 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 21 : if(type==diffusion) {
36 2 : unsigned ncv=paction->getNumberOfArguments();
37 2 : vector<double> average(ncv*(ncv+1)/2);
38 2 : vector<double> variance(ncv*(ncv+1)/2);
39 : }
40 21 : paction->log<<" Limits for sigmas using adaptive hills: \n";
41 61 : for(unsigned i=0; i<paction->getNumberOfArguments(); ++i) {
42 40 : paction->log<<" CV "<<paction->getPntrToArgument(i)->getName()<<":\n";
43 40 : 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 40 : limitmin.push_back(false);
49 40 : paction->log<<" Min No ";
50 : }
51 40 : 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 40 : limitmax.push_back(false);
57 40 : paction->log<<" Max No ";
58 : }
59 40 : paction->log<<" \n";
60 : }
61 :
62 21 : }
63 : /// Update the flexible bin
64 : /// in case of diffusion based: update at every step
65 : /// in case of gradient based: update only when you add the hill
66 778 : void FlexibleBin::update(bool nowAddAHill) {
67 778 : unsigned ncv=paction->getNumberOfArguments();
68 778 : unsigned dimension=ncv*(ncv+1)/2;
69 : // this is done all the times from scratch. It is not an accumulator
70 778 : unsigned k=0;
71 : unsigned i;
72 778 : vector<double> cv;
73 1556 : vector<double> delta;
74 : // if you use this below then the decay is in time units
75 : //double decay=paction->getTimeStep()/sigma;
76 : // to be consistent with the rest of the program: everything is better to be in timesteps
77 778 : double decay=1./sigma;
78 : // here update the flexible bin according to the needs
79 778 : switch (type) {
80 : // This should be called every time
81 : case diffusion:
82 : //
83 : // THE AVERAGE VALUE
84 : //
85 : // beware: the pbc
86 586 : delta.resize(ncv);
87 586 : for(i=0; i<ncv; i++)cv.push_back(paction->getArgument(i));
88 586 : if(average.size()==0) { // initial time: just set the initial vector
89 2 : average.resize(ncv);
90 2 : for(i=0; i<ncv; i++)average[i]=cv[i];
91 : } else { // accumulate
92 1168 : for(i=0; i<ncv; i++) {
93 584 : delta[i]=paction->difference(i,average[i],cv[i]);
94 584 : average[i]+=decay*delta[i];
95 584 : average[i]=paction->bringBackInPbc(i,average[i]); // equation 8 of "Metadynamics with adaptive Gaussians"
96 : }
97 :
98 : }
99 : //
100 : // THE VARIANCE
101 : //
102 586 : if(variance.size()==0) {
103 2 : variance.resize(dimension,0.); // nonredundant members dimension=ncv*(ncv+1)/2;
104 : } else {
105 584 : k=0;
106 1168 : for(i=0; i<ncv; i++) {
107 1168 : for(unsigned j=i; j<ncv; j++) { // upper diagonal loop
108 584 : variance[k]+=decay*(delta[i]*delta[j]-variance[k]);
109 584 : k++;
110 : }
111 : }
112 : }
113 586 : break;
114 : case geometry:
115 : //
116 : //this calculates in variance the \nabla CV_i \dot \nabla CV_j
117 : //
118 192 : variance.resize(dimension);
119 : //cerr<< "Doing geometry "<<endl;
120 : // now the signal for retrieving the gradients should be already given by checkNeedsGradients.
121 : // here just do the projections
122 : // note that the call checkNeedsGradients() in BiasMetaD takes care of switching on the call to gradients
123 192 : if (nowAddAHill) { // geometry is sync with hill deposition
124 : //cerr<< "add a hill "<<endl;
125 83 : k=0;
126 249 : for(unsigned i=0; i<ncv; i++) {
127 415 : for(unsigned j=i; j<ncv; j++) {
128 : // eq 12 of "Metadynamics with adaptive Gaussians"
129 249 : variance[k]=sigma*sigma*(paction->getProjection(i,j));
130 249 : k++;
131 : }
132 : };
133 : };
134 192 : break;
135 : default:
136 0 : cerr<< "This flexible bin is not recognized "<<endl;
137 0 : exit(1) ;
138 778 : }
139 :
140 778 : }
141 :
142 0 : vector<double> FlexibleBin::getMatrix() const {
143 0 : return variance;
144 : }
145 :
146 : ///
147 : /// Calculate the matrix of (dcv_i/dx)*(dcv_j/dx)^-1
148 : /// that is needed for the metrics in metadynamics
149 : ///
150 : ///
151 374 : vector<double> FlexibleBin::getInverseMatrix() const {
152 374 : unsigned ncv=paction->getNumberOfArguments();
153 374 : Matrix<double> matrix(ncv,ncv);
154 : unsigned i,j,k;
155 374 : k=0;
156 : //paction->log<<"------------ GET INVERSE MATRIX ---------------\n";
157 : // place the matrix in a complete matrix for compatibility
158 831 : for (i=0; i<ncv; i++) {
159 997 : for (j=i; j<ncv; j++) {
160 540 : matrix(j,i)=matrix(i,j)=variance[k];
161 540 : k++;
162 : }
163 : }
164 : // paction->log<<"MATRIX\n";
165 : // matrixOut(paction->log,matrix);
166 : #define NEWFLEX
167 : #ifdef NEWFLEX
168 : // diagonalize to impose boundaries (only if boundaries are set)
169 748 : Matrix<double> eigenvecs(ncv,ncv);
170 748 : vector<double> eigenvals(ncv);
171 :
172 : //eigenvecs: first is eigenvec number, second is eigenvec component
173 374 : if(diagMat( matrix, eigenvals, eigenvecs )!=0) {plumed_merror("diagonalization in FlexibleBin failed! This matrix is weird\n");};
174 :
175 : // paction->log<<"EIGENVECS \n";
176 : // matrixOut(paction->log,eigenvecs);
177 :
178 : //for (i=0;i<ncv;i++){//loop on the dimension
179 : // for (j=0;j<ncv;j++){//loop on the dimension
180 : // eigenvecs[i][j]=0.;
181 : // if(i==j)eigenvecs[i][j]=1.;
182 : // }
183 : //}
184 :
185 831 : for (i=0; i<ncv; i++) { //loop on the dimension
186 457 : if( limitmax[i] ) {
187 : //limit every component that is larger
188 0 : for (j=0; j<ncv; j++) { //loop on components
189 0 : if(pow(eigenvals[j]*eigenvecs[j][i],2)>pow(sigmamax[i],2) ) {
190 0 : eigenvals[j]=sqrt(pow(sigmamax[i]/(eigenvecs[j][i]),2))*copysign(1.,eigenvals[j]);
191 : }
192 : }
193 : }
194 : }
195 831 : for (i=0; i<ncv; i++) { //loop on the dimension
196 : // find the largest one: if it is smaller than min then rescale
197 457 : if( limitmin[i] ) {
198 0 : unsigned imax=0;
199 0 : double fmax=-1.e10;
200 0 : for (j=0; j<ncv; j++) { //loop on components
201 0 : double fact=pow(eigenvals[j]*eigenvecs[j][i],2);
202 0 : if(fact>fmax) {
203 0 : fmax=fact; imax=j;
204 : }
205 : }
206 0 : if(fmax<pow(sigmamin[i],2) ) {
207 0 : eigenvals[imax]=sqrt(pow(sigmamin[i]/(eigenvecs[imax][i]),2))*copysign(1.,eigenvals[imax]);
208 : }
209 :
210 : }
211 :
212 : }
213 :
214 : // normalize eigenvecs
215 748 : Matrix<double> newinvmatrix(ncv,ncv);
216 831 : for (i=0; i<ncv; i++) {
217 1080 : for (j=0; j<ncv; j++) {
218 623 : newinvmatrix[j][i]=eigenvecs[j][i]/eigenvals[j];
219 : }
220 : }
221 :
222 374 : vector<double> uppervec(ncv*(ncv+1)/2);
223 374 : k=0;
224 831 : for (i=0; i<ncv; i++) {
225 997 : for (j=i; j<ncv; j++) {
226 540 : double scal=0;
227 1329 : for(unsigned l=0; l<ncv; ++l) {
228 789 : scal+=eigenvecs[l][i]*newinvmatrix[l][j];
229 : }
230 540 : uppervec[k]=scal; k++;
231 : }
232 : }
233 : #else
234 : // get the inverted matrix
235 : Matrix<double> invmatrix(ncv,ncv);
236 : Invert(matrix,invmatrix);
237 : vector<double> uppervec(ncv*(ncv+1)/2);
238 : // upper diagonal of the inverted matrix (that is symmetric)
239 : k=0;
240 : for (i=0; i<ncv; i++) {
241 : for (j=i; j<ncv; j++) {
242 : uppervec[k]=invmatrix(i,j);
243 : //paction->log<<"VV "<<i<<" "<<j<<" "<<uppervec[k]<<"\n";
244 : k++;
245 : }
246 : }
247 : //paction->log<<"------------ END GET INVERSE MATRIX ---------------\n";
248 : #endif
249 :
250 748 : return uppervec;
251 : }
252 :
253 2523 : }
|
