Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2016-2023 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<std::unique_ptr<ReferenceConfiguration>>& pp ):
28 4 : mydpack( 1, pp[0]->getNumberOfReferenceArguments() + 3*pp[0]->getNumberOfReferencePositions() + 9 ),
29 4 : mypack( pp[0]->getNumberOfReferenceArguments(), pp[0]->getNumberOfReferencePositions(), mydpack ),
30 8 : mydir(ReferenceConfigurationOptions("DIRECTION")),
31 4 : pbc(ipbc),
32 4 : args(iargs),
33 4 : mypath(pp),
34 4 : len(pp.size()),
35 4 : sumlen(pp.size()),
36 4 : sfrac(pp.size()),
37 8 : MAXCYCLES(100)
38 : {
39 4 : mypdb.setAtomNumbers( pp[0]->getAbsoluteIndexes() ); mypdb.addBlockEnd( pp[0]->getAbsoluteIndexes().size() );
40 4 : if( pp[0]->getArgumentNames().size()>0 ) mypdb.setArgumentNames( pp[0]->getArgumentNames() );
41 4 : mydir.read( mypdb ); mydir.zeroDirection(); pp[0]->setupPCAStorage( mypack );
42 4 : }
43 :
44 1741 : bool PathReparameterization::loopEnd( const int& index, const int& end, const int& inc ) const {
45 1741 : if( inc>0 && index<end ) return false;
46 270 : else if( inc<0 && index>end ) return false;
47 : return true;
48 : }
49 :
50 71 : void PathReparameterization::calcCurrentPathSpacings( const int& istart, const int& iend ) {
51 : plumed_dbg_assert( istart<len.size() && iend<len.size() );
52 71 : len[istart] = sumlen[istart]=0;
53 : //printf("HELLO PATH SPACINGS ARE CURRENTLY \n");
54 :
55 : // Get the spacings given we can go forward and backwards
56 71 : int incr=1; if( istart>iend ) { incr=-1; }
57 :
58 657 : for(int i=istart+incr; loopEnd(i,iend+incr,incr)==false; i+=incr) {
59 586 : len[i] = mypath[i-incr]->calc( mypath[i]->getReferencePositions(), pbc, args, mypath[i]->getReferenceArguments(), mypack, false );
60 586 : sumlen[i] = sumlen[i-incr] + len[i];
61 : //printf("FRAME %d TO FRAME %d EQUALS %f : %f \n",i-incr,i,len[i],sumlen[i] );
62 : }
63 71 : }
64 :
65 10 : void PathReparameterization::reparameterizePart( const int& istart, const int& iend, const double& target, const double& TOL ) {
66 10 : calcCurrentPathSpacings( istart, iend ); unsigned cfin;
67 : // If a target separation is set we fix where we want the nodes
68 10 : int incr=1; if( istart>iend ) { incr=-1; }
69 :
70 10 : if( target>0 ) {
71 6 : if( iend>istart ) {
72 19 : for(unsigned i=istart; i<iend+1; ++i) sfrac[i] = target*(i-istart);
73 : } else {
74 14 : for(int i=istart-1; i>iend-1; --i) sfrac[i]=target*(istart-i);
75 : }
76 6 : cfin = iend+incr;
77 : } else {
78 4 : cfin = iend;
79 : }
80 :
81 : std::vector<Direction> newpath;
82 170 : for(unsigned i=0; i<mypath.size(); ++i) {
83 320 : newpath.push_back( Direction(ReferenceConfigurationOptions("DIRECTION")) ); newpath[i].read( mypdb );
84 : }
85 :
86 : double prevsum=0.;
87 71 : for(unsigned iter=0; iter<MAXCYCLES; ++iter) {
88 71 : if( std::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 2567 : 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 481 : 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 481 : newpath[i].setDirection( mypath[k]->getReferencePositions(), mypath[k]->getReferenceArguments() );
111 : // Get the displacement of the path
112 481 : 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 amount
116 481 : 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 481 : mypdb.setAtomPositions( newpath[i].getReferencePositions() );
122 1415 : for(unsigned j=0; j<newpath[i].getNumberOfReferenceArguments(); ++j) mypdb.setArgumentValue( mypath[i]->getArgumentNames()[j], newpath[i].getReferenceArgument(j) );
123 481 : mypath[i]->read( mypdb );
124 : }
125 :
126 : // Recompute the separations between frames
127 61 : calcCurrentPathSpacings( istart, iend );
128 : }
129 10 : }
130 :
131 4 : void PathReparameterization::reparameterize( const int& ifix1, const int& ifix2, const double& TOL ) {
132 : plumed_dbg_assert( ifix1<ifix2 );
133 : // First reparameterize the part between the fixed frames
134 4 : reparameterizePart( ifix1, ifix2, -1.0, TOL );
135 :
136 : // Get the separation between frames which we will use to set the remaining frames
137 4 : double target = sumlen[ifix2] / ( ifix2 - ifix1 );
138 :
139 : // And reparameterize the beginning and end of the path
140 4 : if( ifix1>0 ) reparameterizePart( ifix1, 0, target, TOL );
141 4 : if( ifix2<(mypath.size()-1) ) reparameterizePart( ifix2, mypath.size()-1, target, TOL );
142 : // calcCurrentPathSpacings( 0, mypath.size()-1 );
143 4 : }
144 :
145 : }
146 : }
|
