Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2011-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 "PDB.h"
23 : #include "Tools.h"
24 : #include "Log.h"
25 : #include "h36.h"
26 : #include <cstdio>
27 : #include <iostream>
28 : #include "core/SetupMolInfo.h"
29 : #include "Tensor.h"
30 :
31 : using namespace std;
32 :
33 : //+PLUMEDOC INTERNAL pdbreader
34 : /*
35 : PLUMED can use the PDB format in several places
36 :
37 : - To read molecular structure (\ref MOLINFO).
38 : - To read reference conformations (\ref RMSD, but also many other methods in \ref dists, \ref FIT_TO_TEMPLATE, etc).
39 :
40 : The implemented PDB reader expects a file formatted correctly according to the
41 : [PDB standard](http://www.wwpdb.org/documentation/file-format-content/format33/v3.3.html).
42 : In particular, the following columns are read from ATOM records
43 : \verbatim
44 : columns | content
45 : 1-6 | record name (ATOM or HETATM)
46 : 7-11 | serial number of the atom (starting from 1)
47 : 13-16 | atom name
48 : 18-20 | residue name
49 : 22 | chain id
50 : 23-26 | residue number
51 : 31-38 | x coordinate
52 : 39-46 | y coordinate
53 : 47-54 | z coordinate
54 : 55-60 | occupancy
55 : 61-66 | beta factor
56 : \endverbatim
57 : PLUMED parser is slightly more permissive than the official PDB format
58 : in the fact that the format of real numbers is not fixed. In other words,
59 : any real number that can be parsed is OK and the dot can be placed anywhere. However,
60 : __columns are interpret strictly__. A sample PDB should look like the following
61 : \verbatim
62 : ATOM 2 CH3 ACE 1 12.932 -14.718 -6.016 1.00 1.00
63 : ATOM 5 C ACE 1 21.312 -9.928 -5.946 1.00 1.00
64 : ATOM 9 CA ALA 2 19.462 -11.088 -8.986 1.00 1.00
65 : \endverbatim
66 :
67 : Notice that serial numbers need not to be consecutive. In the three-line example above,
68 : only the coordinates of three atoms are provided. This is perfectly legal and indicates PLUMED
69 : that information about these atoms only is available. This could be both for structural
70 : information in \ref MOLINFO, where the other atoms would have no name assigned, and for
71 : reference structures used in \ref RMSD, where only the provided atoms would be used to compute RMSD.
72 :
73 : \par Occupancy and beta factors
74 :
75 : PLUMED reads also occupancy and beta factors that however are given a very special meaning.
76 : In cases where the PDB structure is used as a reference for an alignment (that's the case
77 : for instance in \ref RMSD and in \ref FIT_TO_TEMPLATE), the occupancy column is used
78 : to provide the weight of each atom in the alignment. In cases where, perhaps after alignment,
79 : the displacement between running coordinates and the provided PDB is computed, the beta factors
80 : are used as weight for the displacement.
81 : Since setting the weights to zero is the same as __not__ including an atom in the alignment or
82 : displacement calculation, the two following reference files would be equivalent when used in an \ref RMSD
83 : calculation. First file:
84 : \verbatim
85 : ATOM 2 CH3 ACE 1 12.932 -14.718 -6.016 1.00 1.00
86 : ATOM 5 C ACE 1 21.312 -9.928 -5.946 1.00 1.00
87 : ATOM 9 CA ALA 2 19.462 -11.088 -8.986 0.00 0.00
88 : \endverbatim
89 : Second file:
90 : \verbatim
91 : ATOM 2 CH3 ACE 1 12.932 -14.718 -6.016 1.00 1.00
92 : ATOM 5 C ACE 1 21.312 -9.928 -5.946 1.00 1.00
93 : \endverbatim
94 : However notice that many extra atoms with zero weight might slow down the calculation, so
95 : removing lines is better than setting their weights to zero.
96 : In addition, weights for alignment need not to be equivalent to weights for displacement.
97 :
98 : \par Systems with more than 100k atoms
99 :
100 : Notice that it very likely does not make any sense to compute the \ref RMSD or any other structural
101 : deviation __using__ so many atoms. However, if the protein for which you want to compute \ref RMSD
102 : has atoms with large serial numbers (e.g. because it is located __after__ solvent in the sorted list of atoms)
103 : you might end up with troubles with the limitations of the PDB format. Indeed, since there are 5
104 : columns available for atom serial number, this number cannot be larger than 99999.
105 : In addition, providing \ref MOLINFO with names associated to atoms with a serial larger than 99999 would be impossible.
106 :
107 : Since PLUMED 2.4 we allow [hybrid 36](http://cci.lbl.gov/hybrid_36/) format to be used to specify atom numbers.
108 : This format is not particularly widespread, but has the nice feature that it provides a one-to-one mapping
109 : between numbers up to approximately 80 millions and strings with 5 characters, plus it is backward compatible
110 : for numbers smaller than 100000. This is not true for notations like the hex notation exported by VMD.
111 : Using the hybrid 36 format, the ATOM records for atom ranging from 99997 to 100002 would read like these:
112 : \verbatim
113 : ATOM 99997 Ar X 1 45.349 38.631 15.116 1.00 1.00
114 : ATOM 99998 Ar X 1 46.189 38.631 15.956 1.00 1.00
115 : ATOM 99999 Ar X 1 46.189 39.471 15.116 1.00 1.00
116 : ATOM A0000 Ar X 1 45.349 39.471 15.956 1.00 1.00
117 : ATOM A0000 Ar X 1 45.349 38.631 16.796 1.00 1.00
118 : ATOM A0001 Ar X 1 46.189 38.631 17.636 1.00 1.00
119 : \endverbatim
120 : There are tools that can be found to translate from integers to strings and back using hybrid 36 format
121 : (a simple python script can be found [here](https://sourceforge.net/p/cctbx/code/HEAD/tree/trunk/iotbx/pdb/hybrid_36.py)).
122 : In addition, as of PLUMED 2.5, we provide a \ref pdbrenumber "command line tool" that can be used to renumber atoms in a PDB file.
123 :
124 : */
125 : //+ENDPLUMEDOC
126 :
127 :
128 : namespace PLMD {
129 :
130 27 : void PDB::setAtomNumbers( const std::vector<AtomNumber>& atoms ) {
131 81 : positions.resize( atoms.size() ); occupancy.resize( atoms.size() );
132 81 : beta.resize( atoms.size() ); numbers.resize( atoms.size() );
133 969 : for(unsigned i=0; i<atoms.size(); ++i) { numbers[i]=atoms[i]; beta[i]=1.0; occupancy[i]=1.0; }
134 27 : }
135 :
136 34 : void PDB::setArgumentNames( const std::vector<std::string>& argument_names ) {
137 68 : argnames.resize( argument_names.size() );
138 260 : for(unsigned i=0; i<argument_names.size(); ++i) {
139 : argnames[i]=argument_names[i];
140 64 : arg_data.insert( std::pair<std::string,double>( argnames[i], 0.0 ) );
141 : }
142 34 : }
143 :
144 1504260 : bool PDB::getArgumentValue( const std::string& name, double& value ) const {
145 : std::map<std::string,double>::const_iterator it = arg_data.find(name);
146 1504260 : if( it!=arg_data.end() ) { value = it->second; return true; }
147 : return false;
148 : }
149 :
150 493287 : void PDB::setAtomPositions( const std::vector<Vector>& pos ) {
151 493287 : plumed_assert( pos.size()==positions.size() );
152 1065105 : for(unsigned i=0; i<positions.size(); ++i) positions[i]=pos[i];
153 493287 : }
154 :
155 1502798 : void PDB::setArgumentValue( const std::string& argname, const double& val ) {
156 : // First set the value of the value of the argument in the map
157 1502798 : arg_data.find(argname)->second = val;
158 1502798 : }
159 :
160 : // bool PDB::hasRequiredProperties( const std::vector<std::string>& inproperties ){
161 : // bool hasprop=false;
162 : // for(unsigned i=0;i<remark.size();++i){
163 : // if( remark[i].find("PROPERTIES=")!=std::string::npos){ hasprop=true; break; }
164 : // }
165 : // if( !hasprop ){
166 : // std::string mypropstr="PROPERTIES=" + inproperties[0];
167 : // for(unsigned i=1;i<inproperties.size();++i) mypropstr += "," + inproperties[i];
168 : // remark.push_back( mypropstr );
169 : // }
170 : // // Now check that all required properties are there
171 : // for(unsigned i=0;i<inproperties.size();++i){
172 : // hasprop=false;
173 : // for(unsigned j=0;j<remark.size();++j){
174 : // if( remark[j].find(inproperties[i]+"=")!=std::string::npos){ hasprop=true; break; }
175 : // }
176 : // if( !hasprop ) return false;
177 : // }
178 : // return true;
179 : // }
180 :
181 17 : void PDB::addBlockEnd( const unsigned& end ) {
182 17 : block_ends.push_back( end );
183 17 : }
184 :
185 472 : unsigned PDB::getNumberOfAtomBlocks()const {
186 472 : return block_ends.size();
187 : }
188 :
189 14 : const std::vector<unsigned> & PDB::getAtomBlockEnds()const {
190 14 : return block_ends;
191 : }
192 :
193 23048038 : const std::vector<Vector> & PDB::getPositions()const {
194 23048038 : return positions;
195 : }
196 :
197 0 : void PDB::setPositions(const std::vector<Vector> &v ) {
198 0 : plumed_assert( v.size()==positions.size() );
199 0 : positions=v;
200 0 : }
201 :
202 16508 : const std::vector<double> & PDB::getOccupancy()const {
203 16508 : return occupancy;
204 : }
205 :
206 16508 : const std::vector<double> & PDB::getBeta()const {
207 16508 : return beta;
208 : }
209 :
210 1330 : void PDB::addRemark( std::vector<std::string>& v1 ) {
211 2660 : Tools::parse(v1,"TYPE",mtype);
212 2660 : Tools::parseVector(v1,"ARG",argnames);
213 9827 : for(unsigned i=0; i<v1.size(); ++i) {
214 2389 : if( v1[i].find("=")!=std::string::npos ) {
215 : std::size_t eq=v1[i].find_first_of('=');
216 1864 : std::string name=v1[i].substr(0,eq);
217 1864 : std::string sval=v1[i].substr(eq+1);
218 1864 : double val; Tools::convert( sval, val );
219 1864 : arg_data.insert( std::pair<std::string,double>( name, val ) );
220 : } else {
221 1050 : flags.push_back(v1[i]);
222 : }
223 : }
224 1330 : }
225 :
226 8 : bool PDB::hasFlag( const std::string& fname ) const {
227 16 : for(unsigned i=0; i<flags.size(); ++i) {
228 0 : if( flags[i]==fname ) return true;
229 : }
230 : return false;
231 : }
232 :
233 :
234 701257 : const std::vector<AtomNumber> & PDB::getAtomNumbers()const {
235 701257 : return numbers;
236 : }
237 :
238 0 : const Vector & PDB::getBoxAxs()const {
239 0 : return BoxXYZ;
240 : }
241 :
242 0 : const Vector & PDB::getBoxAng()const {
243 0 : return BoxABG;
244 : }
245 :
246 12 : const Tensor & PDB::getBoxVec()const {
247 12 : return Box;
248 : }
249 :
250 6755656 : std::string PDB::getAtomName(AtomNumber a)const {
251 : const auto p=number2index.find(a);
252 6755656 : if(p==number2index.end()) return "";
253 6755654 : else return atomsymb[p->second];
254 : }
255 :
256 109652 : unsigned PDB::getResidueNumber(AtomNumber a)const {
257 : const auto p=number2index.find(a);
258 109652 : if(p==number2index.end()) return 0;
259 219300 : else return residue[p->second];
260 : }
261 :
262 96007 : std::string PDB::getResidueName(AtomNumber a) const {
263 : const auto p=number2index.find(a);
264 96007 : if(p==number2index.end()) return "";
265 96005 : else return residuenames[p->second];
266 : }
267 :
268 204171636 : unsigned PDB::size()const {
269 204171636 : return positions.size();
270 : }
271 :
272 2007 : bool PDB::readFromFilepointer(FILE *fp,bool naturalUnits,double scale) {
273 : //cerr<<file<<endl;
274 : bool file_is_alive=false;
275 2007 : if(naturalUnits) scale=1.0;
276 : string line;
277 : fpos_t pos; bool between_ters=true;
278 141519 : while(Tools::getline(fp,line)) {
279 : //cerr<<line<<"\n";
280 141333 : fgetpos (fp,&pos);
281 964929 : while(line.length()<80) line.push_back(' ');
282 141333 : string record=line.substr(0,6);
283 141333 : string serial=line.substr(6,5);
284 141333 : string atomname=line.substr(12,4);
285 141333 : string residuename=line.substr(17,3);
286 141333 : string chainID=line.substr(21,1);
287 141333 : string resnum=line.substr(22,4);
288 141333 : string x=line.substr(30,8);
289 141333 : string y=line.substr(38,8);
290 141333 : string z=line.substr(46,8);
291 141333 : string occ=line.substr(54,6);
292 141333 : string bet=line.substr(60,6);
293 141333 : string BoxX=line.substr(6,9);
294 141333 : string BoxY=line.substr(15,9);
295 141333 : string BoxZ=line.substr(24,9);
296 141333 : string BoxA=line.substr(33,7);
297 141333 : string BoxB=line.substr(40,7);
298 141333 : string BoxG=line.substr(47,7);
299 141333 : Tools::trim(record);
300 141693 : if(record=="TER") { between_ters=false; block_ends.push_back( positions.size() ); }
301 141333 : if(record=="END") { file_is_alive=true; break;}
302 139639 : if(record=="ENDMDL") { file_is_alive=true; break;}
303 139512 : if(record=="REMARK") {
304 5320 : vector<string> v1; v1=Tools::getWords(line.substr(6));
305 1330 : addRemark( v1 );
306 : }
307 139512 : if(record=="CRYST1") {
308 57 : Tools::convert(BoxX,BoxXYZ[0]);
309 57 : Tools::convert(BoxY,BoxXYZ[1]);
310 57 : Tools::convert(BoxZ,BoxXYZ[2]);
311 57 : Tools::convert(BoxA,BoxABG[0]);
312 57 : Tools::convert(BoxB,BoxABG[1]);
313 57 : Tools::convert(BoxG,BoxABG[2]);
314 57 : BoxXYZ*=scale;
315 57 : double cosA=cos(BoxABG[0]*pi/180.);
316 57 : double cosB=cos(BoxABG[1]*pi/180.);
317 57 : double cosG=cos(BoxABG[2]*pi/180.);
318 57 : double sinG=sin(BoxABG[2]*pi/180.);
319 228 : for (unsigned i=0; i<3; i++) {Box[i][0]=0.; Box[i][1]=0.; Box[i][2]=0.;}
320 57 : Box[0][0]=BoxXYZ[0];
321 57 : Box[1][0]=BoxXYZ[1]*cosG;
322 57 : Box[1][1]=BoxXYZ[1]*sinG;
323 57 : Box[2][0]=BoxXYZ[2]*cosB;
324 57 : Box[2][1]=(BoxXYZ[2]*BoxXYZ[1]*cosA-Box[2][0]*Box[1][0])/Box[1][1];
325 57 : Box[2][2]=sqrt(BoxXYZ[2]*BoxXYZ[2]-Box[2][0]*Box[2][0]-Box[2][1]*Box[2][1]);
326 : }
327 141476 : if(record=="ATOM" || record=="HETATM") {
328 : between_ters=true;
329 : AtomNumber a; unsigned resno;
330 : double o,b;
331 137548 : Vector p;
332 : {
333 : int result;
334 : auto trimmed=serial;
335 137548 : Tools::trim(trimmed);
336 137641 : while(trimmed.length()<5) trimmed = std::string(" ") + trimmed;
337 275096 : const char* errmsg = h36::hy36decode(5, trimmed.c_str(),trimmed.length(), &result);
338 137548 : if(errmsg) {
339 0 : std::string msg(errmsg);
340 0 : plumed_merror(msg);
341 : }
342 137548 : a.setSerial(result);
343 : }
344 :
345 137548 : Tools::convert(resnum,resno);
346 137548 : Tools::convert(occ,o);
347 137548 : Tools::convert(bet,b);
348 137548 : Tools::convert(x,p[0]);
349 137548 : Tools::convert(y,p[1]);
350 137548 : Tools::convert(z,p[2]);
351 : // scale into nm
352 137548 : p*=scale;
353 137548 : numbers.push_back(a);
354 137548 : number2index[a]=positions.size();
355 : std::size_t startpos=atomname.find_first_not_of(" \t");
356 : std::size_t endpos=atomname.find_last_not_of(" \t");
357 275096 : atomsymb.push_back( atomname.substr(startpos, endpos-startpos+1) );
358 137548 : residue.push_back(resno);
359 137548 : chain.push_back(chainID);
360 137548 : occupancy.push_back(o);
361 137548 : beta.push_back(b);
362 137548 : positions.push_back(p);
363 137548 : residuenames.push_back(residuename);
364 : }
365 : }
366 5715 : if( between_ters ) block_ends.push_back( positions.size() );
367 2007 : return file_is_alive;
368 : }
369 :
370 261 : bool PDB::read(const std::string&file,bool naturalUnits,double scale) {
371 261 : FILE* fp=fopen(file.c_str(),"r");
372 261 : if(!fp) return false;
373 259 : readFromFilepointer(fp,naturalUnits,scale);
374 259 : fclose(fp);
375 259 : return true;
376 : }
377 :
378 213 : void PDB::getChainNames( std::vector<std::string>& chains ) const {
379 213 : chains.resize(0);
380 213 : chains.push_back( chain[0] );
381 811873 : for(unsigned i=1; i<size(); ++i) {
382 1217490 : if( chains[chains.size()-1]!=chain[i] ) chains.push_back( chain[i] );
383 : }
384 213 : }
385 :
386 11557 : void PDB::getResidueRange( const std::string& chainname, unsigned& res_start, unsigned& res_end, std::string& errmsg ) const {
387 : bool inres=false, foundchain=false;
388 59999539 : for(unsigned i=0; i<size(); ++i) {
389 59987982 : if( chain[i]==chainname ) {
390 26545607 : if(!inres) {
391 11557 : if(foundchain) errmsg="found second start of chain named " + chainname;
392 11557 : res_start=residue[i];
393 : }
394 : inres=true; foundchain=true;
395 3448384 : } else if( inres && chain[i]!=chainname ) {
396 : inres=false;
397 22148 : res_end=residue[i-1];
398 : }
399 : }
400 12040 : if(inres) res_end=residue[size()-1];
401 11557 : }
402 :
403 124 : void PDB::getAtomRange( const std::string& chainname, AtomNumber& a_start, AtomNumber& a_end, std::string& errmsg ) const {
404 : bool inres=false, foundchain=false;
405 336610 : for(unsigned i=0; i<size(); ++i) {
406 336486 : if( chain[i]==chainname ) {
407 31175 : if(!inres) {
408 124 : if(foundchain) errmsg="found second start of chain named " + chainname;
409 124 : a_start=numbers[i];
410 : }
411 : inres=true; foundchain=true;
412 137068 : } else if( inres && chain[i]!=chainname ) {
413 : inres=false;
414 116 : a_end=numbers[i-1];
415 : }
416 : }
417 190 : if(inres) a_end=numbers[size()-1];
418 124 : }
419 :
420 5958 : std::string PDB::getResidueName( const unsigned& resnum ) const {
421 14629242 : for(unsigned i=0; i<size(); ++i) {
422 14635200 : if( residue[i]==resnum ) return residuenames[i];
423 : }
424 0 : return "";
425 : }
426 :
427 46586 : std::string PDB::getResidueName(const unsigned& resnum,const std::string& chainid ) const {
428 118178432 : for(unsigned i=0; i<size(); ++i) {
429 118317988 : if( residue[i]==resnum && ( chainid=="*" || chain[i]==chainid) ) return residuenames[i];
430 : }
431 0 : return "";
432 : }
433 :
434 :
435 13260 : AtomNumber PDB::getNamedAtomFromResidue( const std::string& aname, const unsigned& resnum ) const {
436 32681100 : for(unsigned i=0; i<size(); ++i) {
437 32812980 : if( residue[i]==resnum && atomsymb[i]==aname ) return numbers[i];
438 : }
439 0 : std::string num; Tools::convert( resnum, num );
440 0 : plumed_merror("residue " + num + " does not contain an atom named " + aname );
441 : return numbers[0]; // This is to stop compiler errors
442 : }
443 :
444 2192 : AtomNumber PDB::getNamedAtomFromResidueAndChain( const std::string& aname, const unsigned& resnum, const std::string& chainid ) const {
445 2137032 : for(unsigned i=0; i<size(); ++i) {
446 2175604 : if( residue[i]==resnum && atomsymb[i]==aname && ( chainid=="*" || chain[i]==chainid) ) return numbers[i];
447 : }
448 0 : std::string num; Tools::convert( resnum, num );
449 0 : plumed_merror("residue " + num + " from chain " + chainid + " does not contain an atom named " + aname );
450 : return numbers[0]; // This is to stop compiler errors
451 : }
452 :
453 32148 : std::vector<AtomNumber> PDB::getAtomsInResidue(const unsigned& resnum,const std::string& chainid)const {
454 : std::vector<AtomNumber> tmp;
455 167973300 : for(unsigned i=0; i<size(); ++i) {
456 169398720 : if( residue[i]==resnum && ( chainid=="*" || chain[i]==chainid) ) tmp.push_back(numbers[i]);
457 : }
458 32148 : if(tmp.size()==0) {
459 0 : std::string num; Tools::convert( resnum, num );
460 0 : plumed_merror("Cannot find residue " + num + " from chain " + chainid );
461 : }
462 32148 : return tmp;
463 : }
464 :
465 0 : std::vector<AtomNumber> PDB::getAtomsInChain(const std::string& chainid)const {
466 : std::vector<AtomNumber> tmp;
467 0 : for(unsigned i=0; i<size(); ++i) {
468 0 : if( chainid=="*" || chain[i]==chainid ) tmp.push_back(numbers[i]);
469 : }
470 0 : if(tmp.size()==0) {
471 0 : plumed_merror("Cannot find atoms from chain " + chainid );
472 : }
473 0 : return tmp;
474 : }
475 :
476 4638 : std::string PDB::getChainID(const unsigned& resnumber) const {
477 11373706 : for(unsigned i=0; i<size(); ++i) {
478 11382982 : if(resnumber==residue[i]) return chain[i];
479 : }
480 0 : plumed_merror("Not enough residues in pdb input file");
481 : }
482 :
483 0 : bool PDB::checkForResidue( const std::string& name ) const {
484 0 : for(unsigned i=0; i<size(); ++i) {
485 0 : if( residuenames[i]==name ) return true;
486 : }
487 : return false;
488 : }
489 :
490 0 : bool PDB::checkForAtom( const std::string& name ) const {
491 0 : for(unsigned i=0; i<size(); ++i) {
492 0 : if( atomsymb[i]==name ) return true;
493 : }
494 : return false;
495 : }
496 :
497 0 : Log& operator<<(Log& ostr, const PDB& pdb) {
498 : char buffer[1000];
499 0 : for(unsigned i=0; i<pdb.positions.size(); i++) {
500 0 : sprintf(buffer,"ATOM %3d %8.3f %8.3f %8.3f\n",pdb.numbers[i].serial(),pdb.positions[i][0],pdb.positions[i][1],pdb.positions[i][2]);
501 0 : ostr<<buffer;
502 : }
503 0 : return ostr;
504 : }
505 :
506 852 : Vector PDB::getPosition(AtomNumber a)const {
507 : const auto p=number2index.find(a);
508 852 : if(p==number2index.end()) plumed_merror("atom not available");
509 1704 : else return positions[p->second];
510 : }
511 :
512 2539432 : std::vector<std::string> PDB::getArgumentNames()const {
513 2539432 : return argnames;
514 : }
515 :
516 10 : std::string PDB::getMtype() const {
517 10 : return mtype;
518 : }
519 :
520 497 : void PDB::print( const double& lunits, SetupMolInfo* mymoldat, OFile& ofile, const std::string& fmt ) {
521 497 : if( argnames.size()>0 ) {
522 476 : ofile.printf("REMARK ARG=%s", argnames[0].c_str() );
523 4756 : for(unsigned i=1; i<argnames.size(); ++i) ofile.printf(",%s",argnames[i].c_str() );
524 476 : ofile.printf("\n"); ofile.printf("REMARK ");
525 : }
526 : std::string descr2;
527 497 : if(fmt.find("-")!=std::string::npos) {
528 0 : descr2="%s=" + fmt + " ";
529 : } else {
530 : // This ensures numbers are left justified (i.e. next to the equals sign
531 : std::size_t psign=fmt.find("%");
532 497 : plumed_assert( psign!=std::string::npos );
533 1988 : descr2="%s=%-" + fmt.substr(psign+1) + " ";
534 : }
535 3985 : for(std::map<std::string,double>::iterator it=arg_data.begin(); it!=arg_data.end(); ++it) ofile.printf( descr2.c_str(),it->first.c_str(), it->second );
536 497 : if( argnames.size()>0 ) ofile.printf("\n");
537 497 : if( !mymoldat ) {
538 6295 : for(unsigned i=0; i<positions.size(); ++i) {
539 : std::array<char,6> at;
540 1769 : const char* msg = h36::hy36encode(5,numbers[i].serial(),&at[0]);
541 1769 : plumed_assert(msg==nullptr) << msg;
542 1769 : at[5]=0;
543 8845 : ofile.printf("ATOM %s X RES %4u %8.3f%8.3f%8.3f%6.2f%6.2f\n",
544 : &at[0], i,
545 7076 : lunits*positions[i][0], lunits*positions[i][1], lunits*positions[i][2],
546 : occupancy[i], beta[i] );
547 : }
548 : } else {
549 204 : for(unsigned i=0; i<positions.size(); ++i) {
550 : std::array<char,6> at;
551 66 : const char* msg = h36::hy36encode(5,numbers[i].serial(),&at[0]);
552 66 : plumed_assert(msg==nullptr) << msg;
553 66 : at[5]=0;
554 462 : ofile.printf("ATOM %s %-4s %3s %4u %8.3f%8.3f%8.3f%6.2f%6.2f\n",
555 132 : &at[0], mymoldat->getAtomName(numbers[i]).c_str(),
556 132 : mymoldat->getResidueName(numbers[i]).c_str(), mymoldat->getResidueNumber(numbers[i]),
557 264 : lunits*positions[i][0], lunits*positions[i][1], lunits*positions[i][2],
558 : occupancy[i], beta[i] );
559 : }
560 : }
561 497 : ofile.printf("END\n");
562 497 : }
563 :
564 :
565 5517 : }
566 :
|
