Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2011-2019 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 <cstdio>
26 : #include <iostream>
27 :
28 : using namespace std;
29 :
30 : //+PLUMEDOC INTERNAL pdbreader
31 : /*
32 : PLUMED can use the PDB format in several places
33 :
34 : - To read molecular structure (\ref MOLINFO).
35 : - To read reference conformations (\ref RMSD, but also many other methods in \ref dists, \ref FIT_TO_TEMPLATE, etc).
36 :
37 : The implemented PDB reader expects a file formatted correctly according to the
38 : [PDB standard](http://www.wwpdb.org/documentation/file-format-content/format33/v3.3.html).
39 : In particular, the following columns are read from ATOM records
40 : \verbatim
41 : columns | content
42 : 1-6 | record name (ATOM or HETATM)
43 : 7-11 | serial number of the atom (starting from 1)
44 : 13-16 | atom name
45 : 18-20 | residue name
46 : 22 | chain id
47 : 23-26 | residue number
48 : 31-38 | x coordinate
49 : 39-46 | y coordinate
50 : 47-54 | z coordinate
51 : 55-60 | occupancy
52 : 61-66 | beta factor
53 : \endverbatim
54 : PLUMED parser is slightly more permissive than the official PDB format
55 : in the fact that the format of real numbers is not fixed. In other words,
56 : any parsable real number is ok and the dot can be placed anywhere. However,
57 : __columns are interpret strictly__. A sample PDB should look like the following
58 : \verbatim
59 : ATOM 2 CH3 ACE 1 12.932 -14.718 -6.016 1.00 1.00
60 : ATOM 5 C ACE 1 21.312 -9.928 -5.946 1.00 1.00
61 : ATOM 9 CA ALA 2 19.462 -11.088 -8.986 1.00 1.00
62 : \endverbatim
63 :
64 : Notice that serial numbers need not to be consecutive. In the three-line example above,
65 : only the coordinates of three atoms are provided. This is perfectly legal and indicates PLUMED
66 : that information about these atoms only is available. This could be both for structural
67 : information in \ref MOLINFO, where the other atoms would have no name assigned, and for
68 : reference structures used in \ref RMSD, where only the provided atoms would be used to compute RMSD.
69 :
70 : \par Occupancy and beta factors
71 :
72 : PLUMED reads also occupancy and beta factors that however are given a very special meaning.
73 : In cases where the PDB structure is used as a reference for an alignment (that's the case
74 : for instance in \ref RMSD and in \ref FIT_TO_TEMPLATE), the occupancy column is used
75 : to provide the weight of each atom in the alignment. In cases where, perhaps after alignment,
76 : the displacement between running coordinates and the provided PDB is computed, the beta factors
77 : are used as weight for the displacement.
78 : Since setting the weights to zero is the same as __not__ including an atom in the alignement or
79 : displacement calculation, the two following reference files would be equivalent when used in an \ref RMSD
80 : calculation. First file:
81 : \verbatim
82 : ATOM 2 CH3 ACE 1 12.932 -14.718 -6.016 1.00 1.00
83 : ATOM 5 C ACE 1 21.312 -9.928 -5.946 1.00 1.00
84 : ATOM 9 CA ALA 2 19.462 -11.088 -8.986 0.00 0.00
85 : \endverbatim
86 : Second file:
87 : \verbatim
88 : ATOM 2 CH3 ACE 1 12.932 -14.718 -6.016 1.00 1.00
89 : ATOM 5 C ACE 1 21.312 -9.928 -5.946 1.00 1.00
90 : \endverbatim
91 : However notice that many extra atoms with zero weight might slow down the calculation, so
92 : removing lines is better than setting their weights to zero.
93 : In addition, weights for alignment need not to be equivalent to weights for displacement.
94 :
95 : \par Systems with more than 100k atoms
96 :
97 : Notice that it very likely does not make any sense to compute the \ref RMSD or any other structural
98 : deviation __using__ so many atoms. However, if the protein for which you want to compute \ref RMSD
99 : has atoms with large serial numbers (e.g. because it is located __after__ solvent in the sorted list of atoms)
100 : you might end up with troubles with the limitations of the PDB format. Indeed, since there are 5
101 : columns available for atom serial number, this number cannot be larger than 99999.
102 : In addition, providing \ref MOLINFO with names associated to atoms with a serial larger than 99999 would be impossible.
103 :
104 : Since PLUMED 2.4 we allow [hybrid 36](http://cci.lbl.gov/hybrid_36/) format to be used to specify atom numbers.
105 : This format is not particularly widespread, but has the nice feature that it provides a one-to-one mapping
106 : between numbers up to approximately 80 millions and strings with 5 characters, plus it is backward compatible
107 : for numbers smaller than 100000. This is not true for notations like the hex notation exported by VMD.
108 : Using the hybrid 36 format, the ATOM records for atom ranging from 99997 to 100002 would read like these:
109 : \verbatim
110 : ATOM 99997 Ar X 1 45.349 38.631 15.116 1.00 1.00
111 : ATOM 99998 Ar X 1 46.189 38.631 15.956 1.00 1.00
112 : ATOM 99999 Ar X 1 46.189 39.471 15.116 1.00 1.00
113 : ATOM A0000 Ar X 1 45.349 39.471 15.956 1.00 1.00
114 : ATOM A0000 Ar X 1 45.349 38.631 16.796 1.00 1.00
115 : ATOM A0001 Ar X 1 46.189 38.631 17.636 1.00 1.00
116 : \endverbatim
117 : There are tools that can be found to translate from integers to strings and back using hybrid 36 format
118 : (a simple python script can be found [here](https://sourceforge.net/p/cctbx/code/HEAD/tree/trunk/iotbx/pdb/hybrid_36.py)).
119 :
120 : */
121 : //+ENDPLUMEDOC
122 :
123 :
124 : namespace PLMD {
125 :
126 : /// Tiny namespace for hybrid36 format.
127 : /// This namespace includes freely available tools for h36 format.
128 : /// I place them here for usage within PDB class. In case we need them
129 : /// in other places they might be better encapsulated in a c++ class
130 : /// and placed in a separate file.
131 : namespace h36 {
132 :
133 :
134 : /*! C port of the hy36encode() and hy36decode() functions in the
135 : hybrid_36.py Python prototype/reference implementation.
136 : See the Python script for more information.
137 :
138 : This file has no external dependencies, NOT even standard C headers.
139 : Optionally, use hybrid_36_c.h, or simply copy the declarations
140 : into your code.
141 :
142 : This file is unrestricted Open Source (cctbx.sf.net).
143 : Please send corrections and enhancements to cctbx@cci.lbl.gov .
144 :
145 : See also: http://cci.lbl.gov/hybrid_36/
146 :
147 : Ralf W. Grosse-Kunstleve, Feb 2007.
148 : */
149 :
150 : /* The following #include may be commented out.
151 : It is here only to enforce consistency of the declarations
152 : and the definitions.
153 : */
154 : // #include <iotbx/pdb/hybrid_36_c.h>
155 :
156 : /* All static functions below are implementation details
157 : (and not accessible from other translation units).
158 : */
159 :
160 : static
161 : const char*
162 5796 : digits_upper() { return "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ"; }
163 :
164 : static
165 : const char*
166 5796 : digits_lower() { return "0123456789abcdefghijklmnopqrstuvwxyz"; }
167 :
168 : static
169 : const char*
170 0 : value_out_of_range() { return "value out of range."; }
171 :
172 : static
173 0 : const char* invalid_number_literal() { return "invalid number literal."; }
174 :
175 : static
176 0 : const char* unsupported_width() { return "unsupported width."; }
177 :
178 : static
179 : void
180 0 : fill_with_stars(unsigned width, char* result)
181 : {
182 0 : while (width) {
183 0 : *result++ = '*';
184 0 : width--;
185 : }
186 0 : *result = '\0';
187 0 : }
188 :
189 : static
190 : void
191 0 : encode_pure(
192 : const char* digits,
193 : unsigned digits_size,
194 : unsigned width,
195 : int value,
196 : char* result)
197 : {
198 : char buf[16];
199 : int rest;
200 : unsigned i, j;
201 : i = 0;
202 : j = 0;
203 0 : if (value < 0) {
204 : j = 1;
205 0 : value = -value;
206 : }
207 : while (1) {
208 0 : rest = value / digits_size;
209 0 : buf[i++] = digits[value - rest * digits_size];
210 0 : if (rest == 0) break;
211 : value = rest;
212 : }
213 0 : if (j) buf[i++] = '-';
214 0 : for(j=i; j<width; j++) *result++ = ' ';
215 0 : while (i != 0) *result++ = buf[--i];
216 0 : *result = '\0';
217 0 : }
218 :
219 : static
220 : const char*
221 96733 : decode_pure(
222 : const int* digits_values,
223 : unsigned digits_size,
224 : const char* s,
225 : unsigned s_size,
226 : int* result)
227 : {
228 : int si, dv;
229 : int have_minus = 0;
230 : int have_non_blank = 0;
231 : int value = 0;
232 : unsigned i = 0;
233 1064063 : for(; i<s_size; i++) {
234 483665 : si = s[i];
235 483665 : if (si < 0 || si > 127) {
236 0 : *result = 0;
237 0 : return invalid_number_literal();
238 : }
239 483665 : if (si == ' ') {
240 199160 : if (!have_non_blank) continue;
241 0 : value *= digits_size;
242 : }
243 284505 : else if (si == '-') {
244 0 : if (have_non_blank) {
245 0 : *result = 0;
246 0 : return invalid_number_literal();
247 : }
248 : have_non_blank = 1;
249 : have_minus = 1;
250 0 : continue;
251 : }
252 : else {
253 : have_non_blank = 1;
254 284505 : dv = digits_values[si];
255 284505 : if (dv < 0 || dv >= digits_size) {
256 0 : *result = 0;
257 0 : return invalid_number_literal();
258 : }
259 284505 : value *= digits_size;
260 284505 : value += dv;
261 : }
262 : }
263 96733 : if (have_minus) value = -value;
264 96733 : *result = value;
265 96733 : return 0;
266 : }
267 :
268 : /*! hybrid-36 encoder: converts integer value to string result
269 :
270 : width: must be 4 (e.g. for residue sequence numbers)
271 : or 5 (e.g. for atom serial numbers)
272 :
273 : value: integer value to be converted
274 :
275 : result: pointer to char array of size width+1 or greater
276 : on return result is null-terminated
277 :
278 : return value: pointer to error message, if any,
279 : or 0 on success
280 :
281 : Example usage (from C++):
282 : char result[4+1];
283 : const char* errmsg = hy36encode(4, 12345, result);
284 : if (errmsg) throw std::runtime_error(errmsg);
285 : */
286 : const char*
287 0 : hy36encode(unsigned width, int value, char* result)
288 : {
289 : int i = value;
290 0 : if (width == 4U) {
291 0 : if (i >= -999) {
292 0 : if (i < 10000) {
293 0 : encode_pure(digits_upper(), 10U, 4U, i, result);
294 0 : return 0;
295 : }
296 0 : i -= 10000;
297 0 : if (i < 1213056 /* 26*36**3 */) {
298 0 : i += 466560 /* 10*36**3 */;
299 0 : encode_pure(digits_upper(), 36U, 0U, i, result);
300 0 : return 0;
301 : }
302 0 : i -= 1213056;
303 0 : if (i < 1213056) {
304 0 : i += 466560;
305 0 : encode_pure(digits_lower(), 36U, 0U, i, result);
306 0 : return 0;
307 : }
308 : }
309 : }
310 0 : else if (width == 5U) {
311 0 : if (i >= -9999) {
312 0 : if (i < 100000) {
313 0 : encode_pure(digits_upper(), 10U, 5U, i, result);
314 0 : return 0;
315 : }
316 0 : i -= 100000;
317 0 : if (i < 43670016 /* 26*36**4 */) {
318 0 : i += 16796160 /* 10*36**4 */;
319 0 : encode_pure(digits_upper(), 36U, 0U, i, result);
320 0 : return 0;
321 : }
322 0 : i -= 43670016;
323 0 : if (i < 43670016) {
324 0 : i += 16796160;
325 0 : encode_pure(digits_lower(), 36U, 0U, i, result);
326 0 : return 0;
327 : }
328 : }
329 : }
330 : else {
331 0 : fill_with_stars(width, result);
332 0 : return unsupported_width();
333 : }
334 0 : fill_with_stars(width, result);
335 0 : return value_out_of_range();
336 : }
337 :
338 : /*! hybrid-36 decoder: converts string s to integer result
339 :
340 : width: must be 4 (e.g. for residue sequence numbers)
341 : or 5 (e.g. for atom serial numbers)
342 :
343 : s: string to be converted
344 : does not have to be null-terminated
345 :
346 : s_size: size of s
347 : must be equal to width, or an error message is
348 : returned otherwise
349 :
350 : result: integer holding the conversion result
351 :
352 : return value: pointer to error message, if any,
353 : or 0 on success
354 :
355 : Example usage (from C++):
356 : int result;
357 : const char* errmsg = hy36decode(width, "A1T5", 4, &result);
358 : if (errmsg) throw std::runtime_error(errmsg);
359 : */
360 : const char*
361 96733 : hy36decode(unsigned width, const char* s, unsigned s_size, int* result)
362 : {
363 : static int first_call = 1;
364 : static int digits_values_upper[128U];
365 : static int digits_values_lower[128U];
366 : static const char*
367 : ie_range = "internal error hy36decode: integer value out of range.";
368 : unsigned i;
369 : int di;
370 : const char* errmsg;
371 96733 : if (first_call) {
372 161 : first_call = 0;
373 20769 : for(i=0; i<128U; i++) digits_values_upper[i] = -1;
374 20769 : for(i=0; i<128U; i++) digits_values_lower[i] = -1;
375 11753 : for(i=0; i<36U; i++) {
376 5796 : di = digits_upper()[i];
377 5796 : if (di < 0 || di > 127) {
378 0 : *result = 0;
379 0 : return ie_range;
380 : }
381 5796 : digits_values_upper[di] = i;
382 : }
383 11753 : for(i=0; i<36U; i++) {
384 5796 : di = digits_lower()[i];
385 5796 : if (di < 0 || di > 127) {
386 0 : *result = 0;
387 0 : return ie_range;
388 : }
389 5796 : digits_values_lower[di] = i;
390 : }
391 : }
392 96733 : if (s_size == width) {
393 96733 : di = s[0];
394 96733 : if (di >= 0 && di <= 127) {
395 96733 : if (digits_values_upper[di] >= 10) {
396 2 : errmsg = decode_pure(digits_values_upper, 36U, s, s_size, result);
397 2 : if (errmsg == 0) {
398 : /* result - 10*36**(width-1) + 10**width */
399 2 : if (width == 4U) (*result) -= 456560;
400 2 : else if (width == 5U) (*result) -= 16696160;
401 : else {
402 0 : *result = 0;
403 0 : return unsupported_width();
404 : }
405 : return 0;
406 : }
407 : }
408 96731 : else if (digits_values_lower[di] >= 10) {
409 0 : errmsg = decode_pure(digits_values_lower, 36U, s, s_size, result);
410 0 : if (errmsg == 0) {
411 : /* result + 16*36**(width-1) + 10**width */
412 0 : if (width == 4U) (*result) += 756496;
413 0 : else if (width == 5U) (*result) += 26973856;
414 : else {
415 0 : *result = 0;
416 0 : return unsupported_width();
417 : }
418 : return 0;
419 : }
420 : }
421 : else {
422 96731 : errmsg = decode_pure(digits_values_upper, 10U, s, s_size, result);
423 96731 : if (errmsg) return errmsg;
424 96731 : if (!(width == 4U || width == 5U)) {
425 0 : *result = 0;
426 0 : return unsupported_width();
427 : }
428 : return 0;
429 : }
430 : }
431 : }
432 0 : *result = 0;
433 0 : return invalid_number_literal();
434 : }
435 :
436 : }
437 :
438 453 : unsigned PDB::getNumberOfAtomBlocks()const {
439 453 : return block_ends.size();
440 : }
441 :
442 18 : const std::vector<unsigned> & PDB::getAtomBlockEnds()const {
443 18 : return block_ends;
444 : }
445 :
446 8956 : const std::vector<Vector> & PDB::getPositions()const {
447 8956 : return positions;
448 : }
449 :
450 0 : void PDB::setPositions(const std::vector<Vector> &v ) {
451 0 : plumed_assert( v.size()==positions.size() );
452 0 : positions=v;
453 0 : }
454 :
455 14314 : const std::vector<double> & PDB::getOccupancy()const {
456 14314 : return occupancy;
457 : }
458 :
459 14314 : const std::vector<double> & PDB::getBeta()const {
460 14314 : return beta;
461 : }
462 :
463 1465 : const std::vector<std::string> & PDB::getRemark()const {
464 1465 : return remark;
465 : }
466 :
467 914 : void PDB::addRemark( const std::vector<std::string>& v1 ) {
468 1828 : remark.insert(remark.begin(),v1.begin(),v1.end());
469 914 : }
470 :
471 20674 : const std::vector<AtomNumber> & PDB::getAtomNumbers()const {
472 20674 : return numbers;
473 : }
474 :
475 1441030 : std::string PDB::getAtomName(AtomNumber a)const {
476 : const auto p=number2index.find(a);
477 1441030 : if(p==number2index.end()) return "";
478 1441028 : else return atomsymb[p->second];
479 : }
480 :
481 76420 : unsigned PDB::getResidueNumber(AtomNumber a)const {
482 : const auto p=number2index.find(a);
483 76420 : if(p==number2index.end()) return 0;
484 152836 : else return residue[p->second];
485 : }
486 :
487 79380 : std::string PDB::getResidueName(AtomNumber a) const {
488 : const auto p=number2index.find(a);
489 79380 : if(p==number2index.end()) return "";
490 79378 : else return residuenames[p->second];
491 : }
492 :
493 74494064 : unsigned PDB::size()const {
494 74494064 : return positions.size();
495 : }
496 :
497 1790 : bool PDB::readFromFilepointer(FILE *fp,bool naturalUnits,double scale) {
498 : //cerr<<file<<endl;
499 : bool file_is_alive=false;
500 1790 : if(naturalUnits) scale=1.0;
501 : string line;
502 : fpos_t pos; bool between_ters=true;
503 99791 : while(Tools::getline(fp,line)) {
504 : //cerr<<line<<"\n";
505 99621 : fgetpos (fp,&pos);
506 772109 : while(line.length()<80) line.push_back(' ');
507 99621 : string record=line.substr(0,6);
508 99621 : string serial=line.substr(6,5);
509 99621 : string atomname=line.substr(12,4);
510 99621 : string residuename=line.substr(17,3);
511 99621 : string chainID=line.substr(21,1);
512 99621 : string resnum=line.substr(22,4);
513 99621 : string x=line.substr(30,8);
514 99621 : string y=line.substr(38,8);
515 99621 : string z=line.substr(46,8);
516 99621 : string occ=line.substr(54,6);
517 99621 : string bet=line.substr(60,6);
518 99621 : Tools::trim(record);
519 99851 : if(record=="TER") { between_ters=false; block_ends.push_back( positions.size() ); }
520 99621 : if(record=="END") { file_is_alive=true; break;}
521 98075 : if(record=="ENDMDL") { file_is_alive=true; break;}
522 98001 : if(record=="REMARK") {
523 3640 : vector<string> v1; v1=Tools::getWords(line.substr(6));
524 910 : addRemark( v1 );
525 : }
526 99269 : if(record=="ATOM" || record=="HETATM") {
527 : between_ters=true;
528 : AtomNumber a; unsigned resno;
529 : double o,b;
530 96733 : Vector p;
531 : {
532 : int result;
533 : auto trimmed=serial;
534 96733 : Tools::trim(trimmed);
535 96826 : while(trimmed.length()<5) trimmed = std::string(" ") + trimmed;
536 193466 : const char* errmsg = h36::hy36decode(5, trimmed.c_str(),trimmed.length(), &result);
537 96733 : if(errmsg) {
538 0 : std::string msg(errmsg);
539 0 : plumed_merror(msg);
540 : }
541 96733 : a.setSerial(result);
542 : }
543 :
544 96733 : Tools::convert(resnum,resno);
545 96733 : Tools::convert(occ,o);
546 96733 : Tools::convert(bet,b);
547 96733 : Tools::convert(x,p[0]);
548 96733 : Tools::convert(y,p[1]);
549 96733 : Tools::convert(z,p[2]);
550 : // scale into nm
551 96733 : p*=scale;
552 96733 : numbers.push_back(a);
553 96733 : number2index[a]=positions.size();
554 : std::size_t startpos=atomname.find_first_not_of(" \t");
555 : std::size_t endpos=atomname.find_last_not_of(" \t");
556 193466 : atomsymb.push_back( atomname.substr(startpos, endpos-startpos+1) );
557 96733 : residue.push_back(resno);
558 96733 : chain.push_back(chainID);
559 96733 : occupancy.push_back(o);
560 96733 : beta.push_back(b);
561 96733 : positions.push_back(p);
562 96733 : residuenames.push_back(residuename);
563 : }
564 : }
565 5182 : if( between_ters ) block_ends.push_back( positions.size() );
566 1790 : return file_is_alive;
567 : }
568 :
569 80 : void PDB::setArgKeyword( const std::string& new_args ) {
570 : bool replaced=false;
571 880 : for(unsigned i=0; i<remark.size(); ++i) {
572 240 : if( remark[i].find("ARG=")!=std::string::npos) {
573 : remark[i]=new_args; replaced=true;
574 : }
575 : }
576 80 : plumed_assert( replaced );
577 80 : }
578 :
579 232 : bool PDB::read(const std::string&file,bool naturalUnits,double scale) {
580 232 : FILE* fp=fopen(file.c_str(),"r");
581 232 : if(!fp) return false;
582 230 : readFromFilepointer(fp,naturalUnits,scale);
583 230 : fclose(fp);
584 230 : return true;
585 : }
586 :
587 132 : void PDB::getChainNames( std::vector<std::string>& chains ) const {
588 132 : chains.resize(0);
589 132 : chains.push_back( chain[0] );
590 428960 : for(unsigned i=1; i<size(); ++i) {
591 643242 : if( chains[chains.size()-1]!=chain[i] ) chains.push_back( chain[i] );
592 : }
593 132 : }
594 :
595 463 : void PDB::getResidueRange( const std::string& chainname, unsigned& res_start, unsigned& res_end, std::string& errmsg ) const {
596 : bool inres=false, foundchain=false;
597 2083995 : for(unsigned i=0; i<size(); ++i) {
598 2083532 : if( chain[i]==chainname ) {
599 104842 : if(!inres) {
600 463 : if(foundchain) errmsg="found second start of chain named " + chainname;
601 463 : res_start=residue[i];
602 : }
603 : inres=true; foundchain=true;
604 936924 : } else if( inres && chain[i]!=chainname ) {
605 : inres=false;
606 746 : res_end=residue[i-1];
607 : }
608 : }
609 553 : if(inres) res_end=residue[size()-1];
610 463 : }
611 :
612 192 : void PDB::getAtomRange( const std::string& chainname, AtomNumber& a_start, AtomNumber& a_end, std::string& errmsg ) const {
613 : bool inres=false, foundchain=false;
614 742516 : for(unsigned i=0; i<size(); ++i) {
615 742324 : if( chain[i]==chainname ) {
616 132222 : if(!inres) {
617 192 : if(foundchain) errmsg="found second start of chain named " + chainname;
618 192 : a_start=numbers[i];
619 : }
620 : inres=true; foundchain=true;
621 238940 : } else if( inres && chain[i]!=chainname ) {
622 : inres=false;
623 190 : a_end=numbers[i-1];
624 : }
625 : }
626 289 : if(inres) a_end=numbers[size()-1];
627 192 : }
628 :
629 8422 : std::string PDB::getResidueName( const unsigned& resnum ) const {
630 21075150 : for(unsigned i=0; i<size(); ++i) {
631 21083572 : if( residue[i]==resnum ) return residuenames[i];
632 : }
633 0 : return "";
634 : }
635 :
636 10302 : std::string PDB::getResidueName(const unsigned& resnum,const std::string& chainid ) const {
637 26230466 : for(unsigned i=0; i<size(); ++i) {
638 26261244 : if( residue[i]==resnum && ( chainid=="*" || chain[i]==chainid) ) return residuenames[i];
639 : }
640 0 : return "";
641 : }
642 :
643 :
644 13260 : AtomNumber PDB::getNamedAtomFromResidue( const std::string& aname, const unsigned& resnum ) const {
645 32681100 : for(unsigned i=0; i<size(); ++i) {
646 32812980 : if( residue[i]==resnum && atomsymb[i]==aname ) return numbers[i];
647 : }
648 0 : std::string num; Tools::convert( resnum, num );
649 0 : plumed_merror("residue " + num + " does not contain an atom named " + aname );
650 : return numbers[0]; // This is to stop compiler errors
651 : }
652 :
653 1891 : AtomNumber PDB::getNamedAtomFromResidueAndChain( const std::string& aname, const unsigned& resnum, const std::string& chainid ) const {
654 2068917 : for(unsigned i=0; i<size(); ++i) {
655 2104726 : if( residue[i]==resnum && atomsymb[i]==aname && ( chainid=="*" || chain[i]==chainid) ) return numbers[i];
656 : }
657 0 : std::string num; Tools::convert( resnum, num );
658 0 : plumed_merror("residue " + num + " from chain " + chainid + " does not contain an atom named " + aname );
659 : return numbers[0]; // This is to stop compiler errors
660 : }
661 :
662 9954 : std::vector<AtomNumber> PDB::getAtomsInResidue(const unsigned& resnum,const std::string& chainid)const {
663 : std::vector<AtomNumber> tmp;
664 52009650 : for(unsigned i=0; i<size(); ++i) {
665 52447710 : if( residue[i]==resnum && ( chainid=="*" || chain[i]==chainid) ) tmp.push_back(numbers[i]);
666 : }
667 9954 : if(tmp.size()==0) {
668 0 : std::string num; Tools::convert( resnum, num );
669 0 : plumed_merror("Cannot find residue " + num + " from chain " + chainid );
670 : }
671 9954 : return tmp;
672 : }
673 :
674 28 : std::vector<AtomNumber> PDB::getAtomsInChain(const std::string& chainid)const {
675 : std::vector<AtomNumber> tmp;
676 146300 : for(unsigned i=0; i<size(); ++i) {
677 182840 : if( chainid=="*" || chain[i]==chainid ) tmp.push_back(numbers[i]);
678 : }
679 28 : if(tmp.size()==0) {
680 0 : plumed_merror("Cannot find atoms from chain " + chainid );
681 : }
682 28 : return tmp;
683 : }
684 :
685 4638 : std::string PDB::getChainID(const unsigned& resnumber) const {
686 11373706 : for(unsigned i=0; i<size(); ++i) {
687 11382982 : if(resnumber==residue[i]) return chain[i];
688 : }
689 0 : plumed_merror("Not enough residues in pdb input file");
690 : }
691 :
692 0 : bool PDB::checkForResidue( const std::string& name ) const {
693 0 : for(unsigned i=0; i<size(); ++i) {
694 0 : if( residuenames[i]==name ) return true;
695 : }
696 : return false;
697 : }
698 :
699 0 : bool PDB::checkForAtom( const std::string& name ) const {
700 0 : for(unsigned i=0; i<size(); ++i) {
701 0 : if( atomsymb[i]==name ) return true;
702 : }
703 : return false;
704 : }
705 :
706 0 : Log& operator<<(Log& ostr, const PDB& pdb) {
707 : char buffer[1000];
708 0 : for(unsigned i=0; i<pdb.positions.size(); i++) {
709 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]);
710 0 : ostr<<buffer;
711 : }
712 0 : return ostr;
713 : }
714 :
715 852 : Vector PDB::getPosition(AtomNumber a)const {
716 : const auto p=number2index.find(a);
717 852 : if(p==number2index.end()) plumed_merror("atom not available");
718 1704 : else return positions[p->second];
719 : }
720 :
721 :
722 :
723 4839 : }
724 :
|
