Molmodel

RNA.h

Go to the documentation of this file.
00001 #ifndef SimTK_MOLMODEL_RNA_H_
00002 #define SimTK_MOLMODEL_RNA_H_
00003 
00004 /* -------------------------------------------------------------------------- *
00005  *                      SimTK Core: SimTK Molmodel                            *
00006  * -------------------------------------------------------------------------- *
00007  * This is part of the SimTK Core biosimulation toolkit originating from      *
00008  * Simbios, the NIH National Center for Physics-Based Simulation of           *
00009  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
00010  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
00011  *                                                                            *
00012  * Portions copyright (c) 2007 Stanford University and the Authors.           *
00013  * Authors: Michael Sherman, Christopher Bruns                                *
00014  * Contributors: Samuel Flores                                                *
00015  *                                                                            *
00016  * Permission is hereby granted, free of charge, to any person obtaining a    *
00017  * copy of this software and associated documentation files (the "Software"), *
00018  * to deal in the Software without restriction, including without limitation  *
00019  * the rights to use, copy, modify, merge, publish, distribute, sublicense,   *
00020  * and/or sell copies of the Software, and to permit persons to whom the      *
00021  * Software is furnished to do so, subject to the following conditions:       *
00022  *                                                                            *
00023  * The above copyright notice and this permission notice shall be included in *
00024  * all copies or substantial portions of the Software.                        *
00025  *                                                                            *
00026  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
00027  * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,   *
00028  * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL    *
00029  * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,    *
00030  * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR      *
00031  * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE  *
00032  * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
00033  * -------------------------------------------------------------------------- */
00034 
00035 
00036 #include "molmodel/internal/common.h"
00037 #include "molmodel/internal/Compound.h"
00038 
00039 namespace SimTK {
00040 
00041 // RiboseCore leaves 4 bond centers open:
00042 // 1) inboard center at O5* atom, for binding 5' phosphate
00043 // 2) 
00044 class SimTK_MOLMODEL_EXPORT RiboseCore : public BiopolymerResidue {
00045 public:
00046     RiboseCore(const String& name, const String& tlc, char olc) 
00047         : BiopolymerResidue(name, tlc, olc)
00048     {
00049         setBaseAtom(BivalentAtom("O5*", Element::Oxygen(), 120.90*Deg2Rad));
00050 
00051         bondAtom(AliphaticCarbon("C5*"), "O5*/bond2", 0.1423, 180*Deg2Rad);
00052         bondAtom(AliphaticHydrogen("H5*1"), "C5*/bond4");
00053         bondAtom(AliphaticHydrogen("H5*2"), "C5*/bond3");
00054 
00055         bondAtom(AliphaticCarbon("C4*"), "C5*/bond2", 0.1510); // length from (Gelbin et al 1996)
00056         bondAtom(AliphaticHydrogen("H4*"), "C4*/bond4");
00057 
00058         bondAtom(BivalentAtom("O4*", Element::Oxygen(), 112.0*Deg2Rad), "C4*/bond2", 0.1453);
00059 
00060         bondAtom(AliphaticCarbon("C3*"), "C4*/bond3", 0.1524);
00061         bondAtom(AliphaticHydrogen("H3*"), "C3*/bond4");
00062 
00063         bondAtom(BivalentAtom("O3*", Element::Oxygen(), 119.7*Deg2Rad), "C3*/bond2", 0.1423);
00064  
00065         bondAtom(AliphaticCarbon("C2*"), "C3*/bond3", 0.1525);
00066         bondAtom(AliphaticHydrogen("H2*"), "C2*/bond3");
00067 
00068         bondAtom(AliphaticCarbon("C1*"), "C2*/bond2", 0.1528);
00069         bondAtom(AliphaticHydrogen("H1*"), "C1*/bond4");
00070 
00071         addRingClosingBond("C1*/bond2", "O4*/bond2", 0.1414);
00072 
00073         // Ring strain reduces ring bond angles
00074         setDefaultBondAngle(109.6*Deg2Rad, "C1*", "O4*", "C4*"); 
00075         setDefaultBondAngle(105.5*Deg2Rad, "O4*", "C4*", "C3*"); 
00076         setDefaultBondAngle(102.7*Deg2Rad, "C4*", "C3*", "C2*"); 
00077         setDefaultBondAngle(101.5*Deg2Rad, "C3*", "C2*", "C1*"); 
00078         setDefaultBondAngle(106.4*Deg2Rad, "C2*", "C1*", "O4*");
00079 
00080         nameBondCenter("bondO5", "O5*/bond1");
00081         nameBondCenter("bondO3", "O3*/bond2");
00082         nameBondCenter("bondC1", "C1*/bond3");
00083 
00084         // RNA backbone dihedral angles
00085         // defineDihedralAngle("alpha", "bondO5", ""); // TODO requires phosphate
00086         //defineDihedralAngle("beta", "bondO5", "C5*/bond2");
00087         //defineDihedralAngle("gamma", "O5*", "C5*", "C4*", "C3*");
00088         //defineDihedralAngle("delta", "C5*", "C4*", "C3*", "C2*");
00089         //defineDihedralAngle("epsilon", "bondO3", "C3*/bond1");
00090         // TODO zeta requires next phosphate
00091 
00092         // ribose ring dihedral angles
00093         defineDihedralAngle( "nu0", "C4*", "O4*", "C1*", "C2*" );
00094         defineDihedralAngle( "nu1", "O4*", "C1*", "C2*", "C3*" );
00095         defineDihedralAngle( "nu2", "C1*", "C2*", "C3*", "C4*" );
00096         defineDihedralAngle( "nu3", "C2*", "C3*", "C4*", "O4*" );
00097         defineDihedralAngle( "nu4", "C3*", "C4*", "O4*", "C1*" );
00098 
00099         // Guestimated from Figure 2 of 
00100         // Schneider, B; Moravek, Z; and Berman, H.M. (2004)
00101         // "RNA conformational classes"
00102         // Nucleic Acids Research 32(5): 1666-1677
00103         // setDefaultDihedralAngle("alpha", -60*Deg2Rad);
00104         //setDefaultDihedralAngle("beta", 180*Deg2Rad);
00105         //setDefaultDihedralAngle("gamma", 60*Deg2Rad);
00106         //setDefaultDihedralAngle("delta", 80*Deg2Rad);
00107         //setDefaultDihedralAngle("epsilon", 190*Deg2Rad);
00108         // setDefaultDihedralAngle("zeta", -70*Deg2Rad);
00109 
00110         // nu angles selected from resdue U52 of 1EHZ
00111         setDefaultDihedralAngle("nu0", 2.2*Deg2Rad);
00112         setDefaultDihedralAngle("nu1", -24.5*Deg2Rad);
00113         setDefaultDihedralAngle("nu2", 35.8*Deg2Rad);
00114         setDefaultDihedralAngle("nu3", -35.9*Deg2Rad); // same bond as delta, different reference atoms
00115         setDefaultDihedralAngle("nu4", 21.6*Deg2Rad);
00116 
00117         nameBondCenter("bondNext", "bondO3");
00118         nameBondCenter("bondPrevious", "bondO5");
00119 
00120         setDefaultInboardBondLength(0.16100);
00121 
00122         // alternate atom names
00123         nameAtom("O5'", "O5*");
00124         nameAtom("C5'", "C5*");
00125         nameAtom("H5'", "H5*1");
00126         nameAtom("H5''", "H5*2");
00127         nameAtom("C4'", "C4*");
00128         nameAtom("H4'", "H4*");
00129         nameAtom("O4'", "O4*");
00130         nameAtom("C3'", "C3*");
00131         nameAtom("H3'", "H3*");
00132         nameAtom("O3'", "O3*");
00133         nameAtom("C2'", "C2*");
00134         nameAtom("H2'", "H2*");
00135         nameAtom("C1'", "C1*");
00136         nameAtom("H1'", "H1*");
00137     }
00138 
00139 };
00140 
00141 
00143 class SimTK_MOLMODEL_EXPORT FivePrimeRnaHydroxylGroup : public BiopolymerResidue
00144 {
00145 public:
00146     explicit FivePrimeRnaHydroxylGroup(String name, String threeLetterCode = "Unk", char oneLetterCode = '?') 
00147     : BiopolymerResidue(name, threeLetterCode, oneLetterCode)
00148     {
00149         instantiateBiotypes();
00150         
00151         setBaseAtom( UnivalentAtom("H5T", Element::Hydrogen()) );
00152         
00153         nameBondCenter("bondNext", "H5T/bond");
00154 
00155         setBiotypeIndex( "H5T", Biotype::get("Hydroxyl, RNA", "H5T", SimTK::Ordinality::Initial).getIndex() );
00156     }
00157     
00158     static void instantiateBiotypes() 
00159     {
00160         if ( ! Biotype::exists("Hydroxyl, RNA", "H5T", SimTK::Ordinality::Initial) )
00161             Biotype::defineBiotype(Element::Hydrogen(), 1, "Hydroxyl, RNA", "H5T", SimTK::Ordinality::Initial);
00162         
00163         // This oxygen biotype should get applied to the 5' oxygen of the nucleotide this hydroxyl attaches to
00164         if ( ! Biotype::exists("Hydroxyl, RNA", "O5*", SimTK::Ordinality::Initial) )
00165             Biotype::defineBiotype(Element::Oxygen(), 2, "Hydroxyl, RNA", "O5*", SimTK::Ordinality::Initial);
00166     }
00167     
00168 };
00169 
00170 
00172 class SimTK_MOLMODEL_EXPORT FivePrimeRnaPhosphateGroup : public BiopolymerResidue
00173 {
00174 public:
00175     explicit FivePrimeRnaPhosphateGroup(String name, String threeLetterCode = "Unk", char oneLetterCode = '?') 
00176     : BiopolymerResidue(name, threeLetterCode, oneLetterCode)
00177     {
00178         instantiateBiotypes();
00179         
00180         setBaseAtom( QuadrivalentAtom("P", Element::Phosphorus()) );
00181         
00182         bondAtom(UnivalentAtom("OP1", Element::Oxygen()), "P/bond4", 0.14800);
00183         bondAtom(UnivalentAtom("OP2", Element::Oxygen()), "P/bond3", 0.14800);
00184         bondAtom(UnivalentAtom("OP3", Element::Oxygen()), "P/bond2", 0.14800);
00185         
00186         nameBondCenter("bondNext", "P/bond1");
00187 
00188         setBiotypeIndex( "P", Biotype::get("Phosphate, RNA", "P", SimTK::Ordinality::Initial).getIndex() );
00189         setBiotypeIndex( "OP1", Biotype::get("Phosphate, RNA", "OP", SimTK::Ordinality::Initial).getIndex() );
00190         setBiotypeIndex( "OP2", Biotype::get("Phosphate, RNA", "OP", SimTK::Ordinality::Initial).getIndex() );
00191         setBiotypeIndex( "OP3", Biotype::get("Phosphate, RNA", "OP", SimTK::Ordinality::Initial).getIndex() );
00192     }
00193     
00194     static void instantiateBiotypes() 
00195     {
00196         if ( ! Biotype::exists("Phosphate, RNA", "P", SimTK::Ordinality::Initial) )
00197             Biotype::defineBiotype(Element::Phosphorus(), 4, "Phosphate, RNA", "P", SimTK::Ordinality::Initial);
00198         if ( ! Biotype::exists("Phosphate, RNA", "OP", SimTK::Ordinality::Initial) )
00199             Biotype::defineBiotype(Element::Oxygen(), 1, "Phosphate, RNA", "OP", SimTK::Ordinality::Initial);
00200     }
00201     
00202 };
00203 
00204 
00205 // Phosphodiester linkage between two RNA nucleotides
00206 class SimTK_MOLMODEL_EXPORT RnaPhosphodiesterLinkage : public BiopolymerResidue
00207 {
00208 public:
00209     explicit RnaPhosphodiesterLinkage(String name, String threeLetterCode = "Unk", char oneLetterCode = '?') 
00210         : BiopolymerResidue(name, threeLetterCode, oneLetterCode)
00211     {
00212         instantiateBiotypes();
00213 
00214         // TODO - set bond angles: 102.6 between ether oxygens, 119.9 between lone oxygen, 108.23 between mismatched pairs
00215         setBaseAtom( QuadrivalentAtom("P", Element::Phosphorus() 
00216                    ,  104.0*Deg2Rad // 03'-P-O5' from Gelbin et al
00217                    ,  107.9*Deg2Rad // 03'-P-OP? from Gelbin et al
00218                    ,  107.9*Deg2Rad // 03'-P-OP? from Gelbin et al
00219                    , -108.15*Deg2Rad // to make OP1-P-OP2 angle 119.6
00220                    ,  108.15*Deg2Rad // to make OP1-P-OP2 angle 119.6
00221                    ) );
00222 
00223         bondAtom(UnivalentAtom("OP1", Element::Oxygen()), "P/bond3", 0.14800);
00224         bondAtom(UnivalentAtom("OP2", Element::Oxygen()), "P/bond4", 0.14800);
00225 
00226         nameBondCenter("bondPrevious", "P/bond1");
00227         nameBondCenter("bondNext", "P/bond2");
00228 
00229         setDefaultInboardBondLength(0.16100);
00230 
00231         addCompoundSynonym("Phosphodiester, RNA");
00232 
00233         setBiotypeIndex("P", Biotype::get("Phosphodiester, RNA", "P").getIndex() );
00234         setBiotypeIndex("OP1", Biotype::get("Phosphodiester, RNA", "OP").getIndex() );
00235         setBiotypeIndex("OP2", Biotype::get("Phosphodiester, RNA", "OP").getIndex() );
00236     }
00237     
00238     
00239     static void instantiateBiotypes() 
00240     {
00241         if ( ! Biotype::exists("Phosphodiester, RNA", "P") )
00242             Biotype::defineBiotype(Element::Phosphorus(), 4, "Phosphodiester, RNA", "P");
00243         if ( ! Biotype::exists("Phosphodiester, RNA", "OP") )
00244             Biotype::defineBiotype(Element::Oxygen(), 1, "Phosphodiester, RNA", "OP");
00245     }
00246 };
00247 
00248 
00250 class SimTK_MOLMODEL_EXPORT ThreePrimeRnaHydroxylGroup : public BiopolymerResidue
00251 {
00252 public:
00253     explicit ThreePrimeRnaHydroxylGroup(String name, String threeLetterCode = "Unk", char oneLetterCode = '?') 
00254     : BiopolymerResidue(name, threeLetterCode, oneLetterCode)
00255     {
00256         instantiateBiotypes();
00257         
00258         setBaseAtom( UnivalentAtom("H3T", Element::Hydrogen()) );
00259         
00260         nameBondCenter("bondNext", "H3T/bond");
00261 
00262         setBiotypeIndex( "H3T", Biotype::get("Hydroxyl, RNA", "H3T", SimTK::Ordinality::Final).getIndex() );
00263     }
00264     
00265     static void instantiateBiotypes() 
00266     {
00267         if ( ! Biotype::exists("Hydroxyl, RNA", "H3T", SimTK::Ordinality::Final) )
00268             Biotype::defineBiotype(Element::Hydrogen(), 1, "Hydroxyl, RNA", "H3T", SimTK::Ordinality::Final);
00269         
00270         // This oxygen biotype should get applied to the 3' oxygen of the nucleotide this hydroxyl attaches to
00271         if ( ! Biotype::exists("Hydroxyl, RNA", "O3*", SimTK::Ordinality::Final) )
00272             Biotype::defineBiotype(Element::Oxygen(), 2, "Hydroxyl, RNA", "O3*", SimTK::Ordinality::Final);
00273     }
00274     
00275 };
00276 
00278 class SimTK_MOLMODEL_EXPORT ThreePrimeRnaPhosphateGroup : public BiopolymerResidue
00279 {
00280 public:
00281     explicit ThreePrimeRnaPhosphateGroup(String name, String threeLetterCode = "Unk", char oneLetterCode = '?') 
00282     : BiopolymerResidue(name, threeLetterCode, oneLetterCode)
00283     {
00284         instantiateBiotypes();
00285         
00286         setBaseAtom( QuadrivalentAtom("P", Element::Phosphorus()) );
00287         
00288         bondAtom(UnivalentAtom("OP1", Element::Oxygen()), "P/bond4", 0.14800);
00289         bondAtom(UnivalentAtom("OP2", Element::Oxygen()), "P/bond3", 0.14800);
00290         bondAtom(UnivalentAtom("OP3", Element::Oxygen()), "P/bond2", 0.14800);
00291         
00292         nameBondCenter("bondNext", "P/bond1");
00293 
00294         setBiotypeIndex( "P", Biotype::get("Phosphate, RNA", "P", SimTK::Ordinality::Final).getIndex() );
00295         setBiotypeIndex( "OP1", Biotype::get("Phosphate, RNA", "OP", SimTK::Ordinality::Final).getIndex() );
00296         setBiotypeIndex( "OP2", Biotype::get("Phosphate, RNA", "OP", SimTK::Ordinality::Final).getIndex() );
00297         setBiotypeIndex( "OP3", Biotype::get("Phosphate, RNA", "OP", SimTK::Ordinality::Final).getIndex() );
00298     }
00299     
00300     static void instantiateBiotypes() 
00301     {
00302         if ( ! Biotype::exists("Phosphate, RNA", "P", SimTK::Ordinality::Final) )
00303             Biotype::defineBiotype(Element::Phosphorus(), 4, "Phosphate, RNA", "P", SimTK::Ordinality::Final);
00304         if ( ! Biotype::exists("Phosphate, RNA", "OP", SimTK::Ordinality::Final) )
00305             Biotype::defineBiotype(Element::Oxygen(), 1, "Phosphate, RNA", "OP", SimTK::Ordinality::Final);
00306     }
00307     
00308 };
00309 
00310 
00311 // Nucleoside has sugar and base but not phosphate
00312 class SimTK_MOLMODEL_EXPORT RibonucleosideResidue : public RiboseCore {
00313 public:
00314     explicit RibonucleosideResidue(String name, String threeLetterCode = "Unk", char oneLetterCode = '?')
00315         : RiboseCore(name, threeLetterCode, oneLetterCode)
00316     {
00317         // 2' hydroxyl group present in RNA but not in DNA
00318         bondAtom(BivalentAtom("O2*", Element::Oxygen(), 108.50*Deg2Rad), "C2*/bond4", 0.1413, -60*Deg2Rad);
00319         bondAtom(UnivalentAtom("2HO*", Element::Hydrogen()), "O2*/bond2", 0.0960);
00320 
00321         // alternate atom name
00322         nameAtom("O2'", "O2*");
00323         nameAtom("HO2'", "2HO*");
00324     }
00325     
00326     // Create larger residue with phosphodiester prepended
00327     BiopolymerResidue withPhosphodiester() const 
00328     {
00329         const RibonucleosideResidue& residue = *this;
00330         
00331         RnaPhosphodiesterLinkage po2(residue.getCompoundName(), residue.getPdbResidueName(), residue.getOneLetterCode());
00332         po2.bondCompound("nucleoside", residue, "bondNext");
00333         po2.inheritAtomNames("nucleoside");
00334         po2.inheritCompoundSynonyms(residue);
00335         po2.nameBondCenter("bondNext", "nucleoside/bondNext");
00336 
00337         po2.setPdbResidueNumber(residue.getPdbResidueNumber());
00338         po2.setPdbChainId(residue.getPdbChainId());
00339 
00340         po2.defineDihedralAngle("alpha", "bondPrevious", "O5*/bond2");
00341         po2.setDefaultDihedralAngle("alpha", 295.0*Deg2Rad); // A-RNA (Schneider et al 2004)
00342 
00343         return po2;
00344     }
00345 };
00346 
00347 // Nucleotide has sugar and base and phosphate
00348 class SimTK_MOLMODEL_EXPORT RibonucleotideResidue : public RibonucleosideResidue {
00349 public:
00350     // Factory method
00351     static RibonucleotideResidue create(char oneLetterCode);
00352     static RibonucleotideResidue create(const PdbResidue& pdbResidue, char chainId = ' ');
00353 
00354     explicit RibonucleotideResidue(String name, String threeLetterCode = "Unk", char oneLetterCode = '?')
00355         : RibonucleosideResidue(name, threeLetterCode, oneLetterCode)
00356     {
00357         defineDihedralAngle("beta", "O5*/bond1", "C5*/bond2");
00358         setDefaultDihedralAngle("beta", 173.0*Deg2Rad); //  A-RNA (Schneider et al 2004)
00359 
00360         defineDihedralAngle("gamma", "C5*/bond1", "C4*/bond3");
00361         setDefaultDihedralAngle("gamma", 54.0*Deg2Rad); //  A-RNA (Schneider et al 2004)
00362 
00363         defineDihedralAngle("delta", "C4*/bond1", "C3*/bond2");
00364         setDefaultDihedralAngle("delta", 80.0*Deg2Rad); //  A-RNA (Schneider et al 2004)
00365 
00366         defineDihedralAngle("epsilon", "C3*/bond1", "O3*/bond2");
00367         setDefaultDihedralAngle("epsilon", 210.0*Deg2Rad); //  A-RNA (Schneider et al 2004)
00368 
00369         nameBondCenter("bondBase", "C1*/bond3");
00370     }
00371 
00372 
00373     class Adenylate;
00374     class Cytidylate;
00375     class Guanylate;
00376     class Uridylate;
00377     
00378     // Modified residues from tRNA
00379     class TwoNMethylGuanylate;
00380 };
00381 
00382 
00383 // PurineBaseCore represents the atoms in common among purine bases such as Adenine and Guanine
00384 class PurineBaseCore : public Compound {
00385 public:
00386     PurineBaseCore()
00387     {
00388         setBaseAtom( TrivalentAtom("N9",  Element::Nitrogen(), 125.8*Deg2Rad, 128.8*Deg2Rad));
00389 
00390         bondAtom( TrivalentAtom("C8",  Element::Carbon(),   113.9*Deg2Rad, 123.05*Deg2Rad),  "N9/bond3", 0.13710, 180*Deg2Rad,   BondMobility::Rigid );
00391 
00392         bondAtom( UnivalentAtom("H8",  Element::Hydrogen()),                                 "C8/bond3",  0.10800 );
00393         bondAtom(  BivalentAtom("N7",  Element::Nitrogen(), 103.8*Deg2Rad),                  "C8/bond2",  0.13040, 0*Deg2Rad,    BondMobility::Rigid );
00394         bondAtom( TrivalentAtom("C5",  Element::Carbon(),   110.4*Deg2Rad, 132.40*Deg2Rad),  "N7/bond2",  0.13910, 0*Deg2Rad,    BondMobility::Rigid );
00395         
00396         bondAtom( TrivalentAtom("C4",  Element::Carbon(),   126.2*Deg2Rad, 106.2*Deg2Rad),   "N9/bond2",  0.13740, 0*Deg2Rad,    BondMobility::Rigid );
00397 
00398         addRingClosingBond("C5/bond2",                                                       "C4/bond3",  0.13700, 0*Deg2Rad,    BondMobility::Rigid );
00399 
00400         // 114.30 is average between amber angle for guanine(111.30) and adenine(117.30)
00401         bondAtom( TrivalentAtom("C6",  Element::Carbon(), 114.30*Deg2Rad, 122.85*Deg2Rad),   "C5/bond3",  0.14040, 180*Deg2Rad,  BondMobility::Rigid );
00402 
00403         bondAtom(  BivalentAtom("N3",  Element::Nitrogen(), 118.60*Deg2Rad),                 "C4/bond2",  0.13540, 180*Deg2Rad,  BondMobility::Rigid );
00404         bondAtom( TrivalentAtom("C2",  Element::Carbon()),                                   "N3/bond2",  0.13240, 0*Deg2Rad,    BondMobility::Rigid );
00405 
00406         // Two bond centers open, inboard (N9), and C6
00407         nameBondCenter("bondC6", "C6/bond3");
00408         nameBondCenter("bondC2", "C2/bond3");
00409     }
00410 };
00411 
00412 
00413 class AdenineBase : public PurineBaseCore
00414 {
00415 public:
00416     AdenineBase()
00417     {
00418         // Lone pair on N1
00419         bondAtom(  BivalentAtom("N1",  Element::Nitrogen(), 118.60*Deg2Rad), "C6/bond2",  0.13390, 0*Deg2Rad, BondMobility::Rigid );
00420         addRingClosingBond("N1/bond2",                                       "C2/bond2",  0.13240, 0*Deg2Rad, BondMobility::Rigid);
00421 
00422         // Hydrogen on C2
00423         bondAtom( UnivalentAtom("H2",  Element::Hydrogen()),                 "bondC2",    0.10800,            BondMobility::Rigid );
00424 
00425         // Amine on C6
00426         bondAtom( TrivalentAtom("N6",  Element::Nitrogen()),                 "bondC6",    0.13400, 0*Deg2Rad, BondMobility::Rigid );
00427         bondAtom( UnivalentAtom("H61", Element::Hydrogen()),                 "N6/bond2",  0.10100,            BondMobility::Rigid );
00428         bondAtom( UnivalentAtom("H62", Element::Hydrogen()),                 "N6/bond3",  0.10100,            BondMobility::Rigid );
00429     }
00430 };
00431 
00432 class GuanineBase : public PurineBaseCore
00433 {
00434 public:
00435     GuanineBase()
00436     {
00437         // Hydrogen on N1
00438         bondAtom( TrivalentAtom("N1",  Element::Nitrogen(), 125.20*Deg2Rad, 116.8*Deg2Rad), "C6/bond2",  0.13880, 0*Deg2Rad, BondMobility::Rigid );
00439         bondAtom(UnivalentAtom("H1", Element::Hydrogen()),   "N1/bond3", 0.10100, BondMobility::Rigid);
00440         addRingClosingBond("N1/bond2",                       "C2/bond2",  0.13240, 0*Deg2Rad, BondMobility::Rigid);
00441 
00442         // Amine on C2
00443         bondAtom( TrivalentAtom("N2",  Element::Nitrogen()), "bondC2",    0.13400, 0*Deg2Rad, BondMobility::Rigid );
00444         bondAtom( UnivalentAtom("H21", Element::Hydrogen()), "N2/bond2",  0.10100,            BondMobility::Rigid );
00445         bondAtom( UnivalentAtom("H22", Element::Hydrogen()), "N2/bond3",  0.10100,            BondMobility::Rigid );
00446 
00447         // Carbonyl on C6
00448         bondAtom( UnivalentAtom("O6",  Element::Oxygen()),   "bondC6",    0.12290,            BondMobility::Rigid );
00449 
00450     }
00451 };
00452 
00453 class TwoNMethylGuanineBaseGroup : public PurineBaseCore
00454 {
00455 public:
00456     TwoNMethylGuanineBaseGroup() {
00457         addCompoundSynonym("Guanosine"); // for resolving biotypes KLUDGE
00458         // TODO
00459     }
00460 };
00461 
00462 
00463 // PyrimidineBaseCore represents the atoms in common among pyrimidine bases such as Cytosine, Uracil, and Thymine
00464 class PyrimidineBaseCore : public Compound {
00465 public:
00466     PyrimidineBaseCore()
00467     {
00468     }
00469 };
00470 
00471 
00472 class CytosineBase : public PyrimidineBaseCore
00473 {
00474 public:
00475     CytosineBase()
00476     {
00477         setBaseAtom(   TrivalentAtom("N1",  Element::Nitrogen(), 125.8*Deg2Rad, 128.8*Deg2Rad));
00478     //scf added temporarily
00479         //bondAtom( UnivalentAtom("HN1", Element::Hydrogen()), "N1/bond4", 0.1010, 0*Deg2Rad, BondMobility::Rigid);
00480 
00481         bondAtom( TrivalentAtom("C2", Element::Carbon(), 118.6*Deg2Rad, 120.9*Deg2Rad),  "N1/bond2", 0.1383, 180*Deg2Rad,   BondMobility::Rigid );
00482         bondAtom( UnivalentAtom("O2", Element::Oxygen()), "C2/bond3", 0.1229, 0*Deg2Rad, BondMobility::Rigid);
00483 
00484         bondAtom( BivalentAtom("N3", Element::Nitrogen(), 120.5*Deg2Rad), "C2/bond2", 0.1358, 0*Deg2Rad, BondMobility::Rigid);
00485 
00486         bondAtom( TrivalentAtom("C4", Element::Carbon(), 121.5*Deg2Rad, 119.3*Deg2Rad), "N3/bond2", 0.1339, 0*Deg2Rad, BondMobility::Rigid);
00487         bondAtom( TrivalentAtom("N4", Element::Nitrogen(), 120.0*Deg2Rad, 120.0*Deg2Rad), "C4/bond3", 0.1340, 0*Deg2Rad, BondMobility::Rigid);
00488         bondAtom( UnivalentAtom("H41", Element::Hydrogen()), "N4/bond2", 0.1010, 0*Deg2Rad, BondMobility::Rigid);
00489         bondAtom( UnivalentAtom("H42", Element::Hydrogen()), "N4/bond3", 0.1010, 0*Deg2Rad, BondMobility::Rigid);
00490 
00491         bondAtom( TrivalentAtom("C5", Element::Carbon(), 117.0*Deg2Rad, 119.70*Deg2Rad), "C4/bond2", 0.1433, 0*Deg2Rad, BondMobility::Rigid);
00492         bondAtom( UnivalentAtom("H5", Element::Hydrogen()), "C5/bond3", 0.1080, 0*Deg2Rad, BondMobility::Rigid);
00493 
00494         bondAtom( TrivalentAtom("C6", Element::Carbon(), 121.20*Deg2Rad, 119.70*Deg2Rad), "C5/bond2", 0.1350, 0*Deg2Rad, BondMobility::Rigid);
00495         bondAtom( UnivalentAtom("H6", Element::Hydrogen()), "C6/bond3", 0.1080, 0*Deg2Rad, BondMobility::Rigid);
00496 
00497         addRingClosingBond("C6/bond2", "N1/bond3", 0.1365, 0*Deg2Rad, BondMobility::Rigid);
00498     }
00499 };
00500 
00501 class UracilBase : public PyrimidineBaseCore
00502 {
00503 public:
00504     UracilBase()
00505     {
00506         setBaseAtom( TrivalentAtom("N1",  Element::Nitrogen(), 125.8*Deg2Rad, 128.8*Deg2Rad));
00507 
00508         bondAtom( TrivalentAtom("C2", Element::Carbon(), 118.6*Deg2Rad, 120.9*Deg2Rad),  "N1/bond2", 0.1383, 180*Deg2Rad,   BondMobility::Rigid );
00509         bondAtom( UnivalentAtom("O2", Element::Oxygen()), "C2/bond3", 0.1229, 0*Deg2Rad, BondMobility::Rigid);
00510 
00511         bondAtom(TrivalentAtom("N3", Element::Nitrogen(), 126.4*Deg2Rad, 116.8*Deg2Rad), "C2/bond2", 0.1358, 0*Deg2Rad, BondMobility::Rigid);
00512         bondAtom(UnivalentAtom("H3", Element::Hydrogen()), "N3/bond3", 0.1010, 0*Deg2Rad, BondMobility::Rigid);
00513 
00514         bondAtom( TrivalentAtom("C4", Element::Carbon(), 114.0*Deg2Rad, 120.6*Deg2Rad),  "N3/bond2", 0.1388, 0*Deg2Rad,   BondMobility::Rigid );
00515         bondAtom( UnivalentAtom("O4", Element::Oxygen()), "C4/bond3", 0.1229, 0*Deg2Rad, BondMobility::Rigid);
00516 
00517         bondAtom( TrivalentAtom("C5", Element::Carbon(), 120.7*Deg2Rad, 119.70*Deg2Rad), "C4/bond2", 0.1444, 0*Deg2Rad, BondMobility::Rigid);
00518         bondAtom( UnivalentAtom("H5", Element::Hydrogen()), "C5/bond3", 0.1080, 0*Deg2Rad, BondMobility::Rigid);
00519 
00520         bondAtom( TrivalentAtom("C6", Element::Carbon(), 121.20*Deg2Rad, 119.70*Deg2Rad), "C5/bond2", 0.1350, 0*Deg2Rad, BondMobility::Rigid);
00521         bondAtom( UnivalentAtom("H6", Element::Hydrogen()), "C6/bond3", 0.1080, 0*Deg2Rad, BondMobility::Rigid);
00522 
00523         addRingClosingBond("C6/bond2", "N1/bond3", 0.1365, 0*Deg2Rad, BondMobility::Rigid);
00524     }
00525 };
00526 
00527 // TODO - these should not be classes, e.g. Adenylate should be an instance of RibonucleotideResidue
00528 class  RibonucleotideResidue::Adenylate : public RibonucleotideResidue {
00529 public:
00530     Adenylate() : RibonucleotideResidue("adenylate", "A  ", 'A') 
00531     {
00532         bondCompound("base", AdenineBase(), "bondBase", 0.14750, 200*Deg2Rad, BondMobility::Torsion);
00533         inheritAtomNames("base");
00534 
00535         defineDihedralAngle("chi", "O4*", "C1*", "N9", "C4");
00536         setDefaultDihedralAngle("chi", 199.0*Deg2Rad);
00537 
00538         addCompoundSynonym("Adenosine"); // for resolving biotypes
00539     }
00540 };
00541 
00542 class  RibonucleotideResidue::Guanylate : public RibonucleotideResidue {
00543 public:
00544     Guanylate() : RibonucleotideResidue("guanylate", "G  ", 'G') 
00545     {
00546         bondCompound("base", GuanineBase(), "bondBase", 0.14710, 200*Deg2Rad, BondMobility::Torsion);
00547         inheritAtomNames("base");
00548 
00549         defineDihedralAngle("chi", "O4*", "C1*", "N9", "C4");
00550         setDefaultDihedralAngle("chi", 199.0*Deg2Rad);
00551 
00552         addCompoundSynonym("Guanosine"); // for resolving biotypes
00553     }
00554 };
00555 
00556 class  RibonucleotideResidue::Cytidylate : public RibonucleotideResidue {
00557 public:
00558     Cytidylate() : RibonucleotideResidue("cytidylate", "C  ", 'C') 
00559     {
00560         bondCompound("base", CytosineBase(), "bondBase", 0.14710, 200*Deg2Rad, BondMobility::Torsion);
00561         inheritAtomNames("base");
00562 
00563         defineDihedralAngle("chi", "O4*", "C1*", "N1", "C2");
00564         setDefaultDihedralAngle("chi", 199.0*Deg2Rad);
00565 
00566         addCompoundSynonym("Cytidine"); // for resolving biotypes
00567     }
00568 };
00569 
00570 class  RibonucleotideResidue::Uridylate : public RibonucleotideResidue {
00571 public:
00572     Uridylate() : RibonucleotideResidue("uridylate", "U  ", 'U') 
00573     {
00574         bondCompound("base", UracilBase(), "bondBase", 0.14710, 200*Deg2Rad, BondMobility::Torsion);
00575         inheritAtomNames("base");
00576 
00577         defineDihedralAngle("chi", "O4*", "C1*", "N1", "C2");
00578         setDefaultDihedralAngle("chi", 199.0*Deg2Rad);
00579 
00580         addCompoundSynonym("Uridine"); // for resolving biotypes
00581     }
00582 };
00583 
00584 
00585 class TwoNMethylGuanidineGroup : public Compound {
00586 public:
00587     TwoNMethylGuanidineGroup()
00588     {
00589         instantiateBiotypes();
00590 
00591         setPdbResidueName("2MG");
00592 
00593         setCompoundName("TwoNMethylGuanidineGroup");
00594 
00595         setBaseAtom( TrivalentAtom("N9", Element::Nitrogen(), 127.767*Deg2Rad, 112.885*Deg2Rad) );
00596 
00597         bondAtom( TrivalentAtom("C8", Element::Carbon(), 113.564*Deg2Rad, 120.934*Deg2Rad), "N9/bond2", 0.13740, 180.00*Deg2Rad, BondMobility::Rigid );
00598         bondAtom( BivalentAtom("N7", Element::Nitrogen(), 104.422*Deg2Rad), "C8/bond2", 0.12800, 0.02*Deg2Rad, BondMobility::Rigid );
00599         bondAtom( TrivalentAtom("C5", Element::Carbon(), 130.918*Deg2Rad, 112.885*Deg2Rad), "N7/bond2", 0.13760, 179.54*Deg2Rad, BondMobility::Rigid );
00600         bondAtom( TrivalentAtom("C6", Element::Carbon(), 131.213*Deg2Rad, 109.672*Deg2Rad), "C5/bond2", 0.14350, -0.17*Deg2Rad, BondMobility::Rigid );
00601         bondAtom( UnivalentAtom("O6", Element::Oxygen()), "C6/bond2", 0.11940 );
00602         bondAtom( TrivalentAtom("N1", Element::Nitrogen(), 126.495*Deg2Rad, 113.649*Deg2Rad), "C6/bond3", 0.14150, 0.32*Deg2Rad, BondMobility::Rigid );
00603         bondAtom( TrivalentAtom("C2", Element::Carbon(), 116.091*Deg2Rad, 123.417*Deg2Rad), "N1/bond2", 0.13620, 178.30*Deg2Rad, BondMobility::Rigid );
00604         bondAtom( TrivalentAtom("N2", Element::Nitrogen(), 116.703*Deg2Rad, 120.926*Deg2Rad), "C2/bond2", 0.13550, 20.71*Deg2Rad, BondMobility::Rigid );
00605         bondAtom( UnivalentAtom("1H2", Element::Hydrogen()), "N2/bond2", 0.09940 );
00606         bondAtom( QuadrivalentAtom("C10", Element::Carbon()), "N2/bond3", 0.14490, 176.13*Deg2Rad, BondMobility::Torsion );
00607         bondAtom( UnivalentAtom("H20", Element::Hydrogen()), "C10/bond2", 0.10810 );
00608         bondAtom( UnivalentAtom("H21", Element::Hydrogen()), "C10/bond3", 0.10850 );
00609         bondAtom( UnivalentAtom("H22", Element::Hydrogen()), "C10/bond4", 0.10790 );
00610         bondAtom( BivalentAtom("N3", Element::Nitrogen(), 112.885*Deg2Rad), "C2/bond3", 0.12910, 0.73*Deg2Rad, BondMobility::Rigid );
00611         bondAtom( TrivalentAtom("C4", Element::Carbon(), 109.368*Deg2Rad, 104.422*Deg2Rad), "N3/bond2", 0.13560, -119.37*Deg2Rad, BondMobility::Rigid );
00612         bondAtom( UnivalentAtom("H1", Element::Hydrogen()), "N1/bond3", 0.09980 );
00613         bondAtom( UnivalentAtom("H8", Element::Hydrogen()), "C8/bond3", 0.10710 );
00614 
00615         addRingClosingBond( "C4/bond2", "N9/bond3", 0.15);
00616         addRingClosingBond( "C4/bond3", "C5/bond3", 0.15);
00617 
00618         setBiotypeIndex( "N9", Biotype::get("TwoNMethylGuanidineGroup", "N9").getIndex());
00619         setBiotypeIndex( "C8", Biotype::get("TwoNMethylGuanidineGroup", "C8").getIndex());
00620         setBiotypeIndex( "N7", Biotype::get("TwoNMethylGuanidineGroup", "N7").getIndex());
00621         setBiotypeIndex( "C5", Biotype::get("TwoNMethylGuanidineGroup", "C5").getIndex());
00622         setBiotypeIndex( "C6", Biotype::get("TwoNMethylGuanidineGroup", "C6").getIndex());
00623         setBiotypeIndex( "O6", Biotype::get("TwoNMethylGuanidineGroup", "O6").getIndex());
00624         setBiotypeIndex( "N1", Biotype::get("TwoNMethylGuanidineGroup", "N1").getIndex());
00625         setBiotypeIndex( "C2", Biotype::get("TwoNMethylGuanidineGroup", "C2").getIndex());
00626         setBiotypeIndex( "N2", Biotype::get("TwoNMethylGuanidineGroup", "N2").getIndex());
00627         setBiotypeIndex( "1H2", Biotype::get("TwoNMethylGuanidineGroup", "1H2").getIndex());
00628         setBiotypeIndex( "C10", Biotype::get("TwoNMethylGuanidineGroup", "C10").getIndex());
00629         setBiotypeIndex( "H20", Biotype::get("TwoNMethylGuanidineGroup", "H20").getIndex());
00630         setBiotypeIndex( "H21", Biotype::get("TwoNMethylGuanidineGroup", "H21").getIndex());
00631         setBiotypeIndex( "H22", Biotype::get("TwoNMethylGuanidineGroup", "H22").getIndex());
00632         setBiotypeIndex( "N3", Biotype::get("TwoNMethylGuanidineGroup", "N3").getIndex());
00633         setBiotypeIndex( "C4", Biotype::get("TwoNMethylGuanidineGroup", "C4").getIndex());
00634         setBiotypeIndex( "H1", Biotype::get("TwoNMethylGuanidineGroup", "H1").getIndex());
00635         setBiotypeIndex( "H8", Biotype::get("TwoNMethylGuanidineGroup", "H8").getIndex());
00636 
00637         setDefaultDihedralAngle(0.020*Deg2Rad, "C5", "N7", "C8", "N9");
00638         setDefaultDihedralAngle(179.544*Deg2Rad, "C6", "C5", "N7", "C8");
00639         setDefaultDihedralAngle(-0.165*Deg2Rad, "O6", "C6", "C5", "N7");
00640         setDefaultDihedralAngle(-179.886*Deg2Rad, "N1", "C6", "C5", "N7");
00641         setDefaultDihedralAngle(0.318*Deg2Rad, "C2", "N1", "C6", "C5");
00642         setDefaultDihedralAngle(178.303*Deg2Rad, "N2", "C2", "N1", "C6");
00643         setDefaultDihedralAngle(20.715*Deg2Rad, "1H2", "N2", "C2", "N1");
00644         setDefaultDihedralAngle(173.145*Deg2Rad, "C10", "N2", "C2", "N1");
00645         setDefaultDihedralAngle(176.133*Deg2Rad, "H20", "C10", "N2", "C2");
00646         setDefaultDihedralAngle(-63.298*Deg2Rad, "H21", "C10", "N2", "C2");
00647         setDefaultDihedralAngle(57.486*Deg2Rad, "H22", "C10", "N2", "C2");
00648         setDefaultDihedralAngle(-0.550*Deg2Rad, "N3", "C2", "N1", "C6");
00649         setDefaultDihedralAngle(0.727*Deg2Rad, "C4", "N3", "C2", "N1");
00650         setDefaultDihedralAngle(176.522*Deg2Rad, "H1", "N1", "C6", "C5");
00651     } // end constructor TwoNMethylGuanidineGroup 
00652 
00653     static void instantiateBiotypes() {
00654         if (! Biotype::exists("TwoNMethylGuanidineGroup", "N9") )
00655             Biotype::defineBiotype(Element::Nitrogen(), 3, "TwoNMethylGuanidineGroup", "N9"); 
00656         if (! Biotype::exists("TwoNMethylGuanidineGroup", "C8") )
00657             Biotype::defineBiotype(Element::Carbon(), 3, "TwoNMethylGuanidineGroup", "C8"); 
00658         if (! Biotype::exists("TwoNMethylGuanidineGroup", "N7") )
00659             Biotype::defineBiotype(Element::Nitrogen(), 2, "TwoNMethylGuanidineGroup", "N7"); 
00660         if (! Biotype::exists("TwoNMethylGuanidineGroup", "C5") )
00661             Biotype::defineBiotype(Element::Carbon(), 3, "TwoNMethylGuanidineGroup", "C5"); 
00662         if (! Biotype::exists("TwoNMethylGuanidineGroup", "C6") )
00663             Biotype::defineBiotype(Element::Carbon(), 3, "TwoNMethylGuanidineGroup", "C6"); 
00664         if (! Biotype::exists("TwoNMethylGuanidineGroup", "O6") )
00665             Biotype::defineBiotype(Element::Oxygen(), 1, "TwoNMethylGuanidineGroup", "O6"); 
00666         if (! Biotype::exists("TwoNMethylGuanidineGroup", "N1") )
00667             Biotype::defineBiotype(Element::Nitrogen(), 3, "TwoNMethylGuanidineGroup", "N1"); 
00668         if (! Biotype::exists("TwoNMethylGuanidineGroup", "C2") )
00669             Biotype::defineBiotype(Element::Carbon(), 3, "TwoNMethylGuanidineGroup", "C2"); 
00670         if (! Biotype::exists("TwoNMethylGuanidineGroup", "N2") )
00671             Biotype::defineBiotype(Element::Nitrogen(), 3, "TwoNMethylGuanidineGroup", "N2"); 
00672         if (! Biotype::exists("TwoNMethylGuanidineGroup", "1H2") )
00673             Biotype::defineBiotype(Element::Hydrogen(), 1, "TwoNMethylGuanidineGroup", "1H2"); 
00674         if (! Biotype::exists("TwoNMethylGuanidineGroup", "C10") )
00675             Biotype::defineBiotype(Element::Carbon(), 4, "TwoNMethylGuanidineGroup", "C10"); 
00676         if (! Biotype::exists("TwoNMethylGuanidineGroup", "H20") )
00677             Biotype::defineBiotype(Element::Hydrogen(), 1, "TwoNMethylGuanidineGroup", "H20"); 
00678         if (! Biotype::exists("TwoNMethylGuanidineGroup", "H21") )
00679             Biotype::defineBiotype(Element::Hydrogen(), 1, "TwoNMethylGuanidineGroup", "H21"); 
00680         if (! Biotype::exists("TwoNMethylGuanidineGroup", "H22") )
00681             Biotype::defineBiotype(Element::Hydrogen(), 1, "TwoNMethylGuanidineGroup", "H22"); 
00682         if (! Biotype::exists("TwoNMethylGuanidineGroup", "N3") )
00683             Biotype::defineBiotype(Element::Nitrogen(), 2, "TwoNMethylGuanidineGroup", "N3"); 
00684         if (! Biotype::exists("TwoNMethylGuanidineGroup", "C4") )
00685             Biotype::defineBiotype(Element::Carbon(), 3, "TwoNMethylGuanidineGroup", "C4"); 
00686         if (! Biotype::exists("TwoNMethylGuanidineGroup", "H1") )
00687             Biotype::defineBiotype(Element::Hydrogen(), 1, "TwoNMethylGuanidineGroup", "H1"); 
00688         if (! Biotype::exists("TwoNMethylGuanidineGroup", "H8") )
00689             Biotype::defineBiotype(Element::Hydrogen(), 1, "TwoNMethylGuanidineGroup", "H8"); 
00690     } // end instantiateBiotypes
00691 
00692 }; // end class TwoNMethylGuanidineGroup 
00693 
00694 
00695 class  RibonucleotideResidue::TwoNMethylGuanylate : public RibonucleotideResidue {
00696 public:
00697     TwoNMethylGuanylate() : RibonucleotideResidue("2-N-methyl-guanylate", "2MG", 'X') 
00698     {
00699         bondCompound("base", TwoNMethylGuanidineGroup(), "bondBase", 0.14750, 200*Deg2Rad, BondMobility::Torsion);
00700         inheritAtomNames("base");
00701 
00702         defineDihedralAngle("chi", "O4*", "C1*", "N9", "C4");
00703         setDefaultDihedralAngle("chi", 210*Deg2Rad);
00704 
00705         addCompoundSynonym("TwoNMethylGuanidineGroup"); // for resolving biotypes
00706         
00707         setBiotypeIndex( "N9", Biotype::get("TwoNMethylGuanidineGroup", "N9").getIndex());
00708         setBiotypeIndex( "C8", Biotype::get("TwoNMethylGuanidineGroup", "C8").getIndex());
00709         setBiotypeIndex( "N7", Biotype::get("TwoNMethylGuanidineGroup", "N7").getIndex());
00710         setBiotypeIndex( "C5", Biotype::get("TwoNMethylGuanidineGroup", "C5").getIndex());
00711         setBiotypeIndex( "C6", Biotype::get("TwoNMethylGuanidineGroup", "C6").getIndex());
00712         setBiotypeIndex( "O6", Biotype::get("TwoNMethylGuanidineGroup", "O6").getIndex());
00713         setBiotypeIndex( "N1", Biotype::get("TwoNMethylGuanidineGroup", "N1").getIndex());
00714         setBiotypeIndex( "C2", Biotype::get("TwoNMethylGuanidineGroup", "C2").getIndex());
00715         setBiotypeIndex( "N2", Biotype::get("TwoNMethylGuanidineGroup", "N2").getIndex());
00716         setBiotypeIndex( "1H2", Biotype::get("TwoNMethylGuanidineGroup", "1H2").getIndex());
00717         setBiotypeIndex( "C10", Biotype::get("TwoNMethylGuanidineGroup", "C10").getIndex());
00718         setBiotypeIndex( "H20", Biotype::get("TwoNMethylGuanidineGroup", "H20").getIndex());
00719         setBiotypeIndex( "H21", Biotype::get("TwoNMethylGuanidineGroup", "H21").getIndex());
00720         setBiotypeIndex( "H22", Biotype::get("TwoNMethylGuanidineGroup", "H22").getIndex());
00721         setBiotypeIndex( "N3", Biotype::get("TwoNMethylGuanidineGroup", "N3").getIndex());
00722         setBiotypeIndex( "C4", Biotype::get("TwoNMethylGuanidineGroup", "C4").getIndex());
00723         setBiotypeIndex( "H1", Biotype::get("TwoNMethylGuanidineGroup", "H1").getIndex());
00724         setBiotypeIndex( "H8", Biotype::get("TwoNMethylGuanidineGroup", "H8").getIndex());
00725         
00726     }
00727 };
00728 
00729 
00730 
00731 
00732 
00733 
00734 class  RNA : public Biopolymer {
00735 public:
00736 
00737     explicit RNA(const PdbStructure& pdbStructure) 
00738     {
00739         // Assume first chain of first model is wanted
00740         const PdbChain& pdbChain = pdbStructure.getModel(Pdb::ModelIndex(0)).getChain(Pdb::ChainIndex(0));
00741         initializeFromPdbChain(pdbChain);
00742     }
00743 
00744     explicit RNA(const PdbChain& pdbChain)
00745     {
00746         initializeFromPdbChain(pdbChain);
00747     }
00752     explicit RNA(const Sequence& seq, bool useCappingHydroxyls = true) 
00753     {
00754         String previousResidueName;
00755         Biotype::initializePopularBiotypes();
00756 
00757         // TODO - create 5' end cap
00758         for (int resi = 0; resi < (int)seq.size(); ++resi)
00759         {
00760             
00761             RibonucleotideResidue residue = RibonucleotideResidue::create(seq[resi]);
00762             residue.assignBiotypes();
00763 
00764             residue.setPdbResidueNumber(resi + 1);
00765 
00766             // Name residue subcompound after its number in the sequence
00767             // name must be unique within the protein
00768             String residueName(resi);
00769 
00770             // Cap the 3' end with a hydroxyl group
00771             if ((resi == seq.size() - 1) && (useCappingHydroxyls)) {
00772                     residue.bondCompound(
00773                             "3PrimeHydroxyl", 
00774                             ThreePrimeRnaHydroxylGroup(residueName, residue.getPdbResidueName(), 'X'), 
00775                             "bondNext", 
00776                                     0.0960);
00777             residue.inheritAtomNames("3PrimeHydroxyl");
00778                     // Oxygen biotype must be changed for the amber atom types to come out correctly
00779             if ( (residue.hasAtom("O3*")) && ( Biotype::exists("Hydroxyl, RNA", "O3*", SimTK::Ordinality::Final)) )
00780                     residue.setBiotypeIndex( "O3*", Biotype::get("Hydroxyl, RNA", "O3*", SimTK::Ordinality::Final).getIndex() );
00781             }
00782 
00783  if (resi == 0) {
00784                     // first residue needs end cap
00785                     residue.convertInboardBondCenterToOutboard();
00786 
00787                             if (useCappingHydroxyls) { // 5' hydroxyl
00788                             residue.bondCompound(
00789                                             "5PrimeHydroxyl", 
00790                                             FivePrimeRnaHydroxylGroup(residueName, residue.getPdbResidueName(), 'X'), 
00791                                             "bondPrevious", 
00792                                                     0.0960);
00793                             residue.inheritAtomNames("5PrimeHydroxyl");
00794                             if ( (residue.hasAtom("O5*")) && ( Biotype::exists("Hydroxyl, RNA", "O5*", SimTK::Ordinality::Initial)) )
00795                                     residue.setBiotypeIndex( "O5*", Biotype::get("Hydroxyl, RNA", "O5*", SimTK::Ordinality::Initial).getIndex() );
00796                             }
00797 
00798                             else { // 5' phosphate
00799                             residue.bondCompound(
00800                                             "5PrimePhosphate", 
00801                                             FivePrimeRnaPhosphateGroup(residueName, residue.getPdbResidueName(), 'X'), 
00802                                             "bondPrevious");
00803                             residue.inheritAtomNames("5PrimePhosphate");
00804                             }
00805 
00806                     appendResidue( residueName, residue );
00807             }
00808             else {
00809                 // non-first residue needs phosphodiester linkage
00810                 appendResidue( residueName, residue.withPhosphodiester() );
00811 
00812                 // Define zeta angle
00813                 String zetaName = String("zeta") + String(resi - 1);
00814                 defineDihedralAngle(zetaName, previousResidueName + "/O3*/bond1", residueName + "/P/bond2");
00815                 setDefaultDihedralAngle(zetaName, 287.0*Deg2Rad); // (Schneider et al 2004)
00816             }
00817 
00818             previousResidueName = residueName;
00819         }
00820     }
00833     RNA& setRNABondMobility (BondMobility::Mobility  mobility, ResidueInfo::Index startResidue, ResidueInfo::Index endResidue){
00834         for (ResidueInfo::Index q =startResidue; q<=endResidue; q++) {
00835             setResidueBondMobility(q, mobility);
00836             // (updResidue(ResidueInfo::Index(q))).setCompoundBondMobility(mobility);
00837             if (q>startResidue) {
00838                 std::stringstream ss1;
00839                 ss1<<q-1<<"/O3*";
00840                 std::stringstream ss2;
00841                 ss2<<q<<"/P";
00842                 setBondMobility(mobility ,ss1.str() ,ss2.str()  ); 
00843             }
00844         }
00845         return *this;
00846     }
00856     Vec3 calcRNABaseNormal  (const State & state,int residueNumber, int bodyOrGroundFrame, const SimbodyMatterSubsystem & matter) const {
00857         const ResidueInfo& myResidue = getResidue(ResidueInfo::Index(residueNumber));
00858             String myResidueName = myResidue.getPdbResidueName();
00859         //cout<<"[RNA.h:calcRNABaseNormal] myResidueName ="<<myResidueName<<endl;
00860             //String myResidueName = (updResidue(Compound::Index(residueNumber))).getPdbResidueName();
00861         assert((myResidueName.compare("A  ") == 0) || (myResidueName.compare("G  ") == 0)  || (myResidueName.compare("U  ") == 0) || (myResidueName.compare("C  ") == 0));
00862             String atomName1;
00863             String atomName1a;
00864             String atomName1b;
00865         if ((myResidueName.compare("A  ") == 0) || (myResidueName.compare("G  ") == 0)) {
00866             atomName1  = "N9";
00867             atomName1a = "C4";
00868             atomName1b = "C8";
00869         } else {
00870             atomName1  = "N1";
00871             atomName1a = "C2";
00872             atomName1b = "C6";
00873         }
00874             Compound::AtomIndex myAtomIndex1 = myResidue.getAtomIndex(atomName1 );
00875             Compound::AtomIndex myAtomIndex1a= myResidue.getAtomIndex(atomName1a);
00876             Compound::AtomIndex myAtomIndex1b= myResidue.getAtomIndex(atomName1b);
00877         MobilizedBodyIndex myMobilizedBodyIndex1  = getAtomMobilizedBodyIndex(myAtomIndex1 );
00878         MobilizedBodyIndex myMobilizedBodyIndex1a = getAtomMobilizedBodyIndex(myAtomIndex1a);
00879         MobilizedBodyIndex myMobilizedBodyIndex1b = getAtomMobilizedBodyIndex(myAtomIndex1b);
00880             MobilizedBody body1 = matter.getMobilizedBody(myMobilizedBodyIndex1 );
00881             MobilizedBody body1a= matter.getMobilizedBody(myMobilizedBodyIndex1a);
00882             MobilizedBody body1b= matter.getMobilizedBody(myMobilizedBodyIndex1b);
00883         if (bodyOrGroundFrame == 0 ) {
00884                 Vec3 normalInBodyFrame = (
00885                     (~(body1.getBodyTransform(state)) * body1a.getBodyTransform(state) * getAtomLocationInMobilizedBodyFrame(myAtomIndex1a)  
00886                     - getAtomLocationInMobilizedBodyFrame(myAtomIndex1)) %
00887                     (~(body1.getBodyTransform(state)) * body1b.getBodyTransform(state) * getAtomLocationInMobilizedBodyFrame(myAtomIndex1b)
00888                     - getAtomLocationInMobilizedBodyFrame(myAtomIndex1)));
00889         normalInBodyFrame /= normalInBodyFrame.norm();
00890             return normalInBodyFrame;
00891         } else if (bodyOrGroundFrame == 1) {
00892                 Vec3 normalInGroundFrame = 
00893                     (    
00894                     body1a.getBodyTransform(state) * getAtomLocationInMobilizedBodyFrame(myAtomIndex1a)
00895                     - body1.getBodyTransform(state) * getAtomLocationInMobilizedBodyFrame(myAtomIndex1)
00896                     ) % (
00897                     body1b.getBodyTransform(state) * getAtomLocationInMobilizedBodyFrame(myAtomIndex1b)
00898                     - body1.getBodyTransform(state) * getAtomLocationInMobilizedBodyFrame(myAtomIndex1)
00899                     );
00900                 normalInGroundFrame /= normalInGroundFrame.norm();
00901                 return normalInGroundFrame;
00902         }
00903 
00904     }
00905 
00906 private:
00907     void initializeFromPdbChain(
00908             const PdbChain& pdbChain, 
00909             Compound::MatchStratagem matchStratagem = 
00910                     Compound::Match_TopologyOnly)
00911     {
00912         setPdbChainId( pdbChain.getChainId() );
00913         
00914         String previousResidueName("");
00915         
00916         // TODO end caps
00917         for (Pdb::ResidueIndex r(0); r < (Pdb::ResidueIndex)pdbChain.getNumResidues(); ++r) 
00918         {
00919             const PdbResidue& pdbResidue = pdbChain.getResidue(r);
00920             String residueName = pdbResidue.getName() + String(pdbResidue.getPdbResidueNumber());
00921             if (pdbResidue.getInsertionCode() != ' ') // include unusual insertion codes in residue name
00922                     residueName += pdbResidue.getInsertionCode();
00923             
00924             RibonucleotideResidue residue = RibonucleotideResidue::create(pdbResidue, pdbChain.getChainId());
00925             residue.assignBiotypes();
00926             
00927             if (r == 0) { // first residue does not get phosphate
00928                 appendResidue( residueName, residue );
00929             }
00930             else {
00931                 // non-first residue needs phosphodiester linkage
00932                 BiopolymerResidue resWithP = residue.withPhosphodiester();
00933                 resWithP.assignBiotypes();
00934                 appendResidue( residueName, resWithP );
00935 
00936                 // Define zeta angle
00937                 String zetaName = String("zeta") + String(residue.getPdbResidueNumber());
00938                 defineDihedralAngle(zetaName, previousResidueName + "/O3*/bond1", residueName + "/P/bond2");
00939                 setDefaultDihedralAngle(zetaName, 287.0*Deg2Rad); // (Schneider et al 2004)
00940             } 
00941             
00942             previousResidueName = residueName;            
00943         }
00944         
00945         if (matchStratagem != Compound::Match_TopologyOnly) 
00946         {
00947             Compound::AtomTargetLocations atomTargets = 
00948                     createAtomTargets(pdbChain);
00949             
00950             matchDefaultConfiguration(atomTargets, matchStratagem);
00951         }        
00952     }
00953 };
00954 
00955 
00956 } // namespace SimTK
00957 
00958 #endif // SimTK_MOLMODEL_RNA_H_
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines