Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2016-2019 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 "PathReparameterization.h"
23 :
24 : namespace PLMD {
25 : namespace mapping {
26 :
27 4 : PathReparameterization::PathReparameterization( const Pbc& ipbc, const std::vector<Value*>& iargs, std::vector<ReferenceConfiguration*>& pp ):
28 12 : mydpack( 1, pp[0]->getNumberOfReferenceArguments() + 3*pp[0]->getNumberOfReferencePositions() + 9 ),
29 12 : mypack( pp[0]->getNumberOfReferenceArguments(), pp[0]->getNumberOfReferencePositions(), mydpack ),
30 : mypath(pp),
31 : pbc(ipbc),
32 : args(iargs),
33 12 : mydir(ReferenceConfigurationOptions("DIRECTION")),
34 : len(pp.size()),
35 : sumlen(pp.size()),
36 : sfrac(pp.size()),
37 40 : MAXCYCLES(100)
38 : {
39 12 : mydir.setNamesAndAtomNumbers( pp[0]->getAbsoluteIndexes(), pp[0]->getArgumentNames() );
40 8 : mydir.zeroDirection(); pp[0]->setupPCAStorage( mypack );
41 4 : }
42 :
43 1741 : bool PathReparameterization::loopEnd( const int& index, const int& end, const int& inc ) const {
44 1741 : if( inc>0 && index<end ) return false;
45 270 : else if( inc<0 && index>end ) return false;
46 193 : return true;
47 : }
48 :
49 71 : void PathReparameterization::calcCurrentPathSpacings( const int& istart, const int& iend ) {
50 : plumed_dbg_assert( istart<len.size() && iend<len.size() );
51 213 : len[istart] = sumlen[istart]=0;
52 : //printf("HELLO PATH SPACINGS ARE CURRENTLY \n");
53 :
54 : // Get the spacings given we can go forward and backwards
55 71 : int incr=1; if( istart>iend ) { incr=-1; }
56 :
57 657 : for(int i=istart+incr; loopEnd(i,iend+incr,incr)==false; i+=incr) {
58 2930 : len[i] = mypath[i-incr]->calc( mypath[i]->getReferencePositions(), pbc, args, mypath[i]->getReferenceArguments(), mypack, false );
59 2344 : sumlen[i] = sumlen[i-incr] + len[i];
60 : //printf("FRAME %d TO FRAME %d EQUALS %f : %f \n",i-incr,i,len[i],sumlen[i] );
61 : }
62 71 : }
63 :
64 10 : void PathReparameterization::reparameterizePart( const int& istart, const int& iend, const double& target, const double& TOL ) {
65 10 : calcCurrentPathSpacings( istart, iend ); unsigned cfin;
66 : // If a target separation is set we fix where we want the nodes
67 10 : int incr=1; if( istart>iend ) { incr=-1; }
68 :
69 10 : if( target>0 ) {
70 6 : if( iend>istart ) {
71 19 : for(unsigned i=istart; i<iend+1; ++i) sfrac[i] = target*(i-istart);
72 : } else {
73 14 : for(int i=istart-1; i>iend-1; --i) sfrac[i]=target*(istart-i);
74 : }
75 6 : cfin = iend+incr;
76 : } else {
77 4 : cfin = iend;
78 : }
79 :
80 10 : std::vector<Direction> newpath;
81 500 : for(unsigned i=0; i<mypath.size(); ++i) {
82 640 : newpath.push_back( Direction(ReferenceConfigurationOptions("DIRECTION")) );
83 480 : newpath[i].setNamesAndAtomNumbers( mypath[i]->getAbsoluteIndexes(), mypath[i]->getArgumentNames() );
84 : }
85 :
86 : double prevsum=0.;
87 132 : for(unsigned iter=0; iter<MAXCYCLES; ++iter) {
88 142 : if( fabs(sumlen[iend] - prevsum)<=TOL ) break ;
89 : prevsum = sumlen[iend];
90 : // If no target is set we redistribute length
91 61 : if( target<0 ) {
92 49 : plumed_assert( istart<iend );
93 49 : double dr = sumlen[iend] / static_cast<double>( iend - istart );
94 531 : for(unsigned i=istart; i<iend; ++i) sfrac[i] = dr*(i-istart);
95 : }
96 :
97 : // Now compute positions of new nodes in path
98 542 : for(int i=istart+incr; loopEnd(i,cfin,incr)==false; i+=incr) {
99 481 : int k = istart;
100 10268 : while( !((sumlen[k] < sfrac[i]) && (sumlen[k+incr]>=sfrac[i])) ) {
101 2093 : k+=incr;
102 2093 : if( cfin==iend && k>= iend+1 ) plumed_merror("path reparameterization error");
103 2093 : else if( cfin==(iend+1) && k>=iend ) { k=iend-1; break; }
104 2090 : else if( cfin==(iend-1) && k<=iend ) { k=iend+1; break; }
105 : }
106 1443 : double dr = (sfrac[i]-sumlen[k])/len[k+incr];
107 : // Calculate the displacement between the appropriate points
108 : // double dd = mypath[k]->calc( mypath[k+incr]->getReferencePositions(), pbc, args, mypath[k+incr]->getReferenceArguments(), mypack, true );
109 : // Copy the reference configuration from the configuration to a tempory direction
110 1443 : newpath[i].setDirection( mypath[k]->getReferencePositions(), mypath[k]->getReferenceArguments() );
111 : // Get the displacement of the path
112 1924 : mypath[k]->extractDisplacementVector( mypath[k+incr]->getReferencePositions(), args, mypath[k+incr]->getReferenceArguments(), false, mydir );
113 : // Set our direction equal to the displacement
114 : // mydir.setDirection( mypack );
115 : // Shift the reference configuration by this ammount
116 962 : newpath[i].displaceReferenceConfiguration( dr, mydir );
117 : }
118 :
119 : // Copy the positions of the new path to the new paths
120 542 : for(int i=istart+incr; loopEnd(i,cfin,incr)==false; i+=incr) {
121 1924 : mypath[i]->setReferenceConfig( newpath[i].getReferencePositions(), newpath[i].getReferenceArguments(), mypath[i]->getReferenceMetric() );
122 : }
123 :
124 : // Recompute the separations between frames
125 61 : calcCurrentPathSpacings( istart, iend );
126 : }
127 10 : }
128 :
129 4 : void PathReparameterization::reparameterize( const int& ifix1, const int& ifix2, const double& TOL ) {
130 : plumed_dbg_assert( ifix1<ifix2 );
131 : // First reparameterize the part between the fixed frames
132 4 : reparameterizePart( ifix1, ifix2, -1.0, TOL );
133 :
134 : // Get the separation between frames which we will use to set the remaining frames
135 8 : double target = sumlen[ifix2] / ( ifix2 - ifix1 );
136 :
137 : // And reparameterize the begining and end of the path
138 4 : if( ifix1>0 ) reparameterizePart( ifix1, 0, target, TOL );
139 8 : if( ifix2<(mypath.size()-1) ) reparameterizePart( ifix2, mypath.size()-1, target, TOL );
140 : // calcCurrentPathSpacings( 0, mypath.size()-1 );
141 4 : }
142 :
143 : }
144 4839 : }
|
