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