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 4 : mypdb.setAtomNumbers( pp[0]->getAbsoluteIndexes() );
39 4 : mypdb.addBlockEnd( pp[0]->getAbsoluteIndexes().size() );
40 4 : if( pp[0]->getArgumentNames().size()>0 ) {
41 3 : mypdb.setArgumentNames( pp[0]->getArgumentNames() );
42 : }
43 4 : mydir.read( mypdb );
44 4 : mydir.zeroDirection();
45 4 : pp[0]->setupPCAStorage( mypack );
46 4 : }
47 :
48 1741 : bool PathReparameterization::loopEnd( const int& index, const int& end, const int& inc ) const {
49 1741 : if( inc>0 && index<end ) {
50 : return false;
51 270 : } else if( inc<0 && index>end ) {
52 77 : return false;
53 : }
54 : return true;
55 : }
56 :
57 71 : void PathReparameterization::calcCurrentPathSpacings( const int& istart, const int& iend ) {
58 : plumed_dbg_assert( istart<len.size() && iend<len.size() );
59 71 : len[istart] = sumlen[istart]=0;
60 : //printf("HELLO PATH SPACINGS ARE CURRENTLY \n");
61 :
62 : // Get the spacings given we can go forward and backwards
63 71 : int incr=1;
64 71 : if( istart>iend ) {
65 9 : incr=-1;
66 : }
67 :
68 657 : for(int i=istart+incr; loopEnd(i,iend+incr,incr)==false; i+=incr) {
69 586 : len[i] = mypath[i-incr]->calc( mypath[i]->getReferencePositions(), pbc, args, mypath[i]->getReferenceArguments(), mypack, false );
70 586 : sumlen[i] = sumlen[i-incr] + len[i];
71 : //printf("FRAME %d TO FRAME %d EQUALS %f : %f \n",i-incr,i,len[i],sumlen[i] );
72 : }
73 71 : }
74 :
75 10 : void PathReparameterization::reparameterizePart( const int& istart, const int& iend, const double& target, const double& TOL ) {
76 10 : calcCurrentPathSpacings( istart, iend );
77 : unsigned cfin;
78 : // If a target separation is set we fix where we want the nodes
79 10 : int incr=1;
80 10 : if( istart>iend ) {
81 3 : incr=-1;
82 : }
83 :
84 10 : if( target>0 ) {
85 6 : if( iend>istart ) {
86 19 : for(unsigned i=istart; i<iend+1; ++i) {
87 16 : sfrac[i] = target*(i-istart);
88 : }
89 : } else {
90 14 : for(int i=istart-1; i>iend-1; --i) {
91 11 : sfrac[i]=target*(istart-i);
92 : }
93 : }
94 6 : cfin = iend+incr;
95 : } else {
96 4 : cfin = iend;
97 : }
98 :
99 : std::vector<Direction> newpath;
100 170 : for(unsigned i=0; i<mypath.size(); ++i) {
101 320 : newpath.push_back( Direction(ReferenceConfigurationOptions("DIRECTION")) );
102 160 : newpath[i].read( mypdb );
103 : }
104 :
105 : double prevsum=0.;
106 71 : for(unsigned iter=0; iter<MAXCYCLES; ++iter) {
107 71 : if( std::fabs(sumlen[iend] - prevsum)<=TOL ) {
108 : break ;
109 : }
110 : prevsum = sumlen[iend];
111 : // If no target is set we redistribute length
112 61 : if( target<0 ) {
113 49 : plumed_assert( istart<iend );
114 49 : double dr = sumlen[iend] / static_cast<double>( iend - istart );
115 531 : for(unsigned i=istart; i<iend; ++i) {
116 482 : sfrac[i] = dr*(i-istart);
117 : }
118 : }
119 :
120 : // Now compute positions of new nodes in path
121 542 : for(int i=istart+incr; loopEnd(i,cfin,incr)==false; i+=incr) {
122 481 : int k = istart;
123 2567 : while( !((sumlen[k] < sfrac[i]) && (sumlen[k+incr]>=sfrac[i])) ) {
124 2093 : k+=incr;
125 2093 : if( cfin==iend && k>= iend+1 ) {
126 0 : plumed_merror("path reparameterization error");
127 2093 : } else if( cfin==(iend+1) && k>=iend ) {
128 3 : k=iend-1;
129 3 : break;
130 2090 : } else if( cfin==(iend-1) && k<=iend ) {
131 : k=iend+1;
132 : break;
133 : }
134 : }
135 481 : double dr = (sfrac[i]-sumlen[k])/len[k+incr];
136 : // Calculate the displacement between the appropriate points
137 : // double dd = mypath[k]->calc( mypath[k+incr]->getReferencePositions(), pbc, args, mypath[k+incr]->getReferenceArguments(), mypack, true );
138 : // Copy the reference configuration from the configuration to a tempory direction
139 481 : newpath[i].setDirection( mypath[k]->getReferencePositions(), mypath[k]->getReferenceArguments() );
140 : // Get the displacement of the path
141 481 : mypath[k]->extractDisplacementVector( mypath[k+incr]->getReferencePositions(), args, mypath[k+incr]->getReferenceArguments(), false, mydir );
142 : // Set our direction equal to the displacement
143 : // mydir.setDirection( mypack );
144 : // Shift the reference configuration by this amount
145 481 : newpath[i].displaceReferenceConfiguration( dr, mydir );
146 : }
147 :
148 : // Copy the positions of the new path to the new paths
149 542 : for(int i=istart+incr; loopEnd(i,cfin,incr)==false; i+=incr) {
150 481 : mypdb.setAtomPositions( newpath[i].getReferencePositions() );
151 1415 : for(unsigned j=0; j<newpath[i].getNumberOfReferenceArguments(); ++j) {
152 934 : mypdb.setArgumentValue( mypath[i]->getArgumentNames()[j], newpath[i].getReferenceArgument(j) );
153 : }
154 481 : mypath[i]->read( mypdb );
155 : }
156 :
157 : // Recompute the separations between frames
158 61 : calcCurrentPathSpacings( istart, iend );
159 : }
160 10 : }
161 :
162 4 : void PathReparameterization::reparameterize( const int& ifix1, const int& ifix2, const double& TOL ) {
163 : plumed_dbg_assert( ifix1<ifix2 );
164 : // First reparameterize the part between the fixed frames
165 4 : reparameterizePart( ifix1, ifix2, -1.0, TOL );
166 :
167 : // Get the separation between frames which we will use to set the remaining frames
168 4 : double target = sumlen[ifix2] / ( ifix2 - ifix1 );
169 :
170 : // And reparameterize the beginning and end of the path
171 4 : if( ifix1>0 ) {
172 3 : reparameterizePart( ifix1, 0, target, TOL );
173 : }
174 4 : if( ifix2<(mypath.size()-1) ) {
175 3 : reparameterizePart( ifix2, mypath.size()-1, target, TOL );
176 : }
177 : // calcCurrentPathSpacings( 0, mypath.size()-1 );
178 4 : }
179 :
180 : }
181 : }
|
