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