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