Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2016-2020 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 "core/CLTool.h"
23 : #include "core/CLToolRegister.h"
24 : #include "tools/Tools.h"
25 : #include "tools/Pbc.h"
26 : #include "tools/PDB.h"
27 : #include "core/ActionWithValue.h"
28 : #include "core/ActionToPutData.h"
29 : #include "core/ActionSet.h"
30 : #include "core/Value.h"
31 : #include "core/PlumedMain.h"
32 : #include "Path.h"
33 : #include <cstdio>
34 : #include <string>
35 : #include <vector>
36 : #include <iostream>
37 :
38 : namespace PLMD {
39 : namespace mapping {
40 :
41 : //+PLUMEDOC TOOLS pathtools
42 : /*
43 : pathtools can be used to construct paths from pdb data
44 :
45 : The path CVs in PLUMED are curvilinear coordinates through a high dimensional vector space.
46 : Enhanced sampling calculations are often run using the progress along the paths and the distance from the path as CVs
47 : as this provides a convenient way of defining a reaction coordinate for a complicated process. This method is explained
48 : in the documentation for \ref PATH.
49 :
50 : The path itself is an ordered set of equally-spaced, high-dimensional frames the way in which these frames
51 : should be constructed will depend on the problem in hand. In other words, you will need to understand the reaction
52 : you wish to study in order to select a sensible set of frames to use in your path CV. This tool provides two
53 : methods that may be useful when it comes to constructing paths; namely:
54 :
55 : - A tool that takes in an initial guess path in which the frames are not equally spaced. This tool adjusts the positions
56 : of the frames in order to make them equally spaced so that they can be used as the basis for a path CV.
57 :
58 : - A tool that takes two frames as input and that allows you to return a linear path connecting these two frames. The
59 : output from this method may be useful as an initial guess path. It is arguable that a linear path rather defeats the
60 : purpose of the path CV method, however, as the whole purpose is to be able to define non-linear paths.
61 :
62 : Notice that you can use these two methods and take advantage of all the ways of measuring \ref dists that are available within
63 : PLUMED. The way you do this with each of these tools described above is explained in the example below.
64 :
65 : \par Examples
66 :
67 : The example below shows how you can take a set of unequally spaced frames from a pdb file named in_path.pdb
68 : and use pathtools to make them equally spaced so that they can be used as the basis for a path CV. The file
69 : containing this final path is named final_path.pdb.
70 :
71 : \verbatim
72 : plumed pathtools --path in_path.pdb --metric EUCLIDEAN --out final_path.pdb
73 : \endverbatim
74 :
75 : The example below shows how can create an initial linear path connecting the two pdb frames in start.pdb and
76 : end.pdb. In this case the path output to path.pdb will consist of 6 frames: the initial and final frames that
77 : were contained in start.pdb and end.pdb as well as four equally spaced frames along the vector connecting
78 : start.pdb to end.pdb.
79 :
80 : \verbatim
81 : plumed pathtools --start start.pdb --end end.pdb --nframes 4 --metric OPTIMAL --out path.pdb
82 : \endverbatim
83 :
84 : Often the idea with path collective variables is to create a path connecting some initial state A to some final state B. You would
85 : in this case have representative configurations from your A and B states defined in the input files to pathtools
86 : that we have called start.pdb and end.pdb in the example above. Furthermore, it may be useful to have a few frames
87 : before your start frame and after your end frame. You can use path tools to create these extended paths as shown below.
88 : In this case the final path would now consist of 8 frames. Four of these frames would lie on the vector connecting state
89 : 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
90 : frame just after end.pdb. All these frames would be equally spaced.
91 :
92 : \verbatim
93 : plumed pathtools --start start.pdb --end end.pdb --nframes 4 --metric OPTIMAL --out path.pdb --nframes-before-start 2 --nframes-after-end 2
94 : \endverbatim
95 :
96 : Notice also that when you re-parameterize paths you must choose two frames to fix. Generally you chose to fix the states
97 : that are representative of your states A and B. By default pathtools will fix the first and last frames. You can, however,
98 : change the states to fix by taking advantage of the fixed flag as shown below.
99 :
100 : \verbatim
101 : plumed pathtools --path inpath.pdb --metric EUCLIDEAN --out outpath.pdb --fixed 2,12
102 : \endverbatim
103 :
104 : */
105 : //+ENDPLUMEDOC
106 :
107 : class PathTools :
108 : public CLTool {
109 : public:
110 : static void registerKeywords( Keywords& keys );
111 : explicit PathTools(const CLToolOptions& co );
112 : int main(FILE* in, FILE*out,Communicator& pc);
113 : void printLambda( const std::string& mtype, const std::vector<std::string>& argstr, const std::string& ofile );
114 5 : std::string description()const {
115 5 : return "print out a description of the keywords for an action in html";
116 : }
117 : };
118 :
119 16263 : PLUMED_REGISTER_CLTOOL(PathTools,"pathtools")
120 :
121 5418 : void PathTools::registerKeywords( Keywords& keys ) {
122 5418 : CLTool::registerKeywords( keys );
123 5418 : keys.add("atoms","--start","a pdb file that contains the structure for the initial frame of your path");
124 5418 : keys.add("atoms","--end","a pdb file that contains the structure for the final frame of your path");
125 5418 : keys.add("atoms-1","--path","a pdb file that contains an initial path in which the frames are not equally spaced");
126 5418 : keys.add("optional","--arg","the arguments that should be read in from the pdb files");
127 5418 : keys.add("compulsory","--fixed","0","the frames to fix when constructing the path using --path");
128 5418 : keys.add("compulsory","--metric","the measure to use to calculate the distance between frames");
129 5418 : keys.add("compulsory","--out","the name of the file on which to output your path");
130 5418 : keys.add("compulsory","--arg-fmt","%f","the format to use for argument values in your frames");
131 5418 : keys.add("compulsory","--tolerance","1E-4","the tolerance to use for the algorithm that is used to re-parameterize the path");
132 5418 : keys.add("compulsory","--nframes-before-start","1","the number of frames to include in the path before the first frame");
133 5418 : keys.add("compulsory","--nframes","1","the number of frames between the start and end frames in your path");
134 5418 : keys.add("compulsory","--nframes-after-end","1","the number of frames to put after the last frame of your path");
135 5418 : }
136 :
137 9 : PathTools::PathTools(const CLToolOptions& co ):
138 9 : CLTool(co) {
139 9 : inputdata=commandline;
140 9 : }
141 :
142 4 : int PathTools::main(FILE* in, FILE*out,Communicator& pc) {
143 : // Create a PLUMED object
144 4 : PlumedMain plmd;
145 4 : int s=sizeof(double);
146 4 : plmd.cmd("setRealPrecision",&s);
147 4 : plmd.cmd("setMDEngine","pathtools");
148 : // plmd.cmd("setNoVirial");
149 4 : double tstep=1.0;
150 4 : plmd.cmd("setTimestep",&tstep);
151 4 : plmd.cmd("init");
152 4 : int step=1;
153 4 : plmd.cmd("setStep",&step);
154 :
155 : std::string mtype;
156 8 : parse("--metric",mtype);
157 : std::string ifilename;
158 8 : parse("--path",ifilename);
159 : std::string ofmt;
160 8 : parse("--arg-fmt",ofmt);
161 : std::string ofilename;
162 8 : parse("--out",ofilename);
163 : std::vector<std::string> argstr;
164 8 : parseVector("--arg",argstr);
165 4 : if( ifilename.length()>0 ) {
166 : fprintf(out,"Reparameterising path in file named %s so that all frames are equally spaced \n",ifilename.c_str() );
167 2 : FILE* fp=fopen(ifilename.c_str(),"r");
168 2 : PDB mypdb;
169 2 : bool do_read=mypdb.readFromFilepointer(fp,false,0.1);
170 2 : std::vector<double> alig( mypdb.getOccupancy() ), disp( mypdb.getBeta() );
171 2 : std::vector<AtomNumber> indices( mypdb.getAtomNumbers() );
172 2 : std::vector<unsigned> residue_indices( indices.size() );
173 15 : for(unsigned i=0; i<residue_indices.size(); ++i) {
174 13 : residue_indices[i] = mypdb.getResidueNumber( indices[i] );
175 : }
176 : // Count the number of frames in the input file
177 : unsigned nfram=1;
178 20 : while ( do_read ) {
179 20 : if( !mypdb.readFromFilepointer(fp,false,0.1) ) {
180 : break;
181 : }
182 18 : nfram++;
183 : }
184 2 : if( argstr.size()>0 ) {
185 3 : for(unsigned i=0; i<argstr.size(); ++i ) {
186 4 : std::string input = argstr[i] + ": PDB2CONSTANT NOARGS REFERENCE=" + ifilename + " ARG=" + argstr[i];
187 : const char* inpt = input.c_str();
188 2 : plmd.cmd("readInputLine", inpt);
189 : }
190 : } else {
191 1 : std::string input = "data: PDB2CONSTANT REFERENCE=" + ifilename;
192 : const char* inpt = input.c_str();
193 1 : plmd.cmd("readInputLine", inpt );
194 : }
195 2 : std::string reparam_str = "REPARAMETERIZE_PATH REFERENCE=";
196 2 : if( argstr.size()>0 ) {
197 : reparam_str += argstr[0];
198 3 : for(unsigned i=0; i<argstr.size(); ++i ) {
199 4 : reparam_str += "," + argstr[i];
200 : }
201 : } else {
202 : reparam_str += "data";
203 : }
204 : std::vector<unsigned> fixed;
205 4 : parseVector("--fixed",fixed);
206 2 : if( fixed.size()==1 ) {
207 1 : if( fixed[0]!=0 ) {
208 0 : error("input to --fixed should be two integers");
209 : }
210 1 : fixed.resize(2);
211 1 : fixed[0]=1;
212 1 : fixed[1]=nfram;
213 1 : } else if( fixed.size()==2 ) {
214 1 : if( fixed[0]<1 || fixed[1]<1 || fixed[0]>nfram || fixed[1]>nfram ) {
215 0 : error("input to --fixed should be two numbers between 0 and the number of frames-1");
216 : }
217 : } else {
218 0 : error("input to --fixed should be two integers");
219 : }
220 : std::string fix1, fix2;
221 2 : Tools::convert( fixed[0], fix1 );
222 2 : Tools::convert( fixed[1], fix2 );
223 4 : reparam_str += " FIXED=" + fix1 + "," + fix2;
224 : std::string tol;
225 2 : parse("--tolerance",tol);
226 4 : reparam_str += " TOL=" + tol;
227 : // Now create the metric object
228 : reparam_str += " METRIC=";
229 5 : if( mtype=="OPTIMAL-FAST" || mtype=="OPTIMAL" || mtype=="SIMPLE" ) {
230 2 : reparam_str += "{RMSD_VECTOR DISPLACEMENT SQUARED UNORMALIZED TYPE=" + mtype;
231 : // Get the align values
232 : std::string anum;
233 1 : Tools::convert( alig[0], anum );
234 1 : reparam_str += " ALIGN=" + anum;
235 13 : for(unsigned i=1; i<alig.size(); ++i) {
236 12 : Tools::convert( alig[i], anum );
237 24 : reparam_str += "," + anum;
238 : }
239 : // Get the displace values
240 : std::string dnum;
241 1 : Tools::convert( disp[0], dnum );
242 1 : reparam_str += " DISPLACE=" + dnum;
243 13 : for(unsigned i=1; i<disp.size(); ++i) {
244 12 : Tools::convert( disp[i], dnum );
245 24 : reparam_str += "," + dnum;
246 : }
247 : reparam_str += "} METRIC_COMPONENT=disp";
248 1 : } else if( mtype=="EUCLIDEAN" ) {
249 : reparam_str += "DIFFERENCE";
250 : } else {
251 : // Add functionality to read plumed input here
252 0 : plumed_merror("metric type " + mtype + " has not been implemented");
253 : }
254 :
255 : // Now do the reparameterization
256 : const char* icinp= reparam_str.c_str();
257 2 : plmd.cmd("readInputLine",icinp);
258 : Action* raction = plmd.getActionSet()[plmd.getActionSet().size()-1].get();
259 2 : raction->update();
260 :
261 : // And print the final reference configurations
262 4 : std::string pinput="DUMPPDB STRIDE=1 DESCRIPTION=PATH FILE=" + ofilename + " FMT=" + ofmt;
263 2 : if( argstr.size()>0 ) {
264 1 : pinput += " ARG=" + argstr[0];
265 2 : for(unsigned i=1; i<argstr.size(); ++i ) {
266 2 : pinput += "," + argstr[i];
267 : }
268 : } else {
269 : std::string num;
270 1 : Tools::convert( indices[0].serial(), num );
271 1 : pinput += " ATOMS=data ATOM_INDICES=" + num;
272 13 : for(unsigned i=1; i<indices.size(); ++i ) {
273 12 : Tools::convert( indices[i].serial(), num );
274 24 : pinput += "," + num;
275 : }
276 1 : Tools::convert( residue_indices[0], num );
277 1 : pinput += " RESIDUE_INDICES=" + num;
278 13 : for(unsigned i=1; i<residue_indices.size(); ++i ) {
279 12 : Tools::convert( residue_indices[i], num );
280 24 : pinput += "," + num;
281 : }
282 : std::string anum, dnum;
283 1 : Tools::convert( alig[0], anum );
284 1 : Tools::convert( disp[0], dnum );
285 1 : pinput += " OCCUPANCY=" + anum;
286 13 : for(unsigned i=1; i<alig.size(); ++i) {
287 12 : Tools::convert( alig[i], anum );
288 24 : pinput += "," + anum;
289 : }
290 1 : pinput += " BETA=" + dnum;
291 13 : for(unsigned i=1; i<disp.size(); ++i) {
292 12 : Tools::convert( disp[i], dnum );
293 24 : pinput += "," + dnum;
294 : }
295 : }
296 : const char* pcinp=pinput.c_str();
297 2 : plmd.cmd("readInputLine",pcinp);
298 : Action* paction = plmd.getActionSet()[plmd.getActionSet().size()-1].get();
299 2 : paction->update();
300 :
301 : // Ouput data on spacings
302 2 : printLambda( mtype, argstr, ofilename );
303 : return 0;
304 2 : }
305 4 : for(unsigned i=0; i<argstr.size(); ++i) {
306 2 : std::string input = argstr[i] + ": CONSTANT VALUE=1";
307 : const char* inpt = input.c_str();
308 2 : plmd.cmd("readInputLine", inpt);
309 : }
310 :
311 : // Read in the instructions
312 : unsigned nbefore, nbetween, nafter;
313 : std::string istart;
314 4 : parse("--start",istart);
315 : std::string iend;
316 2 : parse("--end",iend);
317 2 : parse("--nframes-before-start",nbefore);
318 2 : parse("--nframes",nbetween);
319 2 : parse("--nframes-after-end",nafter);
320 2 : nbetween++;
321 : fprintf(out,"Generating linear path connecting structure in file named %s to structure in file named %s \n",istart.c_str(),iend.c_str() );
322 2 : fprintf(out,"A path consisting of %u equally-spaced frames before the initial structure, %u frames between the initial and final structures "
323 : "and %u frames after the final structure will be created \n",nbefore,nbetween,nafter);
324 :
325 2 : std::vector<double> start_pos( argstr.size() ), end_pos( argstr.size() );
326 : std::vector<AtomNumber> indices;
327 2 : if( argstr.size()>0 ) {
328 3 : for(unsigned i=0; i<argstr.size(); ++i ) {
329 4 : std::string input = argstr[i] + "_start: PDB2CONSTANT REFERENCE=" + istart + " ARG=" + argstr[i];
330 : const char* inpt = input.c_str();
331 2 : plmd.cmd("readInputLine", inpt );
332 : long srank;
333 4 : plmd.cmd("getDataRank " + argstr[i] + "_start", &srank );
334 2 : if( srank!=1 ) {
335 0 : error("should only be one frame in input pdb");
336 : }
337 2 : std::vector<long> sshape( 1 );
338 4 : plmd.cmd("getDataShape " + argstr[i] + "_start", &sshape[0] );
339 2 : if( sshape[0]!=1 ) {
340 0 : error("should only be one frame in input pdb");
341 : }
342 4 : plmd.cmd("setMemoryForData " + argstr[i] + "_start", &start_pos[i] );
343 4 : std::string input2 = argstr[i] + "_end: PDB2CONSTANT REFERENCE=" + iend + " ARG=" + argstr[i];
344 : const char* inpt2 = input2.c_str();
345 2 : plmd.cmd("readInputLine", inpt2 );
346 : long erank;
347 4 : plmd.cmd("getDataRank " + argstr[i] + "_end", &erank );
348 2 : if( erank!=1 ) {
349 0 : error("should only be one frame in input pdb");
350 : }
351 2 : std::vector<long> eshape( 1 );
352 4 : plmd.cmd("getDataShape " + argstr[i] + "_end", &eshape[0] );
353 2 : if( eshape[0]!=1 ) {
354 0 : error("should only be one frame in input pdb");
355 : }
356 4 : plmd.cmd("setMemoryForData " + argstr[i] + "_end", &end_pos[i] );
357 : }
358 : } else {
359 1 : std::string input = "start: PDB2CONSTANT REFERENCE=" + istart;
360 : const char* inpt = input.c_str();
361 1 : plmd.cmd("readInputLine", inpt );
362 1 : FILE* fp=fopen(istart.c_str(),"r");
363 1 : PDB mypdb;
364 1 : bool do_read=mypdb.readFromFilepointer(fp,false,0.1);
365 1 : indices.resize( mypdb.getAtomNumbers().size() );
366 14 : for(unsigned i=0; i<indices.size(); ++i) {
367 13 : indices[i] = mypdb.getAtomNumbers()[i];
368 : }
369 : long srank;
370 1 : plmd.cmd("getDataRank start", &srank );
371 1 : if( srank!=2 ) {
372 0 : error("should only be one frame in input pdb:" + istart);
373 : }
374 1 : std::vector<long> sshape( srank );
375 1 : plmd.cmd("getDataShape start", &sshape[0] );
376 1 : start_pos.resize( sshape[0]*sshape[1] );
377 1 : if( sshape[0]!=1 ) {
378 0 : error("should only be one frame in input pdb: " + istart);
379 : }
380 1 : plmd.cmd("setMemoryForData start", &start_pos[0] );
381 1 : std::string input2 = "end: PDB2CONSTANT REFERENCE=" + iend;
382 : const char* inpt2 = input2.c_str();
383 1 : plmd.cmd("readInputLine", inpt2 );
384 : long erank;
385 1 : plmd.cmd("getDataRank end", &erank );
386 1 : if( erank!=2 ) {
387 0 : error("should only be one frame in input pdb: " + iend);
388 : }
389 1 : std::vector<long> eshape( erank );
390 1 : plmd.cmd("getDataShape end", &eshape[0] );
391 1 : end_pos.resize( eshape[0]*eshape[1] );
392 1 : if( eshape[0]!=1 ) {
393 0 : error("should only be one frame in input pdb: " + iend );
394 : }
395 1 : plmd.cmd("setMemoryForData end", &end_pos[0] );
396 1 : }
397 :
398 : // Now create the metric object
399 : std::vector<double> alig, disp;
400 : std::vector<unsigned> residue_indices;
401 5 : if( mtype=="OPTIMAL-FAST" || mtype=="OPTIMAL" || mtype=="SIMPLE" ) {
402 1 : std::string trinput = "endT: TRANSPOSE ARG=end";
403 : const char* tinp=trinput.c_str();
404 1 : plmd.cmd("readInputLine",tinp);
405 1 : PDB pdb;
406 1 : pdb.read(istart,false,0.1);
407 1 : alig.resize( pdb.getOccupancy().size() );
408 1 : disp.resize( pdb.getBeta().size() );
409 14 : for(unsigned i=0; i<alig.size(); ++i) {
410 13 : alig[i] = pdb.getOccupancy()[i];
411 13 : disp[i]= pdb.getBeta()[i];
412 : }
413 2 : std::string minput = "d: RMSD_VECTOR DISPLACEMENT SQUARED UNORMALIZED TYPE=" + mtype + " ARG=endT,start";
414 : // Get the align values
415 : std::string anum;
416 1 : Tools::convert( alig[0], anum );
417 1 : minput += " ALIGN=" + anum;
418 13 : for(unsigned i=1; i<alig.size(); ++i) {
419 12 : Tools::convert( alig[i], anum );
420 24 : minput += "," + anum;
421 : }
422 : // Get the displace values
423 : std::string dnum;
424 1 : Tools::convert( disp[0], dnum );
425 1 : minput += " DISPLACE=" + dnum;
426 13 : for(unsigned i=1; i<disp.size(); ++i) {
427 12 : Tools::convert( disp[i], dnum );
428 24 : minput += "," + dnum;
429 : }
430 : const char* mcinp=minput.c_str();
431 1 : plmd.cmd("readInputLine",mcinp);
432 1 : std::string tinput = "CUSTOM ARG=d.disp FUNC=x PERIODIC=NO";
433 : const char* tcinp=tinput.c_str();
434 1 : plmd.cmd("readInputLine",tcinp);
435 : // Get the residue numbers
436 1 : residue_indices.resize( pdb.size() );
437 14 : for(unsigned i=0; i<residue_indices.size(); ++i) {
438 13 : residue_indices[i] = pdb.getResidueNumber( indices[i] );
439 : }
440 2 : } else if( mtype=="EUCLIDEAN" ) {
441 2 : std::string end_args="ARG1=" + argstr[0] + "_end", start_args="ARG2=" + argstr[0] + "_start";
442 2 : for(unsigned i=1; i<argstr.size(); ++i ) {
443 2 : end_args += "," + argstr[i] + "_end";
444 2 : start_args += "," + argstr[i] + "_start";
445 : }
446 2 : std::string minput = "d: DISPLACEMENT " + end_args + " " + start_args;
447 : const char* mcinp=minput.c_str();
448 1 : plmd.cmd("readInputLine",mcinp);
449 1 : std::string tinput = "TRANSPOSE ARG=d";
450 : const char* tcinp=tinput.c_str();
451 1 : plmd.cmd("readInputLine",tcinp);
452 : } else {
453 : // Add functionality to read plumed input here
454 0 : plumed_merror("metric type " + mtype + " has not been implemented");
455 : }
456 :
457 : // Retrieve final displacement vector
458 2 : unsigned aind = plmd.getActionSet().size()-1;
459 : while( true ) {
460 3 : const ActionShortcut* as=dynamic_cast<const ActionShortcut*>( plmd.getActionSet()[aind].get() );
461 3 : if( !as ) {
462 : break ;
463 : }
464 1 : aind = aind - 1;
465 : plumed_assert( aind>=0 );
466 1 : }
467 2 : ActionWithValue* av = dynamic_cast<ActionWithValue*>( plmd.getActionSet()[aind].get() );
468 2 : if( !av ) {
469 0 : error("invalid input for metric" );
470 : }
471 2 : if( av->getNumberOfComponents()!=1 && av->getName()!="RMSD_VECTOR" ) {
472 0 : error("cannot use multi component actions as metric");
473 : }
474 2 : std::string mydisp = av->copyOutput(0)->getName();
475 : // Now add calls so we can grab the data from plumed
476 : long rank;
477 4 : plmd.cmd("getDataRank " + mydisp, &rank );
478 2 : if( rank!=1 ) {
479 0 : error("displacement must be a vector quantity");
480 : }
481 2 : std::vector<long> ishape( rank );
482 4 : plmd.cmd("getDataShape " + mydisp, &ishape[0] );
483 2 : std::vector<double> displacement( ishape[0] );
484 4 : plmd.cmd("setMemoryForData " + mydisp, &displacement[0] );
485 : // And calculate the displacement
486 2 : plmd.cmd("calc");
487 :
488 : // Now read in the initial frame
489 :
490 : // Now create frames
491 2 : double delr = 1.0 / static_cast<double>( nbetween );
492 : std::vector<std::vector<double> > allframes;
493 2 : std::vector<double> frame( displacement.size() );
494 4 : for(int i=0; i<nbefore; ++i) {
495 43 : for(unsigned j=0; j<frame.size(); ++j) {
496 41 : frame[j] = start_pos[j] - i*delr*displacement[j];
497 : }
498 2 : allframes.push_back( frame );
499 : }
500 8 : for(unsigned i=1; i<nbetween; ++i) {
501 55 : for(unsigned j=0; j<frame.size(); ++j) {
502 49 : frame[j] = start_pos[j] + i*delr*displacement[j];
503 : }
504 6 : allframes.push_back( frame );
505 : }
506 7 : for(unsigned i=0; i<nafter; ++i) {
507 52 : for(unsigned j=0; j<frame.size(); ++j) {
508 47 : frame[j] = end_pos[j] + i*delr*displacement[j];
509 : }
510 5 : allframes.push_back( frame );
511 : }
512 : // This prints out our final reference configurations
513 2 : plmd.cmd("clear");
514 4 : plmd.readInputLine("timestep: PUT UNIT=time PERIODIC=NO CONSTANT");
515 2 : ActionToPutData* ts = plmd.getActionSet().selectWithLabel<ActionToPutData*>("timestep");
516 4 : ts->setValuePointer("timestep", &tstep );
517 4 : std::string pinput="DUMPPDB STRIDE=1 DESCRIPTION=PATH FILE=" + ofilename + " FMT=" + ofmt;
518 2 : if( argstr.size()>0 ) {
519 3 : for(unsigned i=0; i<argstr.size(); ++i) {
520 : std::string rnum;
521 2 : Tools::convert( allframes[0][i], rnum );
522 2 : std::string inpt = argstr[i] + ": CONSTANT VALUES=" + rnum;
523 20 : for(unsigned j=1; j<allframes.size(); ++j) {
524 18 : Tools::convert( allframes[j][i], rnum );
525 36 : inpt += "," + rnum;
526 : }
527 : const char* icinp=inpt.c_str();
528 2 : plmd.cmd("readInputLine",icinp);
529 2 : if( i==0 ) {
530 2 : pinput += " ARG=" + argstr[i];
531 : } else {
532 2 : pinput += "," + argstr[i];
533 : }
534 : }
535 : } else {
536 : std::string nc;
537 1 : Tools::convert( frame.size(), nc );
538 : std::string nr;
539 1 : Tools::convert( allframes.size(), nr );
540 : std::string rnum;
541 1 : Tools::convert( allframes[0][0], rnum );
542 2 : std::string inpt = "atom_data: CONSTANT NROWS=" + nr + " NCOLS=" + nc + " VALUES=" + rnum;
543 4 : for(unsigned i=0; i<allframes.size(); ++i) {
544 120 : for(unsigned j=0; j<allframes[i].size(); ++j) {
545 117 : if( i==0 && j==0 ) {
546 1 : continue;
547 : }
548 116 : Tools::convert( allframes[i][j], rnum );
549 232 : inpt += "," + rnum;
550 : }
551 : }
552 : const char* icinp=inpt.c_str();
553 1 : plmd.cmd("readInputLine",icinp);
554 : std::string num;
555 1 : Tools::convert( indices[0].serial(), num );
556 1 : pinput += " ATOMS=atom_data ATOM_INDICES=" + num;
557 13 : for(unsigned i=1; i<indices.size(); ++i ) {
558 12 : Tools::convert( indices[i].serial(), num );
559 24 : pinput += "," + num;
560 : }
561 1 : Tools::convert( residue_indices[0], num );
562 1 : pinput += " RESIDUE_INDICES=" + num;
563 13 : for(unsigned i=1; i<residue_indices.size(); ++i ) {
564 12 : Tools::convert( residue_indices[i], num );
565 24 : pinput += "," + num;
566 : }
567 : std::string anum, dnum;
568 1 : Tools::convert( alig[0], anum );
569 1 : Tools::convert( disp[0], dnum );
570 1 : pinput += " OCCUPANCY=" + anum;
571 13 : for(unsigned i=1; i<alig.size(); ++i) {
572 12 : Tools::convert( alig[i], anum );
573 24 : pinput += "," + anum;
574 : }
575 1 : pinput += " BETA=" + dnum;
576 13 : for(unsigned i=1; i<disp.size(); ++i) {
577 12 : Tools::convert( disp[i], dnum );
578 24 : pinput += "," + dnum;
579 : }
580 : }
581 : const char* pcinp=pinput.c_str();
582 2 : plmd.cmd("readInputLine",pcinp);
583 : Action* paction = plmd.getActionSet()[plmd.getActionSet().size()-1].get();
584 2 : paction->update();
585 : // And output suggestions on the value of Lambda
586 2 : printLambda( mtype, argstr, ofilename );
587 : return 0;
588 10 : }
589 :
590 4 : void PathTools::printLambda( const std::string& mtype, const std::vector<std::string>& argstr, const std::string& ofile ) {
591 : // Create a PLUMED object
592 4 : PlumedMain plmd;
593 4 : int s=sizeof(double);
594 4 : plmd.cmd("setRealPrecision",&s);
595 4 : plmd.cmd("setMDEngine","pathtools");
596 4 : double tstep=1.0;
597 4 : plmd.cmd("setTimestep",&tstep);
598 4 : plmd.cmd("init");
599 4 : int step=1;
600 4 : plmd.cmd("setStep",&step);
601 :
602 4 : FILE* fp=fopen(ofile.c_str(),"r");
603 : bool do_read=true;
604 : unsigned nfram=0;
605 : std::vector<double> alig, disp;
606 : // Read in the argument names
607 8 : for(unsigned i=0; i<argstr.size(); ++i ) {
608 4 : std::string input2 = argstr[i] + ": CONSTANT VALUE=1";
609 : const char* inpt2 = input2.c_str();
610 4 : plmd.cmd("readInputLine", inpt2);
611 : }
612 37 : while (do_read ) {
613 37 : PDB mypdb;
614 : // Read the pdb file
615 37 : do_read=mypdb.readFromFilepointer(fp,false,0.1);
616 37 : if( !do_read ) {
617 : break;
618 : }
619 : std::string num;
620 33 : Tools::convert( nfram+1, num );
621 : nfram++;
622 : std::string iinput;
623 33 : if( argstr.size()>0 ) {
624 60 : for(unsigned i=0; i<argstr.size(); ++i ) {
625 80 : std::string input = "ref_" + num + "_" + argstr[i] + ": PDB2CONSTANT REFERENCE=" + ofile + " ARG=" + argstr[i] + " NUMBER=" + num;
626 : const char* inpt = input.c_str();
627 40 : plmd.cmd("readInputLine", inpt );
628 : }
629 : } else {
630 26 : std::string input = "ref_" + num + ": PDB2CONSTANT REFERENCE=" + ofile + " NUMBER=" + num;
631 : const char* inpt = input.c_str();
632 13 : plmd.cmd("readInputLine", inpt );
633 : }
634 :
635 33 : if( nfram==1 ) {
636 4 : alig.resize( mypdb.getOccupancy().size() );
637 30 : for(unsigned i=0; i<alig.size(); ++i) {
638 26 : alig[i]=mypdb.getOccupancy()[i];
639 : }
640 4 : disp.resize( mypdb.getBeta().size() );
641 30 : for(unsigned i=0; i<disp.size(); ++i) {
642 26 : disp[i]=mypdb.getBeta()[i];
643 : }
644 : }
645 37 : }
646 : // Now create the objects to measure the distances between the frames
647 4 : std::vector<double> data( nfram );
648 33 : for(unsigned j=1; j<nfram; ++j) {
649 : std::string istr, jstr;
650 29 : Tools::convert( j, istr);
651 29 : Tools::convert( j+1, jstr );
652 76 : if( mtype=="OPTIMAL-FAST" || mtype=="OPTIMAL" || mtype=="SIMPLE" ) {
653 22 : std::string strv = "ref_" + jstr + "T: TRANSPOSE ARG=ref_" + jstr;
654 : const char* cstrv = strv.c_str();
655 11 : plmd.cmd("readInputLine",cstrv);
656 22 : std::string dstring = "d" + istr + ": RMSD_VECTOR TYPE=" + mtype + " ARG=ref_" + jstr + "T,ref_" + istr;
657 : // Get the align values
658 : std::string anum;
659 11 : Tools::convert( alig[0], anum );
660 11 : dstring += " ALIGN=" + anum;
661 143 : for(unsigned i=1; i<alig.size(); ++i) {
662 132 : Tools::convert( alig[i], anum );
663 264 : dstring += "," + anum;
664 : }
665 : // Get the displace values
666 : std::string dnum;
667 11 : Tools::convert( disp[0], dnum );
668 11 : dstring += " DISPLACE=" + dnum;
669 143 : for(unsigned i=1; i<disp.size(); ++i) {
670 132 : Tools::convert( disp[i], dnum );
671 264 : dstring += "," + dnum;
672 : }
673 : const char* icinp=dstring.c_str();
674 11 : plmd.cmd("readInputLine",icinp);
675 18 : } else if( mtype=="EUCLIDEAN" ) {
676 54 : std::string end_args="ARG1=ref_" + istr + "_" + argstr[0], start_args="ARG2=ref_" + jstr + "_" + argstr[0];
677 36 : for(unsigned i=1; i<argstr.size(); ++i ) {
678 36 : end_args += ",ref_" + istr + "_" + argstr[i];
679 36 : start_args += ",ref_" + jstr + "_" + argstr[i];
680 : }
681 36 : std::string fstr = "d" + istr + ": EUCLIDEAN_DISTANCE " + end_args + " " + start_args;
682 : const char* icinp=fstr.c_str();
683 18 : plmd.cmd("readInputLine",icinp);
684 : } else {
685 0 : plumed_merror("metric type " + mtype + " has not been implemented");
686 : }
687 : long rank;
688 58 : plmd.cmd("getDataRank d" + istr, &rank );
689 29 : if( rank!=1 ) {
690 0 : error("distance should be of rank 1");
691 : }
692 29 : std::vector<long> ishape(1);
693 58 : plmd.cmd("getDataShape d" + istr, &ishape[0] );
694 29 : if( ishape[0]!=1 ) {
695 0 : error("distance should be of rank 1");
696 : }
697 58 : plmd.cmd("setMemoryForData d" + istr, &data[j-1] );
698 : }
699 4 : plmd.cmd("calc");
700 : double mean=0;
701 : printf("N.B. THIS CODE ALWAYS AIMS TO CREATE EQUALLY SPACED FRAMES \n");
702 : printf("THERE MAY BE SMALL DESCREPENCIES IN THE NUMBERS BELOW, HOWEVER, BECAUSE OF ROUNDING ERRORS \n");
703 33 : for(unsigned i=1; i<nfram; ++i) {
704 29 : printf("FINAL DISTANCE BETWEEN FRAME %u AND %u IS %f \n",i-1,i,data[i-1]);
705 29 : mean += data[i-1];
706 : }
707 4 : printf("SUGGESTED LAMBDA PARAMETER IS THUS %f \n",2.3/mean/static_cast<double>( nfram-1 ) );
708 4 : }
709 :
710 : } // End of namespace
711 : }
|
