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 ); ActionPilot::registerKeywords( keys );
77 8 : PathProjectionCalculator::registerKeywords(keys);
78 16 : keys.add("compulsory","STRIDE","1","the frequency with which to reparameterize the path");
79 16 : keys.add("compulsory","FIXED","0","the frames in the path to fix");
80 16 : keys.add("compulsory","MAXCYLES","100","number of cycles of the algorithm to run");
81 16 : keys.add("compulsory","TOL","1E-4","the tolerance to use for the path reparameterization algorithm");
82 16 : keys.add("optional","DISPLACE_FRAMES","label of an action that tells us how to displace the frames. These displacements are applied before "
83 : "running the reparameterization algorith");
84 8 : }
85 :
86 4 : PathReparameterization::PathReparameterization(const ActionOptions&ao):
87 : Action(ao),
88 : ActionPilot(ao),
89 4 : displace_value(NULL),
90 4 : path_projector(this)
91 : {
92 8 : parse("MAXCYLES",maxcycles); parse("TOL",TOL);
93 4 : log.printf(" running till change is less than %f or until there have been %d optimization cycles \n", TOL, maxcycles);
94 4 : len.resize( path_projector.getNumberOfFrames() ); sumlen.resize( path_projector.getNumberOfFrames() ); sfrac.resize( path_projector.getNumberOfFrames() );
95 8 : std::vector<unsigned> fixed; parseVector("FIXED",fixed);
96 4 : if( fixed.size()==1 ) {
97 0 : if( fixed[0]!=0 ) error("input to FIXED should be two integers");
98 0 : ifix1=0; ifix2=path_projector.getNumberOfFrames()-1;
99 4 : } else if( fixed.size()==2 ) {
100 4 : if( fixed[0]<1 || fixed[1]<1 || fixed[0]>path_projector.getNumberOfFrames() || fixed[1]>path_projector.getNumberOfFrames() ) {
101 0 : error("input to FIXED should be two numbers between 1 and the number of frames");
102 : }
103 4 : ifix1=fixed[0]-1; ifix2=fixed[1]-1;
104 0 : } else error("input to FIXED should be two integers");
105 4 : log.printf(" fixing frames %d and %d when reparameterizing \n", ifix1, ifix2 );
106 8 : std::string dframe; parse("DISPLACE_FRAMES",dframe);
107 4 : if( dframe.length()>0 ) {
108 2 : ActionWithValue* av = plumed.getActionSet().selectWithLabel<ActionWithValue*>( dframe );
109 2 : if( !av ) error("could not find action with label " + dframe + " specified to DISPLACE_FRAMES keyword in input file");
110 2 : if( av->getName()!="AVERAGE_PATH_DISPLACEMENT" ) error("displace object is not of correct type");
111 2 : displace_value = av->copyOutput(0);
112 : }
113 4 : }
114 :
115 11542 : bool PathReparameterization::loopEnd( const int& index, const int& end, const int& inc ) const {
116 11542 : if( inc>0 && index<end ) return false;
117 428 : else if( inc<0 && index>end ) return false;
118 : return true;
119 : }
120 :
121 4459 : double PathReparameterization::computeSpacing( const unsigned& ifrom, const unsigned& ito ) {
122 4459 : path_projector.getDisplaceVector( ifrom, ito, data );
123 170534 : double length=0; for(unsigned i=0; i<data.size(); ++i) length += data[i]*data[i];
124 4459 : return sqrt( length );
125 : }
126 :
127 139 : void PathReparameterization::calcCurrentPathSpacings( const int& istart, const int& iend ) {
128 : plumed_dbg_assert( istart<len.size() && iend<len.size() );
129 139 : len[istart] = sumlen[istart]=0;
130 : //printf("HELLO PATH SPACINGS ARE CURRENTLY \n");
131 :
132 : // Get the spacings given we can go forward and backwards
133 139 : int incr=1; if( istart>iend ) { incr=-1; }
134 :
135 4598 : for(int i=istart+incr; loopEnd(i,iend+incr,incr)==false; i+=incr) {
136 4459 : len[i] = computeSpacing( i-incr, i ); sumlen[i] = sumlen[i-incr] + len[i];
137 : //printf("FRAME %d TO FRAME %d EQUALS %f : %f \n",i-incr,i,len[i],sumlen[i] );
138 : }
139 139 : }
140 :
141 33 : void PathReparameterization::reparameterizePart( const int& istart, const int& iend, const double& target ) {
142 33 : calcCurrentPathSpacings( istart, iend ); unsigned cfin;
143 : // If a target separation is set we fix where we want the nodes
144 33 : int incr=1; if( istart>iend ) { incr=-1; }
145 :
146 33 : if( target>0 ) {
147 6 : if( iend>istart ) {
148 19 : for(unsigned i=istart; i<iend+1; ++i) sfrac[i] = target*(i-istart);
149 : } else {
150 14 : for(int i=istart-1; i>iend-1; --i) sfrac[i]=target*(istart-i);
151 : }
152 6 : cfin = iend+incr;
153 : } else {
154 27 : cfin = iend;
155 : }
156 :
157 33 : double prevsum=0.; Matrix<double> newmatrix( path_projector.getNumberOfFrames(), data.size() );
158 139 : for(unsigned iter=0; iter<maxcycles; ++iter) {
159 139 : if( fabs(sumlen[iend] - prevsum)<=TOL ) break ;
160 : prevsum = sumlen[iend];
161 : // If no target is set we redistribute length
162 106 : if( target<0 ) {
163 94 : plumed_assert( istart<iend );
164 94 : double dr = sumlen[iend] / static_cast<double>( iend - istart );
165 3506 : for(unsigned i=istart; i<iend; ++i) sfrac[i] = dr*(i-istart);
166 : }
167 :
168 : // Now compute positions of new nodes in path
169 3472 : for(int i=istart+incr; loopEnd(i,cfin,incr)==false; i+=incr) {
170 3366 : int k = istart;
171 69429 : while( !((sumlen[k] < sfrac[i]) && (sumlen[k+incr]>=sfrac[i])) ) {
172 66069 : k+=incr;
173 66069 : if( cfin==iend && k>= iend+1 ) plumed_merror("path reparameterization error");
174 66069 : else if( cfin==(iend+1) && k>=iend ) { k=iend-1; break; }
175 66067 : else if( cfin==(iend-1) && k<=iend ) { k=iend+1; break; }
176 : }
177 3366 : double dr = (sfrac[i]-sumlen[k])/len[k+incr];
178 : // Copy the reference configuration to the row of a matrix
179 3366 : path_projector.getReferenceConfiguration( k, data );
180 129024 : for(unsigned j=0; j<data.size(); ++j) newmatrix(i,j) = data[j];
181 3366 : path_projector.getDisplaceVector( k, k+incr, data );
182 : // Shift the reference configuration by this ammount
183 129024 : for(unsigned j=0; j<data.size(); ++j) newmatrix(i,j) += dr*data[j];
184 : }
185 :
186 : // Copy the positions of the new path to the new paths
187 3472 : for(int i=istart+incr; loopEnd(i,cfin,incr)==false; i+=incr) {
188 129024 : for(unsigned j=0; j<data.size(); ++j) data[j] = newmatrix(i,j);
189 3366 : path_projector.setReferenceConfiguration( i, data );
190 : }
191 :
192 : // Recompute the separations between frames
193 106 : calcCurrentPathSpacings( istart, iend );
194 : }
195 33 : }
196 :
197 29 : void PathReparameterization::update() {
198 : // We never run this on the first step
199 29 : if( getStep()==0 ) return ;
200 :
201 : // Shift the frames using the displacements
202 27 : if( displace_value ) {
203 1031 : for(unsigned i=0; i<path_projector.getNumberOfFrames(); ++i) {
204 1006 : if( i==ifix1 || i==ifix2 ) continue ;
205 : // Retrieve the current position of the frame
206 956 : path_projector.getReferenceConfiguration( i, data );
207 : // Shift using the averages accumulated in the action that accumulates the displacements
208 956 : unsigned kstart = i*data.size();
209 36908 : for(unsigned j=0; j<data.size(); ++j) data[j] += displace_value->get( kstart + j );
210 : // And now set the new position of the refernce frame
211 956 : path_projector.setReferenceConfiguration( i, data );
212 : }
213 : }
214 :
215 : // First reparameterize the part between the fixed frames
216 27 : reparameterizePart( ifix1, ifix2, -1.0 );
217 :
218 : // Get the separation between frames which we will use to set the remaining frames
219 27 : double target = sumlen[ifix2] / ( ifix2 - ifix1 );
220 :
221 : // And reparameterize the begining and end of the path
222 27 : if( ifix1>0 ) reparameterizePart( ifix1, 0, target );
223 27 : if( ifix2<(path_projector.getNumberOfFrames()-1) ) reparameterizePart( ifix2, path_projector.getNumberOfFrames()-1, target );
224 : // And update any RMSD objects that depend on input values
225 27 : path_projector.updateDepedentRMSDObjects();
226 : }
227 :
228 : }
229 : }
|
