Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2016-2020 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/ActionRegister.h"
25 : #include "core/PlumedMain.h"
26 : #include "core/ActionSet.h"
27 : #include "tools/Pbc.h"
28 : #include "tools/Matrix.h"
29 : #include "PathProjectionCalculator.h"
30 :
31 : //+PLUMEDOC ANALYSIS REPARAMETERIZE_PATH
32 : /*
33 : Take an input path with frames that are not equally spaced and make the frames equally spaced
34 :
35 : \par Examples
36 :
37 : */
38 : //+ENDPLUMEDOC
39 :
40 : namespace PLMD {
41 : namespace mapping {
42 :
43 : class PathReparameterization : public ActionPilot {
44 : private:
45 : /// Number of cycles of the optimization algorithm to run
46 : unsigned maxcycles;
47 : /// The points on the path to fix
48 : unsigned ifix1, ifix2;
49 : /// Tolerance for the minimization algorithm
50 : double TOL;
51 : /// Value containing ammount to displace each reference configuration by
52 : Value* displace_value;
53 : /// The action for calculating the distances between the frames
54 : PathProjectionCalculator path_projector;
55 : /// Used to store current spacing between frames in path
56 : std::vector<double> data, len, sumlen, sfrac;
57 : ///
58 : bool loopEnd( const int& index, const int& end, const int& inc ) const ;
59 : ///
60 : double computeSpacing( const unsigned& ifrom, const unsigned& ito );
61 : ///
62 : void calcCurrentPathSpacings( const int& istart, const int& iend );
63 : ///
64 : void reparameterizePart( const int& istart, const int& iend, const double& target );
65 : public:
66 : static void registerKeywords( Keywords& keys );
67 : PathReparameterization(const ActionOptions&);
68 27 : void calculate() {}
69 27 : void apply() {}
70 : void update();
71 : };
72 :
73 : PLUMED_REGISTER_ACTION(PathReparameterization,"REPARAMETERIZE_PATH")
74 :
75 8 : void PathReparameterization::registerKeywords( Keywords& keys ) {
76 8 : Action::registerKeywords( keys );
77 8 : ActionPilot::registerKeywords( keys );
78 8 : PathProjectionCalculator::registerKeywords(keys);
79 8 : keys.add("compulsory","STRIDE","1","the frequency with which to reparameterize the path");
80 8 : keys.add("compulsory","FIXED","0","the frames in the path to fix");
81 8 : keys.add("compulsory","MAXCYLES","100","number of cycles of the algorithm to run");
82 8 : keys.add("compulsory","TOL","1E-4","the tolerance to use for the path reparameterization algorithm");
83 8 : keys.add("optional","DISPLACE_FRAMES","label of an action that tells us how to displace the frames. These displacements are applied before "
84 : "running the reparameterization algorith");
85 8 : }
86 :
87 4 : PathReparameterization::PathReparameterization(const ActionOptions&ao):
88 : Action(ao),
89 : ActionPilot(ao),
90 4 : displace_value(NULL),
91 4 : path_projector(this) {
92 4 : parse("MAXCYLES",maxcycles);
93 4 : parse("TOL",TOL);
94 4 : log.printf(" running till change is less than %f or until there have been %d optimization cycles \n", TOL, maxcycles);
95 4 : len.resize( path_projector.getNumberOfFrames() );
96 4 : sumlen.resize( path_projector.getNumberOfFrames() );
97 4 : sfrac.resize( path_projector.getNumberOfFrames() );
98 : std::vector<unsigned> fixed;
99 8 : parseVector("FIXED",fixed);
100 4 : if( fixed.size()==1 ) {
101 0 : if( fixed[0]!=0 ) {
102 0 : error("input to FIXED should be two integers");
103 : }
104 0 : ifix1=0;
105 0 : ifix2=path_projector.getNumberOfFrames()-1;
106 4 : } else if( fixed.size()==2 ) {
107 4 : if( fixed[0]<1 || fixed[1]<1 || fixed[0]>path_projector.getNumberOfFrames() || fixed[1]>path_projector.getNumberOfFrames() ) {
108 0 : error("input to FIXED should be two numbers between 1 and the number of frames");
109 : }
110 4 : ifix1=fixed[0]-1;
111 4 : ifix2=fixed[1]-1;
112 : } else {
113 0 : error("input to FIXED should be two integers");
114 : }
115 4 : log.printf(" fixing frames %d and %d when reparameterizing \n", ifix1, ifix2 );
116 : std::string dframe;
117 8 : parse("DISPLACE_FRAMES",dframe);
118 4 : if( dframe.length()>0 ) {
119 2 : ActionWithValue* av = plumed.getActionSet().selectWithLabel<ActionWithValue*>( dframe );
120 2 : if( !av ) {
121 0 : error("could not find action with label " + dframe + " specified to DISPLACE_FRAMES keyword in input file");
122 : }
123 2 : if( av->getName()!="AVERAGE_PATH_DISPLACEMENT" ) {
124 0 : error("displace object is not of correct type");
125 : }
126 2 : displace_value = av->copyOutput(0);
127 : }
128 4 : }
129 :
130 11542 : bool PathReparameterization::loopEnd( const int& index, const int& end, const int& inc ) const {
131 11542 : if( inc>0 && index<end ) {
132 : return false;
133 428 : } else if( inc<0 && index>end ) {
134 77 : return false;
135 : }
136 : return true;
137 : }
138 :
139 4459 : double PathReparameterization::computeSpacing( const unsigned& ifrom, const unsigned& ito ) {
140 4459 : path_projector.getDisplaceVector( ifrom, ito, data );
141 : double length=0;
142 170534 : for(unsigned i=0; i<data.size(); ++i) {
143 166075 : length += data[i]*data[i];
144 : }
145 4459 : return sqrt( length );
146 : }
147 :
148 139 : void PathReparameterization::calcCurrentPathSpacings( const int& istart, const int& iend ) {
149 : plumed_dbg_assert( istart<len.size() && iend<len.size() );
150 139 : len[istart] = sumlen[istart]=0;
151 : //printf("HELLO PATH SPACINGS ARE CURRENTLY \n");
152 :
153 : // Get the spacings given we can go forward and backwards
154 139 : int incr=1;
155 139 : if( istart>iend ) {
156 9 : incr=-1;
157 : }
158 :
159 4598 : for(int i=istart+incr; loopEnd(i,iend+incr,incr)==false; i+=incr) {
160 4459 : len[i] = computeSpacing( i-incr, i );
161 4459 : sumlen[i] = sumlen[i-incr] + len[i];
162 : //printf("FRAME %d TO FRAME %d EQUALS %f : %f \n",i-incr,i,len[i],sumlen[i] );
163 : }
164 139 : }
165 :
166 33 : void PathReparameterization::reparameterizePart( const int& istart, const int& iend, const double& target ) {
167 33 : calcCurrentPathSpacings( istart, iend );
168 : unsigned cfin;
169 : // If a target separation is set we fix where we want the nodes
170 33 : int incr=1;
171 33 : if( istart>iend ) {
172 3 : incr=-1;
173 : }
174 :
175 33 : if( target>0 ) {
176 6 : if( iend>istart ) {
177 19 : for(unsigned i=istart; i<iend+1; ++i) {
178 16 : sfrac[i] = target*(i-istart);
179 : }
180 : } else {
181 14 : for(int i=istart-1; i>iend-1; --i) {
182 11 : sfrac[i]=target*(istart-i);
183 : }
184 : }
185 6 : cfin = iend+incr;
186 : } else {
187 27 : cfin = iend;
188 : }
189 :
190 : double prevsum=0.;
191 33 : Matrix<double> newmatrix( path_projector.getNumberOfFrames(), data.size() );
192 139 : for(unsigned iter=0; iter<maxcycles; ++iter) {
193 139 : if( fabs(sumlen[iend] - prevsum)<=TOL ) {
194 : break ;
195 : }
196 : prevsum = sumlen[iend];
197 : // If no target is set we redistribute length
198 106 : if( target<0 ) {
199 94 : plumed_assert( istart<iend );
200 94 : double dr = sumlen[iend] / static_cast<double>( iend - istart );
201 3506 : for(unsigned i=istart; i<iend; ++i) {
202 3412 : sfrac[i] = dr*(i-istart);
203 : }
204 : }
205 :
206 : // Now compute positions of new nodes in path
207 3472 : for(int i=istart+incr; loopEnd(i,cfin,incr)==false; i+=incr) {
208 3366 : int k = istart;
209 69429 : while( !((sumlen[k] < sfrac[i]) && (sumlen[k+incr]>=sfrac[i])) ) {
210 66069 : k+=incr;
211 66069 : if( cfin==iend && k>= iend+1 ) {
212 0 : plumed_merror("path reparameterization error");
213 66069 : } else if( cfin==(iend+1) && k>=iend ) {
214 2 : k=iend-1;
215 2 : break;
216 66067 : } else if( cfin==(iend-1) && k<=iend ) {
217 : k=iend+1;
218 : break;
219 : }
220 : }
221 3366 : double dr = (sfrac[i]-sumlen[k])/len[k+incr];
222 : // Copy the reference configuration to the row of a matrix
223 3366 : path_projector.getReferenceConfiguration( k, data );
224 129024 : for(unsigned j=0; j<data.size(); ++j) {
225 125658 : newmatrix(i,j) = data[j];
226 : }
227 3366 : path_projector.getDisplaceVector( k, k+incr, data );
228 : // Shift the reference configuration by this ammount
229 129024 : for(unsigned j=0; j<data.size(); ++j) {
230 125658 : newmatrix(i,j) += dr*data[j];
231 : }
232 : }
233 :
234 : // Copy the positions of the new path to the new paths
235 3472 : for(int i=istart+incr; loopEnd(i,cfin,incr)==false; i+=incr) {
236 129024 : for(unsigned j=0; j<data.size(); ++j) {
237 125658 : data[j] = newmatrix(i,j);
238 : }
239 3366 : path_projector.setReferenceConfiguration( i, data );
240 : }
241 :
242 : // Recompute the separations between frames
243 106 : calcCurrentPathSpacings( istart, iend );
244 : }
245 33 : }
246 :
247 29 : void PathReparameterization::update() {
248 : // We never run this on the first step
249 29 : if( getStep()==0 ) {
250 2 : return ;
251 : }
252 :
253 : // Shift the frames using the displacements
254 27 : if( displace_value ) {
255 1031 : for(unsigned i=0; i<path_projector.getNumberOfFrames(); ++i) {
256 1006 : if( i==ifix1 || i==ifix2 ) {
257 50 : continue ;
258 : }
259 : // Retrieve the current position of the frame
260 956 : path_projector.getReferenceConfiguration( i, data );
261 : // Shift using the averages accumulated in the action that accumulates the displacements
262 956 : unsigned kstart = i*data.size();
263 36908 : for(unsigned j=0; j<data.size(); ++j) {
264 35952 : data[j] += displace_value->get( kstart + j );
265 : }
266 : // And now set the new position of the refernce frame
267 956 : path_projector.setReferenceConfiguration( i, data );
268 : }
269 : }
270 :
271 : // First reparameterize the part between the fixed frames
272 27 : reparameterizePart( ifix1, ifix2, -1.0 );
273 :
274 : // Get the separation between frames which we will use to set the remaining frames
275 27 : double target = sumlen[ifix2] / ( ifix2 - ifix1 );
276 :
277 : // And reparameterize the begining and end of the path
278 27 : if( ifix1>0 ) {
279 3 : reparameterizePart( ifix1, 0, target );
280 : }
281 27 : if( ifix2<(path_projector.getNumberOfFrames()-1) ) {
282 3 : reparameterizePart( ifix2, path_projector.getNumberOfFrames()-1, target );
283 : }
284 : // And update any RMSD objects that depend on input values
285 27 : path_projector.updateDepedentRMSDObjects();
286 : }
287 :
288 : }
289 : }
|
