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 "Mapping.h"
23 : #include "TrigonometricPathVessel.h"
24 : #include "PathReparameterization.h"
25 : #include "reference/Direction.h"
26 : #include "core/ActionRegister.h"
27 : #include "core/PlumedMain.h"
28 :
29 : //+PLUMEDOC COLVAR ADAPTIVE_PATH
30 : /*
31 : Compute path collective variables that adapt to the lowest free energy path connecting states A and B.
32 :
33 : The Path Collective Variables developed by Branduardi and co-workers \cite brand07 allow one
34 : to compute the progress along a high-dimensional path and the distance from the high-dimensional
35 : path. The progress along the path (s) is computed using:
36 :
37 : \f[
38 : s = i_2 + \textrm{sign}(i_2-i_1) \frac{ \sqrt{( \mathbf{v}_1\cdot\mathbf{v}_2 )^2 - |\mathbf{v}_3|^2(|\mathbf{v}_1|^2 - |\mathbf{v}_2|^2) } }{2|\mathbf{v}_3|^2} - \frac{\mathbf{v}_1\cdot\mathbf{v}_3 - |\mathbf{v}_3|^2}{2|\mathbf{v}_3|^2}
39 : \f]
40 :
41 : In this expression \f$\mathbf{v}_1\f$ and \f$\mathbf{v}_3\f$ are the vectors connecting the current position to the closest and second closest node of the path,
42 : respectfully and \f$i_1\f$ and \f$i_2\f$ are the projections of the closest and second closest frames of the path. \f$\mathbf{v}_2\f$, meanwhile, is the
43 : vector connecting the closest frame to the second closest frame. The distance from the path, \f$z\f$ is calculated using:
44 :
45 : \f[
46 : z = \sqrt{ \left[ |\mathbf{v}_1|^2 - |\mathbf{v}_2| \left( \frac{ \sqrt{( \mathbf{v}_1\cdot\mathbf{v}_2 )^2 - |\mathbf{v}_3|^2(|\mathbf{v}_1|^2 - |\mathbf{v}_2|^2) } }{2|\mathbf{v}_3|^2} - \frac{\mathbf{v}_1\cdot\mathbf{v}_3 - |\mathbf{v}_3|^2}{2|\mathbf{v}_3|^2} \right) \right]^2 }
47 : \f]
48 :
49 : Notice that these are the definitions of \f$s\f$ and \f$z\f$ that are used by \ref PATH when the GPATH option is employed. The reason for this is that
50 : the adaptive path method implemented in this action was inspired by the work of Diaz and Ensing in which these formula were used \cite BerndAdaptivePath.
51 : To learn more about how the path is adapted we strongly recommend reading this paper.
52 :
53 : \par Examples
54 :
55 : The input below provides an example of how the adaptive path works in practise. The path is updated every 50 steps of
56 : MD based on the data accumulated during the preceding 50 time steps.
57 :
58 : \plumedfile
59 : d1: DISTANCE ATOMS=1,2 COMPONENTS
60 : pp: ADAPTIVE_PATH TYPE=EUCLIDEAN FIXED=5,15 UPDATE=50 WFILE=out-path.pdb WSTRIDE=50 REFERENCE=mypath.pdb
61 : PRINT ARG=d1.x,d1.y,pp.* FILE=colvar
62 : \endplumedfile
63 :
64 : In the case above the distance between frames is calculated based on the \f$x\f$ and \f$y\f$ components of the vector connecting
65 : atoms 1 and 2. As such an extract from the input reference path (mypath.pdb) would look as follows:
66 :
67 : \verbatim
68 : REMARK ARG=d1.x,d1.y d1.x=1.12 d1.y=-.60
69 : END
70 : REMARK ARG=d1.x,d1.y d1.x=.99 d1.y=-.45
71 : END
72 : \endverbatim
73 :
74 : Notice that one can also use RMSD frames in place of arguments like those above.
75 :
76 : */
77 : //+ENDPLUMEDOC
78 :
79 : namespace PLMD {
80 : namespace mapping {
81 :
82 3 : class AdaptivePath : public Mapping {
83 : private:
84 : OFile pathfile;
85 : std::string ofmt;
86 : double fadefact, tolerance;
87 : unsigned update_str, wstride;
88 : std::vector<unsigned> fixedn;
89 : TrigonometricPathVessel* mypathv;
90 : std::vector<double> wsum;
91 : Direction displacement,displacement2;
92 : std::vector<Direction> pdisplacements;
93 : public:
94 : static void registerKeywords( Keywords& keys );
95 : explicit AdaptivePath(const ActionOptions&);
96 : void calculate();
97 : void performTask( const unsigned&, const unsigned&, MultiValue& ) const ;
98 101 : double getLambda() { return 0.0; }
99 : double transformHD( const double& dist, double& df ) const ;
100 : void update();
101 : };
102 :
103 6453 : PLUMED_REGISTER_ACTION(AdaptivePath,"ADAPTIVE_PATH")
104 :
105 2 : void AdaptivePath::registerKeywords( Keywords& keys ) {
106 4 : Mapping::registerKeywords( keys ); keys.remove("PROPERTY");
107 8 : keys.add("compulsory","FIXED","the positions in the list of input frames of the two path nodes whose positions remain fixed during the path optimization");
108 10 : keys.add("compulsory","HALFLIFE","-1","the number of MD steps after which a previously measured path distance weighs only 50% in the average. This option may increase convergence by allowing to \"forget\" the memory of a bad initial guess path. The default is to set this to infinity");
109 8 : keys.add("compulsory","UPDATE","the frequency with which the path should be updated");
110 10 : keys.add("compulsory","TOLERANCE","1E-6","the tolerance to use for the path updating algorithm that makes all frames equidistant");
111 8 : keys.add("optional","WFILE","file on which to write out the path");
112 10 : keys.add("compulsory","FMT","%f","the format to use for output files");
113 8 : keys.add("optional","WSTRIDE","frequency with which to write out the path");
114 2 : }
115 :
116 1 : AdaptivePath::AdaptivePath(const ActionOptions& ao):
117 : Action(ao),
118 : Mapping(ao),
119 : fixedn(2),
120 3 : displacement(ReferenceConfigurationOptions("DIRECTION")),
121 7 : displacement2(ReferenceConfigurationOptions("DIRECTION"))
122 : {
123 2 : setLowMemOption( true ); parseVector("FIXED",fixedn);
124 2 : if( fixedn[0]<1 || fixedn[1]>getNumberOfReferencePoints() ) error("fixed nodes must be in range from 0 to number of nodes");
125 1 : if( fixedn[0]>=fixedn[1] ) error("invalid selection for fixed nodes first index provided must be smaller than second index");
126 2 : log.printf(" fixing position of frames numbered %u and %u \n",fixedn[0],fixedn[1]);
127 2 : fixedn[0]--; fixedn[1]--; // Set fixed notes with c++ indexing starting from zero
128 2 : parse("UPDATE",update_str); if( update_str<1 ) error("update frequency for path should be greater than or equal to one");
129 1 : log.printf(" updating path every %u MD steps \n",update_str);
130 :
131 2 : double halflife; parse("HALFLIFE",halflife);
132 1 : log.printf(" weight of contribution to frame halves every %f steps \n",halflife);
133 1 : if( halflife<0 ) fadefact=1.0;
134 0 : else fadefact = exp( -0.693147180559945 / static_cast<double>(halflife) );
135 :
136 : // Create the list of tasks (and reset projections of frames)
137 2 : std::vector<std::string> argument_names( getNumberOfArguments() );
138 5 : for(unsigned i=0; i<getNumberOfArguments(); ++i) argument_names[i] = getPntrToArgument(i)->getName();
139 2 : displacement.setNamesAndAtomNumbers( getAbsoluteIndexes(), argument_names );
140 1 : displacement2.setNamesAndAtomNumbers( getAbsoluteIndexes(), argument_names );
141 41 : for(int i=0; i<getNumberOfReferencePoints(); ++i) {
142 80 : addTaskToList( i ); pdisplacements.push_back( Direction(ReferenceConfigurationOptions("DIRECTION")) );
143 20 : setPropertyValue( i, 0, static_cast<double>( i - static_cast<int>(fixedn[0]) ) / static_cast<double>( fixedn[1] - fixedn[0] ) );
144 60 : pdisplacements[i].setNamesAndAtomNumbers( getAbsoluteIndexes(), argument_names ); wsum.push_back( 0.0 );
145 : }
146 2 : plumed_assert( getPropertyValue( fixedn[0], 0 )==0.0 && getPropertyValue( fixedn[1], 0 )==1.0 );
147 : // And activate them all
148 1 : deactivateAllTasks();
149 41 : for(unsigned i=0; i<getFullNumberOfTasks(); ++i) taskFlags[i]=1;
150 1 : lockContributors();
151 :
152 : // Setup the vessel to hold the trig path
153 2 : std::string input; addVessel("GPATH", input, -1 );
154 1 : readVesselKeywords();
155 : // Check that there is only one vessel - the one holding the trig path
156 : plumed_dbg_assert( getNumberOfVessels()==1 );
157 : // Retrieve the path vessel
158 1 : mypathv = dynamic_cast<TrigonometricPathVessel*>( getPntrToVessel(0) );
159 1 : plumed_assert( mypathv );
160 :
161 : // Information for write out
162 2 : std::string wfilename; parse("WFILE",wfilename);
163 1 : if( wfilename.length()>0 ) {
164 3 : wstride=0; parse("WSTRIDE",wstride); parse("FMT",ofmt);
165 1 : pathfile.link(*this); pathfile.open( wfilename ); pathfile.setHeavyFlush();
166 1 : if( wstride<update_str ) error("makes no sense to write out path more frequently than update stride");
167 2 : log.printf(" writing path out every %u steps to file named %s with format %s \n",wstride,wfilename.c_str(),ofmt.c_str());
168 : }
169 3 : log<<" Bibliography "<<plumed.cite("Diaz Leines and Ensing, Phys. Rev. Lett. 109, 020601 (2012)")<<"\n";
170 1 : }
171 :
172 101 : void AdaptivePath::calculate() {
173 101 : runAllTasks();
174 101 : }
175 :
176 2020 : void AdaptivePath::performTask( const unsigned& task_index, const unsigned& current, MultiValue& myvals ) const {
177 : // This builds a pack to hold the derivatives
178 4040 : ReferenceValuePack mypack( getNumberOfArguments(), getNumberOfAtoms(), myvals );
179 2020 : finishPackSetup( current, mypack );
180 : // Calculate the distance from the frame
181 2020 : double val=calculateDistanceFunction( current, mypack, true );
182 : // Put the element value in element zero
183 : myvals.setValue( 0, val ); myvals.setValue( 1, 1.0 );
184 2020 : return;
185 : }
186 :
187 2020 : double AdaptivePath::transformHD( const double& dist, double& df ) const {
188 2020 : df=1.0; return dist;
189 : }
190 :
191 101 : void AdaptivePath::update() {
192 101 : double weight2 = -1.*mypathv->dx;
193 101 : double weight1 = 1.0 + mypathv->dx;
194 101 : if( weight1>1.0 ) {
195 0 : weight1=1.0; weight2=0.0;
196 101 : } else if( weight2>1.0 ) {
197 0 : weight1=0.0; weight2=1.0;
198 : }
199 : // Add projections to dispalcement accumulators
200 101 : ReferenceConfiguration* myref = getReferenceConfiguration( mypathv->iclose1 );
201 202 : myref->extractDisplacementVector( getPositions(), getArguments(), mypathv->cargs, false, displacement );
202 101 : getReferenceConfiguration( mypathv->iclose2 )->extractDisplacementVector( myref->getReferencePositions(), getArguments(), myref->getReferenceArguments(), false, displacement2 );
203 101 : displacement.addDirection( -mypathv->dx, displacement2 );
204 202 : pdisplacements[mypathv->iclose1].addDirection( weight1, displacement );
205 202 : pdisplacements[mypathv->iclose2].addDirection( weight2, displacement );
206 : // Update weight accumulators
207 202 : wsum[mypathv->iclose1] *= fadefact;
208 202 : wsum[mypathv->iclose2] *= fadefact;
209 202 : wsum[mypathv->iclose1] += weight1;
210 202 : wsum[mypathv->iclose2] += weight2;
211 :
212 : // This does the update of the path if it is time to
213 101 : if( (getStep()>0) && (getStep()%update_str==0) ) {
214 8 : wsum[fixedn[0]]=wsum[fixedn[1]]=0.;
215 84 : for(unsigned inode=0; inode<getNumberOfReferencePoints(); ++inode) {
216 80 : if( wsum[inode]>0 ) {
217 : // First displace the node by the weighted direction
218 12 : getReferenceConfiguration( inode )->displaceReferenceConfiguration( 1./wsum[inode], pdisplacements[inode] );
219 : // Reset the displacement
220 12 : pdisplacements[inode].zeroDirection();
221 : }
222 : }
223 : // Now ensure all the nodes of the path are equally spaced
224 4 : PathReparameterization myspacings( getPbc(), getArguments(), getAllReferenceConfigurations() );
225 4 : myspacings.reparameterize( fixedn[0], fixedn[1], tolerance );
226 : }
227 101 : if( (getStep()>0) && (getStep()%wstride==0) ) {
228 2 : pathfile.printf("# PATH AT STEP %d TIME %f \n", getStep(), getTime() );
229 : std::vector<ReferenceConfiguration*>& myconfs=getAllReferenceConfigurations();
230 124 : for(unsigned i=0; i<myconfs.size(); ++i) myconfs[i]->print( pathfile, ofmt, atoms.getUnits().getLength()/0.1 );
231 2 : pathfile.flush();
232 : }
233 101 : }
234 :
235 : }
236 4839 : }
|
