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 : {
110 : public:
111 : static void registerKeywords( Keywords& keys );
112 : explicit PathTools(const CLToolOptions& co );
113 : int main(FILE* in, FILE*out,Communicator& pc);
114 : void printLambda( const std::string& mtype, const std::vector<std::string>& argstr, const std::string& ofile );
115 4 : std::string description()const {
116 4 : return "print out a description of the keywords for an action in html";
117 : }
118 : };
119 :
120 15956 : PLUMED_REGISTER_CLTOOL(PathTools,"pathtools")
121 :
122 5316 : void PathTools::registerKeywords( Keywords& keys ) {
123 5316 : CLTool::registerKeywords( keys );
124 10632 : keys.add("atoms","--start","a pdb file that contains the structure for the initial frame of your path");
125 10632 : keys.add("atoms","--end","a pdb file that contains the structure for the final frame of your path");
126 10632 : keys.add("atoms-1","--path","a pdb file that contains an initial path in which the frames are not equally spaced");
127 10632 : keys.add("optional","--arg","the arguments that should be read in from the pdb files");
128 10632 : keys.add("compulsory","--fixed","0","the frames to fix when constructing the path using --path");
129 10632 : keys.add("compulsory","--metric","the measure to use to calculate the distance between frames");
130 10632 : keys.add("compulsory","--out","the name of the file on which to output your path");
131 10632 : keys.add("compulsory","--arg-fmt","%f","the format to use for argument values in your frames");
132 10632 : keys.add("compulsory","--tolerance","1E-4","the tolerance to use for the algorithm that is used to re-parameterize the path");
133 10632 : keys.add("compulsory","--nframes-before-start","1","the number of frames to include in the path before the first frame");
134 10632 : keys.add("compulsory","--nframes","1","the number of frames between the start and end frames in your path");
135 10632 : keys.add("compulsory","--nframes-after-end","1","the number of frames to put after the last frame of your path");
136 5316 : }
137 :
138 8 : PathTools::PathTools(const CLToolOptions& co ):
139 8 : CLTool(co)
140 : {
141 8 : inputdata=commandline;
142 8 : }
143 :
144 4 : int PathTools::main(FILE* in, FILE*out,Communicator& pc) {
145 : // Create a PLUMED object
146 4 : PlumedMain plmd; int s=sizeof(double);
147 4 : plmd.cmd("setRealPrecision",&s);
148 4 : plmd.cmd("setMDEngine","pathtools");
149 : // plmd.cmd("setNoVirial");
150 4 : double tstep=1.0; plmd.cmd("setTimestep",&tstep);
151 4 : plmd.cmd("init");
152 4 : int step=1; plmd.cmd("setStep",&step);
153 :
154 8 : std::string mtype; parse("--metric",mtype);
155 8 : std::string ifilename; parse("--path",ifilename);
156 8 : std::string ofmt; parse("--arg-fmt",ofmt);
157 8 : std::string ofilename; parse("--out",ofilename);
158 8 : std::vector<std::string> argstr; parseVector("--arg",argstr);
159 4 : if( ifilename.length()>0 ) {
160 : fprintf(out,"Reparameterising path in file named %s so that all frames are equally spaced \n",ifilename.c_str() );
161 2 : FILE* fp=fopen(ifilename.c_str(),"r"); PDB mypdb; bool do_read=mypdb.readFromFilepointer(fp,false,0.1);
162 2 : std::vector<double> alig( mypdb.getOccupancy() ), disp( mypdb.getBeta() ); std::vector<AtomNumber> indices( mypdb.getAtomNumbers() );
163 15 : std::vector<unsigned> residue_indices( indices.size() ); for(unsigned i=0; i<residue_indices.size(); ++i) residue_indices[i] = mypdb.getResidueNumber( indices[i] );
164 : // Count the number of frames in the input file
165 20 : unsigned nfram=1; while ( do_read ) { if( !mypdb.readFromFilepointer(fp,false,0.1) ) break; nfram++; }
166 2 : if( argstr.size()>0 ) {
167 3 : for(unsigned i=0; i<argstr.size(); ++i ) {
168 4 : std::string input = argstr[i] + ": PDB2CONSTANT NOARGS REFERENCE=" + ifilename + " ARG=" + argstr[i];
169 2 : const char* inpt = input.c_str(); plmd.cmd("readInputLine", inpt);
170 : }
171 : } else {
172 1 : std::string input = "data: PDB2CONSTANT REFERENCE=" + ifilename;
173 1 : const char* inpt = input.c_str(); plmd.cmd("readInputLine", inpt );
174 : }
175 2 : std::string reparam_str = "REPARAMETERIZE_PATH REFERENCE=";
176 5 : if( argstr.size()>0 ) { reparam_str += argstr[0]; for(unsigned i=0; i<argstr.size(); ++i ) reparam_str += "," + argstr[i]; }
177 4 : else { reparam_str += "data"; } std::vector<unsigned> fixed; parseVector("--fixed",fixed);
178 2 : if( fixed.size()==1 ) {
179 1 : if( fixed[0]!=0 ) error("input to --fixed should be two integers");
180 1 : fixed.resize(2); fixed[0]=1; fixed[1]=nfram;
181 1 : } else if( fixed.size()==2 ) {
182 1 : if( fixed[0]<1 || fixed[1]<1 || fixed[0]>nfram || fixed[1]>nfram ) {
183 0 : error("input to --fixed should be two numbers between 0 and the number of frames-1");
184 : }
185 0 : } else error("input to --fixed should be two integers");
186 2 : std::string fix1, fix2; Tools::convert( fixed[0], fix1 ); Tools::convert( fixed[1], fix2 );
187 4 : reparam_str += " FIXED=" + fix1 + "," + fix2;
188 6 : std::string tol; parse("--tolerance",tol); reparam_str += " TOL=" + tol;
189 : // Now create the metric object
190 : reparam_str += " METRIC=";
191 5 : if( mtype=="OPTIMAL-FAST" || mtype=="OPTIMAL" || mtype=="SIMPLE" ) {
192 2 : reparam_str += "{RMSD_VECTOR DISPLACEMENT SQUARED UNORMALIZED TYPE=" + mtype;
193 : // Get the align values
194 1 : std::string anum; Tools::convert( alig[0], anum ); reparam_str += " ALIGN=" + anum;
195 13 : for(unsigned i=1; i<alig.size(); ++i) { Tools::convert( alig[i], anum ); reparam_str += "," + anum; }
196 : // Get the displace values
197 1 : std::string dnum; Tools::convert( disp[0], dnum ); reparam_str += " DISPLACE=" + dnum;
198 13 : for(unsigned i=1; i<disp.size(); ++i) { Tools::convert( disp[i], dnum ); reparam_str += "," + dnum; }
199 : reparam_str += "} METRIC_COMPONENT=disp";
200 1 : } else if( mtype=="EUCLIDEAN" ) {
201 : reparam_str += "DIFFERENCE";
202 : } else {
203 : // Add functionality to read plumed input here
204 0 : plumed_merror("metric type " + mtype + " has not been implemented");
205 : }
206 :
207 : // Now do the reparameterization
208 2 : const char* icinp= reparam_str.c_str(); plmd.cmd("readInputLine",icinp);
209 : Action* raction = plmd.getActionSet()[plmd.getActionSet().size()-1].get();
210 2 : raction->update();
211 :
212 : // And print the final reference configurations
213 4 : std::string pinput="DUMPPDB STRIDE=1 DESCRIPTION=PATH FILE=" + ofilename + " FMT=" + ofmt;
214 2 : if( argstr.size()>0 ) {
215 3 : pinput += " ARG=" + argstr[0]; for(unsigned i=1; i<argstr.size(); ++i ) pinput += "," + argstr[i];
216 : } else {
217 14 : std::string num; Tools::convert( indices[0].serial(), num ); pinput += " ATOMS=data ATOM_INDICES=" + num; for(unsigned i=1; i<indices.size(); ++i ) { Tools::convert( indices[i].serial(), num ); pinput += "," + num; }
218 14 : Tools::convert( residue_indices[0], num ); pinput += " RESIDUE_INDICES=" + num; for(unsigned i=1; i<residue_indices.size(); ++i ) { Tools::convert( residue_indices[i], num ); pinput += "," + num; }
219 1 : std::string anum, dnum; Tools::convert( alig[0], anum ); Tools::convert( disp[0], dnum );
220 14 : pinput += " OCCUPANCY=" + anum; for(unsigned i=1; i<alig.size(); ++i) { Tools::convert( alig[i], anum ); pinput += "," + anum; }
221 14 : pinput += " BETA=" + dnum; for(unsigned i=1; i<disp.size(); ++i) { Tools::convert( disp[i], dnum ); pinput += "," + dnum; }
222 : }
223 2 : const char* pcinp=pinput.c_str(); plmd.cmd("readInputLine",pcinp);
224 : Action* paction = plmd.getActionSet()[plmd.getActionSet().size()-1].get();
225 2 : paction->update();
226 :
227 : // Ouput data on spacings
228 2 : printLambda( mtype, argstr, ofilename );
229 : return 0;
230 2 : }
231 4 : for(unsigned i=0; i<argstr.size(); ++i) {
232 2 : std::string input = argstr[i] + ": CONSTANT VALUE=1";
233 2 : const char* inpt = input.c_str(); plmd.cmd("readInputLine", inpt);
234 : }
235 :
236 : // Read in the instructions
237 : unsigned nbefore, nbetween, nafter;
238 4 : std::string istart; parse("--start",istart); std::string iend; parse("--end",iend);
239 6 : parse("--nframes-before-start",nbefore); parse("--nframes",nbetween); parse("--nframes-after-end",nafter);
240 2 : nbetween++;
241 : fprintf(out,"Generating linear path connecting structure in file named %s to structure in file named %s \n",istart.c_str(),iend.c_str() );
242 2 : fprintf(out,"A path consisting of %u equally-spaced frames before the initial structure, %u frames between the initial and final structures "
243 : "and %u frames after the final structure will be created \n",nbefore,nbetween,nafter);
244 :
245 2 : std::vector<double> start_pos( argstr.size() ), end_pos( argstr.size() ); std::vector<AtomNumber> indices;
246 2 : if( argstr.size()>0 ) {
247 3 : for(unsigned i=0; i<argstr.size(); ++i ) {
248 4 : std::string input = argstr[i] + "_start: PDB2CONSTANT REFERENCE=" + istart + " ARG=" + argstr[i];
249 2 : const char* inpt = input.c_str(); plmd.cmd("readInputLine", inpt );
250 4 : long srank; plmd.cmd("getDataRank " + argstr[i] + "_start", &srank );
251 2 : if( srank!=1 ) error("should only be one frame in input pdb");
252 4 : std::vector<long> sshape( 1 ); plmd.cmd("getDataShape " + argstr[i] + "_start", &sshape[0] );
253 2 : if( sshape[0]!=1 ) error("should only be one frame in input pdb");
254 4 : plmd.cmd("setMemoryForData " + argstr[i] + "_start", &start_pos[i] );
255 4 : std::string input2 = argstr[i] + "_end: PDB2CONSTANT REFERENCE=" + iend + " ARG=" + argstr[i];
256 2 : const char* inpt2 = input2.c_str(); plmd.cmd("readInputLine", inpt2 );
257 4 : long erank; plmd.cmd("getDataRank " + argstr[i] + "_end", &erank );
258 2 : if( erank!=1 ) error("should only be one frame in input pdb");
259 4 : std::vector<long> eshape( 1 ); plmd.cmd("getDataShape " + argstr[i] + "_end", &eshape[0] );
260 2 : if( eshape[0]!=1 ) error("should only be one frame in input pdb");
261 4 : plmd.cmd("setMemoryForData " + argstr[i] + "_end", &end_pos[i] );
262 : }
263 : } else {
264 1 : std::string input = "start: PDB2CONSTANT REFERENCE=" + istart;
265 1 : const char* inpt = input.c_str(); plmd.cmd("readInputLine", inpt );
266 1 : FILE* fp=fopen(istart.c_str(),"r"); PDB mypdb; bool do_read=mypdb.readFromFilepointer(fp,false,0.1);
267 14 : indices.resize( mypdb.getAtomNumbers().size() ); for(unsigned i=0; i<indices.size(); ++i) indices[i] = mypdb.getAtomNumbers()[i];
268 1 : long srank; plmd.cmd("getDataRank start", &srank ); if( srank!=2 ) error("should only be one frame in input pdb:" + istart);
269 1 : std::vector<long> sshape( srank ); plmd.cmd("getDataShape start", &sshape[0] ); start_pos.resize( sshape[0]*sshape[1] );
270 1 : if( sshape[0]!=1 ) error("should only be one frame in input pdb: " + istart);
271 1 : plmd.cmd("setMemoryForData start", &start_pos[0] );
272 1 : std::string input2 = "end: PDB2CONSTANT REFERENCE=" + iend;
273 1 : const char* inpt2 = input2.c_str(); plmd.cmd("readInputLine", inpt2 );
274 1 : long erank; plmd.cmd("getDataRank end", &erank ); if( erank!=2 ) error("should only be one frame in input pdb: " + iend);
275 1 : std::vector<long> eshape( erank ); plmd.cmd("getDataShape end", &eshape[0] ); end_pos.resize( eshape[0]*eshape[1] );
276 1 : if( eshape[0]!=1 ) error("should only be one frame in input pdb: " + iend );
277 1 : plmd.cmd("setMemoryForData end", &end_pos[0] );
278 1 : }
279 :
280 : // Now create the metric object
281 : std::vector<double> alig, disp; std::vector<unsigned> residue_indices;
282 5 : if( mtype=="OPTIMAL-FAST" || mtype=="OPTIMAL" || mtype=="SIMPLE" ) {
283 1 : std::string trinput = "endT: TRANSPOSE ARG=end";
284 1 : const char* tinp=trinput.c_str(); plmd.cmd("readInputLine",tinp);
285 1 : PDB pdb; pdb.read(istart,false,0.1);
286 1 : alig.resize( pdb.getOccupancy().size() ); disp.resize( pdb.getBeta().size() );
287 14 : for(unsigned i=0; i<alig.size(); ++i) { alig[i] = pdb.getOccupancy()[i]; disp[i]= pdb.getBeta()[i]; }
288 2 : std::string minput = "d: RMSD_VECTOR DISPLACEMENT SQUARED UNORMALIZED TYPE=" + mtype + " ARG=endT,start";
289 : // Get the align values
290 1 : std::string anum; Tools::convert( alig[0], anum ); minput += " ALIGN=" + anum;
291 13 : for(unsigned i=1; i<alig.size(); ++i) { Tools::convert( alig[i], anum ); minput += "," + anum; }
292 : // Get the displace values
293 1 : std::string dnum; Tools::convert( disp[0], dnum ); minput += " DISPLACE=" + dnum;
294 13 : for(unsigned i=1; i<disp.size(); ++i) { Tools::convert( disp[i], dnum ); minput += "," + dnum; }
295 1 : const char* mcinp=minput.c_str(); plmd.cmd("readInputLine",mcinp);
296 1 : std::string tinput = "CUSTOM ARG=d.disp FUNC=x PERIODIC=NO"; const char* tcinp=tinput.c_str(); plmd.cmd("readInputLine",tcinp);
297 : // Get the residue numbers
298 1 : residue_indices.resize( pdb.size() );
299 14 : for(unsigned i=0; i<residue_indices.size(); ++i) residue_indices[i] = pdb.getResidueNumber( indices[i] );
300 2 : } else if( mtype=="EUCLIDEAN" ) {
301 2 : std::string end_args="ARG1=" + argstr[0] + "_end", start_args="ARG2=" + argstr[0] + "_start";
302 3 : for(unsigned i=1; i<argstr.size(); ++i ) { end_args += "," + argstr[i] + "_end"; start_args += "," + argstr[i] + "_start"; }
303 2 : std::string minput = "d: DISPLACEMENT " + end_args + " " + start_args; const char* mcinp=minput.c_str(); plmd.cmd("readInputLine",mcinp);
304 1 : std::string tinput = "TRANSPOSE ARG=d"; const char* tcinp=tinput.c_str(); plmd.cmd("readInputLine",tcinp);
305 : } else {
306 : // Add functionality to read plumed input here
307 0 : plumed_merror("metric type " + mtype + " has not been implemented");
308 : }
309 :
310 : // Retrieve final displacement vector
311 2 : unsigned aind = plmd.getActionSet().size()-1;
312 : while( true ) {
313 3 : const ActionShortcut* as=dynamic_cast<const ActionShortcut*>( plmd.getActionSet()[aind].get() );
314 3 : if( !as ) break ; aind = aind - 1; plumed_assert( aind>=0 );
315 1 : }
316 2 : ActionWithValue* av = dynamic_cast<ActionWithValue*>( plmd.getActionSet()[aind].get() );
317 2 : if( !av ) error("invalid input for metric" );
318 2 : if( av->getNumberOfComponents()!=1 && av->getName()!="RMSD_VECTOR" ) error("cannot use multi component actions as metric");
319 2 : std::string mydisp = av->copyOutput(0)->getName();
320 : // Now add calls so we can grab the data from plumed
321 4 : long rank; plmd.cmd("getDataRank " + mydisp, &rank );
322 2 : if( rank!=1 ) error("displacement must be a vector quantity");
323 4 : std::vector<long> ishape( rank ); plmd.cmd("getDataShape " + mydisp, &ishape[0] );
324 4 : std::vector<double> displacement( ishape[0] ); plmd.cmd("setMemoryForData " + mydisp, &displacement[0] );
325 : // And calculate the displacement
326 2 : plmd.cmd("calc");
327 :
328 : // Now read in the initial frame
329 :
330 : // Now create frames
331 2 : double delr = 1.0 / static_cast<double>( nbetween ); std::vector<std::vector<double> > allframes; std::vector<double> frame( displacement.size() );
332 4 : for(int i=0; i<nbefore; ++i) {
333 43 : for(unsigned j=0; j<frame.size(); ++j) frame[j] = start_pos[j] - i*delr*displacement[j];
334 2 : allframes.push_back( frame );
335 : }
336 8 : for(unsigned i=1; i<nbetween; ++i) {
337 55 : for(unsigned j=0; j<frame.size(); ++j) frame[j] = start_pos[j] + i*delr*displacement[j];
338 6 : allframes.push_back( frame );
339 : }
340 7 : for(unsigned i=0; i<nafter; ++i) {
341 52 : for(unsigned j=0; j<frame.size(); ++j) frame[j] = end_pos[j] + i*delr*displacement[j];
342 5 : allframes.push_back( frame );
343 : }
344 : // This prints out our final reference configurations
345 4 : plmd.cmd("clear"); plmd.readInputLine("timestep: PUT UNIT=time PERIODIC=NO CONSTANT");
346 2 : ActionToPutData* ts = plmd.getActionSet().selectWithLabel<ActionToPutData*>("timestep");
347 4 : ts->setValuePointer("timestep", &tstep );
348 4 : std::string pinput="DUMPPDB STRIDE=1 DESCRIPTION=PATH FILE=" + ofilename + " FMT=" + ofmt;
349 2 : if( argstr.size()>0 ) {
350 3 : for(unsigned i=0; i<argstr.size(); ++i) {
351 2 : std::string rnum; Tools::convert( allframes[0][i], rnum );
352 2 : std::string inpt = argstr[i] + ": CONSTANT VALUES=" + rnum;
353 20 : for(unsigned j=1; j<allframes.size(); ++j) { Tools::convert( allframes[j][i], rnum ); inpt += "," + rnum; }
354 2 : const char* icinp=inpt.c_str(); plmd.cmd("readInputLine",icinp);
355 4 : if( i==0 ) pinput += " ARG=" + argstr[i]; else pinput += "," + argstr[i];
356 : }
357 : } else {
358 1 : std::string nc; Tools::convert( frame.size(), nc );
359 1 : std::string nr; Tools::convert( allframes.size(), nr );
360 1 : std::string rnum; Tools::convert( allframes[0][0], rnum );
361 2 : std::string inpt = "atom_data: CONSTANT NROWS=" + nr + " NCOLS=" + nc + " VALUES=" + rnum;
362 4 : for(unsigned i=0; i<allframes.size(); ++i) {
363 120 : for(unsigned j=0; j<allframes[i].size(); ++j) {
364 233 : if( i==0 && j==0 ) continue; Tools::convert( allframes[i][j], rnum ); inpt += "," + rnum;
365 : }
366 : }
367 1 : const char* icinp=inpt.c_str(); plmd.cmd("readInputLine",icinp); std::string num; Tools::convert( indices[0].serial(), num );
368 14 : pinput += " ATOMS=atom_data ATOM_INDICES=" + num; for(unsigned i=1; i<indices.size(); ++i ) { Tools::convert( indices[i].serial(), num ); pinput += "," + num; }
369 14 : Tools::convert( residue_indices[0], num ); pinput += " RESIDUE_INDICES=" + num; for(unsigned i=1; i<residue_indices.size(); ++i ) { Tools::convert( residue_indices[i], num ); pinput += "," + num; }
370 1 : std::string anum, dnum; Tools::convert( alig[0], anum ); Tools::convert( disp[0], dnum );
371 14 : pinput += " OCCUPANCY=" + anum; for(unsigned i=1; i<alig.size(); ++i) { Tools::convert( alig[i], anum ); pinput += "," + anum; }
372 14 : pinput += " BETA=" + dnum; for(unsigned i=1; i<disp.size(); ++i) { Tools::convert( disp[i], dnum ); pinput += "," + dnum; }
373 : }
374 2 : const char* pcinp=pinput.c_str(); plmd.cmd("readInputLine",pcinp);
375 : Action* paction = plmd.getActionSet()[plmd.getActionSet().size()-1].get();
376 2 : paction->update();
377 : // And output suggestions on the value of Lambda
378 2 : printLambda( mtype, argstr, ofilename );
379 : return 0;
380 10 : }
381 :
382 4 : void PathTools::printLambda( const std::string& mtype, const std::vector<std::string>& argstr, const std::string& ofile ) {
383 : // Create a PLUMED object
384 4 : PlumedMain plmd; int s=sizeof(double);
385 4 : plmd.cmd("setRealPrecision",&s);
386 4 : plmd.cmd("setMDEngine","pathtools");
387 4 : double tstep=1.0; plmd.cmd("setTimestep",&tstep);
388 4 : plmd.cmd("init");
389 4 : int step=1; plmd.cmd("setStep",&step);
390 :
391 4 : FILE* fp=fopen(ofile.c_str(),"r");
392 : bool do_read=true; unsigned nfram=0;
393 : std::vector<double> alig, disp;
394 : // Read in the argument names
395 8 : for(unsigned i=0; i<argstr.size(); ++i ) {
396 4 : std::string input2 = argstr[i] + ": CONSTANT VALUE=1";
397 4 : const char* inpt2 = input2.c_str(); plmd.cmd("readInputLine", inpt2);
398 : }
399 37 : while (do_read ) {
400 37 : PDB mypdb;
401 : // Read the pdb file
402 37 : do_read=mypdb.readFromFilepointer(fp,false,0.1); if( !do_read ) break;
403 33 : std::string num; Tools::convert( nfram+1, num ); nfram++; std::string iinput;
404 33 : if( argstr.size()>0 ) {
405 60 : for(unsigned i=0; i<argstr.size(); ++i ) {
406 80 : std::string input = "ref_" + num + "_" + argstr[i] + ": PDB2CONSTANT REFERENCE=" + ofile + " ARG=" + argstr[i] + " NUMBER=" + num;
407 40 : const char* inpt = input.c_str(); plmd.cmd("readInputLine", inpt );
408 : }
409 : } else {
410 26 : std::string input = "ref_" + num + ": PDB2CONSTANT REFERENCE=" + ofile + " NUMBER=" + num;
411 13 : const char* inpt = input.c_str(); plmd.cmd("readInputLine", inpt );
412 : }
413 :
414 33 : if( nfram==1 ) {
415 4 : alig.resize( mypdb.getOccupancy().size() );
416 30 : for(unsigned i=0; i<alig.size(); ++i) alig[i]=mypdb.getOccupancy()[i];
417 4 : disp.resize( mypdb.getBeta().size() );
418 30 : for(unsigned i=0; i<disp.size(); ++i) disp[i]=mypdb.getBeta()[i];
419 : }
420 37 : }
421 : // Now create the objects to measure the distances between the frames
422 4 : std::vector<double> data( nfram );
423 33 : for(unsigned j=1; j<nfram; ++j) {
424 29 : std::string istr, jstr; Tools::convert( j, istr); Tools::convert( j+1, jstr );
425 76 : if( mtype=="OPTIMAL-FAST" || mtype=="OPTIMAL" || mtype=="SIMPLE" ) {
426 22 : std::string strv = "ref_" + jstr + "T: TRANSPOSE ARG=ref_" + jstr;
427 11 : const char* cstrv = strv.c_str(); plmd.cmd("readInputLine",cstrv);
428 22 : std::string dstring = "d" + istr + ": RMSD_VECTOR TYPE=" + mtype + " ARG=ref_" + jstr + "T,ref_" + istr;
429 : // Get the align values
430 11 : std::string anum; Tools::convert( alig[0], anum ); dstring += " ALIGN=" + anum;
431 143 : for(unsigned i=1; i<alig.size(); ++i) { Tools::convert( alig[i], anum ); dstring += "," + anum; }
432 : // Get the displace values
433 11 : std::string dnum; Tools::convert( disp[0], dnum ); dstring += " DISPLACE=" + dnum;
434 143 : for(unsigned i=1; i<disp.size(); ++i) { Tools::convert( disp[i], dnum ); dstring += "," + dnum; }
435 11 : const char* icinp=dstring.c_str(); plmd.cmd("readInputLine",icinp);
436 18 : } else if( mtype=="EUCLIDEAN" ) {
437 54 : std::string end_args="ARG1=ref_" + istr + "_" + argstr[0], start_args="ARG2=ref_" + jstr + "_" + argstr[0];
438 54 : for(unsigned i=1; i<argstr.size(); ++i ) { end_args += ",ref_" + istr + "_" + argstr[i]; start_args += ",ref_" + jstr + "_" + argstr[i]; }
439 36 : std::string fstr = "d" + istr + ": EUCLIDEAN_DISTANCE " + end_args + " " + start_args;
440 18 : const char* icinp=fstr.c_str(); plmd.cmd("readInputLine",icinp);
441 : } else {
442 0 : plumed_merror("metric type " + mtype + " has not been implemented");
443 : }
444 58 : long rank; plmd.cmd("getDataRank d" + istr, &rank );
445 29 : if( rank!=1 ) error("distance should be of rank 1");
446 58 : std::vector<long> ishape(1); plmd.cmd("getDataShape d" + istr, &ishape[0] );
447 29 : if( ishape[0]!=1 ) error("distance should be of rank 1");
448 58 : plmd.cmd("setMemoryForData d" + istr, &data[j-1] );
449 : }
450 4 : plmd.cmd("calc"); double mean=0;
451 : printf("N.B. THIS CODE ALWAYS AIMS TO CREATE EQUALLY SPACED FRAMES \n");
452 : printf("THERE MAY BE SMALL DESCREPENCIES IN THE NUMBERS BELOW, HOWEVER, BECAUSE OF ROUNDING ERRORS \n");
453 33 : for(unsigned i=1; i<nfram; ++i) {
454 29 : printf("FINAL DISTANCE BETWEEN FRAME %u AND %u IS %f \n",i-1,i,data[i-1]); mean += data[i-1];
455 : }
456 4 : printf("SUGGESTED LAMBDA PARAMETER IS THUS %f \n",2.3/mean/static_cast<double>( nfram-1 ) );
457 4 : }
458 :
459 : } // End of namespace
460 : }
|
