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 {
109 101 : return 0.0;
110 : }
111 : double transformHD( const double& dist, double& df ) const override;
112 : void update() override;
113 : };
114 :
115 13787 : PLUMED_REGISTER_ACTION(AdaptivePath,"ADAPTIVE_PATH")
116 :
117 5 : void AdaptivePath::registerKeywords( Keywords& keys ) {
118 5 : Mapping::registerKeywords( keys );
119 5 : keys.remove("PROPERTY");
120 10 : 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");
121 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");
122 10 : keys.add("compulsory","UPDATE","the frequency with which the path should be updated");
123 10 : keys.add("compulsory","TOLERANCE","1E-6","the tolerance to use for the path updating algorithm that makes all frames equidistant");
124 10 : keys.add("optional","WFILE","file on which to write out the path");
125 10 : keys.add("compulsory","FMT","%f","the format to use for output files");
126 10 : keys.add("optional","WSTRIDE","frequency with which to write out the path");
127 5 : }
128 :
129 1 : AdaptivePath::AdaptivePath(const ActionOptions& ao):
130 : Action(ao),
131 : Mapping(ao),
132 1 : fixedn(2),
133 2 : displacement(ReferenceConfigurationOptions("DIRECTION")),
134 3 : displacement2(ReferenceConfigurationOptions("DIRECTION")) {
135 : setLowMemOption( true );
136 2 : parseVector("FIXED",fixedn);
137 1 : if( fixedn[0]<1 || fixedn[1]>getNumberOfReferencePoints() ) {
138 0 : error("fixed nodes must be in range from 0 to number of nodes");
139 : }
140 1 : if( fixedn[0]>=fixedn[1] ) {
141 0 : error("invalid selection for fixed nodes first index provided must be smaller than second index");
142 : }
143 1 : log.printf(" fixing position of frames numbered %u and %u \n",fixedn[0],fixedn[1]);
144 1 : fixedn[0]--;
145 1 : fixedn[1]--; // Set fixed notes with c++ indexing starting from zero
146 1 : parse("UPDATE",update_str);
147 1 : if( update_str<1 ) {
148 0 : error("update frequency for path should be greater than or equal to one");
149 : }
150 1 : log.printf(" updating path every %u MD steps \n",update_str);
151 :
152 : double halflife;
153 1 : parse("HALFLIFE",halflife);
154 1 : if( halflife<0 ) {
155 1 : fadefact=1.0;
156 : } else {
157 0 : fadefact = exp( -0.693147180559945 / static_cast<double>(halflife) );
158 0 : log.printf(" weight of contribution to frame halves every %f steps \n",halflife);
159 : }
160 :
161 : // Create the list of tasks (and reset projections of frames)
162 1 : PDB mypdb;
163 1 : mypdb.setAtomNumbers( getAbsoluteIndexes() );
164 1 : mypdb.addBlockEnd( getAbsoluteIndexes().size() );
165 1 : std::vector<std::string> argument_names( getNumberOfArguments() );
166 3 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
167 : argument_names[i] = getPntrToArgument(i)->getName();
168 : }
169 1 : if( argument_names.size()>0 ) {
170 1 : mypdb.setArgumentNames( argument_names );
171 : }
172 1 : displacement.read( mypdb );
173 1 : displacement2.read( mypdb );
174 21 : for(int i=0; i<getNumberOfReferencePoints(); ++i) {
175 20 : addTaskToList( i );
176 40 : pdisplacements.push_back( Direction(ReferenceConfigurationOptions("DIRECTION")) );
177 40 : property.find("spath")->second[i] = static_cast<double>( i - static_cast<int>(fixedn[0]) ) / static_cast<double>( fixedn[1] - fixedn[0] );
178 20 : pdisplacements[i].read( mypdb );
179 20 : wsum.push_back( 0.0 );
180 : }
181 2 : plumed_assert( getPropertyValue( fixedn[0], "spath" )==0.0 && getPropertyValue( fixedn[1], "spath" )==1.0 );
182 : // And activate them all
183 1 : deactivateAllTasks();
184 21 : for(unsigned i=0; i<getFullNumberOfTasks(); ++i) {
185 20 : taskFlags[i]=1;
186 : }
187 1 : lockContributors();
188 :
189 : // Setup the vessel to hold the trig path
190 : std::string input;
191 1 : addVessel("GPATH", input, -1 );
192 1 : readVesselKeywords();
193 : // Check that there is only one vessel - the one holding the trig path
194 : plumed_dbg_assert( getNumberOfVessels()==1 );
195 : // Retrieve the path vessel
196 1 : mypathv = dynamic_cast<TrigonometricPathVessel*>( getPntrToVessel(0) );
197 1 : plumed_assert( mypathv );
198 :
199 : // Information for write out
200 : std::string wfilename;
201 2 : parse("WFILE",wfilename);
202 1 : if( wfilename.length()>0 ) {
203 1 : wstride=0;
204 1 : parse("WSTRIDE",wstride);
205 1 : parse("FMT",ofmt);
206 1 : pathfile.link(*this);
207 1 : pathfile.open( wfilename );
208 : pathfile.setHeavyFlush();
209 1 : if( wstride<update_str ) {
210 0 : error("makes no sense to write out path more frequently than update stride");
211 : }
212 1 : log.printf(" writing path out every %u steps to file named %s with format %s \n",wstride,wfilename.c_str(),ofmt.c_str());
213 : }
214 2 : log<<" Bibliography "<<plumed.cite("Diaz Leines and Ensing, Phys. Rev. Lett. 109, 020601 (2012)")<<"\n";
215 1 : }
216 :
217 101 : void AdaptivePath::calculate() {
218 101 : runAllTasks();
219 101 : }
220 :
221 2020 : void AdaptivePath::performTask( const unsigned& task_index, const unsigned& current, MultiValue& myvals ) const {
222 : // This builds a pack to hold the derivatives
223 2020 : ReferenceValuePack mypack( getNumberOfArguments(), getNumberOfAtoms(), myvals );
224 2020 : finishPackSetup( current, mypack );
225 : // Calculate the distance from the frame
226 2020 : double val=calculateDistanceFunction( current, mypack, true );
227 : // Put the element value in element zero
228 : myvals.setValue( 0, val );
229 : myvals.setValue( 1, 1.0 );
230 2020 : return;
231 2020 : }
232 :
233 2020 : double AdaptivePath::transformHD( const double& dist, double& df ) const {
234 2020 : df=1.0;
235 2020 : return dist;
236 : }
237 :
238 101 : void AdaptivePath::update() {
239 101 : double weight2 = -1.*mypathv->dx;
240 101 : double weight1 = 1.0 + mypathv->dx;
241 101 : if( weight1>1.0 ) {
242 0 : weight1=1.0;
243 0 : weight2=0.0;
244 101 : } else if( weight2>1.0 ) {
245 0 : weight1=0.0;
246 0 : weight2=1.0;
247 : }
248 : // Add projections to dispalcement accumulators
249 101 : ReferenceConfiguration* myref = getReferenceConfiguration( mypathv->iclose1 );
250 101 : myref->extractDisplacementVector( getPositions(), getArguments(), mypathv->cargs, false, displacement );
251 101 : getReferenceConfiguration( mypathv->iclose2 )->extractDisplacementVector( myref->getReferencePositions(), getArguments(), myref->getReferenceArguments(), false, displacement2 );
252 101 : displacement.addDirection( -mypathv->dx, displacement2 );
253 101 : pdisplacements[mypathv->iclose1].addDirection( weight1, displacement );
254 101 : pdisplacements[mypathv->iclose2].addDirection( weight2, displacement );
255 : // Update weight accumulators
256 101 : wsum[mypathv->iclose1] *= fadefact;
257 101 : wsum[mypathv->iclose2] *= fadefact;
258 101 : wsum[mypathv->iclose1] += weight1;
259 101 : wsum[mypathv->iclose2] += weight2;
260 :
261 : // This does the update of the path if it is time to
262 101 : if( (getStep()>0) && (getStep()%update_str==0) ) {
263 2 : wsum[fixedn[0]]=wsum[fixedn[1]]=0.;
264 42 : for(unsigned inode=0; inode<getNumberOfReferencePoints(); ++inode) {
265 40 : if( wsum[inode]>0 ) {
266 : // First displace the node by the weighted direction
267 6 : getReferenceConfiguration( inode )->displaceReferenceConfiguration( 1./wsum[inode], pdisplacements[inode] );
268 : // Reset the displacement
269 6 : pdisplacements[inode].zeroDirection();
270 : }
271 : }
272 : // Now ensure all the nodes of the path are equally spaced
273 2 : PathReparameterization myspacings( getPbc(), getArguments(), getAllReferenceConfigurations() );
274 2 : myspacings.reparameterize( fixedn[0], fixedn[1], tolerance );
275 2 : }
276 101 : if( (getStep()>0) && (getStep()%wstride==0) ) {
277 2 : pathfile<<"# PATH AT STEP "<<getStep();
278 2 : pathfile.printf(" TIME %f \n",getTime());
279 : std::vector<std::unique_ptr<ReferenceConfiguration>>& myconfs=getAllReferenceConfigurations();
280 2 : auto* mymoldat=plumed.getActionSet().selectLatest<GenericMolInfo*>(this);
281 2 : std::vector<std::string> argument_names( getNumberOfArguments() );
282 6 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
283 : argument_names[i] = getPntrToArgument(i)->getName();
284 : }
285 2 : PDB mypdb;
286 2 : mypdb.setArgumentNames( argument_names );
287 42 : for(unsigned i=0; i<myconfs.size(); ++i) {
288 80 : pathfile.printf("REMARK TYPE=%s\n", myconfs[i]->getName().c_str() );
289 40 : mypdb.setAtomPositions( myconfs[i]->getReferencePositions() );
290 120 : for(unsigned j=0; j<getNumberOfArguments(); ++j) {
291 80 : mypdb.setArgumentValue( getPntrToArgument(j)->getName(), myconfs[i]->getReferenceArgument(j) );
292 : }
293 40 : mypdb.print( atoms.getUnits().getLength()/0.1, mymoldat, pathfile, ofmt );
294 : }
295 2 : pathfile.flush();
296 2 : }
297 101 : }
298 :
299 : }
300 : }
|
