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 "cltools/CLTool.h"
23 : #include "cltools/CLToolRegister.h"
24 : #include "tools/Tools.h"
25 : #include "tools/Pbc.h"
26 : #include "core/Value.h"
27 : #include "reference/ReferenceConfiguration.h"
28 : #include "PathReparameterization.h"
29 : #include "reference/MetricRegister.h"
30 : #include <cstdio>
31 : #include <string>
32 : #include <vector>
33 : #include <iostream>
34 :
35 : using namespace std;
36 :
37 : namespace PLMD {
38 : namespace mapping {
39 :
40 : //+PLUMEDOC TOOLS pathtools
41 : /*
42 : pathtools can be used to construct paths from pdb data
43 :
44 : The path CVs in PLUMED are curvilinear coordinates through a high dimensional vector space.
45 : Enhanced sampling calculations are often run using the progress along the paths and the distance from the path as CVs
46 : as this provides a convenient way of defining a reaction coordinate for a complicated process. This method is explained
47 : in the documentation for \ref PATH.
48 :
49 : The path itself is an ordered set of equally-spaced, high-dimensional frames the way in which these frames
50 : should be constructed will depend on the problem in hand. In other words, you will need to understand the reaction
51 : you wish to study in order to select a sensible set of frames to use in your path CV. This tool provides two
52 : methods that may be useful when it comes to constructing paths; namely:
53 :
54 : - A tool that takes in an initial guess path in which the frames are not equally spaced. This tool adjusts the positions
55 : of the frames in order to make them equally spaced so that they can be used as the basis for a path CV.
56 :
57 : - A tool that takes two frames as input and that allows you to return a linear path connecting these two frames. The
58 : output from this method may be useful as an initial guess path. It is arguable that a linear path rather defeats the
59 : purpose of the path CV method, however, as the whole purpose is to be able to define non-linear paths.
60 :
61 : Notice that you can use these two methods and take advantage of all the ways of measuring \ref dists that are available within
62 : PLUMED. The way you do this with each of these tools described above is explained in the example below.
63 :
64 : \par Examples
65 :
66 : The example below shows how you can take a set of unequally spaced frames from a pdb file named in_path.pdb
67 : and use pathtools to make them equally spaced so that they can be used as the basis for a path CV. The file
68 : containing this final path is named final_path.pdb.
69 :
70 : \verbatim
71 : plumed pathtools --path in_path.pdb --metric EUCLIDEAN --out final_path.pdb
72 : \endverbatim
73 :
74 : The example below shows how can create an initial linear path connecting the two pdb frames in start.pdb and
75 : end.pdb. In this case the path output to path.pdb will consist of 6 frames: the initial and final frames that
76 : were contained in start.pdb and end.pdb as well as four equally spaced frames along the vector connecting
77 : start.pdb to end.pdb.
78 :
79 : \verbatim
80 : plumed pathtools --start start.pdb --end end.pdb --nframes 4 --metric OPTIMAL --out path.pdb
81 : \endverbatim
82 :
83 : Often the idea with path collective variables is to create a path connecting some initial state A to some final state B. You would
84 : in this case have representative configurations from your A and B states defined in the input files to pathtools
85 : that we have called start.pdb and end.pdb in the example above. Furthermore, it may be useful to have a few frames
86 : before your start frame and after your end frame. You can use path tools to create these extended paths as shown below.
87 : In this case the final path would now consist of 8 frames. Four of these frames would lie on the vector connecting state
88 : A to state B, there would be one frame each at start.pdb and end.pdb as well as one frame just before start.pdb and one
89 : frame just after end.pdb. All these frames would be equally spaced.
90 :
91 : \verbatim
92 : plumed pathtools --start start.pdb --end end.pdb --nframes 4 --metric OPTIMAL --out path.pdb --nframes-before-start 2 --nframes-after-end 2
93 : \endverbatim
94 :
95 : Notice also that when you re-parameterize paths you must choose two frames to fix. Generally you chose to fix the states
96 : that are representative of your states A and B. By default pathtools will fix the first and last frames. You can, however,
97 : change the states to fix by taking advantage of the fixed flag as shown below.
98 :
99 : \verbatim
100 : plumed pathtools --path inpath.pdb --metric EUCLIDEAN --out outpath.pdb --fixed 2,12
101 : \endverbatim
102 :
103 : */
104 : //+ENDPLUMEDOC
105 :
106 4 : class PathTools :
107 : public CLTool
108 : {
109 : public:
110 : static void registerKeywords( Keywords& keys );
111 : explicit PathTools(const CLToolOptions& co );
112 : int main(FILE* in, FILE*out,Communicator& pc);
113 0 : string description()const {
114 0 : return "print out a description of the keywords for an action in html";
115 : }
116 : };
117 :
118 7840 : PLUMED_REGISTER_CLTOOL(PathTools,"pathtools")
119 :
120 1958 : void PathTools::registerKeywords( Keywords& keys ) {
121 1958 : CLTool::registerKeywords( keys );
122 7832 : keys.add("atoms","--start","a pdb file that contains the structure for the initial frame of your path");
123 7832 : keys.add("atoms","--end","a pdb file that contains the structure for the final frame of your path");
124 7832 : keys.add("atoms-1","--path","a pdb file that contains an initial path in which the frames are not equally spaced");
125 9790 : keys.add("compulsory","--fixed","0","the frames to fix when constructing the path using --path");
126 7832 : keys.add("compulsory","--metric","the measure to use to calculate the distance between frames");
127 7832 : keys.add("compulsory","--out","the name of the file on which to output your path");
128 9790 : keys.add("compulsory","--arg-fmt","%f","the format to use for argument values in your frames");
129 9790 : keys.add("compulsory","--tolerance","1E-4","the tolerance to use for the algorithm that is used to re-parameterize the path");
130 9790 : keys.add("compulsory","--nframes-before-start","1","the number of frames to include in the path before the first frame");
131 9790 : keys.add("compulsory","--nframes","1","the number of frames between the start and end frames in your path");
132 9790 : keys.add("compulsory","--nframes-after-end","1","the number of frames to put after the last frame of your path");
133 1958 : }
134 :
135 4 : PathTools::PathTools(const CLToolOptions& co ):
136 4 : CLTool(co)
137 : {
138 4 : inputdata=commandline;
139 4 : }
140 :
141 4 : int PathTools::main(FILE* in, FILE*out,Communicator& pc) {
142 8 : std::string mtype; parse("--metric",mtype);
143 8 : std::string ifilename; parse("--path",ifilename);
144 8 : std::string ofmt; parse("--arg-fmt",ofmt);
145 8 : std::string ofilename; parse("--out",ofilename);
146 4 : if( ifilename.length()>0 ) {
147 : fprintf(out,"Reparameterising path in file named %s so that all frames are equally spaced \n",ifilename.c_str() );
148 2 : FILE* fp=fopen(ifilename.c_str(),"r");
149 : bool do_read=true; std::vector<std::unique_ptr<ReferenceConfiguration>> frames;
150 24 : while (do_read) {
151 22 : PDB mypdb;
152 : // Read the pdb file
153 22 : do_read=mypdb.readFromFilepointer(fp,false,0.1);
154 22 : if( do_read ) {
155 20 : auto mymsd(metricRegister().create<ReferenceConfiguration>( mtype, mypdb ));
156 20 : frames.emplace_back( std::move(mymsd) );
157 : }
158 22 : }
159 4 : std::vector<unsigned> fixed; parseVector("--fixed",fixed);
160 2 : if( fixed.size()==1 ) {
161 1 : if( fixed[0]!=0 ) error("input to --fixed should be two integers");
162 3 : fixed.resize(2); fixed[0]=0; fixed[1]=frames.size()-1;
163 1 : } else if( fixed.size()==2 ) {
164 2 : if( fixed[0]>(frames.size()-1) || fixed[1]>(frames.size()-1) ) {
165 0 : error("input to --fixed should be two numbers between 0 and the number of frames-1");
166 : }
167 : } else {
168 0 : error("input to --fixed should be two integers");
169 : }
170 2 : std::vector<AtomNumber> atoms; std::vector<std::string> arg_names;
171 44 : for(unsigned i=0; i<frames.size(); ++i) {
172 20 : frames[i]->getAtomRequests( atoms);
173 20 : frames[i]->getArgumentRequests( arg_names );
174 : }
175 : // Generate stuff to reparameterize
176 4 : Pbc fake_pbc; std::vector<std::unique_ptr<Value>> vals;
177 8 : for(unsigned i=0; i<frames[0]->getNumberOfReferenceArguments(); ++i) {
178 6 : vals.emplace_back(new Value()); vals[vals.size()-1]->setNotPeriodic();
179 : }
180 :
181 : // temporary pointes used to make the conversion once
182 :
183 2 : auto vals_ptr=Tools::unique2raw(vals);
184 : // And reparameterize
185 4 : PathReparameterization myparam( fake_pbc, vals_ptr, frames );
186 : // And make all points equally spaced
187 6 : double tol; parse("--tolerance",tol); myparam.reparameterize( fixed[0], fixed[1], tol );
188 :
189 : // Ouput data on spacings
190 : double mean=0;
191 6 : MultiValue myvpack( 1, frames[0]->getNumberOfReferenceArguments() + 3*frames[0]->getNumberOfReferencePositions() + 9 );
192 6 : ReferenceValuePack mypack( frames[0]->getNumberOfReferenceArguments(), frames[0]->getNumberOfReferencePositions(), myvpack );
193 40 : for(unsigned i=1; i<frames.size(); ++i) {
194 54 : double len = frames[i]->calc( frames[i-1]->getReferencePositions(), fake_pbc, vals_ptr, frames[i-1]->getReferenceArguments(), mypack, false );
195 : printf("FINAL DISTANCE BETWEEN FRAME %u AND %u IS %f \n",i-1,i,len );
196 18 : mean+=len;
197 : }
198 2 : printf("SUGGESTED LAMBDA PARAMETER IS THUS %f \n",2.3/mean/static_cast<double>( frames.size()-1 ) );
199 :
200 : // Delete all the frames
201 4 : OFile ofile; ofile.open(ofilename);
202 4 : std::vector<std::string> argnames; frames[0]->getArgumentRequests( argnames );
203 2 : std::vector<AtomNumber> atindices; frames[0]->getAtomRequests( atindices );
204 4 : PDB mypdb; mypdb.setAtomNumbers( atindices ); mypdb.setArgumentNames( argnames );
205 42 : for(unsigned i=0; i<frames.size(); ++i) {
206 20 : mypdb.setAtomPositions( frames[i]->getReferencePositions() );
207 80 : for(unsigned j=0; j<argnames.size(); ++j) mypdb.setArgumentValue( argnames[j], frames[i]->getReferenceArguments()[j] );
208 20 : ofile.printf("REMARK TYPE=%s\n",mtype.c_str() );
209 20 : mypdb.print( 10, NULL, ofile, ofmt );
210 : }
211 : // Delete the vals as we don't need them
212 : // for(unsigned i=0; i<vals.size(); ++i) delete vals[i];
213 : // Return as we are done
214 2 : return 0;
215 : }
216 :
217 : // Read initial frame
218 8 : std::string istart; parse("--start",istart); FILE* fp2=fopen(istart.c_str(),"r"); PDB mystartpdb;
219 2 : if( istart.length()==0 ) error("input is missing use --istart + --iend or --path");
220 2 : if( !mystartpdb.readFromFilepointer(fp2,false,0.1) ) error("could not read fila " + istart);
221 2 : auto sframe=metricRegister().create<ReferenceConfiguration>( mtype, mystartpdb );
222 2 : fclose(fp2);
223 :
224 : // Read final frame
225 8 : std::string iend; parse("--end",iend); FILE* fp1=fopen(iend.c_str(),"r"); PDB myendpdb;
226 2 : if( iend.length()==0 ) error("input is missing using --istart + --iend or --path");
227 2 : if( !myendpdb.readFromFilepointer(fp1,false,0.1) ) error("could not read fila " + iend);
228 2 : auto eframe=metricRegister().create<ReferenceConfiguration>( mtype, myendpdb );
229 2 : fclose(fp1);
230 : // Get atoms and arg requests
231 2 : std::vector<AtomNumber> atoms; std::vector<std::string> arg_names;
232 4 : sframe->getAtomRequests( atoms); eframe->getAtomRequests( atoms);
233 4 : sframe->getArgumentRequests( arg_names ); eframe->getArgumentRequests( arg_names );
234 :
235 : // Now read in the rest of the instructions
236 : unsigned nbefore, nbetween, nafter;
237 8 : parse("--nframes-before-start",nbefore); parse("--nframes",nbetween); parse("--nframes-after-end",nafter);
238 2 : nbetween++;
239 : fprintf(out,"Generating linear path connecting structure in file named %s to structure in file named %s \n",istart.c_str(),iend.c_str() );
240 : fprintf(out,"A path consisting of %u equally-spaced frames before the initial structure, %u frames between the intial and final structures "
241 2 : "and %u frames after the final structure will be created \n",nbefore,nbetween,nafter);
242 :
243 : // Create a vector of arguments to use for calculating displacements
244 4 : Pbc fpbc;
245 2 : std::vector<std::unique_ptr<Value>> args;
246 8 : for(unsigned i=0; i<eframe->getNumberOfReferenceArguments(); ++i) {
247 6 : args.emplace_back(new Value()); args[args.size()-1]->setNotPeriodic();
248 : }
249 :
250 : // convert pointer once:
251 2 : auto args_ptr=Tools::unique2raw(args);
252 :
253 : // Calculate the distance between the start and the end
254 6 : MultiValue myvpack( 1, sframe->getNumberOfReferenceArguments() + 3*sframe->getNumberOfReferencePositions() + 9);
255 6 : ReferenceValuePack mypack( sframe->getNumberOfReferenceArguments(), sframe->getNumberOfReferencePositions(), myvpack );
256 6 : sframe->calc( eframe->getReferencePositions(), fpbc, args_ptr, eframe->getReferenceArguments(), mypack, false );
257 : // And the spacing between frames
258 2 : double delr = 1.0 / static_cast<double>( nbetween );
259 : // Calculate the vector connecting the start to the end
260 10 : PDB mypdb; mypdb.setAtomNumbers( sframe->getAbsoluteIndexes() ); mypdb.addBlockEnd( sframe->getAbsoluteIndexes().size() );
261 5 : if( sframe->getArgumentNames().size()>0 ) mypdb.setArgumentNames( sframe->getArgumentNames() );
262 10 : Direction mydir(ReferenceConfigurationOptions("DIRECTION")); sframe->setupPCAStorage( mypack ); mydir.read( mypdb ); mydir.zeroDirection();
263 8 : sframe->extractDisplacementVector( eframe->getReferencePositions(), args_ptr, eframe->getReferenceArguments(), false, mydir );
264 :
265 : // Now create frames
266 4 : OFile ofile; ofile.open(ofilename); unsigned nframes=0;
267 8 : Direction pos(ReferenceConfigurationOptions("DIRECTION")); pos.read( mypdb );
268 2 : for(int i=0; i<nbefore; ++i) {
269 4 : pos.setDirection( sframe->getReferencePositions(), sframe->getReferenceArguments() );
270 2 : pos.displaceReferenceConfiguration( -i*delr, mydir );
271 2 : mypdb.setAtomPositions( pos.getReferencePositions() );
272 10 : for(unsigned j=0; j<pos.getReferenceArguments().size(); ++j) mypdb.setArgumentValue( sframe->getArgumentNames()[j], pos.getReferenceArgument(j) );
273 2 : ofile.printf("REMARK TYPE=%s\n",mtype.c_str() );
274 2 : mypdb.print( 10, NULL, ofile, ofmt ); nframes++;
275 : }
276 6 : for(unsigned i=1; i<nbetween; ++i) {
277 12 : pos.setDirection( sframe->getReferencePositions(), sframe->getReferenceArguments() );
278 6 : pos.displaceReferenceConfiguration( i*delr, mydir );
279 6 : mypdb.setAtomPositions( pos.getReferencePositions() );
280 46 : for(unsigned j=0; j<pos.getReferenceArguments().size(); ++j) mypdb.setArgumentValue( sframe->getArgumentNames()[j], pos.getReferenceArgument(j) );
281 6 : ofile.printf("REMARK TYPE=%s\n",mtype.c_str() );
282 6 : mypdb.print( 10, NULL, ofile, ofmt ); nframes++;
283 : }
284 5 : for(unsigned i=0; i<nafter; ++i) {
285 10 : pos.setDirection( eframe->getReferencePositions(), eframe->getReferenceArguments() );
286 5 : pos.displaceReferenceConfiguration( i*delr, mydir );
287 5 : mypdb.setAtomPositions( pos.getReferencePositions() );
288 37 : for(unsigned j=0; j<pos.getReferenceArguments().size(); ++j) mypdb.setArgumentValue( sframe->getArgumentNames()[j], pos.getReferenceArgument(j) );
289 5 : ofile.printf("REMARK TYPE=%s\n",mtype.c_str() );
290 5 : mypdb.print( 10, NULL, ofile, ofmt ); nframes++;
291 : }
292 :
293 2 : ofile.close(); return 0;
294 : }
295 :
296 : } // End of namespace
297 5874 : }
|
