1 /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 Copyright (c) 2011-2021 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/GenericMolInfo.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 Starting with PLUMED 2.7, if all the weights are set to zero they will be normalized to be equal to the
98 inverse of the number of involved atoms. This means that it will be possible to use files with
99 the weight columns set to zero obtaining a meaningful result. In previous PLUMED versions,
100 setting all weights to zero was resulting in an error instead.
101
102
103 \par Systems with more than 100k atoms
104
105 Notice that it very likely does not make any sense to compute the \ref RMSD or any other structural
106 deviation __using__ so many atoms. However, if the protein for which you want to compute \ref RMSD
107 has atoms with large serial numbers (e.g. because it is located __after__ solvent in the sorted list of atoms)
108 you might end up with troubles with the limitations of the PDB format. Indeed, since there are 5
109 columns available for atom serial number, this number cannot be larger than 99999.
110 In addition, providing \ref MOLINFO with names associated to atoms with a serial larger than 99999 would be impossible.
111
112 Since PLUMED 2.4 we allow [hybrid 36](http://cci.lbl.gov/hybrid_36/) format to be used to specify atom numbers.
113 This format is not particularly widespread, but has the nice feature that it provides a one-to-one mapping
114 between numbers up to approximately 80 millions and strings with 5 characters, plus it is backward compatible
115 for numbers smaller than 100000. This is not true for notations like the hex notation exported by VMD.
116 Using the hybrid 36 format, the ATOM records for atom ranging from 99997 to 100002 would read like these:
117 \verbatim
118 ATOM 99997 Ar X 1 45.349 38.631 15.116 1.00 1.00
119 ATOM 99998 Ar X 1 46.189 38.631 15.956 1.00 1.00
120 ATOM 99999 Ar X 1 46.189 39.471 15.116 1.00 1.00
121 ATOM A0000 Ar X 1 45.349 39.471 15.956 1.00 1.00
122 ATOM A0000 Ar X 1 45.349 38.631 16.796 1.00 1.00
123 ATOM A0001 Ar X 1 46.189 38.631 17.636 1.00 1.00
124 \endverbatim
125 There are tools that can be found to translate from integers to strings and back using hybrid 36 format
126 (a simple python script can be found [here](https://sourceforge.net/p/cctbx/code/HEAD/tree/trunk/iotbx/pdb/hybrid_36.py)).
127 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.
128
129 */
130 //+ENDPLUMEDOC
131
132
133 namespace PLMD {
134
setAtomNumbers(const std::vector<AtomNumber> & atoms)135 void PDB::setAtomNumbers( const std::vector<AtomNumber>& atoms ) {
136 positions.resize( atoms.size() ); occupancy.resize( atoms.size() );
137 beta.resize( atoms.size() ); numbers.resize( atoms.size() );
138 for(unsigned i=0; i<atoms.size(); ++i) { numbers[i]=atoms[i]; beta[i]=1.0; occupancy[i]=1.0; }
139 }
140
setArgumentNames(const std::vector<std::string> & argument_names)141 void PDB::setArgumentNames( const std::vector<std::string>& argument_names ) {
142 argnames.resize( argument_names.size() );
143 for(unsigned i=0; i<argument_names.size(); ++i) {
144 argnames[i]=argument_names[i];
145 arg_data.insert( std::pair<std::string,double>( argnames[i], 0.0 ) );
146 }
147 }
148
getArgumentValue(const std::string & name,double & value) const149 bool PDB::getArgumentValue( const std::string& name, double& value ) const {
150 std::map<std::string,double>::const_iterator it = arg_data.find(name);
151 if( it!=arg_data.end() ) { value = it->second; return true; }
152 return false;
153 }
154
setAtomPositions(const std::vector<Vector> & pos)155 void PDB::setAtomPositions( const std::vector<Vector>& pos ) {
156 plumed_assert( pos.size()==positions.size() );
157 for(unsigned i=0; i<positions.size(); ++i) positions[i]=pos[i];
158 }
159
setArgumentValue(const std::string & argname,const double & val)160 void PDB::setArgumentValue( const std::string& argname, const double& val ) {
161 // First set the value of the value of the argument in the map
162 arg_data.find(argname)->second = val;
163 }
164
165 // bool PDB::hasRequiredProperties( const std::vector<std::string>& inproperties ){
166 // bool hasprop=false;
167 // for(unsigned i=0;i<remark.size();++i){
168 // if( remark[i].find("PROPERTIES=")!=std::string::npos){ hasprop=true; break; }
169 // }
170 // if( !hasprop ){
171 // std::string mypropstr="PROPERTIES=" + inproperties[0];
172 // for(unsigned i=1;i<inproperties.size();++i) mypropstr += "," + inproperties[i];
173 // remark.push_back( mypropstr );
174 // }
175 // // Now check that all required properties are there
176 // for(unsigned i=0;i<inproperties.size();++i){
177 // hasprop=false;
178 // for(unsigned j=0;j<remark.size();++j){
179 // if( remark[j].find(inproperties[i]+"=")!=std::string::npos){ hasprop=true; break; }
180 // }
181 // if( !hasprop ) return false;
182 // }
183 // return true;
184 // }
185
addBlockEnd(const unsigned & end)186 void PDB::addBlockEnd( const unsigned& end ) {
187 block_ends.push_back( end );
188 }
189
getNumberOfAtomBlocks() const190 unsigned PDB::getNumberOfAtomBlocks()const {
191 return block_ends.size();
192 }
193
getAtomBlockEnds() const194 const std::vector<unsigned> & PDB::getAtomBlockEnds()const {
195 return block_ends;
196 }
197
getPositions() const198 const std::vector<Vector> & PDB::getPositions()const {
199 return positions;
200 }
201
setPositions(const std::vector<Vector> & v)202 void PDB::setPositions(const std::vector<Vector> &v ) {
203 plumed_assert( v.size()==positions.size() );
204 positions=v;
205 }
206
getOccupancy() const207 const std::vector<double> & PDB::getOccupancy()const {
208 return occupancy;
209 }
210
getBeta() const211 const std::vector<double> & PDB::getBeta()const {
212 return beta;
213 }
214
addRemark(std::vector<std::string> & v1)215 void PDB::addRemark( std::vector<std::string>& v1 ) {
216 Tools::parse(v1,"TYPE",mtype);
217 Tools::parseVector(v1,"ARG",argnames);
218 for(unsigned i=0; i<v1.size(); ++i) {
219 if( v1[i].find("=")!=std::string::npos ) {
220 std::size_t eq=v1[i].find_first_of('=');
221 std::string name=v1[i].substr(0,eq);
222 std::string sval=v1[i].substr(eq+1);
223 double val; Tools::convert( sval, val );
224 arg_data.insert( std::pair<std::string,double>( name, val ) );
225 } else {
226 flags.push_back(v1[i]);
227 }
228 }
229 }
230
hasFlag(const std::string & fname) const231 bool PDB::hasFlag( const std::string& fname ) const {
232 for(unsigned i=0; i<flags.size(); ++i) {
233 if( flags[i]==fname ) return true;
234 }
235 return false;
236 }
237
238
getAtomNumbers() const239 const std::vector<AtomNumber> & PDB::getAtomNumbers()const {
240 return numbers;
241 }
242
getBoxAxs() const243 const Vector & PDB::getBoxAxs()const {
244 return BoxXYZ;
245 }
246
getBoxAng() const247 const Vector & PDB::getBoxAng()const {
248 return BoxABG;
249 }
250
getBoxVec() const251 const Tensor & PDB::getBoxVec()const {
252 return Box;
253 }
254
getAtomName(AtomNumber a) const255 std::string PDB::getAtomName(AtomNumber a)const {
256 const auto p=number2index.find(a);
257 if(p==number2index.end()) {
258 std::string num; Tools::convert( a.serial(), num );
259 plumed_merror("Name of atom " + num + " not found" );
260 return "";
261 } else return atomsymb[p->second];
262 }
263
getResidueNumber(AtomNumber a) const264 unsigned PDB::getResidueNumber(AtomNumber a)const {
265 const auto p=number2index.find(a);
266 if(p==number2index.end()) {
267 std::string num; Tools::convert( a.serial(), num );
268 plumed_merror("Residue for atom " + num + " not found" );
269 return 0;
270 } else return residue[p->second];
271 }
272
getResidueName(AtomNumber a) const273 std::string PDB::getResidueName(AtomNumber a) const {
274 const auto p=number2index.find(a);
275 if(p==number2index.end()) {
276 std::string num; Tools::convert( a.serial(), num );
277 plumed_merror("Residue for atom " + num + " not found" );
278 return "";
279 } else return residuenames[p->second];
280 }
281
size() const282 unsigned PDB::size()const {
283 return positions.size();
284 }
285
readFromFilepointer(FILE * fp,bool naturalUnits,double scale)286 bool PDB::readFromFilepointer(FILE *fp,bool naturalUnits,double scale) {
287 //cerr<<file<<endl;
288 bool file_is_alive=false;
289 if(naturalUnits) scale=1.0;
290 string line;
291 fpos_t pos; bool between_ters=true;
292 while(Tools::getline(fp,line)) {
293 //cerr<<line<<"\n";
294 fgetpos (fp,&pos);
295 while(line.length()<80) line.push_back(' ');
296 string record=line.substr(0,6);
297 string serial=line.substr(6,5);
298 string atomname=line.substr(12,4);
299 string residuename=line.substr(17,3);
300 string chainID=line.substr(21,1);
301 string resnum=line.substr(22,4);
302 string x=line.substr(30,8);
303 string y=line.substr(38,8);
304 string z=line.substr(46,8);
305 string occ=line.substr(54,6);
306 string bet=line.substr(60,6);
307 string BoxX=line.substr(6,9);
308 string BoxY=line.substr(15,9);
309 string BoxZ=line.substr(24,9);
310 string BoxA=line.substr(33,7);
311 string BoxB=line.substr(40,7);
312 string BoxG=line.substr(47,7);
313 Tools::trim(record);
314 if(record=="TER") { between_ters=false; block_ends.push_back( positions.size() ); }
315 if(record=="END") { file_is_alive=true; break;}
316 if(record=="ENDMDL") { file_is_alive=true; break;}
317 if(record=="REMARK") {
318 vector<string> v1; v1=Tools::getWords(line.substr(6));
319 addRemark( v1 );
320 }
321 if(record=="CRYST1") {
322 Tools::convert(BoxX,BoxXYZ[0]);
323 Tools::convert(BoxY,BoxXYZ[1]);
324 Tools::convert(BoxZ,BoxXYZ[2]);
325 Tools::convert(BoxA,BoxABG[0]);
326 Tools::convert(BoxB,BoxABG[1]);
327 Tools::convert(BoxG,BoxABG[2]);
328 BoxXYZ*=scale;
329 double cosA=cos(BoxABG[0]*pi/180.);
330 double cosB=cos(BoxABG[1]*pi/180.);
331 double cosG=cos(BoxABG[2]*pi/180.);
332 double sinG=sin(BoxABG[2]*pi/180.);
333 for (unsigned i=0; i<3; i++) {Box[i][0]=0.; Box[i][1]=0.; Box[i][2]=0.;}
334 Box[0][0]=BoxXYZ[0];
335 Box[1][0]=BoxXYZ[1]*cosG;
336 Box[1][1]=BoxXYZ[1]*sinG;
337 Box[2][0]=BoxXYZ[2]*cosB;
338 Box[2][1]=(BoxXYZ[2]*BoxXYZ[1]*cosA-Box[2][0]*Box[1][0])/Box[1][1];
339 Box[2][2]=sqrt(BoxXYZ[2]*BoxXYZ[2]-Box[2][0]*Box[2][0]-Box[2][1]*Box[2][1]);
340 }
341 if(record=="ATOM" || record=="HETATM") {
342 between_ters=true;
343 AtomNumber a; unsigned resno;
344 double o,b;
345 Vector p;
346 {
347 int result;
348 auto trimmed=serial;
349 Tools::trim(trimmed);
350 while(trimmed.length()<5) trimmed = std::string(" ") + trimmed;
351 const char* errmsg = h36::hy36decode(5, trimmed.c_str(),trimmed.length(), &result);
352 if(errmsg) {
353 std::string msg(errmsg);
354 plumed_merror(msg);
355 }
356 a.setSerial(result);
357 }
358
359 Tools::convert(resnum,resno);
360 Tools::convert(occ,o);
361 Tools::convert(bet,b);
362 Tools::convert(x,p[0]);
363 Tools::convert(y,p[1]);
364 Tools::convert(z,p[2]);
365 // scale into nm
366 p*=scale;
367 numbers.push_back(a);
368 number2index[a]=positions.size();
369 std::size_t startpos=atomname.find_first_not_of(" \t");
370 std::size_t endpos=atomname.find_last_not_of(" \t");
371 atomsymb.push_back( atomname.substr(startpos, endpos-startpos+1) );
372 residue.push_back(resno);
373 chain.push_back(chainID);
374 occupancy.push_back(o);
375 beta.push_back(b);
376 positions.push_back(p);
377 residuenames.push_back(residuename);
378 }
379 }
380 if( between_ters ) block_ends.push_back( positions.size() );
381 return file_is_alive;
382 }
383
read(const std::string & file,bool naturalUnits,double scale)384 bool PDB::read(const std::string&file,bool naturalUnits,double scale) {
385 FILE* fp=fopen(file.c_str(),"r");
386 if(!fp) return false;
387 readFromFilepointer(fp,naturalUnits,scale);
388 fclose(fp);
389 return true;
390 }
391
getChainNames(std::vector<std::string> & chains) const392 void PDB::getChainNames( std::vector<std::string>& chains ) const {
393 chains.resize(0);
394 chains.push_back( chain[0] );
395 for(unsigned i=1; i<size(); ++i) {
396 if( chains[chains.size()-1]!=chain[i] ) chains.push_back( chain[i] );
397 }
398 }
399
getResidueRange(const std::string & chainname,unsigned & res_start,unsigned & res_end,std::string & errmsg) const400 void PDB::getResidueRange( const std::string& chainname, unsigned& res_start, unsigned& res_end, std::string& errmsg ) const {
401 bool inres=false, foundchain=false;
402 for(unsigned i=0; i<size(); ++i) {
403 if( chain[i]==chainname ) {
404 if(!inres) {
405 if(foundchain) errmsg="found second start of chain named " + chainname;
406 res_start=residue[i];
407 }
408 inres=true; foundchain=true;
409 } else if( inres && chain[i]!=chainname ) {
410 inres=false;
411 res_end=residue[i-1];
412 }
413 }
414 if(inres) res_end=residue[size()-1];
415 }
416
getAtomRange(const std::string & chainname,AtomNumber & a_start,AtomNumber & a_end,std::string & errmsg) const417 void PDB::getAtomRange( const std::string& chainname, AtomNumber& a_start, AtomNumber& a_end, std::string& errmsg ) const {
418 bool inres=false, foundchain=false;
419 for(unsigned i=0; i<size(); ++i) {
420 if( chain[i]==chainname ) {
421 if(!inres) {
422 if(foundchain) errmsg="found second start of chain named " + chainname;
423 a_start=numbers[i];
424 }
425 inres=true; foundchain=true;
426 } else if( inres && chain[i]!=chainname ) {
427 inres=false;
428 a_end=numbers[i-1];
429 }
430 }
431 if(inres) a_end=numbers[size()-1];
432 }
433
getResidueName(const unsigned & resnum) const434 std::string PDB::getResidueName( const unsigned& resnum ) const {
435 for(unsigned i=0; i<size(); ++i) {
436 if( residue[i]==resnum ) return residuenames[i];
437 }
438 std::string num; Tools::convert( resnum, num );
439 plumed_merror("residue " + num + " not found" );
440 return "";
441 }
442
getResidueName(const unsigned & resnum,const std::string & chainid) const443 std::string PDB::getResidueName(const unsigned& resnum,const std::string& chainid ) const {
444 for(unsigned i=0; i<size(); ++i) {
445 if( residue[i]==resnum && ( chainid=="*" || chain[i]==chainid) ) return residuenames[i];
446 }
447 std::string num; Tools::convert( resnum, num );
448 plumed_merror("residue " + num + " not found in chain " + chainid );
449 return "";
450 }
451
452
getNamedAtomFromResidue(const std::string & aname,const unsigned & resnum) const453 AtomNumber PDB::getNamedAtomFromResidue( const std::string& aname, const unsigned& resnum ) const {
454 for(unsigned i=0; i<size(); ++i) {
455 if( residue[i]==resnum && atomsymb[i]==aname ) return numbers[i];
456 }
457 std::string num; Tools::convert( resnum, num );
458 plumed_merror("residue " + num + " does not contain an atom named " + aname );
459 return numbers[0]; // This is to stop compiler errors
460 }
461
getNamedAtomFromResidueAndChain(const std::string & aname,const unsigned & resnum,const std::string & chainid) const462 AtomNumber PDB::getNamedAtomFromResidueAndChain( const std::string& aname, const unsigned& resnum, const std::string& chainid ) const {
463 for(unsigned i=0; i<size(); ++i) {
464 if( residue[i]==resnum && atomsymb[i]==aname && ( chainid=="*" || chain[i]==chainid) ) return numbers[i];
465 }
466 std::string num; Tools::convert( resnum, num );
467 plumed_merror("residue " + num + " from chain " + chainid + " does not contain an atom named " + aname );
468 return numbers[0]; // This is to stop compiler errors
469 }
470
getAtomsInResidue(const unsigned & resnum,const std::string & chainid) const471 std::vector<AtomNumber> PDB::getAtomsInResidue(const unsigned& resnum,const std::string& chainid)const {
472 std::vector<AtomNumber> tmp;
473 for(unsigned i=0; i<size(); ++i) {
474 if( residue[i]==resnum && ( chainid=="*" || chain[i]==chainid) ) tmp.push_back(numbers[i]);
475 }
476 if(tmp.size()==0) {
477 std::string num; Tools::convert( resnum, num );
478 plumed_merror("Cannot find residue " + num + " from chain " + chainid );
479 }
480 return tmp;
481 }
482
getAtomsInChain(const std::string & chainid) const483 std::vector<AtomNumber> PDB::getAtomsInChain(const std::string& chainid)const {
484 std::vector<AtomNumber> tmp;
485 for(unsigned i=0; i<size(); ++i) {
486 if( chainid=="*" || chain[i]==chainid ) tmp.push_back(numbers[i]);
487 }
488 if(tmp.size()==0) {
489 plumed_merror("Cannot find atoms from chain " + chainid );
490 }
491 return tmp;
492 }
493
getChainID(const unsigned & resnumber) const494 std::string PDB::getChainID(const unsigned& resnumber) const {
495 for(unsigned i=0; i<size(); ++i) {
496 if(resnumber==residue[i]) return chain[i];
497 }
498 plumed_merror("Not enough residues in pdb input file");
499 }
500
checkForResidue(const std::string & name) const501 bool PDB::checkForResidue( const std::string& name ) const {
502 for(unsigned i=0; i<size(); ++i) {
503 if( residuenames[i]==name ) return true;
504 }
505 return false;
506 }
507
checkForAtom(const std::string & name) const508 bool PDB::checkForAtom( const std::string& name ) const {
509 for(unsigned i=0; i<size(); ++i) {
510 if( atomsymb[i]==name ) return true;
511 }
512 return false;
513 }
514
checkForAtom(AtomNumber a) const515 bool PDB::checkForAtom( AtomNumber a ) const {
516 const auto p=number2index.find(a);
517 return (p!=number2index.end());
518 }
519
operator <<(Log & ostr,const PDB & pdb)520 Log& operator<<(Log& ostr, const PDB& pdb) {
521 char buffer[1000];
522 for(unsigned i=0; i<pdb.positions.size(); i++) {
523 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]);
524 ostr<<buffer;
525 }
526 return ostr;
527 }
528
getPosition(AtomNumber a) const529 Vector PDB::getPosition(AtomNumber a)const {
530 const auto p=number2index.find(a);
531 if(p==number2index.end()) plumed_merror("atom not available");
532 else return positions[p->second];
533 }
534
getArgumentNames() const535 std::vector<std::string> PDB::getArgumentNames()const {
536 return argnames;
537 }
538
getMtype() const539 std::string PDB::getMtype() const {
540 return mtype;
541 }
542
print(const double & lunits,GenericMolInfo * mymoldat,OFile & ofile,const std::string & fmt)543 void PDB::print( const double& lunits, GenericMolInfo* mymoldat, OFile& ofile, const std::string& fmt ) {
544 if( argnames.size()>0 ) {
545 ofile.printf("REMARK ARG=%s", argnames[0].c_str() );
546 for(unsigned i=1; i<argnames.size(); ++i) ofile.printf(",%s",argnames[i].c_str() );
547 ofile.printf("\n"); ofile.printf("REMARK ");
548 }
549 std::string descr2;
550 if(fmt.find("-")!=std::string::npos) {
551 descr2="%s=" + fmt + " ";
552 } else {
553 // This ensures numbers are left justified (i.e. next to the equals sign
554 std::size_t psign=fmt.find("%");
555 plumed_assert( psign!=std::string::npos );
556 descr2="%s=%-" + fmt.substr(psign+1) + " ";
557 }
558 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 );
559 if( argnames.size()>0 ) ofile.printf("\n");
560 if( !mymoldat ) {
561 for(unsigned i=0; i<positions.size(); ++i) {
562 std::array<char,6> at;
563 const char* msg = h36::hy36encode(5,numbers[i].serial(),&at[0]);
564 plumed_assert(msg==nullptr) << msg;
565 at[5]=0;
566 ofile.printf("ATOM %s X RES %4u %8.3f%8.3f%8.3f%6.2f%6.2f\n",
567 &at[0], i,
568 lunits*positions[i][0], lunits*positions[i][1], lunits*positions[i][2],
569 occupancy[i], beta[i] );
570 }
571 } else {
572 for(unsigned i=0; i<positions.size(); ++i) {
573 std::array<char,6> at;
574 const char* msg = h36::hy36encode(5,numbers[i].serial(),&at[0]);
575 plumed_assert(msg==nullptr) << msg;
576 at[5]=0;
577 ofile.printf("ATOM %s %-4s %3s %4u %8.3f%8.3f%8.3f%6.2f%6.2f\n",
578 &at[0], mymoldat->getAtomName(numbers[i]).c_str(),
579 mymoldat->getResidueName(numbers[i]).c_str(), mymoldat->getResidueNumber(numbers[i]),
580 lunits*positions[i][0], lunits*positions[i][1], lunits*positions[i][2],
581 occupancy[i], beta[i] );
582 }
583 }
584 ofile.printf("END\n");
585 }
586
587
588 }
589
590