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