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 : This action is used by [ADAPTIVE_PATH](ADAPTIVE_PATH.md) and [pathtools](pathtools.md). The algorithm in this action takes
36 : a set of trajectory frames along a [PATH](PATH.md) in input and adjusts them until they are all are equally spaced. To call
37 : this algorithm from [pathtools](pathtools.md) you would use an input like this one:
38 :
39 : ```plumed
40 : plumed pathtools --path in_path.pdb --metric EUCLIDEAN --out final_path.pdb
41 : ```
42 :
43 : If you are using this action directly and not through [pathtools](pathtools.md) or [ADAPTIVE_PATH](ADAPTIVE_PATH.md) you will
44 : probably be using it in an input something like this if your path is defined using atomic coordinates:
45 :
46 : ```plumed
47 : #SETTINGS INPUTFILES=regtest/trajectories/path_msd/all.pdb
48 : rmsd: RMSD REFERENCE=regtest/trajectories/path_msd/all.pdb DISPLACEMENT TYPE=OPTIMAL
49 : # Accumulate the average displacement between the reference path and the trajectories that have sampled the transition
50 : disp: AVERAGE_PATH_DISPLACEMENT ...
51 : ARG=rmsd.disp STRIDE=1
52 : METRIC={RMSD DISPLACEMENT TYPE=OPTIMAL ALIGN=1,1,1,1,1,1,1,1,1,1,1,1,1 DISPLACE=1,1,1,1,1,1,1,1,1,1,1,1,1}
53 : METRIC_COMPONENT=disp REFERENCE=rmsd_ref
54 : ...
55 : # Now displace the original path by the accumulated displacement and reparameterize so that all frames are equally spaced
56 : REPARAMETERIZE_PATH ...
57 : DISPLACE_FRAMES=disp FIXED=1,42
58 : METRIC={RMSD DISPLACEMENT TYPE=OPTIMAL ALIGN=1,1,1,1,1,1,1,1,1,1,1,1,1 DISPLACE=1,1,1,1,1,1,1,1,1,1,1,1,1}
59 : METRIC_COMPONENT=disp REFERENCE=rmsd_ref
60 : MAXCYLES=100 TOL=1E-4 STRIDE=0
61 : ...
62 : # And output the final reparameterized path at the end of the simulation
63 : DUMPPDB DESCRIPTION=PATH STRIDE=0 FILE=outpatb.pdb ATOMS=rmsd_ref ATOM_INDICES=1-13
64 : ```
65 :
66 : The method for reparameterizing paths that is implemented in this input is discussed in the example documentation for
67 : [AVERAGE_PATH_DISPLACEMENT](AVERAGE_PATH_DISPLACEMENT.md).
68 :
69 : */
70 : //+ENDPLUMEDOC
71 :
72 : namespace PLMD {
73 : namespace mapping {
74 :
75 : class PathReparameterization : public ActionPilot {
76 : private:
77 : /// Number of cycles of the optimization algorithm to run
78 : unsigned maxcycles;
79 : /// The points on the path to fix
80 : unsigned ifix1, ifix2;
81 : /// Tolerance for the minimization algorithm
82 : double TOL;
83 : /// Value containing amount to displace each reference configuration by
84 : Value* displace_value;
85 : /// The action for calculating the distances between the frames
86 : PathProjectionCalculator path_projector;
87 : /// Used to store current spacing between frames in path
88 : std::vector<double> data, len, sumlen, sfrac;
89 : ///
90 : bool loopEnd( int index, int end, int inc ) const;
91 : ///
92 : double computeSpacing( unsigned ifrom, unsigned ito );
93 : ///
94 : void calcCurrentPathSpacings( int istart, int iend );
95 : ///
96 : void reparameterizePart( int istart, int iend, double target );
97 : public:
98 : static void registerKeywords( Keywords& keys );
99 : PathReparameterization(const ActionOptions&);
100 110 : void calculate() {}
101 110 : void apply() {}
102 : void update();
103 : };
104 :
105 : PLUMED_REGISTER_ACTION(PathReparameterization,"REPARAMETERIZE_PATH")
106 :
107 14 : void PathReparameterization::registerKeywords( Keywords& keys ) {
108 14 : Action::registerKeywords( keys );
109 14 : ActionPilot::registerKeywords( keys );
110 14 : PathProjectionCalculator::registerKeywords(keys);
111 28 : keys.remove("ARG");
112 14 : keys.add("compulsory","STRIDE","1","the frequency with which to reparameterize the path");
113 14 : keys.add("compulsory","FIXED","0","the frames in the path to fix");
114 14 : keys.add("compulsory","MAXCYLES","100","number of cycles of the algorithm to run");
115 14 : keys.add("compulsory","TOL","1E-4","the tolerance to use for the path reparameterization algorithm");
116 14 : keys.add("optional","DISPLACE_FRAMES","label of an action that tells us how to displace the frames. These displacements are applied before "
117 : "running the reparameterization algorith");
118 14 : }
119 :
120 6 : PathReparameterization::PathReparameterization(const ActionOptions&ao):
121 : Action(ao),
122 : ActionPilot(ao),
123 6 : displace_value(NULL),
124 6 : path_projector(this) {
125 6 : parse("MAXCYLES",maxcycles);
126 6 : parse("TOL",TOL);
127 6 : log.printf(" running till change is less than %f or until there have been %d optimization cycles \n", TOL, maxcycles);
128 6 : len.resize( path_projector.getNumberOfFrames() );
129 6 : sumlen.resize( path_projector.getNumberOfFrames() );
130 6 : sfrac.resize( path_projector.getNumberOfFrames() );
131 6 : std::vector<unsigned> fixed;
132 12 : parseVector("FIXED",fixed);
133 6 : if( fixed.size()==1 ) {
134 0 : if( fixed[0]!=0 ) {
135 0 : error("input to FIXED should be two integers");
136 : }
137 0 : ifix1=0;
138 0 : ifix2=path_projector.getNumberOfFrames()-1;
139 6 : } else if( fixed.size()==2 ) {
140 6 : if( fixed[0]<1 || fixed[1]<1 || fixed[0]>path_projector.getNumberOfFrames() || fixed[1]>path_projector.getNumberOfFrames() ) {
141 0 : error("input to FIXED should be two numbers between 1 and the number of frames");
142 : }
143 6 : ifix1=fixed[0]-1;
144 6 : ifix2=fixed[1]-1;
145 : } else {
146 0 : error("input to FIXED should be two integers");
147 : }
148 6 : log.printf(" fixing frames %d and %d when reparameterizing \n", ifix1, ifix2 );
149 : std::string dframe;
150 12 : parse("DISPLACE_FRAMES",dframe);
151 6 : if( dframe.length()>0 ) {
152 4 : ActionWithValue* av = plumed.getActionSet().selectWithLabel<ActionWithValue*>( dframe );
153 4 : if( !av ) {
154 0 : error("could not find action with label " + dframe + " specified to DISPLACE_FRAMES keyword in input file");
155 : }
156 4 : if( av->getName()!="AVERAGE_PATH_DISPLACEMENT" ) {
157 0 : error("displace object is not of correct type");
158 : }
159 4 : displace_value = av->copyOutput(0);
160 : }
161 6 : }
162 :
163 25374 : bool PathReparameterization::loopEnd( const int index, const int end, const int inc ) const {
164 25374 : if( inc>0 && index<end ) {
165 : return false;
166 2792 : } else if( inc<0 && index>end ) {
167 311 : return false;
168 : }
169 : return true;
170 : }
171 :
172 9139 : double PathReparameterization::computeSpacing( const unsigned ifrom, const unsigned ito ) {
173 9139 : path_projector.getDisplaceVector( ifrom, ito, data );
174 : double length=0;
175 357734 : for(unsigned i=0; i<data.size(); ++i) {
176 348595 : length += data[i]*data[i];
177 : }
178 9139 : return sqrt( length );
179 : }
180 :
181 975 : void PathReparameterization::calcCurrentPathSpacings( const int istart, const int iend ) {
182 : plumed_dbg_assert( static_cast<unsigned>(istart)<len.size() && static_cast<unsigned>(iend)<len.size() );
183 975 : len[istart] = sumlen[istart]=0;
184 :
185 : // Get the spacings given we can go forward and backwards
186 : int incr=1;
187 975 : if( istart>iend ) {
188 : incr=-1;
189 : }
190 :
191 10114 : for(int i=istart+incr; loopEnd(i,iend+incr,incr)==false; i+=incr) {
192 9139 : len[i] = computeSpacing( i-incr, i );
193 9139 : sumlen[i] = sumlen[i-incr] + len[i];
194 : }
195 975 : }
196 :
197 222 : void PathReparameterization::reparameterizePart( const int istart, const int iend, const double target ) {
198 222 : calcCurrentPathSpacings( istart, iend );
199 : int cfin;
200 : // If a target separation is set we fix where we want the nodes
201 : int incr=1;
202 222 : if( istart>iend ) {
203 : incr=-1;
204 : }
205 :
206 222 : if( target>0 ) {
207 114 : if( iend>istart ) {
208 181 : for(int i=istart; i<iend+1; ++i) {
209 124 : sfrac[i] = target*(i-istart);
210 : }
211 : } else {
212 122 : for(int i=istart-1; i>iend-1; --i) {
213 65 : sfrac[i]=target*(istart-i);
214 : }
215 : }
216 114 : cfin = iend+incr;
217 : } else {
218 : cfin = iend;
219 : }
220 :
221 : double prevsum=0.;
222 222 : Matrix<double> newmatrix( path_projector.getNumberOfFrames(), data.size() );
223 975 : for(unsigned iter=0; iter<maxcycles; ++iter) {
224 975 : if( fabs(sumlen[iend] - prevsum)<=TOL ) {
225 : break ;
226 : }
227 : prevsum = sumlen[iend];
228 : // If no target is set we redistribute length
229 753 : if( target<0 ) {
230 615 : plumed_assert( istart<iend );
231 615 : double dr = sumlen[iend] / static_cast<double>( iend - istart );
232 7933 : for(int i=istart; i<iend; ++i) {
233 7318 : sfrac[i] = dr*(i-istart);
234 : }
235 : }
236 :
237 : // Now compute positions of new nodes in path
238 7630 : for(int i=istart+incr; loopEnd(i,cfin,incr)==false; i+=incr) {
239 : int k = istart;
240 85696 : while( !((sumlen[k] < sfrac[i]) && (sumlen[k+incr]>=sfrac[i])) ) {
241 78949 : k+=incr;
242 78949 : if( cfin==iend && k>= iend+1 ) {
243 0 : plumed_merror("path reparameterization error");
244 78949 : } else if( cfin==(iend+1) && k>=iend ) {
245 67 : k=iend-1;
246 67 : break;
247 78882 : } else if( cfin==(iend-1) && k<=iend ) {
248 : k=iend+1;
249 : break;
250 : }
251 : }
252 6877 : double dr = (sfrac[i]-sumlen[k])/len[k+incr];
253 : // Copy the reference configuration to the row of a matrix
254 6877 : path_projector.getReferenceConfiguration( k, data );
255 269464 : for(unsigned j=0; j<data.size(); ++j) {
256 262587 : newmatrix(i,j) = data[j];
257 : }
258 6877 : path_projector.getDisplaceVector( k, k+incr, data );
259 : // Shift the reference configuration by this amount
260 269464 : for(unsigned j=0; j<data.size(); ++j) {
261 262587 : newmatrix(i,j) += dr*data[j];
262 : }
263 : }
264 :
265 : // Copy the positions of the new path to the new paths
266 7630 : for(int i=istart+incr; loopEnd(i,cfin,incr)==false; i+=incr) {
267 269464 : for(unsigned j=0; j<data.size(); ++j) {
268 262587 : data[j] = newmatrix(i,j);
269 : }
270 6877 : path_projector.setReferenceConfiguration( i, data );
271 : }
272 :
273 : // Recompute the separations between frames
274 753 : calcCurrentPathSpacings( istart, iend );
275 : }
276 222 : }
277 :
278 112 : void PathReparameterization::update() {
279 : // We never run this on the first step
280 112 : if( getStep()==0 ) {
281 : return ;
282 : }
283 :
284 : // Shift the frames using the displacements
285 108 : if( displace_value ) {
286 1841 : for(unsigned i=0; i<path_projector.getNumberOfFrames(); ++i) {
287 1735 : if( i==ifix1 || i==ifix2 ) {
288 212 : continue ;
289 : }
290 : // Retrieve the current position of the frame
291 1523 : path_projector.getReferenceConfiguration( i, data );
292 : // Shift using the averages accumulated in the action that accumulates the displacements
293 1523 : unsigned kstart = i*data.size();
294 59588 : for(unsigned j=0; j<data.size(); ++j) {
295 58065 : data[j] += displace_value->get( kstart + j );
296 : }
297 : // And now set the new position of the refernce frame
298 1523 : path_projector.setReferenceConfiguration( i, data );
299 : }
300 : }
301 :
302 : // First reparameterize the part between the fixed frames
303 108 : reparameterizePart( ifix1, ifix2, -1.0 );
304 :
305 : // Get the separation between frames which we will use to set the remaining frames
306 108 : double target = sumlen[ifix2] / ( ifix2 - ifix1 );
307 :
308 : // And reparameterize the begining and end of the path
309 108 : if( ifix1>0 ) {
310 57 : reparameterizePart( ifix1, 0, target );
311 : }
312 108 : if( ifix2<(path_projector.getNumberOfFrames()-1) ) {
313 57 : reparameterizePart( ifix2, path_projector.getNumberOfFrames()-1, target );
314 : }
315 : // And update any RMSD objects that depend on input values
316 108 : path_projector.updateDepedentRMSDObjects();
317 : }
318 :
319 : }
320 : }
|
