Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2011-2023 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 <set>
28 : #include "core/GenericMolInfo.h"
29 : #include "Tensor.h"
30 :
31 : //+PLUMEDOC INTERNAL pdbreader
32 : /*
33 : PLUMED can use the PDB format in several places
34 :
35 : - To read molecular structure (\ref MOLINFO).
36 : - To read reference conformations (\ref RMSD, but also many other methods in \ref dists, \ref FIT_TO_TEMPLATE, etc).
37 :
38 : The implemented PDB reader expects a file formatted correctly according to the
39 : [PDB standard](http://www.wwpdb.org/documentation/file-format-content/format33/v3.3.html).
40 : In particular, the following columns are read from ATOM records
41 : \verbatim
42 : columns | content
43 : 1-6 | record name (ATOM or HETATM)
44 : 7-11 | serial number of the atom (starting from 1)
45 : 13-16 | atom name
46 : 18-20 | residue name
47 : 22 | chain id
48 : 23-26 | residue number
49 : 31-38 | x coordinate
50 : 39-46 | y coordinate
51 : 47-54 | z coordinate
52 : 55-60 | occupancy
53 : 61-66 | beta factor
54 : \endverbatim
55 : PLUMED parser is slightly more permissive than the official PDB format
56 : in the fact that the format of real numbers is not fixed. In other words,
57 : any real number that can be parsed is OK and the dot can be placed anywhere. However,
58 : __columns are interpret strictly__. A sample PDB should look like the following
59 : \verbatim
60 : ATOM 2 CH3 ACE 1 12.932 -14.718 -6.016 1.00 1.00
61 : ATOM 5 C ACE 1 21.312 -9.928 -5.946 1.00 1.00
62 : ATOM 9 CA ALA 2 19.462 -11.088 -8.986 1.00 1.00
63 : \endverbatim
64 :
65 : Notice that serial numbers need not to be consecutive. In the three-line example above,
66 : only the coordinates of three atoms are provided. This is perfectly legal and indicates PLUMED
67 : that information about these atoms only is available. This could be both for structural
68 : information in \ref MOLINFO, where the other atoms would have no name assigned, and for
69 : reference structures used in \ref RMSD, where only the provided atoms would be used to compute RMSD.
70 :
71 : \par Occupancy and beta factors
72 :
73 : PLUMED reads also occupancy and beta factors that however are given a very special meaning.
74 : In cases where the PDB structure is used as a reference for an alignment (that's the case
75 : for instance in \ref RMSD and in \ref FIT_TO_TEMPLATE), the occupancy column is used
76 : to provide the weight of each atom in the alignment. In cases where, perhaps after alignment,
77 : the displacement between running coordinates and the provided PDB is computed, the beta factors
78 : are used as weight for the displacement.
79 : Since setting the weights to zero is the same as __not__ including an atom in the alignment or
80 : displacement calculation, the two following reference files would be equivalent when used in an \ref RMSD
81 : calculation. First file:
82 : \verbatim
83 : ATOM 2 CH3 ACE 1 12.932 -14.718 -6.016 1.00 1.00
84 : ATOM 5 C ACE 1 21.312 -9.928 -5.946 1.00 1.00
85 : ATOM 9 CA ALA 2 19.462 -11.088 -8.986 0.00 0.00
86 : \endverbatim
87 : Second file:
88 : \verbatim
89 : ATOM 2 CH3 ACE 1 12.932 -14.718 -6.016 1.00 1.00
90 : ATOM 5 C ACE 1 21.312 -9.928 -5.946 1.00 1.00
91 : \endverbatim
92 : However notice that many extra atoms with zero weight might slow down the calculation, so
93 : removing lines is better than setting their weights to zero.
94 : In addition, weights for alignment need not to be equivalent to weights for displacement.
95 : Starting with PLUMED 2.7, if all the weights are set to zero they will be normalized to be equal to the
96 : inverse of the number of involved atoms. This means that it will be possible to use files with
97 : the weight columns set to zero obtaining a meaningful result. In previous PLUMED versions,
98 : setting all weights to zero was resulting in an error instead.
99 :
100 :
101 : \par Systems with more than 100k atoms
102 :
103 : Notice that it very likely does not make any sense to compute the \ref RMSD or any other structural
104 : deviation __using__ so many atoms. However, if the protein for which you want to compute \ref RMSD
105 : has atoms with large serial numbers (e.g. because it is located __after__ solvent in the sorted list of atoms)
106 : you might end up with troubles with the limitations of the PDB format. Indeed, since there are 5
107 : columns available for atom serial number, this number cannot be larger than 99999.
108 : In addition, providing \ref MOLINFO with names associated to atoms with a serial larger than 99999 would be impossible.
109 :
110 : Since PLUMED 2.4 we allow [hybrid 36](http://cci.lbl.gov/hybrid_36/) format to be used to specify atom numbers.
111 : This format is not particularly widespread, but has the nice feature that it provides a one-to-one mapping
112 : between numbers up to approximately 80 millions and strings with 5 characters, plus it is backward compatible
113 : for numbers smaller than 100000. This is not true for notations like the hex notation exported by VMD.
114 : Using the hybrid 36 format, the ATOM records for atom ranging from 99997 to 100002 would read like these:
115 : \verbatim
116 : ATOM 99997 Ar X 1 45.349 38.631 15.116 1.00 1.00
117 : ATOM 99998 Ar X 1 46.189 38.631 15.956 1.00 1.00
118 : ATOM 99999 Ar X 1 46.189 39.471 15.116 1.00 1.00
119 : ATOM A0000 Ar X 1 45.349 39.471 15.956 1.00 1.00
120 : ATOM A0000 Ar X 1 45.349 38.631 16.796 1.00 1.00
121 : ATOM A0001 Ar X 1 46.189 38.631 17.636 1.00 1.00
122 : \endverbatim
123 : There are tools that can be found to translate from integers to strings and back using hybrid 36 format
124 : (a simple python script can be found [here](https://sourceforge.net/p/cctbx/code/HEAD/tree/trunk/iotbx/pdb/hybrid_36.py)).
125 : 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.
126 :
127 : */
128 : //+ENDPLUMEDOC
129 :
130 :
131 : namespace PLMD {
132 :
133 0 : void PDB::setAtomNumbers( const std::vector<AtomNumber>& atoms ) {
134 0 : positions.resize( atoms.size() );
135 0 : occupancy.resize( atoms.size() );
136 0 : beta.resize( atoms.size() );
137 0 : numbers.resize( atoms.size() );
138 0 : for(unsigned i=0; i<atoms.size(); ++i) {
139 0 : numbers[i]=atoms[i];
140 0 : beta[i]=1.0;
141 0 : occupancy[i]=1.0;
142 : }
143 0 : }
144 :
145 0 : void PDB::setArgumentNames( const std::vector<std::string>& argument_names ) {
146 0 : argnames.resize( argument_names.size() );
147 0 : std::vector<double> tmp(1,0);
148 0 : for(unsigned i=0; i<argument_names.size(); ++i) {
149 : argnames[i]=argument_names[i];
150 0 : arg_data.insert( std::pair<std::string,std::vector<double> >( argnames[i], tmp ) );
151 : }
152 0 : }
153 :
154 2148 : bool PDB::getArgumentValue( const std::string& name, std::vector<double>& value ) const {
155 : std::map<std::string,std::vector<double> >::const_iterator it = arg_data.find(name);
156 2148 : if( it!=arg_data.end() ) {
157 2148 : if( value.size()!=it->second.size() ) {
158 : return false;
159 : }
160 4300 : for(unsigned i=0; i<value.size(); ++i) {
161 2152 : value[i] = it->second[i];
162 : }
163 : return true;
164 : }
165 : return false;
166 : }
167 :
168 0 : void PDB::setAtomPositions( const std::vector<Vector>& pos ) {
169 0 : plumed_assert( pos.size()==positions.size() );
170 0 : for(unsigned i=0; i<positions.size(); ++i) {
171 0 : positions[i]=pos[i];
172 : }
173 0 : }
174 :
175 0 : void PDB::setArgumentValue( const std::string& argname, const double& val ) {
176 : // First set the value of the value of the argument in the map
177 0 : arg_data.find(argname)->second.resize(1);
178 0 : arg_data.find(argname)->second[0] = val;
179 0 : }
180 :
181 : // bool PDB::hasRequiredProperties( const std::vector<std::string>& inproperties ){
182 : // bool hasprop=false;
183 : // for(unsigned i=0;i<remark.size();++i){
184 : // if( remark[i].find("PROPERTIES=")!=std::string::npos){ hasprop=true; break; }
185 : // }
186 : // if( !hasprop ){
187 : // std::string mypropstr="PROPERTIES=" + inproperties[0];
188 : // for(unsigned i=1;i<inproperties.size();++i) mypropstr += "," + inproperties[i];
189 : // remark.push_back( mypropstr );
190 : // }
191 : // // Now check that all required properties are there
192 : // for(unsigned i=0;i<inproperties.size();++i){
193 : // hasprop=false;
194 : // for(unsigned j=0;j<remark.size();++j){
195 : // if( remark[j].find(inproperties[i]+"=")!=std::string::npos){ hasprop=true; break; }
196 : // }
197 : // if( !hasprop ) return false;
198 : // }
199 : // return true;
200 : // }
201 :
202 0 : void PDB::addBlockEnd( const unsigned& end ) {
203 0 : block_ends.push_back( end );
204 0 : }
205 :
206 3 : unsigned PDB::getNumberOfAtomBlocks()const {
207 3 : return block_ends.size();
208 : }
209 :
210 6 : const std::vector<unsigned> & PDB::getAtomBlockEnds()const {
211 6 : return block_ends;
212 : }
213 :
214 23595961 : const std::vector<Vector> & PDB::getPositions()const {
215 23595961 : return positions;
216 : }
217 :
218 0 : void PDB::setPositions(const std::vector<Vector> &v ) {
219 0 : plumed_assert( v.size()==positions.size() );
220 0 : positions=v;
221 0 : }
222 :
223 154221 : const std::vector<double> & PDB::getOccupancy()const {
224 154221 : return occupancy;
225 : }
226 :
227 152365 : const std::vector<double> & PDB::getBeta()const {
228 152365 : return beta;
229 : }
230 :
231 4620 : void PDB::addRemark( std::vector<std::string>& v1 ) {
232 4620 : Tools::parse(v1,"TYPE",mtype);
233 4620 : Tools::parseVector(v1,"ARG",argnames);
234 12565 : for(unsigned i=0; i<v1.size(); ++i) {
235 7945 : if( v1[i].find("=")!=std::string::npos ) {
236 : std::size_t eq=v1[i].find_first_of('=');
237 7286 : std::string name=v1[i].substr(0,eq);
238 7286 : std::string sval=v1[i].substr(eq+1);
239 7286 : std::vector<std::string> words=Tools::getWords(sval,"\t\n ,");
240 7286 : std::vector<double> val( words.size() );
241 14580 : for(unsigned ii=0; ii<val.size(); ++ii) {
242 7294 : Tools::convert( words[ii], val[ii] );
243 : }
244 7286 : arg_data.insert( std::pair<std::string,std::vector<double> >( name, val ) );
245 7286 : } else {
246 659 : flags.push_back(v1[i]);
247 : }
248 : }
249 4620 : }
250 :
251 0 : bool PDB::hasFlag( const std::string& fname ) const {
252 0 : for(unsigned i=0; i<flags.size(); ++i) {
253 0 : if( flags[i]==fname ) {
254 : return true;
255 : }
256 : }
257 : return false;
258 : }
259 :
260 :
261 318050 : const std::vector<AtomNumber> & PDB::getAtomNumbers()const {
262 318050 : return numbers;
263 : }
264 :
265 18 : const Vector & PDB::getBoxAxs()const {
266 18 : return BoxXYZ;
267 : }
268 :
269 3 : const Vector & PDB::getBoxAng()const {
270 3 : return BoxABG;
271 : }
272 :
273 28 : const Tensor & PDB::getBoxVec()const {
274 28 : return Box;
275 : }
276 :
277 16 : const std::string & PDB::getSpaceGroup()const {
278 16 : return spaceGroup;
279 : }
280 :
281 8575206 : std::string PDB::getAtomName(AtomNumber a)const {
282 : const auto p=number2index.find(a);
283 8575206 : if(p==number2index.end()) {
284 : std::string num;
285 0 : Tools::convert( a.serial(), num );
286 0 : plumed_merror("Name of atom " + num + " not found" );
287 : } else {
288 8575206 : return atomsymb[p->second];
289 : }
290 : }
291 :
292 175769 : int PDB::getResidueNumber(AtomNumber a)const {
293 : const auto p=number2index.find(a);
294 175769 : if(p==number2index.end()) {
295 : std::string num;
296 0 : Tools::convert( a.serial(), num );
297 0 : plumed_merror("Residue for atom " + num + " not found" );
298 : } else {
299 175769 : return residue[p->second];
300 : }
301 : }
302 :
303 252980 : std::string PDB::getResidueName(AtomNumber a) const {
304 : const auto p=number2index.find(a);
305 252980 : if(p==number2index.end()) {
306 : std::string num;
307 0 : Tools::convert( a.serial(), num );
308 0 : plumed_merror("Residue for atom " + num + " not found" );
309 : } else {
310 252980 : return residuenames[p->second];
311 : }
312 : }
313 :
314 248924144 : unsigned PDB::size()const {
315 248924144 : return positions.size();
316 : }
317 :
318 5087 : bool PDB::readFromFilepointer(FILE *fp,bool naturalUnits,double scale) {
319 : bool file_is_alive=false;
320 5087 : if(naturalUnits) {
321 : scale=1.0;
322 : }
323 : std::string line;
324 : fpos_t pos;
325 : bool between_ters=true;
326 516785 : while(Tools::getline(fp,line)) {
327 516260 : fgetpos (fp,&pos);
328 3929523 : while(line.length()<80) {
329 3413263 : line.push_back(' ');
330 : }
331 516260 : std::string record=line.substr(0,6);
332 516260 : std::string serial=line.substr(6,5);
333 516260 : std::string atomname=line.substr(12,4);
334 516260 : std::string residuename=line.substr(17,3);
335 516260 : std::string chainID=line.substr(21,1);
336 516260 : std::string resnum=line.substr(22,4);
337 516260 : std::string x=line.substr(30,8);
338 516260 : std::string y=line.substr(38,8);
339 516260 : std::string z=line.substr(46,8);
340 516260 : std::string occ=line.substr(54,6);
341 516260 : std::string bet=line.substr(60,6);
342 516260 : std::string BoxX=line.substr(6,9);
343 516260 : std::string BoxY=line.substr(15,9);
344 516260 : std::string BoxZ=line.substr(24,9);
345 516260 : std::string BoxA=line.substr(33,7);
346 516260 : std::string BoxB=line.substr(40,7);
347 516260 : std::string BoxG=line.substr(47,7);
348 516260 : Tools::trim(record);
349 516260 : if(record=="TER") {
350 : between_ters=false;
351 311 : block_ends.push_back( positions.size() );
352 : }
353 516260 : if(record=="END") {
354 : file_is_alive=true;
355 : break;
356 : }
357 511864 : if(record=="ENDMDL") {
358 : file_is_alive=true;
359 : break;
360 : }
361 511699 : if(record=="REMARK") {
362 4620 : std::vector<std::string> v1;
363 9240 : v1=Tools::getWords(line.substr(6));
364 4620 : addRemark( v1 );
365 4620 : }
366 511699 : if(record=="CRYST1") {
367 112 : Tools::convert(BoxX,BoxXYZ[0]);
368 112 : Tools::convert(BoxY,BoxXYZ[1]);
369 112 : Tools::convert(BoxZ,BoxXYZ[2]);
370 112 : Tools::convert(BoxA,BoxABG[0]);
371 112 : Tools::convert(BoxB,BoxABG[1]);
372 112 : Tools::convert(BoxG,BoxABG[2]);
373 112 : spaceGroup=line.substr(55,11);
374 112 : Tools::trim(spaceGroup);
375 : BoxXYZ*=scale;
376 112 : double cosA=std::cos(BoxABG[0]*pi/180.);
377 112 : double cosB=std::cos(BoxABG[1]*pi/180.);
378 112 : double cosG=std::cos(BoxABG[2]*pi/180.);
379 112 : double sinG=std::sin(BoxABG[2]*pi/180.);
380 448 : for (unsigned i=0; i<3; i++) {
381 336 : Box[i][0]=0.;
382 336 : Box[i][1]=0.;
383 336 : Box[i][2]=0.;
384 : }
385 112 : Box[0][0]=BoxXYZ[0];
386 112 : Box[1][0]=BoxXYZ[1]*cosG;
387 112 : Box[1][1]=BoxXYZ[1]*sinG;
388 112 : Box[2][0]=BoxXYZ[2]*cosB;
389 112 : Box[2][1]=(BoxXYZ[2]*BoxXYZ[1]*cosA-Box[2][0]*Box[1][0])/Box[1][1];
390 112 : Box[2][2]=std::sqrt(BoxXYZ[2]*BoxXYZ[2]-Box[2][0]*Box[2][0]-Box[2][1]*Box[2][1]);
391 : }
392 511699 : if(record=="ATOM" || record=="HETATM") {
393 : between_ters=true;
394 : AtomNumber a;
395 504907 : int resno=0; // GB: when resnum string is not present, we set res number to zero
396 : double o,b;
397 504907 : Vector p;
398 : {
399 : int result;
400 504907 : auto trimmed=serial;
401 504907 : Tools::trim(trimmed);
402 504969 : while(trimmed.length()<5) {
403 62 : trimmed = std::string(" ") + trimmed;
404 : }
405 504907 : const char* errmsg = h36::hy36decode(5, trimmed.c_str(),trimmed.length(), &result);
406 504907 : if(errmsg) {
407 0 : std::string msg(errmsg);
408 0 : plumed_merror(msg);
409 : }
410 504907 : a.setSerial(result);
411 : }
412 :
413 : // allow skipping residue number
414 : {
415 504907 : auto trimmed=resnum;
416 504907 : Tools::trim(trimmed);
417 504907 : if(trimmed.length()>0) {
418 : int result;
419 504907 : while(trimmed.length()<4) {
420 0 : trimmed = std::string(" ") + trimmed;
421 : }
422 504907 : const char* errmsg = h36::hy36decode(4, trimmed.c_str(),trimmed.length(), &result);
423 504907 : if(errmsg) {
424 0 : std::string msg(errmsg);
425 0 : plumed_merror(msg);
426 : }
427 504907 : resno=result;
428 : }
429 : }
430 :
431 : // a field that cannot be read is an error that shows the line (it was an assertion that showed only the text)
432 2524531 : auto readField=[&line](const std::string& field, double& val, const std::string& what) {
433 2524531 : if( !Tools::convertNoexcept(field,val) ) {
434 4 : plumed_merror("cannot read the " + what + " from the following line of a PDB file: in a PDB file each field must be in its own columns\n"
435 : + line.substr(0,line.find_last_not_of(' ')+1));
436 : }
437 2524530 : };
438 504907 : readField(occ,o,"occupancy (columns 55-60)");
439 504906 : readField(bet,b,"beta (columns 61-66)");
440 504906 : readField(x,p[0],"x coordinate (columns 31-38)");
441 504906 : readField(y,p[1],"y coordinate (columns 39-46)");
442 1009812 : readField(z,p[2],"z coordinate (columns 47-54)");
443 : // scale into nm
444 : p*=scale;
445 504906 : numbers.push_back(a);
446 504906 : number2index[a]=positions.size();
447 504906 : std::size_t startpos=atomname.find_first_not_of(" \t");
448 504906 : std::size_t endpos=atomname.find_last_not_of(" \t");
449 504907 : atomsymb.push_back( atomname.substr(startpos, endpos-startpos+1) );
450 504906 : residue.push_back(resno);
451 504906 : chain.push_back(chainID);
452 504906 : occupancy.push_back(o);
453 504906 : beta.push_back(b);
454 504906 : positions.push_back(p);
455 504906 : residuenames.push_back(residuename);
456 : }
457 : }
458 5086 : if( between_ters ) {
459 4844 : block_ends.push_back( positions.size() );
460 : }
461 5086 : return file_is_alive;
462 124 : }
463 :
464 549 : bool PDB::read(const std::string&file,bool naturalUnits,double scale) {
465 549 : FILE* fp=std::fopen(file.c_str(),"r");
466 549 : if(!fp) {
467 : return false;
468 : }
469 : // call fclose when exiting this function
470 : auto deleter=[](auto f) {
471 548 : std::fclose(f);
472 : };
473 : std::unique_ptr<FILE,decltype(deleter)> fp_deleter(fp,deleter);
474 548 : readFromFilepointer(fp,naturalUnits,scale);
475 : return true;
476 : }
477 :
478 429 : void PDB::getChainNames( std::vector<std::string>& chains ) const {
479 429 : chains.resize(0);
480 429 : if(chain.empty()) {
481 : return;
482 : }
483 429 : chains.push_back( chain[0] );
484 648368 : for(unsigned i=1; i<size(); ++i) {
485 647939 : if( chains[chains.size()-1]!=chain[i] ) {
486 979 : chains.push_back( chain[i] );
487 : }
488 : }
489 : }
490 :
491 11720 : void PDB::getResidueRange( const std::string& chainname, int& res_start, int& res_end, std::string& errmsg ) const {
492 : bool inres=false, foundchain=false;
493 29683704 : for(unsigned i=0; i<size(); ++i) {
494 29671984 : if( chain[i]==chainname ) {
495 26606314 : if(!inres) {
496 11720 : if(foundchain) {
497 0 : errmsg="found second start of chain named " + chainname;
498 : }
499 11720 : res_start=residue[i];
500 : }
501 : inres=true;
502 : foundchain=true;
503 3065670 : } else if( inres && chain[i]!=chainname ) {
504 : inres=false;
505 11010 : res_end=residue[i-1];
506 : }
507 : }
508 11720 : if(inres) {
509 710 : res_end=residue[size()-1];
510 : }
511 11720 : }
512 :
513 782 : std::vector<int> PDB::getResiduesInChain( const std::string& chainname ) const {
514 782 : std::vector<int> res;
515 : std::set<int> seen;
516 1889718 : for(unsigned i=0; i<size(); ++i) {
517 1888936 : if( chain[i]==chainname && seen.insert(residue[i]).second ) {
518 14546 : res.push_back(residue[i]);
519 : }
520 : }
521 782 : return res;
522 0 : }
523 :
524 314 : void PDB::getAtomRange( const std::string& chainname, AtomNumber& a_start, AtomNumber& a_end, std::string& errmsg ) const {
525 : bool inres=false, foundchain=false;
526 659346 : for(unsigned i=0; i<size(); ++i) {
527 659032 : if( chain[i]==chainname ) {
528 135598 : if(!inres) {
529 314 : if(foundchain) {
530 0 : errmsg="found second start of chain named " + chainname;
531 : }
532 314 : a_start=numbers[i];
533 : }
534 : inres=true;
535 : foundchain=true;
536 523434 : } else if( inres && chain[i]!=chainname ) {
537 : inres=false;
538 150 : a_end=numbers[i-1];
539 : }
540 : }
541 314 : if(inres) {
542 164 : a_end=numbers[size()-1];
543 : }
544 314 : }
545 :
546 0 : std::string PDB::getResidueName( int resnum ) const {
547 0 : for(unsigned i=0; i<size(); ++i) {
548 0 : if( residue[i]==resnum ) {
549 0 : return residuenames[i];
550 : }
551 : }
552 : std::string num;
553 0 : Tools::convert( resnum, num );
554 0 : plumed_merror("residue " + num + " not found" );
555 : }
556 :
557 65098 : std::string PDB::getResidueName(int resnum,const std::string& chainid ) const {
558 81655662 : for(unsigned i=0; i<size(); ++i) {
559 81655662 : if( residue[i]==resnum && ( chainid=="*" || chain[i]==chainid) ) {
560 65098 : return residuenames[i];
561 : }
562 : }
563 : std::string num;
564 0 : Tools::convert( resnum, num );
565 0 : plumed_merror("residue " + num + " not found in chain " + chainid );
566 : }
567 :
568 :
569 0 : AtomNumber PDB::getNamedAtomFromResidue( const std::string& aname, int resnum ) const {
570 0 : for(unsigned i=0; i<size(); ++i) {
571 0 : if( residue[i]==resnum && atomsymb[i]==aname ) {
572 0 : return numbers[i];
573 : }
574 : }
575 : std::string num;
576 0 : Tools::convert( resnum, num );
577 0 : plumed_merror("residue " + num + " does not contain an atom named " + aname );
578 : }
579 :
580 31164 : AtomNumber PDB::getNamedAtomFromResidueAndChain( const std::string& aname, int resnum, const std::string& chainid ) const {
581 34828884 : for(unsigned i=0; i<size(); ++i) {
582 34828884 : if( residue[i]==resnum && atomsymb[i]==aname && ( chainid=="*" || chain[i]==chainid) ) {
583 31164 : return numbers[i];
584 : }
585 : }
586 : std::string num;
587 0 : Tools::convert( resnum, num );
588 0 : plumed_merror("residue " + num + " from chain " + chainid + " does not contain an atom named " + aname );
589 : }
590 :
591 37322 : std::vector<AtomNumber> PDB::getAtomsInResidue(int resnum,const std::string& chainid)const {
592 37322 : std::vector<AtomNumber> tmp;
593 98889080 : for(unsigned i=0; i<size(); ++i) {
594 98851758 : if( residue[i]==resnum && ( chainid=="*" || chain[i]==chainid) ) {
595 571162 : tmp.push_back(numbers[i]);
596 : }
597 : }
598 37322 : if(tmp.size()==0) {
599 : std::string num;
600 0 : Tools::convert( resnum, num );
601 0 : plumed_merror("Cannot find residue " + num + " from chain " + chainid );
602 : }
603 37322 : return tmp;
604 0 : }
605 :
606 0 : std::vector<AtomNumber> PDB::getAtomsInChain(const std::string& chainid)const {
607 0 : std::vector<AtomNumber> tmp;
608 0 : for(unsigned i=0; i<size(); ++i) {
609 0 : if( chainid=="*" || chain[i]==chainid ) {
610 0 : tmp.push_back(numbers[i]);
611 : }
612 : }
613 0 : if(tmp.size()==0) {
614 0 : plumed_merror("Cannot find atoms from chain " + chainid );
615 : }
616 0 : return tmp;
617 0 : }
618 :
619 0 : std::string PDB::getChainID(int resnumber) const {
620 0 : for(unsigned i=0; i<size(); ++i) {
621 0 : if(resnumber==residue[i]) {
622 0 : return chain[i];
623 : }
624 : }
625 0 : plumed_merror("Not enough residues in pdb input file");
626 : }
627 :
628 54282 : std::string PDB::getChainID(AtomNumber a) const {
629 : const auto p=number2index.find(a);
630 54282 : if(p==number2index.end()) {
631 : std::string num;
632 0 : Tools::convert( a.serial(), num );
633 0 : plumed_merror("Chain for atom " + num + " not found" );
634 : }
635 54282 : return chain[p->second];
636 : }
637 :
638 0 : bool PDB::checkForResidue( const std::string& name ) const {
639 0 : for(unsigned i=0; i<size(); ++i) {
640 0 : if( residuenames[i]==name ) {
641 : return true;
642 : }
643 : }
644 : return false;
645 : }
646 :
647 0 : bool PDB::checkForAtom( const std::string& name ) const {
648 0 : for(unsigned i=0; i<size(); ++i) {
649 0 : if( atomsymb[i]==name ) {
650 : return true;
651 : }
652 : }
653 : return false;
654 : }
655 :
656 3846 : bool PDB::checkForAtom( AtomNumber a ) const {
657 : const auto p=number2index.find(a);
658 3846 : return (p!=number2index.end());
659 : }
660 :
661 0 : Log& operator<<(Log& ostr, const PDB& pdb) {
662 : const std::size_t bufferlen=1000;
663 : char buffer[bufferlen];
664 0 : for(unsigned i=0; i<pdb.positions.size(); i++) {
665 0 : std::snprintf(buffer,bufferlen,"ATOM %3u %8.3f %8.3f %8.3f\n",pdb.numbers[i].serial(),pdb.positions[i][0],pdb.positions[i][1],pdb.positions[i][2]);
666 0 : ostr<<buffer;
667 : }
668 0 : return ostr;
669 : }
670 :
671 8645 : Vector PDB::getPosition(AtomNumber a)const {
672 : const auto p=number2index.find(a);
673 8645 : if(p==number2index.end()) {
674 0 : plumed_merror("atom not available");
675 : } else {
676 8645 : return positions[p->second];
677 : }
678 : }
679 :
680 0 : std::vector<std::string> PDB::getArgumentNames()const {
681 0 : return argnames;
682 : }
683 :
684 0 : std::string PDB::getMtype() const {
685 0 : return mtype;
686 : }
687 :
688 0 : void PDB::print( const double& lunits, GenericMolInfo* mymoldat, OFile& ofile, const std::string& fmt ) {
689 0 : if( argnames.size()>0 ) {
690 0 : ofile.printf("REMARK ARG=%s", argnames[0].c_str() );
691 0 : for(unsigned i=1; i<argnames.size(); ++i) {
692 0 : ofile.printf(",%s",argnames[i].c_str() );
693 : }
694 0 : ofile.printf("\n");
695 0 : ofile.printf("REMARK ");
696 : }
697 : std::string descr2;
698 0 : if(fmt.find("-")!=std::string::npos) {
699 0 : descr2="%s=" + fmt + " ";
700 : } else {
701 : // This ensures numbers are left justified (i.e. next to the equals sign
702 0 : std::size_t psign=fmt.find("%");
703 0 : plumed_assert( psign!=std::string::npos );
704 0 : descr2="%s=%-" + fmt.substr(psign+1) + " ";
705 : }
706 0 : for(std::map<std::string,std::vector<double> >::iterator it=arg_data.begin(); it!=arg_data.end(); ++it) {
707 0 : ofile.printf( descr2.c_str(),it->first.c_str(), it->second[0] );
708 : }
709 0 : if( argnames.size()>0 ) {
710 0 : ofile.printf("\n");
711 : }
712 0 : if( !mymoldat ) {
713 0 : for(unsigned i=0; i<positions.size(); ++i) {
714 : std::array<char,6> at;
715 : {
716 0 : const char* msg = h36::hy36encode(5,numbers[i].serial(),&at[0]);
717 0 : plumed_assert(msg==nullptr) << msg;
718 0 : at[5]=0;
719 : }
720 : std::array<char,5> res;
721 : {
722 0 : const char* msg = h36::hy36encode(4,i,&res[0]);
723 0 : plumed_assert(msg==nullptr) << msg;
724 0 : res[4]=0;
725 : }
726 0 : ofile.printf("ATOM %s X RES %s %8.3f%8.3f%8.3f%6.2f%6.2f\n",
727 : &at[0], &res[0],
728 0 : lunits*positions[i][0], lunits*positions[i][1], lunits*positions[i][2],
729 : occupancy[i], beta[i] );
730 : }
731 : } else {
732 0 : for(unsigned i=0; i<positions.size(); ++i) {
733 : std::array<char,6> at;
734 : {
735 0 : const char* msg = h36::hy36encode(5,numbers[i].serial(),&at[0]);
736 0 : plumed_assert(msg==nullptr) << msg;
737 0 : at[5]=0;
738 : }
739 : std::array<char,5> res;
740 : {
741 0 : const char* msg = h36::hy36encode(4,mymoldat->getResidueNumber(numbers[i]),&res[0]);
742 0 : plumed_assert(msg==nullptr) << msg;
743 0 : res[4]=0;
744 : }
745 0 : ofile.printf("ATOM %s %-4s %3s %s %8.3f%8.3f%8.3f%6.2f%6.2f\n",
746 0 : &at[0], mymoldat->getAtomName(numbers[i]).c_str(),
747 0 : mymoldat->getResidueName(numbers[i]).c_str(), &res[0],
748 0 : lunits*positions[i][0], lunits*positions[i][1], lunits*positions[i][2],
749 : occupancy[i], beta[i] );
750 : }
751 : }
752 0 : ofile.printf("END\n");
753 0 : }
754 :
755 74492 : bool PDB::allowedResidue( const std::string& type, const std::string& residuename ) const {
756 74492 : if( type=="protein" ) {
757 74492 : if(residuename=="ALA") {
758 : return true;
759 70532 : } else if(residuename=="ARG") {
760 : return true;
761 66052 : } else if(residuename=="ASN") {
762 : return true;
763 62212 : } else if(residuename=="ASP") {
764 : return true;
765 58852 : } else if(residuename=="CYS") {
766 : return true;
767 58592 : } else if(residuename=="GLN") {
768 : return true;
769 53112 : } else if(residuename=="GLU") {
770 : return true;
771 50592 : } else if(residuename=="GLY") {
772 : return true;
773 47072 : } else if(residuename=="HIS") {
774 : return true;
775 47072 : } else if(residuename=="ILE") {
776 : return true;
777 42572 : } else if(residuename=="LEU") {
778 : return true;
779 34592 : } else if(residuename=="LYS") {
780 : return true;
781 29192 : } else if(residuename=="MET") {
782 : return true;
783 26232 : } else if(residuename=="PHE") {
784 : return true;
785 20952 : } else if(residuename=="PRO") {
786 : return true;
787 16992 : } else if(residuename=="SER") {
788 : return true;
789 14210 : } else if(residuename=="THR") {
790 : return true;
791 11600 : } else if(residuename=="TRP") {
792 : return true;
793 11600 : } else if(residuename=="TYR") {
794 : return true;
795 8840 : } else if(residuename=="VAL") {
796 : return true;
797 : }
798 : // Terminal groups
799 3162 : else if(residuename=="ACE") {
800 : return true;
801 3162 : } else if(residuename=="NME") {
802 : return true;
803 3162 : } else if(residuename=="NH2") {
804 : return true;
805 : }
806 : // Alternative residue names in common force fields
807 3162 : else if(residuename=="GLH") {
808 : return true; // neutral GLU
809 3162 : } else if(residuename=="ASH") {
810 : return true; // neutral ASP
811 3162 : } else if(residuename=="HID") {
812 : return true; // HIS-D amber
813 3162 : } else if(residuename=="HSD") {
814 : return true; // HIS-D charmm
815 3162 : } else if(residuename=="HIE") {
816 : return true; // HIS-E amber
817 2402 : } else if(residuename=="HSE") {
818 : return true; // HIS-E charmm
819 2402 : } else if(residuename=="HIP") {
820 : return true; // HIS-P amber
821 2402 : } else if(residuename=="HSP") {
822 : return true; // HIS-P charmm
823 2402 : } else if(residuename=="CYX") {
824 : return true; // disulfide bridge CYS
825 : }
826 : // Weird amino acids
827 2402 : else if(residuename=="NLE") {
828 : return true;
829 2402 : } else if(residuename=="SFO") {
830 : return true;
831 : } else {
832 2402 : return false;
833 : }
834 0 : } else if( type=="dna" ) {
835 0 : if(residuename=="A") {
836 : return true;
837 0 : } else if(residuename=="A5") {
838 : return true;
839 0 : } else if(residuename=="A3") {
840 : return true;
841 0 : } else if(residuename=="AN") {
842 : return true;
843 0 : } else if(residuename=="G") {
844 : return true;
845 0 : } else if(residuename=="G5") {
846 : return true;
847 0 : } else if(residuename=="G3") {
848 : return true;
849 0 : } else if(residuename=="GN") {
850 : return true;
851 0 : } else if(residuename=="T") {
852 : return true;
853 0 : } else if(residuename=="T5") {
854 : return true;
855 0 : } else if(residuename=="T3") {
856 : return true;
857 0 : } else if(residuename=="TN") {
858 : return true;
859 0 : } else if(residuename=="C") {
860 : return true;
861 0 : } else if(residuename=="C5") {
862 : return true;
863 0 : } else if(residuename=="C3") {
864 : return true;
865 0 : } else if(residuename=="CN") {
866 : return true;
867 0 : } else if(residuename=="DA") {
868 : return true;
869 0 : } else if(residuename=="DA5") {
870 : return true;
871 0 : } else if(residuename=="DA3") {
872 : return true;
873 0 : } else if(residuename=="DAN") {
874 : return true;
875 0 : } else if(residuename=="DG") {
876 : return true;
877 0 : } else if(residuename=="DG5") {
878 : return true;
879 0 : } else if(residuename=="DG3") {
880 : return true;
881 0 : } else if(residuename=="DGN") {
882 : return true;
883 0 : } else if(residuename=="DT") {
884 : return true;
885 0 : } else if(residuename=="DT5") {
886 : return true;
887 0 : } else if(residuename=="DT3") {
888 : return true;
889 0 : } else if(residuename=="DTN") {
890 : return true;
891 0 : } else if(residuename=="DC") {
892 : return true;
893 0 : } else if(residuename=="DC5") {
894 : return true;
895 0 : } else if(residuename=="DC3") {
896 : return true;
897 0 : } else if(residuename=="DCN") {
898 : return true;
899 : } else {
900 0 : return false;
901 : }
902 0 : } else if( type=="rna" ) {
903 0 : if(residuename=="A") {
904 : return true;
905 0 : } else if(residuename=="A5") {
906 : return true;
907 0 : } else if(residuename=="A3") {
908 : return true;
909 0 : } else if(residuename=="AN") {
910 : return true;
911 0 : } else if(residuename=="G") {
912 : return true;
913 0 : } else if(residuename=="G5") {
914 : return true;
915 0 : } else if(residuename=="G3") {
916 : return true;
917 0 : } else if(residuename=="GN") {
918 : return true;
919 0 : } else if(residuename=="U") {
920 : return true;
921 0 : } else if(residuename=="U5") {
922 : return true;
923 0 : } else if(residuename=="U3") {
924 : return true;
925 0 : } else if(residuename=="UN") {
926 : return true;
927 0 : } else if(residuename=="C") {
928 : return true;
929 0 : } else if(residuename=="C5") {
930 : return true;
931 0 : } else if(residuename=="C3") {
932 : return true;
933 0 : } else if(residuename=="CN") {
934 : return true;
935 0 : } else if(residuename=="RA") {
936 : return true;
937 0 : } else if(residuename=="RA5") {
938 : return true;
939 0 : } else if(residuename=="RA3") {
940 : return true;
941 0 : } else if(residuename=="RAN") {
942 : return true;
943 0 : } else if(residuename=="RG") {
944 : return true;
945 0 : } else if(residuename=="RG5") {
946 : return true;
947 0 : } else if(residuename=="RG3") {
948 : return true;
949 0 : } else if(residuename=="RGN") {
950 : return true;
951 0 : } else if(residuename=="RU") {
952 : return true;
953 0 : } else if(residuename=="RU5") {
954 : return true;
955 0 : } else if(residuename=="RU3") {
956 : return true;
957 0 : } else if(residuename=="RUN") {
958 : return true;
959 0 : } else if(residuename=="RC") {
960 : return true;
961 0 : } else if(residuename=="RC5") {
962 : return true;
963 0 : } else if(residuename=="RC3") {
964 : return true;
965 0 : } else if(residuename=="RCN") {
966 : return true;
967 : } else {
968 0 : return false;
969 : }
970 0 : } else if( type=="water" ) {
971 0 : if(residuename=="SOL") {
972 : return true;
973 : }
974 0 : if(residuename=="WAT") {
975 : return true;
976 : }
977 0 : return false;
978 0 : } else if( type=="ion" ) {
979 0 : if(residuename=="IB+") {
980 : return true;
981 : }
982 0 : if(residuename=="CA") {
983 : return true;
984 : }
985 0 : if(residuename=="CL") {
986 : return true;
987 : }
988 0 : if(residuename=="NA") {
989 : return true;
990 : }
991 0 : if(residuename=="MG") {
992 : return true;
993 : }
994 0 : if(residuename=="K") {
995 : return true;
996 : }
997 0 : if(residuename=="RB") {
998 : return true;
999 : }
1000 0 : if(residuename=="CS") {
1001 : return true;
1002 : }
1003 0 : if(residuename=="LI") {
1004 : return true;
1005 : }
1006 0 : if(residuename=="ZN") {
1007 : return true;
1008 : }
1009 0 : return false;
1010 : }
1011 : return false;
1012 : }
1013 :
1014 : }
1015 :
|
