Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2013-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 "MolDataClass.h"
23 : #include "Exception.h"
24 : #include "Tools.h"
25 : #include "PDB.h"
26 :
27 : namespace PLMD {
28 :
29 60 : unsigned MolDataClass::numberOfAtomsPerResidueInBackbone( const std::string& type ) {
30 60 : if( type=="protein" ) {
31 : return 5;
32 0 : } else if( type=="dna" ) {
33 : return 6;
34 0 : } else if( type=="rna" ) {
35 : return 6;
36 : } else {
37 0 : return 0;
38 : }
39 : }
40 :
41 14679 : bool MolDataClass::allowedResidue( const std::string& type, const std::string& residuename ) {
42 14679 : if( type=="protein" ) {
43 12454 : if(residuename=="ALA") {
44 : return true;
45 11594 : } else if(residuename=="ARG") {
46 : return true;
47 11584 : } else if(residuename=="ASN") {
48 : return true;
49 11579 : } else if(residuename=="ASP") {
50 : return true;
51 11557 : } else if(residuename=="CYS") {
52 : return true;
53 11557 : } else if(residuename=="GLN") {
54 : return true;
55 11547 : } else if(residuename=="GLU") {
56 : return true;
57 11535 : } else if(residuename=="GLY") {
58 : return true;
59 11475 : } else if(residuename=="HIS") {
60 : return true;
61 11475 : } else if(residuename=="ILE") {
62 : return true;
63 11461 : } else if(residuename=="LEU") {
64 : return true;
65 11441 : } else if(residuename=="LYS") {
66 : return true;
67 11425 : } else if(residuename=="MET") {
68 : return true;
69 11425 : } else if(residuename=="PHE") {
70 : return true;
71 11419 : } else if(residuename=="PRO") {
72 : return true;
73 11372 : } else if(residuename=="SER") {
74 : return true;
75 11325 : } else if(residuename=="THR") {
76 : return true;
77 11295 : } else if(residuename=="TRP") {
78 : return true;
79 11253 : } else if(residuename=="TYR") {
80 : return true;
81 11237 : } else if(residuename=="VAL") {
82 : return true;
83 : }
84 : // Terminal groups
85 575 : else if(residuename=="ACE") {
86 : return true;
87 569 : } else if(residuename=="NME") {
88 : return true;
89 563 : } else if(residuename=="NH2") {
90 : return true;
91 : }
92 : // Alternative residue names in common force fields
93 563 : else if(residuename=="GLH") {
94 : return true; // neutral GLU
95 563 : } else if(residuename=="ASH") {
96 : return true; // neutral ASP
97 563 : } else if(residuename=="HID") {
98 : return true; // HIS-D amber
99 563 : } else if(residuename=="HSD") {
100 : return true; // HIS-D charmm
101 563 : } else if(residuename=="HIE") {
102 : return true; // HIS-E amber
103 563 : } else if(residuename=="HSE") {
104 : return true; // HIS-E charmm
105 563 : } else if(residuename=="HIP") {
106 : return true; // HIS-P amber
107 563 : } else if(residuename=="HSP") {
108 : return true; // HIS-P charmm
109 : }
110 : // Weird amino acids
111 563 : else if(residuename=="NLE") {
112 : return true;
113 559 : } else if(residuename=="SFO") {
114 : return true;
115 : } else {
116 559 : return false;
117 : }
118 2225 : } else if( type=="dna" ) {
119 540 : if(residuename=="A") {
120 : return true;
121 475 : } else if(residuename=="A5") {
122 : return true;
123 475 : } else if(residuename=="A3") {
124 : return true;
125 475 : } else if(residuename=="AN") {
126 : return true;
127 475 : } else if(residuename=="G") {
128 : return true;
129 475 : } else if(residuename=="G5") {
130 : return true;
131 475 : } else if(residuename=="G3") {
132 : return true;
133 475 : } else if(residuename=="GN") {
134 : return true;
135 475 : } else if(residuename=="T") {
136 : return true;
137 471 : } else if(residuename=="T5") {
138 : return true;
139 471 : } else if(residuename=="T3") {
140 : return true;
141 471 : } else if(residuename=="TN") {
142 : return true;
143 471 : } else if(residuename=="C") {
144 : return true;
145 471 : } else if(residuename=="C5") {
146 : return true;
147 471 : } else if(residuename=="C3") {
148 : return true;
149 471 : } else if(residuename=="CN") {
150 : return true;
151 471 : } else if(residuename=="DA") {
152 : return true;
153 463 : } else if(residuename=="DA5") {
154 : return true;
155 463 : } else if(residuename=="DA3") {
156 : return true;
157 463 : } else if(residuename=="DAN") {
158 : return true;
159 463 : } else if(residuename=="DG") {
160 : return true;
161 448 : } else if(residuename=="DG5") {
162 : return true;
163 448 : } else if(residuename=="DG3") {
164 : return true;
165 448 : } else if(residuename=="DGN") {
166 : return true;
167 448 : } else if(residuename=="DT") {
168 : return true;
169 440 : } else if(residuename=="DT5") {
170 : return true;
171 440 : } else if(residuename=="DT3") {
172 : return true;
173 440 : } else if(residuename=="DTN") {
174 : return true;
175 440 : } else if(residuename=="DC") {
176 : return true;
177 425 : } else if(residuename=="DC5") {
178 : return true;
179 425 : } else if(residuename=="DC3") {
180 : return true;
181 425 : } else if(residuename=="DCN") {
182 : return true;
183 : } else {
184 425 : return false;
185 : }
186 1685 : } else if( type=="rna" ) {
187 905 : if(residuename=="A") {
188 : return true;
189 842 : } else if(residuename=="A5") {
190 : return true;
191 842 : } else if(residuename=="A3") {
192 : return true;
193 842 : } else if(residuename=="AN") {
194 : return true;
195 842 : } else if(residuename=="G") {
196 : return true;
197 820 : } else if(residuename=="G5") {
198 : return true;
199 816 : } else if(residuename=="G3") {
200 : return true;
201 816 : } else if(residuename=="GN") {
202 : return true;
203 816 : } else if(residuename=="U") {
204 : return true;
205 793 : } else if(residuename=="U5") {
206 : return true;
207 793 : } else if(residuename=="U3") {
208 : return true;
209 793 : } else if(residuename=="UN") {
210 : return true;
211 793 : } else if(residuename=="C") {
212 : return true;
213 778 : } else if(residuename=="C5") {
214 : return true;
215 778 : } else if(residuename=="C3") {
216 : return true;
217 774 : } else if(residuename=="CN") {
218 : return true;
219 774 : } else if(residuename=="RA") {
220 : return true;
221 666 : } else if(residuename=="RA5") {
222 : return true;
223 635 : } else if(residuename=="RA3") {
224 : return true;
225 635 : } else if(residuename=="RAN") {
226 : return true;
227 635 : } else if(residuename=="RG") {
228 : return true;
229 515 : } else if(residuename=="RG5") {
230 : return true;
231 515 : } else if(residuename=="RG3") {
232 : return true;
233 476 : } else if(residuename=="RGN") {
234 : return true;
235 476 : } else if(residuename=="RU") {
236 : return true;
237 376 : } else if(residuename=="RU5") {
238 : return true;
239 376 : } else if(residuename=="RU3") {
240 : return true;
241 345 : } else if(residuename=="RUN") {
242 : return true;
243 345 : } else if(residuename=="RC") {
244 : return true;
245 211 : } else if(residuename=="RC5") {
246 : return true;
247 178 : } else if(residuename=="RC3") {
248 : return true;
249 178 : } else if(residuename=="RCN") {
250 : return true;
251 : } else {
252 175 : return false;
253 : }
254 780 : } else if( type=="water" ) {
255 390 : if(residuename=="SOL") {
256 : return true;
257 : }
258 274 : if(residuename=="WAT") {
259 : return true;
260 : }
261 274 : return false;
262 390 : } else if( type=="ion" ) {
263 390 : if(residuename=="IB+") {
264 : return true;
265 : }
266 390 : if(residuename=="CA") {
267 : return true;
268 : }
269 390 : if(residuename=="CL") {
270 : return true;
271 : }
272 384 : if(residuename=="NA") {
273 : return true;
274 : }
275 372 : if(residuename=="MG") {
276 : return true;
277 : }
278 372 : if(residuename=="K") {
279 : return true;
280 : }
281 372 : if(residuename=="RB") {
282 : return true;
283 : }
284 372 : if(residuename=="CS") {
285 : return true;
286 : }
287 372 : if(residuename=="LI") {
288 : return true;
289 : }
290 372 : if(residuename=="ZN") {
291 : return true;
292 : }
293 372 : return false;
294 : }
295 : return false;
296 : }
297 :
298 5642 : void MolDataClass::getBackboneForResidue( const std::string& type, int residuenum, const std::string& chain, const PDB& mypdb, std::vector<AtomNumber>& atoms ) {
299 5642 : std::string residuename=mypdb.getResidueName( residuenum, chain );
300 5642 : plumed_massert( MolDataClass::allowedResidue( type, residuename ), "residue " + residuename + " unrecognized for molecule type " + type );
301 5642 : if( type=="protein" ) {
302 5642 : if( residuename=="GLY") {
303 0 : atoms.resize(5);
304 0 : atoms[0]=mypdb.getNamedAtomFromResidueAndChain("N",residuenum,chain);
305 0 : atoms[1]=mypdb.getNamedAtomFromResidueAndChain("CA",residuenum,chain);
306 0 : atoms[2]=mypdb.getNamedAtomFromResidueAndChain("HA2",residuenum,chain);
307 0 : atoms[3]=mypdb.getNamedAtomFromResidueAndChain("C",residuenum,chain);
308 0 : atoms[4]=mypdb.getNamedAtomFromResidueAndChain("O",residuenum,chain);
309 5642 : } else if( residuename=="ACE") {
310 0 : atoms.resize(1);
311 0 : atoms[0]=mypdb.getNamedAtomFromResidueAndChain("C",residuenum,chain);
312 5642 : } else if( residuename=="NME"||residuename=="NH2") {
313 0 : atoms.resize(1);
314 0 : atoms[0]=mypdb.getNamedAtomFromResidueAndChain("N",residuenum,chain);
315 : } else {
316 5642 : atoms.resize(5);
317 5642 : atoms[0]=mypdb.getNamedAtomFromResidueAndChain("N",residuenum,chain);
318 5642 : atoms[1]=mypdb.getNamedAtomFromResidueAndChain("CA",residuenum,chain);
319 5642 : atoms[2]=mypdb.getNamedAtomFromResidueAndChain("CB",residuenum,chain);
320 5642 : atoms[3]=mypdb.getNamedAtomFromResidueAndChain("C",residuenum,chain);
321 5642 : atoms[4]=mypdb.getNamedAtomFromResidueAndChain("O",residuenum,chain);
322 : }
323 0 : } else if( type=="dna" || type=="rna" ) {
324 0 : atoms.resize(6);
325 0 : atoms[0]=mypdb.getNamedAtomFromResidueAndChain("P",residuenum,chain);
326 0 : atoms[1]=mypdb.getNamedAtomFromResidueAndChain("O5\'",residuenum,chain);
327 0 : atoms[2]=mypdb.getNamedAtomFromResidueAndChain("C5\'",residuenum,chain);
328 0 : atoms[3]=mypdb.getNamedAtomFromResidueAndChain("C4\'",residuenum,chain);
329 0 : atoms[4]=mypdb.getNamedAtomFromResidueAndChain("C3\'",residuenum,chain);
330 0 : atoms[5]=mypdb.getNamedAtomFromResidueAndChain("O3\'",residuenum,chain);
331 : } else {
332 0 : plumed_merror(type + " is not a valid molecule type");
333 : }
334 5642 : }
335 :
336 7315 : bool MolDataClass::isTerminalGroup( const std::string& type, const std::string& residuename ) {
337 7315 : if( type=="protein" ) {
338 7315 : if( residuename=="ACE" ) {
339 : return true;
340 6625 : } else if( residuename=="NME" ) {
341 : return true;
342 5935 : } else if( residuename=="NH2" ) {
343 : return true;
344 : } else {
345 5935 : return false;
346 : }
347 : } else {
348 0 : plumed_merror(type + " is not a valid molecule type");
349 : }
350 : return false;
351 : }
352 :
353 777 : void MolDataClass::specialSymbol( const std::string& type, const std::string& symbol, const PDB& mypdb, std::vector<AtomNumber>& numbers ) {
354 777 : if(type=="protein" || type=="rna" || type=="dna") {
355 : // symbol should be something like
356 : // phi-123 i.e. phi torsion of residue 123 of first chain
357 : // psi-A321 i.e. psi torsion of residue 321 of chain A
358 : // psi-4_321 is psi torsion of residue 321 of chain 4
359 : // psi-A_321 is equivalent to psi-A321
360 777 : numbers.resize(0);
361 :
362 : // special cases:
363 777 : if(symbol=="protein") {
364 1 : const auto & all = mypdb.getAtomNumbers();
365 133 : for(auto a : all) {
366 132 : auto resname=mypdb.getResidueName(a);
367 132 : Tools::stripLeadingAndTrailingBlanks(resname);
368 264 : if(allowedResidue("protein",resname)) {
369 132 : numbers.push_back(a);
370 : }
371 : }
372 7 : return;
373 : }
374 :
375 776 : if(symbol=="nucleic") {
376 2 : const auto & all = mypdb.getAtomNumbers();
377 457 : for(auto a : all) {
378 455 : auto resname=mypdb.getResidueName(a);
379 455 : Tools::stripLeadingAndTrailingBlanks(resname);
380 845 : if(allowedResidue("dna",resname) || allowedResidue("rna",resname)) {
381 321 : numbers.push_back(a);
382 : }
383 : }
384 2 : return;
385 : }
386 :
387 774 : if(symbol=="ions") {
388 1 : const auto & all = mypdb.getAtomNumbers();
389 391 : for(auto a : all) {
390 390 : auto resname=mypdb.getResidueName(a);
391 390 : Tools::stripLeadingAndTrailingBlanks(resname);
392 780 : if(allowedResidue("ion",resname)) {
393 18 : numbers.push_back(a);
394 : }
395 : }
396 1 : return;
397 : }
398 :
399 773 : if(symbol=="water") {
400 1 : const auto & all = mypdb.getAtomNumbers();
401 391 : for(auto a : all) {
402 390 : auto resname=mypdb.getResidueName(a);
403 390 : Tools::stripLeadingAndTrailingBlanks(resname);
404 780 : if(allowedResidue("water",resname)) {
405 116 : numbers.push_back(a);
406 : }
407 : }
408 1 : return;
409 : }
410 :
411 772 : if(symbol=="hydrogens") {
412 1 : const auto & all = mypdb.getAtomNumbers();
413 391 : for(auto a : all) {
414 390 : auto atomname=mypdb.getAtomName(a);
415 390 : Tools::stripLeadingAndTrailingBlanks(atomname);
416 390 : auto notnumber=atomname.find_first_not_of("0123456789");
417 390 : if(notnumber!=std::string::npos && atomname[notnumber]=='H') {
418 149 : numbers.push_back(a);
419 : }
420 : }
421 1 : return;
422 : }
423 :
424 771 : if(symbol=="nonhydrogens") {
425 1 : const auto & all = mypdb.getAtomNumbers();
426 391 : for(auto a : all) {
427 390 : auto atomname=mypdb.getAtomName(a);
428 390 : Tools::stripLeadingAndTrailingBlanks(atomname);
429 390 : auto notnumber=atomname.find_first_not_of("0123456789");
430 390 : if(!(notnumber!=std::string::npos && atomname[notnumber]=='H')) {
431 241 : numbers.push_back(a);
432 : }
433 : }
434 1 : return;
435 : }
436 :
437 :
438 : std::size_t dash=symbol.find_first_of('-');
439 770 : if(dash==std::string::npos) {
440 0 : plumed_error() << "Unrecognized symbol "<<symbol;
441 : }
442 :
443 770 : std::size_t firstunderscore=symbol.find_first_of('_',dash+1);
444 770 : std::size_t firstnum=symbol.find_first_of("0123456789",dash+1);
445 : // a minus sign right before the residue number is part of the number (negative residue number),
446 : // e.g. @phi--3, @phi-A-3 (chain A, residue -3) and @phi-A_-3
447 770 : if(firstnum!=std::string::npos && firstnum>dash+1 && symbol[firstnum-1]=='-') {
448 : firstnum--;
449 : }
450 770 : std::string name=symbol.substr(0,dash);
451 : int resnum;
452 : std::string resname;
453 : std::string chainid;
454 770 : if(firstunderscore != std::string::npos) {
455 7 : Tools::convert( symbol.substr(firstunderscore+1), resnum );
456 14 : chainid=symbol.substr(dash+1,firstunderscore-(dash+1));
457 763 : } else if(firstnum==dash+1) {
458 1504 : Tools::convert( symbol.substr(dash+1), resnum );
459 : chainid="*"; // this is going to match the first chain
460 : } else {
461 : // if chain id is provided:
462 11 : Tools::convert( symbol.substr(firstnum), resnum );
463 22 : chainid=symbol.substr(dash+1,firstnum-(dash+1)); // this is going to match the proper chain
464 : }
465 770 : resname= mypdb.getResidueName(resnum,chainid);
466 770 : Tools::stripLeadingAndTrailingBlanks(resname);
467 1540 : if(allowedResidue("protein",resname)) {
468 352 : if( name=="phi" && !isTerminalGroup("protein",resname) ) {
469 63 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C",resnum-1,chainid));
470 63 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("N",resnum,chainid));
471 63 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("CA",resnum,chainid));
472 126 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C",resnum,chainid));
473 446 : } else if( name=="psi" && !isTerminalGroup("protein",resname) ) {
474 220 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("N",resnum,chainid));
475 220 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("CA",resnum,chainid));
476 220 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C",resnum,chainid));
477 440 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("N",resnum+1,chainid));
478 8 : } else if( name=="omega" && !isTerminalGroup("protein",resname) ) {
479 2 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("CA",resnum-1,chainid));
480 2 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C",resnum-1,chainid));
481 2 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("N",resnum,chainid));
482 4 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("CA",resnum,chainid));
483 5 : } else if( name=="chi1" && !isTerminalGroup("protein",resname) ) {
484 1 : if ( resname=="GLY" || resname=="ALA" || resname=="SFO" ) {
485 0 : plumed_merror("chi-1 is not defined for ALA, GLY and SFO");
486 : }
487 1 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("N",resnum,chainid));
488 1 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("CA",resnum,chainid));
489 1 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("CB",resnum,chainid));
490 1 : if(resname=="ILE"||resname=="VAL") {
491 0 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("CG1",resnum,chainid));
492 1 : } else if(resname=="CYS") {
493 0 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("SG",resnum,chainid));
494 1 : } else if(resname=="THR") {
495 0 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("OG1",resnum,chainid));
496 1 : } else if(resname=="SER") {
497 0 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("OG",resnum,chainid));
498 : } else {
499 2 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("CG",resnum,chainid));
500 : }
501 4 : } else if( name=="chi2" && !isTerminalGroup("protein",resname) ) {
502 1 : if ( resname=="GLY" || resname=="ALA" || resname=="SFO" || resname=="CYS" || resname=="SER" ||
503 2 : resname=="THR" || resname=="VAL" ) {
504 0 : plumed_merror("chi-2 is not defined for ALA, GLY, CYS, SER, THR, VAL and SFO");
505 : }
506 1 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("CA",resnum,chainid));
507 1 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("CB",resnum,chainid));
508 :
509 1 : if(resname=="ILE") {
510 0 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("CG1",resnum,chainid));
511 : } else {
512 2 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("CG",resnum,chainid));
513 : }
514 1 : if(resname=="ASN" || resname=="ASP") {
515 0 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("OD1",resnum,chainid));
516 1 : } else if(resname=="HIS") {
517 0 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("ND1",resnum,chainid));
518 1 : } else if(resname=="LEU" || resname=="PHE" || resname=="TRP" || resname=="TYR") {
519 0 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("CD1",resnum,chainid));
520 1 : } else if(resname=="MET") {
521 0 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("SD",resnum,chainid));
522 : } else {
523 2 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("CD",resnum,chainid));
524 : }
525 2 : } else if( name=="chi3" && !isTerminalGroup("protein",resname) ) {
526 0 : if (!( resname=="ARG" || resname=="GLN" || resname=="GLU" || resname=="LYS" ||
527 0 : resname=="MET" )) {
528 0 : plumed_merror("chi-3 is defined only for ARG, GLN, GLU, LYS and MET");
529 : }
530 0 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("CB",resnum,chainid));
531 0 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("CG",resnum,chainid));
532 :
533 0 : if(resname=="MET") {
534 0 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("SD",resnum,chainid));
535 : } else {
536 0 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("CD",resnum,chainid));
537 : }
538 0 : if(resname=="GLN" || resname=="GLU") {
539 0 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("OE1",resnum,chainid));
540 0 : } else if(resname=="LYS" || resname=="MET") {
541 0 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("CE",resnum,chainid));
542 : } else {
543 0 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("NE",resnum,chainid));
544 : }
545 2 : } else if( name=="chi4" && !isTerminalGroup("protein",resname) ) {
546 0 : if (!( resname=="ARG" || resname=="LYS" )) {
547 0 : plumed_merror("chi-4 is defined only for ARG and LYS");
548 : }
549 0 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("CG",resnum,chainid));
550 0 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("CD",resnum,chainid));
551 :
552 0 : if(resname=="ARG") {
553 0 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("NE",resnum,chainid));
554 0 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("CZ",resnum,chainid));
555 : } else {
556 0 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("CE",resnum,chainid));
557 0 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("NZ",resnum,chainid));
558 : }
559 2 : } else if( name=="chi5" && !isTerminalGroup("protein",resname) ) {
560 0 : if (!( resname=="ARG" )) {
561 0 : plumed_merror("chi-5 is defined only for ARG");
562 : }
563 0 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("CD",resnum,chainid));
564 0 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("NE",resnum,chainid));
565 0 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("CZ",resnum,chainid));
566 0 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("NH1",resnum,chainid));
567 3 : } else if( name=="sidechain" && !isTerminalGroup("protein",resname) ) {
568 1 : if(resname=="GLY") {
569 0 : plumed_merror("sidechain is not defined for GLY");
570 : }
571 1 : std::vector<AtomNumber> tmpnum1(mypdb.getAtomsInResidue(resnum,chainid));
572 15 : for(unsigned i=0; i<tmpnum1.size(); i++) {
573 41 : if(mypdb.getAtomName(tmpnum1[i])!="N"&&mypdb.getAtomName(tmpnum1[i])!="CA"&&
574 38 : mypdb.getAtomName(tmpnum1[i])!="C"&&mypdb.getAtomName(tmpnum1[i])!="O"&&
575 34 : mypdb.getAtomName(tmpnum1[i])!="HA"&&mypdb.getAtomName(tmpnum1[i])!="H"&&
576 32 : mypdb.getAtomName(tmpnum1[i])!="HN"&&mypdb.getAtomName(tmpnum1[i])!="H1"&&
577 32 : mypdb.getAtomName(tmpnum1[i])!="H2"&&mypdb.getAtomName(tmpnum1[i])!="H3"&&
578 32 : mypdb.getAtomName(tmpnum1[i])!="OXT"&&mypdb.getAtomName(tmpnum1[i])!="OT1"&&
579 46 : mypdb.getAtomName(tmpnum1[i])!="OT2"&&mypdb.getAtomName(tmpnum1[i])!="OC1"&&
580 23 : mypdb.getAtomName(tmpnum1[i])!="OC2") {
581 9 : numbers.push_back( tmpnum1[i] );
582 : }
583 : }
584 2 : } else if( name=="back" && !isTerminalGroup("protein",resname) ) {
585 1 : std::vector<AtomNumber> tmpnum1(mypdb.getAtomsInResidue(resnum,chainid));
586 15 : for(unsigned i=0; i<tmpnum1.size(); i++) {
587 41 : if(mypdb.getAtomName(tmpnum1[i])=="N"||mypdb.getAtomName(tmpnum1[i])=="CA"||
588 38 : mypdb.getAtomName(tmpnum1[i])=="C"||mypdb.getAtomName(tmpnum1[i])=="O"||
589 34 : mypdb.getAtomName(tmpnum1[i])=="HA"||mypdb.getAtomName(tmpnum1[i])=="H"||
590 32 : mypdb.getAtomName(tmpnum1[i])=="HN"||mypdb.getAtomName(tmpnum1[i])=="H1"||
591 32 : mypdb.getAtomName(tmpnum1[i])=="H2"||mypdb.getAtomName(tmpnum1[i])=="H3"||
592 32 : mypdb.getAtomName(tmpnum1[i])=="OXT"||mypdb.getAtomName(tmpnum1[i])=="OT1"||
593 32 : mypdb.getAtomName(tmpnum1[i])=="OT2"||mypdb.getAtomName(tmpnum1[i])=="OC1"||
594 32 : mypdb.getAtomName(tmpnum1[i])=="OC2"||mypdb.getAtomName(tmpnum1[i])=="HA1"||
595 45 : mypdb.getAtomName(tmpnum1[i])=="HA2"||mypdb.getAtomName(tmpnum1[i])=="HA3") {
596 5 : numbers.push_back( tmpnum1[i] );
597 : }
598 : }
599 : } else {
600 0 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain(name,resnum,chainid));
601 : }
602 :
603 488 : } else if( allowedResidue("rna",resname) || allowedResidue("dna",resname)) {
604 : std::string basetype;
605 480 : if(resname.find_first_of("A")!=std::string::npos) {
606 : basetype+="A";
607 : }
608 480 : if(resname.find_first_of("U")!=std::string::npos) {
609 : basetype+="U";
610 : }
611 480 : if(resname.find_first_of("T")!=std::string::npos) {
612 : basetype+="T";
613 : }
614 480 : if(resname.find_first_of("C")!=std::string::npos) {
615 : basetype+="C";
616 : }
617 480 : if(resname.find_first_of("G")!=std::string::npos) {
618 : basetype+="G";
619 : }
620 480 : plumed_massert(basetype.length()==1,"cannot find type of rna/dna residue "+resname+" "+basetype);
621 480 : if( name=="chi" ) {
622 36 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("O4\'",resnum,chainid));
623 36 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C1\'",resnum,chainid));
624 36 : if(basetype=="T" || basetype=="U" || basetype=="C") {
625 9 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("N1",resnum,chainid));
626 18 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C2",resnum,chainid));
627 27 : } else if(basetype=="G" || basetype=="A") {
628 27 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("N9",resnum,chainid));
629 54 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C4",resnum,chainid));
630 : } else {
631 0 : plumed_error();
632 : }
633 444 : } else if( name=="alpha" ) {
634 7 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("O3\'",resnum-1,chainid));
635 7 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("P",resnum,chainid));
636 7 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("O5\'",resnum,chainid));
637 14 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C5\'",resnum,chainid));
638 437 : } else if( name=="beta" ) {
639 7 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("P",resnum,chainid));
640 7 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("O5\'",resnum,chainid));
641 7 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C5\'",resnum,chainid));
642 14 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C4\'",resnum,chainid));
643 430 : } else if( name=="gamma" ) {
644 29 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("O5\'",resnum,chainid));
645 29 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C5\'",resnum,chainid));
646 29 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C4\'",resnum,chainid));
647 58 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C3\'",resnum,chainid));
648 401 : } else if( name=="delta" ) {
649 1 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C5\'",resnum,chainid));
650 1 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C4\'",resnum,chainid));
651 1 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C3\'",resnum,chainid));
652 2 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("O3\'",resnum,chainid));
653 400 : } else if( name=="epsilon" ) {
654 8 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C4\'",resnum,chainid));
655 8 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C3\'",resnum,chainid));
656 8 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("O3\'",resnum,chainid));
657 16 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("P",resnum+1,chainid));
658 392 : } else if( name=="zeta" ) {
659 7 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C3\'",resnum,chainid));
660 7 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("O3\'",resnum,chainid));
661 7 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("P",resnum+1,chainid));
662 14 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("O5\'",resnum+1,chainid));
663 385 : } else if( name=="v0" ) {
664 3 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C4\'",resnum,chainid));
665 3 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("O4\'",resnum,chainid));
666 3 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C1\'",resnum,chainid));
667 6 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C2\'",resnum,chainid));
668 382 : } else if( name=="v1" ) {
669 3 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("O4\'",resnum,chainid));
670 3 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C1\'",resnum,chainid));
671 3 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C2\'",resnum,chainid));
672 6 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C3\'",resnum,chainid));
673 379 : } else if( name=="v2" ) {
674 3 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C1\'",resnum,chainid));
675 3 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C2\'",resnum,chainid));
676 3 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C3\'",resnum,chainid));
677 6 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C4\'",resnum,chainid));
678 376 : } else if( name=="v3" ) {
679 3 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C2\'",resnum,chainid));
680 3 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C3\'",resnum,chainid));
681 3 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C4\'",resnum,chainid));
682 6 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("O4\'",resnum,chainid));
683 373 : } else if( name=="v4" ) {
684 3 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C3\'",resnum,chainid));
685 3 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C4\'",resnum,chainid));
686 3 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("O4\'",resnum,chainid));
687 6 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C1\'",resnum,chainid));
688 370 : } else if( name=="back" ) {
689 1 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("P",resnum,chainid));
690 1 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("O5\'",resnum,chainid));
691 1 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C5\'",resnum,chainid));
692 1 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C4\'",resnum,chainid));
693 1 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C3\'",resnum,chainid));
694 2 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("O3\'",resnum,chainid));
695 369 : } else if( name=="sugar" ) {
696 54 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C4\'",resnum,chainid));
697 54 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("O4\'",resnum,chainid));
698 54 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C1\'",resnum,chainid));
699 54 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C2\'",resnum,chainid));
700 108 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C3\'",resnum,chainid));
701 315 : } else if( name=="base" ) {
702 25 : if(basetype=="C") {
703 9 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("N1",resnum,chainid));
704 9 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C2",resnum,chainid));
705 9 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("O2",resnum,chainid));
706 9 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("N3",resnum,chainid));
707 9 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C4",resnum,chainid));
708 9 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("N4",resnum,chainid));
709 9 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C5",resnum,chainid));
710 18 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C6",resnum,chainid));
711 16 : } else if(basetype=="U") {
712 2 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("N1",resnum,chainid));
713 2 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C2",resnum,chainid));
714 2 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("O2",resnum,chainid));
715 2 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("N3",resnum,chainid));
716 2 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C4",resnum,chainid));
717 2 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("O4",resnum,chainid));
718 2 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C5",resnum,chainid));
719 4 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C6",resnum,chainid));
720 14 : } else if(basetype=="T") {
721 0 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("N1",resnum,chainid));
722 0 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C2",resnum,chainid));
723 0 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("O2",resnum,chainid));
724 0 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("N3",resnum,chainid));
725 0 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C4",resnum,chainid));
726 0 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("O4",resnum,chainid));
727 0 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C5",resnum,chainid));
728 0 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C7",resnum,chainid));
729 0 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C6",resnum,chainid));
730 14 : } else if(basetype=="G") {
731 3 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("N9",resnum,chainid));
732 3 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C4",resnum,chainid));
733 3 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("N3",resnum,chainid));
734 3 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C2",resnum,chainid));
735 3 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("N2",resnum,chainid));
736 3 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("N1",resnum,chainid));
737 3 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C6",resnum,chainid));
738 3 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("O6",resnum,chainid));
739 3 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C5",resnum,chainid));
740 3 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("N7",resnum,chainid));
741 6 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C8",resnum,chainid));
742 11 : } else if(basetype=="A") {
743 11 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("N9",resnum,chainid));
744 11 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C4",resnum,chainid));
745 11 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("N1",resnum,chainid));
746 11 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C2",resnum,chainid));
747 11 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("N3",resnum,chainid));
748 11 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C6",resnum,chainid));
749 11 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("N6",resnum,chainid));
750 11 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C5",resnum,chainid));
751 11 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("N7",resnum,chainid));
752 22 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C8",resnum,chainid));
753 : } else {
754 0 : plumed_error();
755 : }
756 290 : } else if( name=="lcs" ) {
757 284 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C2",resnum,chainid));
758 284 : if(basetype=="T" || basetype=="U" || basetype=="C") {
759 148 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C4",resnum,chainid));
760 296 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C6",resnum,chainid));
761 136 : } else if(basetype=="G" || basetype=="A") {
762 136 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C6",resnum,chainid));
763 272 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C4",resnum,chainid));
764 : } else {
765 0 : plumed_error();
766 : }
767 : } else {
768 12 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain(name,resnum,chainid));
769 : }
770 : } else {
771 2 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain(name,resnum,chainid));
772 : }
773 : } else {
774 0 : plumed_merror(type + " is not a valid molecule type");
775 : }
776 : }
777 :
778 : }
|
