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 32 : unsigned MolDataClass::numberOfAtomsPerResidueInBackbone( const std::string& type ) {
30 32 : if( type=="protein" ) return 5;
31 0 : else if( type=="dna" ) return 6;
32 0 : else if( type=="rna" ) return 6;
33 0 : else return 0;
34 : }
35 :
36 9955 : bool MolDataClass::allowedResidue( const std::string& type, const std::string& residuename ) {
37 9955 : if( type=="protein" ) {
38 7842 : if(residuename=="ALA") return true;
39 7526 : else if(residuename=="ARG") return true;
40 7526 : else if(residuename=="ASN") return true;
41 7526 : else if(residuename=="ASP") return true;
42 7514 : else if(residuename=="CYS") return true;
43 7514 : else if(residuename=="GLN") return true;
44 7514 : else if(residuename=="GLU") return true;
45 7502 : else if(residuename=="GLY") return true;
46 7492 : else if(residuename=="HIS") return true;
47 7492 : else if(residuename=="ILE") return true;
48 7488 : else if(residuename=="LEU") return true;
49 7488 : else if(residuename=="LYS") return true;
50 7482 : else if(residuename=="MET") return true;
51 7482 : else if(residuename=="PHE") return true;
52 7476 : else if(residuename=="PRO") return true;
53 7469 : else if(residuename=="SER") return true;
54 7465 : else if(residuename=="THR") return true;
55 7435 : else if(residuename=="TRP") return true;
56 7421 : else if(residuename=="TYR") return true;
57 7415 : else if(residuename=="VAL") return true;
58 : // Terminal groups
59 497 : else if(residuename=="ACE") return true;
60 491 : else if(residuename=="NME") return true;
61 485 : else if(residuename=="NH2") return true;
62 : // Alternative residue names in common force fields
63 485 : else if(residuename=="GLH") return true; // neutral GLU
64 485 : else if(residuename=="ASH") return true; // neutral ASP
65 485 : else if(residuename=="HID") return true; // HIS-D amber
66 485 : else if(residuename=="HSD") return true; // HIS-D charmm
67 485 : else if(residuename=="HIE") return true; // HIS-E amber
68 485 : else if(residuename=="HSE") return true; // HIS-E charmm
69 485 : else if(residuename=="HIP") return true; // HIS-P amber
70 485 : else if(residuename=="HSP") return true; // HIS-P charmm
71 : // Weird amino acids
72 485 : else if(residuename=="NLE") return true;
73 481 : else if(residuename=="SFO") return true;
74 481 : else return false;
75 2113 : } else if( type=="dna" ) {
76 462 : if(residuename=="A") return true;
77 397 : else if(residuename=="A5") return true;
78 397 : else if(residuename=="A3") return true;
79 397 : else if(residuename=="AN") return true;
80 397 : else if(residuename=="G") return true;
81 397 : else if(residuename=="G5") return true;
82 397 : else if(residuename=="G3") return true;
83 397 : else if(residuename=="GN") return true;
84 397 : else if(residuename=="T") return true;
85 393 : else if(residuename=="T5") return true;
86 393 : else if(residuename=="T3") return true;
87 393 : else if(residuename=="TN") return true;
88 393 : else if(residuename=="C") return true;
89 393 : else if(residuename=="C5") return true;
90 393 : else if(residuename=="C3") return true;
91 393 : else if(residuename=="CN") return true;
92 393 : else if(residuename=="DA") return true;
93 393 : else if(residuename=="DA5") return true;
94 393 : else if(residuename=="DA3") return true;
95 393 : else if(residuename=="DAN") return true;
96 393 : else if(residuename=="DG") return true;
97 392 : else if(residuename=="DG5") return true;
98 392 : else if(residuename=="DG3") return true;
99 392 : else if(residuename=="DGN") return true;
100 392 : else if(residuename=="DT") return true;
101 392 : else if(residuename=="DT5") return true;
102 392 : else if(residuename=="DT3") return true;
103 392 : else if(residuename=="DTN") return true;
104 392 : else if(residuename=="DC") return true;
105 391 : else if(residuename=="DC5") return true;
106 391 : else if(residuename=="DC3") return true;
107 391 : else if(residuename=="DCN") return true;
108 391 : else return false;
109 1651 : } else if( type=="rna" ) {
110 871 : if(residuename=="A") return true;
111 808 : else if(residuename=="A5") return true;
112 808 : else if(residuename=="A3") return true;
113 808 : else if(residuename=="AN") return true;
114 808 : else if(residuename=="G") return true;
115 786 : else if(residuename=="G5") return true;
116 782 : else if(residuename=="G3") return true;
117 782 : else if(residuename=="GN") return true;
118 782 : else if(residuename=="U") return true;
119 759 : else if(residuename=="U5") return true;
120 759 : else if(residuename=="U3") return true;
121 759 : else if(residuename=="UN") return true;
122 759 : else if(residuename=="C") return true;
123 744 : else if(residuename=="C5") return true;
124 744 : else if(residuename=="C3") return true;
125 740 : else if(residuename=="CN") return true;
126 740 : else if(residuename=="RA") return true;
127 632 : else if(residuename=="RA5") return true;
128 601 : else if(residuename=="RA3") return true;
129 601 : else if(residuename=="RAN") return true;
130 601 : else if(residuename=="RG") return true;
131 481 : else if(residuename=="RG5") return true;
132 481 : else if(residuename=="RG3") return true;
133 442 : else if(residuename=="RGN") return true;
134 442 : else if(residuename=="RU") return true;
135 342 : else if(residuename=="RU5") return true;
136 342 : else if(residuename=="RU3") return true;
137 311 : else if(residuename=="RUN") return true;
138 311 : else if(residuename=="RC") return true;
139 177 : else if(residuename=="RC5") return true;
140 144 : else if(residuename=="RC3") return true;
141 144 : else if(residuename=="RCN") return true;
142 141 : else return false;
143 780 : } else if( type=="water" ) {
144 390 : if(residuename=="SOL") return true;
145 274 : if(residuename=="WAT") return true;
146 274 : return false;
147 390 : } else if( type=="ion" ) {
148 390 : if(residuename=="IB+") return true;
149 390 : if(residuename=="CA") return true;
150 390 : if(residuename=="CL") return true;
151 384 : if(residuename=="NA") return true;
152 372 : if(residuename=="MG") return true;
153 372 : if(residuename=="K") return true;
154 372 : if(residuename=="RB") return true;
155 372 : if(residuename=="CS") return true;
156 372 : if(residuename=="LI") return true;
157 372 : if(residuename=="ZN") return true;
158 372 : return false;
159 : }
160 : return false;
161 : }
162 :
163 3540 : void MolDataClass::getBackboneForResidue( const std::string& type, const unsigned& residuenum, const PDB& mypdb, std::vector<AtomNumber>& atoms ) {
164 3540 : std::string residuename=mypdb.getResidueName( residuenum );
165 3540 : plumed_massert( MolDataClass::allowedResidue( type, residuename ), "residue " + residuename + " unrecognized for molecule type " + type );
166 3540 : if( type=="protein" ) {
167 3540 : if( residuename=="GLY") {
168 0 : atoms.resize(5);
169 0 : atoms[0]=mypdb.getNamedAtomFromResidue("N",residuenum);
170 0 : atoms[1]=mypdb.getNamedAtomFromResidue("CA",residuenum);
171 0 : atoms[2]=mypdb.getNamedAtomFromResidue("HA2",residuenum);
172 0 : atoms[3]=mypdb.getNamedAtomFromResidue("C",residuenum);
173 0 : atoms[4]=mypdb.getNamedAtomFromResidue("O",residuenum);
174 3540 : } else if( residuename=="ACE") {
175 0 : atoms.resize(1);
176 0 : atoms[0]=mypdb.getNamedAtomFromResidue("C",residuenum);
177 7080 : } else if( residuename=="NME"||residuename=="NH2") {
178 0 : atoms.resize(1);
179 0 : atoms[0]=mypdb.getNamedAtomFromResidue("N",residuenum);
180 : } else {
181 3540 : atoms.resize(5);
182 3540 : atoms[0]=mypdb.getNamedAtomFromResidue("N",residuenum);
183 3540 : atoms[1]=mypdb.getNamedAtomFromResidue("CA",residuenum);
184 3540 : atoms[2]=mypdb.getNamedAtomFromResidue("CB",residuenum);
185 3540 : atoms[3]=mypdb.getNamedAtomFromResidue("C",residuenum);
186 3540 : atoms[4]=mypdb.getNamedAtomFromResidue("O",residuenum);
187 : }
188 0 : } else if( type=="dna" || type=="rna" ) {
189 0 : atoms.resize(6);
190 0 : atoms[0]=mypdb.getNamedAtomFromResidue("P",residuenum);
191 0 : atoms[1]=mypdb.getNamedAtomFromResidue("O5\'",residuenum);
192 0 : atoms[2]=mypdb.getNamedAtomFromResidue("C5\'",residuenum);
193 0 : atoms[3]=mypdb.getNamedAtomFromResidue("C4\'",residuenum);
194 0 : atoms[4]=mypdb.getNamedAtomFromResidue("C3\'",residuenum);
195 0 : atoms[5]=mypdb.getNamedAtomFromResidue("O3\'",residuenum);
196 : }
197 : else {
198 0 : plumed_merror(type + " is not a valid molecule type");
199 : }
200 3540 : }
201 :
202 4565 : bool MolDataClass::isTerminalGroup( const std::string& type, const std::string& residuename ) {
203 4565 : if( type=="protein" ) {
204 4565 : if( residuename=="ACE" ) return true;
205 4127 : else if( residuename=="NME" ) return true;
206 3689 : else if( residuename=="NH2" ) return true;
207 3689 : else return false;
208 : } else {
209 0 : plumed_merror(type + " is not a valid molecule type");
210 : }
211 : return false;
212 : }
213 :
214 637 : void MolDataClass::specialSymbol( const std::string& type, const std::string& symbol, const PDB& mypdb, std::vector<AtomNumber>& numbers ) {
215 873 : if(type=="protein" || type=="rna" || type=="dna") {
216 : // symbol should be something like
217 : // phi-123 i.e. phi torsion of residue 123 of first chain
218 : // psi-A321 i.e. psi torsion of residue 321 of chain A
219 : // psi-4_321 is psi torsion of residue 321 of chain 4
220 : // psi-A_321 is equivalent to psi-A321
221 637 : numbers.resize(0);
222 :
223 : // special cases:
224 637 : if(symbol=="protein") {
225 1 : const auto & all = mypdb.getAtomNumbers();
226 133 : for(auto a : all) {
227 132 : auto resname=mypdb.getResidueName(a);
228 132 : Tools::stripLeadingAndTrailingBlanks(resname);
229 264 : if(allowedResidue("protein",resname)) numbers.push_back(a);
230 : }
231 7 : return;
232 : }
233 :
234 636 : if(symbol=="nucleic") {
235 2 : const auto & all = mypdb.getAtomNumbers();
236 457 : for(auto a : all) {
237 455 : auto resname=mypdb.getResidueName(a);
238 455 : Tools::stripLeadingAndTrailingBlanks(resname);
239 1101 : if(allowedResidue("dna",resname) || allowedResidue("rna",resname)) numbers.push_back(a);
240 : }
241 2 : return;
242 : }
243 :
244 634 : if(symbol=="ions") {
245 1 : const auto & all = mypdb.getAtomNumbers();
246 391 : for(auto a : all) {
247 390 : auto resname=mypdb.getResidueName(a);
248 390 : Tools::stripLeadingAndTrailingBlanks(resname);
249 780 : if(allowedResidue("ion",resname)) numbers.push_back(a);
250 : }
251 1 : return;
252 : }
253 :
254 633 : if(symbol=="water") {
255 1 : const auto & all = mypdb.getAtomNumbers();
256 391 : for(auto a : all) {
257 390 : auto resname=mypdb.getResidueName(a);
258 390 : Tools::stripLeadingAndTrailingBlanks(resname);
259 780 : if(allowedResidue("water",resname)) numbers.push_back(a);
260 : }
261 1 : return;
262 : }
263 :
264 632 : if(symbol=="hydrogens") {
265 1 : const auto & all = mypdb.getAtomNumbers();
266 391 : for(auto a : all) {
267 390 : auto atomname=mypdb.getAtomName(a);
268 390 : Tools::stripLeadingAndTrailingBlanks(atomname);
269 390 : auto notnumber=atomname.find_first_not_of("0123456789");
270 390 : if(notnumber!=std::string::npos && atomname[notnumber]=='H') numbers.push_back(a);
271 : }
272 1 : return;
273 : }
274 :
275 631 : if(symbol=="nonhydrogens") {
276 1 : const auto & all = mypdb.getAtomNumbers();
277 391 : for(auto a : all) {
278 390 : auto atomname=mypdb.getAtomName(a);
279 390 : Tools::stripLeadingAndTrailingBlanks(atomname);
280 390 : auto notnumber=atomname.find_first_not_of("0123456789");
281 390 : if(!(notnumber!=std::string::npos && atomname[notnumber]=='H')) numbers.push_back(a);
282 : }
283 1 : return;
284 : }
285 :
286 :
287 : std::size_t dash=symbol.find_first_of('-');
288 630 : if(dash==std::string::npos) plumed_error() << "Unrecognized symbol "<<symbol;
289 :
290 630 : std::size_t firstunderscore=symbol.find_first_of('_',dash+1);
291 630 : std::size_t firstnum=symbol.find_first_of("0123456789",dash+1);
292 630 : std::string name=symbol.substr(0,dash);
293 : unsigned resnum;
294 : std::string resname;
295 : std::string chainid;
296 630 : if(firstunderscore != std::string::npos) {
297 6 : Tools::convert( symbol.substr(firstunderscore+1), resnum );
298 12 : chainid=symbol.substr(dash+1,firstunderscore-(dash+1));
299 624 : } else if(firstnum==dash+1) {
300 1228 : Tools::convert( symbol.substr(dash+1), resnum );
301 : chainid="*"; // this is going to match the first chain
302 : } else {
303 : // if chain id is provided:
304 10 : Tools::convert( symbol.substr(firstnum), resnum );
305 20 : chainid=symbol.substr(dash+1,firstnum-(dash+1)); // this is going to match the proper chain
306 : }
307 630 : resname= mypdb.getResidueName(resnum,chainid);
308 630 : Tools::stripLeadingAndTrailingBlanks(resname);
309 1260 : if(allowedResidue("protein",resname)) {
310 208 : if( name=="phi" && !isTerminalGroup("protein",resname) ) {
311 118 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C",resnum-1,chainid));
312 118 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("N",resnum,chainid));
313 118 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("CA",resnum,chainid));
314 118 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C",resnum,chainid));
315 175 : } else if( name=="psi" && !isTerminalGroup("protein",resname) ) {
316 170 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("N",resnum,chainid));
317 170 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("CA",resnum,chainid));
318 170 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C",resnum,chainid));
319 170 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("N",resnum+1,chainid));
320 6 : } else if( name=="omega" && !isTerminalGroup("protein",resname) ) {
321 2 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("CA",resnum-1,chainid));
322 2 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C",resnum-1,chainid));
323 2 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("N",resnum,chainid));
324 2 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("CA",resnum,chainid));
325 5 : } else if( name=="chi1" && !isTerminalGroup("protein",resname) ) {
326 3 : if ( resname=="GLY" || resname=="ALA" || resname=="SFO" ) plumed_merror("chi-1 is not defined for ALA, GLY and SFO");
327 2 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("N",resnum,chainid));
328 2 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("CA",resnum,chainid));
329 2 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("CB",resnum,chainid));
330 2 : if(resname=="ILE"||resname=="VAL") {
331 0 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("CG1",resnum,chainid));
332 1 : } else if(resname=="CYS") {
333 0 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("SG",resnum,chainid));
334 1 : } else if(resname=="THR") {
335 0 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("OG1",resnum,chainid));
336 1 : } else if(resname=="SER") {
337 0 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("OG",resnum,chainid));
338 : } else {
339 2 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("CG",resnum,chainid));
340 : }
341 4 : } else if( name=="chi2" && !isTerminalGroup("protein",resname) ) {
342 5 : if ( resname=="GLY" || resname=="ALA" || resname=="SFO" || resname=="CYS" || resname=="SER" ||
343 2 : resname=="THR" || resname=="VAL" ) plumed_merror("chi-2 is not defined for ALA, GLY, CYS, SER, THR, VAL and SFO");
344 2 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("CA",resnum,chainid));
345 2 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("CB",resnum,chainid));
346 :
347 1 : if(resname=="ILE") {
348 0 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("CG1",resnum,chainid));
349 : } else {
350 2 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("CG",resnum,chainid));
351 : }
352 2 : if(resname=="ASN" || resname=="ASP") {
353 0 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("OD1",resnum,chainid));
354 1 : } else if(resname=="HIS") {
355 0 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("ND1",resnum,chainid));
356 4 : } else if(resname=="LEU" || resname=="PHE" || resname=="TRP" || resname=="TYR") {
357 0 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("CD1",resnum,chainid));
358 1 : } else if(resname=="MET") {
359 0 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("SD",resnum,chainid));
360 : } else {
361 2 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("CD",resnum,chainid));
362 : }
363 2 : } else if( name=="chi3" && !isTerminalGroup("protein",resname) ) {
364 0 : if (!( resname=="ARG" || resname=="GLN" || resname=="GLU" || resname=="LYS" ||
365 0 : resname=="MET" )) plumed_merror("chi-3 is defined only for ARG, GLN, GLU, LYS and MET");
366 0 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("CB",resnum,chainid));
367 0 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("CG",resnum,chainid));
368 :
369 0 : if(resname=="MET") {
370 0 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("SD",resnum,chainid));
371 : } else {
372 0 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("CD",resnum,chainid));
373 : }
374 0 : if(resname=="GLN" || resname=="GLU") {
375 0 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("OE1",resnum,chainid));
376 0 : } else if(resname=="LYS" || resname=="MET") {
377 0 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("CE",resnum,chainid));
378 : } else {
379 0 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("NE",resnum,chainid));
380 : }
381 2 : } else if( name=="chi4" && !isTerminalGroup("protein",resname) ) {
382 0 : if (!( resname=="ARG" || resname=="LYS" )) plumed_merror("chi-4 is defined only for ARG and LYS");
383 0 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("CG",resnum,chainid));
384 0 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("CD",resnum,chainid));
385 :
386 0 : if(resname=="ARG") {
387 0 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("NE",resnum,chainid));
388 0 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("CZ",resnum,chainid));
389 : } else {
390 0 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("CE",resnum,chainid));
391 0 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("NZ",resnum,chainid));
392 : }
393 2 : } else if( name=="chi5" && !isTerminalGroup("protein",resname) ) {
394 0 : if (!( resname=="ARG" )) plumed_merror("chi-5 is defined only for ARG");
395 0 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("CD",resnum,chainid));
396 0 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("NE",resnum,chainid));
397 0 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("CZ",resnum,chainid));
398 0 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("NH1",resnum,chainid));
399 3 : } else if( name=="sidechain" && !isTerminalGroup("protein",resname) ) {
400 1 : if(resname=="GLY") plumed_merror("sidechain is not defined for GLY");
401 1 : std::vector<AtomNumber> tmpnum1(mypdb.getAtomsInResidue(resnum,chainid));
402 15 : for(unsigned i=0; i<tmpnum1.size(); i++) {
403 66 : if(mypdb.getAtomName(tmpnum1[i])!="N"&&mypdb.getAtomName(tmpnum1[i])!="CA"&&
404 70 : mypdb.getAtomName(tmpnum1[i])!="C"&&mypdb.getAtomName(tmpnum1[i])!="O"&&
405 61 : mypdb.getAtomName(tmpnum1[i])!="HA"&&mypdb.getAtomName(tmpnum1[i])!="H"&&
406 59 : mypdb.getAtomName(tmpnum1[i])!="HN"&&mypdb.getAtomName(tmpnum1[i])!="H1"&&
407 59 : mypdb.getAtomName(tmpnum1[i])!="H2"&&mypdb.getAtomName(tmpnum1[i])!="H3"&&
408 59 : mypdb.getAtomName(tmpnum1[i])!="OXT"&&mypdb.getAtomName(tmpnum1[i])!="OT1"&&
409 73 : mypdb.getAtomName(tmpnum1[i])!="OT2"&&mypdb.getAtomName(tmpnum1[i])!="OC1"&&
410 32 : mypdb.getAtomName(tmpnum1[i])!="OC2") {
411 9 : numbers.push_back( tmpnum1[i] );
412 : }
413 : }
414 2 : } else if( name=="back" && !isTerminalGroup("protein",resname) ) {
415 1 : std::vector<AtomNumber> tmpnum1(mypdb.getAtomsInResidue(resnum,chainid));
416 15 : for(unsigned i=0; i<tmpnum1.size(); i++) {
417 66 : if(mypdb.getAtomName(tmpnum1[i])=="N"||mypdb.getAtomName(tmpnum1[i])=="CA"||
418 70 : mypdb.getAtomName(tmpnum1[i])=="C"||mypdb.getAtomName(tmpnum1[i])=="O"||
419 61 : mypdb.getAtomName(tmpnum1[i])=="HA"||mypdb.getAtomName(tmpnum1[i])=="H"||
420 59 : mypdb.getAtomName(tmpnum1[i])=="HN"||mypdb.getAtomName(tmpnum1[i])=="H1"||
421 59 : mypdb.getAtomName(tmpnum1[i])=="H2"||mypdb.getAtomName(tmpnum1[i])=="H3"||
422 59 : mypdb.getAtomName(tmpnum1[i])=="OXT"||mypdb.getAtomName(tmpnum1[i])=="OT1"||
423 59 : mypdb.getAtomName(tmpnum1[i])=="OT2"||mypdb.getAtomName(tmpnum1[i])=="OC1"||
424 59 : mypdb.getAtomName(tmpnum1[i])=="OC2"||mypdb.getAtomName(tmpnum1[i])=="HA1"||
425 63 : mypdb.getAtomName(tmpnum1[i])=="HA2"||mypdb.getAtomName(tmpnum1[i])=="HA3") {
426 5 : numbers.push_back( tmpnum1[i] );
427 : }
428 : }
429 0 : } else numbers.push_back(mypdb.getNamedAtomFromResidueAndChain(name,resnum,chainid));
430 :
431 494 : } else if( allowedResidue("rna",resname) || allowedResidue("dna",resname)) {
432 : std::string basetype;
433 480 : if(resname.find_first_of("A")!=std::string::npos) basetype+="A";
434 480 : if(resname.find_first_of("U")!=std::string::npos) basetype+="U";
435 480 : if(resname.find_first_of("T")!=std::string::npos) basetype+="T";
436 480 : if(resname.find_first_of("C")!=std::string::npos) basetype+="C";
437 480 : if(resname.find_first_of("G")!=std::string::npos) basetype+="G";
438 480 : plumed_massert(basetype.length()==1,"cannot find type of rna/dna residue "+resname+" "+basetype);
439 480 : if( name=="chi" ) {
440 72 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("O4\'",resnum,chainid));
441 72 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C1\'",resnum,chainid));
442 105 : if(basetype=="T" || basetype=="U" || basetype=="C") {
443 18 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("N1",resnum,chainid));
444 18 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C2",resnum,chainid));
445 51 : } else if(basetype=="G" || basetype=="A") {
446 54 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("N9",resnum,chainid));
447 54 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C4",resnum,chainid));
448 0 : } else plumed_error();
449 444 : } else if( name=="alpha" ) {
450 14 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("O3\'",resnum-1,chainid));
451 14 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("P",resnum,chainid));
452 14 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("O5\'",resnum,chainid));
453 14 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C5\'",resnum,chainid));
454 437 : } else if( name=="beta" ) {
455 14 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("P",resnum,chainid));
456 14 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("O5\'",resnum,chainid));
457 14 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C5\'",resnum,chainid));
458 14 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C4\'",resnum,chainid));
459 430 : } else if( name=="gamma" ) {
460 58 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("O5\'",resnum,chainid));
461 58 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C5\'",resnum,chainid));
462 58 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C4\'",resnum,chainid));
463 58 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C3\'",resnum,chainid));
464 401 : } else if( name=="delta" ) {
465 2 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C5\'",resnum,chainid));
466 2 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C4\'",resnum,chainid));
467 2 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C3\'",resnum,chainid));
468 2 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("O3\'",resnum,chainid));
469 400 : } else if( name=="epsilon" ) {
470 16 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C4\'",resnum,chainid));
471 16 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C3\'",resnum,chainid));
472 16 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("O3\'",resnum,chainid));
473 16 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("P",resnum+1,chainid));
474 392 : } else if( name=="zeta" ) {
475 14 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C3\'",resnum,chainid));
476 14 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("O3\'",resnum,chainid));
477 14 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("P",resnum+1,chainid));
478 14 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("O5\'",resnum+1,chainid));
479 385 : } else if( name=="v0" ) {
480 6 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C4\'",resnum,chainid));
481 6 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("O4\'",resnum,chainid));
482 6 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C1\'",resnum,chainid));
483 6 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C2\'",resnum,chainid));
484 382 : } else if( name=="v1" ) {
485 6 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("O4\'",resnum,chainid));
486 6 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C1\'",resnum,chainid));
487 6 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C2\'",resnum,chainid));
488 6 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C3\'",resnum,chainid));
489 379 : } else if( name=="v2" ) {
490 6 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C1\'",resnum,chainid));
491 6 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C2\'",resnum,chainid));
492 6 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C3\'",resnum,chainid));
493 6 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C4\'",resnum,chainid));
494 376 : } else if( name=="v3" ) {
495 6 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C2\'",resnum,chainid));
496 6 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C3\'",resnum,chainid));
497 6 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C4\'",resnum,chainid));
498 6 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("O4\'",resnum,chainid));
499 373 : } else if( name=="v4" ) {
500 6 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C3\'",resnum,chainid));
501 6 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C4\'",resnum,chainid));
502 6 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("O4\'",resnum,chainid));
503 6 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C1\'",resnum,chainid));
504 370 : } else if( name=="back" ) {
505 2 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("P",resnum,chainid));
506 2 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("O5\'",resnum,chainid));
507 2 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C5\'",resnum,chainid));
508 2 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C4\'",resnum,chainid));
509 2 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C3\'",resnum,chainid));
510 2 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("O3\'",resnum,chainid));
511 369 : } else if( name=="sugar" ) {
512 108 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C4\'",resnum,chainid));
513 108 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("O4\'",resnum,chainid));
514 108 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C1\'",resnum,chainid));
515 108 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C2\'",resnum,chainid));
516 108 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C3\'",resnum,chainid));
517 315 : } else if( name=="base" ) {
518 25 : if(basetype=="C") {
519 18 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("N1",resnum,chainid));
520 18 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C2",resnum,chainid));
521 18 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("O2",resnum,chainid));
522 18 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("N3",resnum,chainid));
523 18 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C4",resnum,chainid));
524 18 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("N4",resnum,chainid));
525 18 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C5",resnum,chainid));
526 18 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C6",resnum,chainid));
527 16 : } else if(basetype=="U") {
528 4 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("N1",resnum,chainid));
529 4 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C2",resnum,chainid));
530 4 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("O2",resnum,chainid));
531 4 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("N3",resnum,chainid));
532 4 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C4",resnum,chainid));
533 4 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("O4",resnum,chainid));
534 4 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C5",resnum,chainid));
535 4 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C6",resnum,chainid));
536 14 : } else if(basetype=="T") {
537 0 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("N1",resnum,chainid));
538 0 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C2",resnum,chainid));
539 0 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("O2",resnum,chainid));
540 0 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("N3",resnum,chainid));
541 0 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C4",resnum,chainid));
542 0 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("O4",resnum,chainid));
543 0 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C5",resnum,chainid));
544 0 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C7",resnum,chainid));
545 0 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C6",resnum,chainid));
546 14 : } else if(basetype=="G") {
547 6 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("N9",resnum,chainid));
548 6 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C4",resnum,chainid));
549 6 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("N3",resnum,chainid));
550 6 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C2",resnum,chainid));
551 6 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("N2",resnum,chainid));
552 6 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("N1",resnum,chainid));
553 6 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C6",resnum,chainid));
554 6 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("O6",resnum,chainid));
555 6 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C5",resnum,chainid));
556 6 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("N7",resnum,chainid));
557 6 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C8",resnum,chainid));
558 11 : } else if(basetype=="A") {
559 22 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("N9",resnum,chainid));
560 22 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C4",resnum,chainid));
561 22 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("N1",resnum,chainid));
562 22 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C2",resnum,chainid));
563 22 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("N3",resnum,chainid));
564 22 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C6",resnum,chainid));
565 22 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("N6",resnum,chainid));
566 22 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C5",resnum,chainid));
567 22 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("N7",resnum,chainid));
568 22 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C8",resnum,chainid));
569 0 : } else plumed_error();
570 290 : } else if( name=="lcs" ) {
571 568 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C2",resnum,chainid));
572 752 : if(basetype=="T" || basetype=="U" || basetype=="C") {
573 296 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C4",resnum,chainid));
574 296 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C6",resnum,chainid));
575 216 : } else if(basetype=="G" || basetype=="A") {
576 272 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C6",resnum,chainid));
577 272 : numbers.push_back(mypdb.getNamedAtomFromResidueAndChain("C4",resnum,chainid));
578 0 : } else plumed_error();
579 12 : } else numbers.push_back(mypdb.getNamedAtomFromResidueAndChain(name,resnum,chainid));
580 2 : } else numbers.push_back(mypdb.getNamedAtomFromResidueAndChain(name,resnum,chainid));
581 : }
582 : else {
583 0 : plumed_merror(type + " is not a valid molecule type");
584 : }
585 : }
586 :
587 : }
|
