Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2016-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 "core/ActionWithValue.h"
23 : #include "core/ActionPilot.h"
24 : #include "core/ActionWithArguments.h"
25 : #include "core/ActionRegister.h"
26 : #include "tools/Matrix.h"
27 : #include "PathProjectionCalculator.h"
28 :
29 : //+PLUMEDOC ANALYSIS AVERAGE_PATH_DISPLACEMENT
30 : /*
31 : Accumulate the distances between the reference frames in the paths and the configurations visited
32 :
33 : \par Examples
34 :
35 : */
36 : //+ENDPLUMEDOC
37 :
38 : namespace PLMD {
39 : namespace mapping {
40 :
41 : class PathDisplacements : public ActionWithValue, public ActionPilot, public ActionWithArguments {
42 : private:
43 : bool clearnextstep;
44 : unsigned clearstride;
45 : double fadefact;
46 : std::vector<double> wsum, displace_v;
47 : Matrix<double> displacements;
48 : PathProjectionCalculator path_projector;
49 : public:
50 : static void registerKeywords( Keywords& keys );
51 : explicit PathDisplacements(const ActionOptions&);
52 : unsigned getNumberOfDerivatives();
53 573 : void clearDerivatives( const bool& force=false ) {}
54 573 : void calculate() {}
55 573 : void apply() {}
56 : void update();
57 : };
58 :
59 : PLUMED_REGISTER_ACTION(PathDisplacements,"AVERAGE_PATH_DISPLACEMENT")
60 :
61 6 : void PathDisplacements::registerKeywords( Keywords& keys ) {
62 6 : Action::registerKeywords( keys );
63 6 : ActionWithValue::registerKeywords( keys );
64 6 : ActionPilot::registerKeywords( keys );
65 6 : ActionWithArguments::registerKeywords( keys );
66 6 : PathProjectionCalculator::registerKeywords( keys );
67 6 : keys.add("compulsory","STRIDE","1","the frequency with which the average displacements should be collected and added to the average displacements");
68 6 : keys.add("compulsory","HALFLIFE","-1","the number of MD steps after which a previously measured path distance weighs only 50 percent in the average. This option may increase convergence by allowing to forget the memory of a bad initial guess path. The default is to set this to infinity");
69 6 : keys.add("compulsory","CLEAR","0","the frequency with which to clear all the accumulated data. The default value "
70 : "of 0 implies that all the data will be used and that the grid will never be cleared");
71 12 : keys.setValueDescription("matrix","matrix containing the average displacement between the trajectory and each of the landmarks that makes up the path");
72 6 : }
73 :
74 2 : PathDisplacements::PathDisplacements(const ActionOptions& ao):
75 : Action(ao),
76 : ActionWithValue(ao),
77 : ActionPilot(ao),
78 : ActionWithArguments(ao),
79 2 : clearnextstep(false),
80 2 : path_projector(this) {
81 : // Read in clear instructions
82 2 : parse("CLEAR",clearstride);
83 2 : if( clearstride>0 ) {
84 2 : if( clearstride%getStride()!=0 ) {
85 0 : error("CLEAR parameter must be a multiple of STRIDE");
86 : }
87 2 : log.printf(" clearing average every %u steps \n",clearstride);
88 : }
89 : double halflife;
90 2 : parse("HALFLIFE",halflife);
91 2 : log.printf(" weight of contribution to frame halves every %f steps \n",halflife);
92 2 : if( halflife<0 ) {
93 2 : fadefact=1.0;
94 : } else {
95 0 : fadefact = exp( -0.693147180559945 / static_cast<double>(halflife) );
96 : }
97 : // Now create the weights vector and displacements matrix
98 2 : unsigned nrows = getPntrToArgument(0)->getShape()[0];
99 2 : unsigned ncols = getPntrToArgument(0)->getShape()[1];
100 2 : wsum.resize( nrows );
101 : displacements.resize( nrows, ncols );
102 64 : for(unsigned i=0; i<nrows; ++i) {
103 62 : wsum[i]=0;
104 1740 : for(unsigned j=0; j<ncols; ++j) {
105 1678 : displacements(i,j)=0;
106 : }
107 : }
108 : // Add bibliography
109 4 : log<<" Bibliography "<<plumed.cite("Diaz Leines and Ensing, Phys. Rev. Lett. 109, 020601 (2012)")<<"\n";
110 : // And create a value to hold the displacements
111 2 : std::vector<unsigned> shape(2);
112 2 : shape[0]=nrows;
113 2 : shape[1]=ncols;
114 2 : addValue( shape );
115 2 : setNotPeriodic();
116 2 : getPntrToComponent(0)->buildDataStore();
117 2 : getPntrToComponent(0)->reshapeMatrixStore( shape[1] );
118 2 : }
119 :
120 0 : unsigned PathDisplacements::getNumberOfDerivatives() {
121 0 : return 0;
122 : }
123 :
124 573 : void PathDisplacements::update() {
125 573 : unsigned nrows = getPntrToArgument(0)->getShape()[0];
126 573 : unsigned ncols = getPntrToArgument(0)->getShape()[1];
127 :
128 573 : if( clearnextstep ) {
129 : unsigned k=0;
130 1010 : for(unsigned i=0; i<nrows; ++i) {
131 38700 : for(unsigned j=0; j<ncols; ++j) {
132 37714 : displacements(i,j)=0;
133 37714 : getPntrToComponent(0)->set(k,0);
134 37714 : k++;
135 : }
136 : }
137 24 : clearnextstep=false;
138 : }
139 :
140 : unsigned k=0, iclose1=0, iclose2=0;
141 : double v1v1=0, v3v3=0;
142 22417 : for(unsigned i=0; i<nrows; ++i) {
143 : double dist = 0;
144 799020 : for(unsigned j=0; j<ncols; ++j) {
145 777176 : double tmp = getPntrToArgument(0)->get(k);
146 777176 : dist += tmp*tmp;
147 777176 : k++;
148 : }
149 21844 : if( i==0 ) {
150 : v1v1 = dist;
151 : iclose1 = 0;
152 21271 : } else if( dist<v1v1 ) {
153 : v3v3=v1v1;
154 : v1v1=dist;
155 : iclose2=iclose1;
156 : iclose1=i;
157 10991 : } else if( i==1 ) {
158 : v3v3=dist;
159 : iclose2=1;
160 10991 : } else if( dist<v3v3 ) {
161 : v3v3=dist;
162 : iclose2=i;
163 : }
164 : }
165 : // And find third closest point
166 573 : int isign = iclose1 - iclose2;
167 : if( isign>1 ) {
168 : isign=1;
169 : } else if( isign<-1 ) {
170 : isign=-1;
171 : }
172 573 : int iclose3 = iclose1 + isign;
173 573 : unsigned ifrom=iclose1, ito=iclose3;
174 573 : if( iclose3<0 || iclose3>=nrows ) {
175 0 : ifrom=iclose2;
176 0 : ito=iclose1;
177 : }
178 :
179 : // Calculate the dot product of v1 with v2
180 573 : path_projector.getDisplaceVector( ifrom, ito, displace_v );
181 : double v2v2=0, v1v2=0;
182 573 : unsigned kclose1 = iclose1*ncols;
183 19183 : for(unsigned i=0; i<displace_v.size(); ++i) {
184 18610 : v2v2 += displace_v[i]*displace_v[i];
185 18610 : v1v2 += displace_v[i]*getPntrToArgument(0)->get(kclose1+i);
186 : }
187 :
188 573 : double root = sqrt( v1v2*v1v2 - v2v2 * ( v1v1 - v3v3) );
189 573 : double dx = 0.5 * ( (root + v1v2) / v2v2 - 1.);
190 573 : double weight2 = -1.* dx;
191 573 : double weight1 = 1.0 + dx;
192 573 : if( weight1>1.0 ) {
193 : weight1=1.0;
194 : weight2=0.0;
195 573 : } else if( weight2>1.0 ) {
196 : weight1=0.0;
197 : weight2=1.0;
198 : }
199 :
200 : // Accumulate displacements for path
201 19183 : for(unsigned i=0; i<ncols; ++i) {
202 18610 : double displace = getPntrToArgument(0)->get(kclose1+i) - dx*displace_v[i];
203 18610 : displacements(iclose1,i) += weight1 * displace;
204 18610 : displacements(iclose2,i) += weight2 * displace;
205 : }
206 :
207 : // Update weight accumulators
208 573 : wsum[iclose1] *= fadefact;
209 573 : wsum[iclose2] *= fadefact;
210 573 : wsum[iclose1] += weight1;
211 573 : wsum[iclose2] += weight2;
212 :
213 : // Update numbers in values
214 573 : if( wsum[iclose1] > epsilon ) {
215 19183 : for(unsigned i=0; i<ncols; ++i) {
216 18610 : getPntrToComponent(0)->set( kclose1+i, displacements(iclose1,i) / wsum[iclose1] );
217 : }
218 : }
219 573 : if( wsum[iclose2] > epsilon ) {
220 573 : unsigned kclose2 = iclose2*ncols;
221 19183 : for(unsigned i=0; i<ncols; ++i) {
222 18610 : getPntrToComponent(0)->set( kclose2+i, displacements(iclose2,i) / wsum[iclose2] );
223 : }
224 : }
225 :
226 : // Clear if required
227 573 : if( (getStep()>0 && clearstride>0 && getStep()%clearstride==0) ) {
228 25 : clearnextstep=true;
229 : }
230 573 : }
231 :
232 : }
233 : }
|
