Compound.h

Go to the documentation of this file.

#ifndef SimTK_MOLMODEL_COMPOUND_H_
#define SimTK_MOLMODEL_COMPOUND_H_

/* -------------------------------------------------------------------------- *
 *                      SimTK Core: SimTK Molmodel                            *
 * -------------------------------------------------------------------------- *
 * This is part of the SimTK Core biosimulation toolkit originating from      *
 * Simbios, the NIH National Center for Physics-Based Simulation of           *
 * Biological Structures at Stanford, funded under the NIH Roadmap for        *
 * Medical Research, grant U54 GM072970. See https://simtk.org.               *
 *                                                                            *
 * Portions copyright (c) 2007-2008 Stanford University and the Authors.      *
 * Authors: Michael Sherman, Christopher Bruns                                *
 * Contributors:                                                              *
 *                                                                            *
 * Permission is hereby granted, free of charge, to any person obtaining a    *
 * copy of this software and associated documentation files (the "Software"), *
 * to deal in the Software without restriction, including without limitation  *
 * the rights to use, copy, modify, merge, publish, distribute, sublicense,   *
 * and/or sell copies of the Software, and to permit persons to whom the      *
 * Software is furnished to do so, subject to the following conditions:       *
 *                                                                            *
 * The above copyright notice and this permission notice shall be included in *
 * all copies or substantial portions of the Software.                        *
 *                                                                            *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,   *
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL    *
 * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,    *
 * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR      *
 * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE  *
 * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
 * -------------------------------------------------------------------------- */


#include "SimTKsimbody.h"

#include "molmodel/internal/common.h"
#include "DuMMForceFieldSubsystem.h"
#include "molmodel/internal/Pdb.h"
#include "molmodel/internal/Superpose.h"
#include "molmodel/internal/units.h"
#include <map>

#include <iosfwd> // declare ostream without all the definitions

namespace SimTK {
class CompoundSystem;

// Because of private data, put Compound implementation behind the curtain
class Compound;
class CompoundRep;

// Make sure class is instantiated just once in the library
// #ifndef DO_INSTANTIATE_COMPOUND_PIMPL_HANDLE
//    template class SimTK_MOLMODEL_EXPORT PIMPLHandle<Compound, CompoundRep>;
// #endif

namespace BondMobility {
    enum Mobility {
        Free = 1, 
        Torsion = 2, 
        Rigid = 3 
    };
    static Mobility Default = Torsion;
}

class SimTK_MOLMODEL_EXPORT Compound : public PIMPLHandle<Compound,CompoundRep> {
public:

    enum MatchStratagem {
        Match_Exact, //< Try to match input atom positions precisely
        Match_Idealized, //< Try to match atom positions with idealized geometry
        Match_TopologyOnly //< Don't match atom positions, only topology
    };
    

    typedef String Name;

    typedef String AtomName;

    typedef String BondCenterName;

    typedef String DihedralName;

    typedef String BondCenterPathName;

    typedef String AtomPathName;

    // Unique classes which behave like type-safe ints. These are local
    // to Compound.

    // SimTK_DEFINE_UNIQUE_LOCAL_INDEX_TYPE(Compound,Index);


    SimTK_DEFINE_UNIQUE_LOCAL_INDEX_TYPE(Compound,AtomIndex);

    SimTK_DEFINE_UNIQUE_LOCAL_INDEX_TYPE(Compound,LocalAtomIndex);

    SimTK_DEFINE_UNIQUE_LOCAL_INDEX_TYPE(Compound,BondCenterIndex);

    SimTK_DEFINE_UNIQUE_LOCAL_INDEX_TYPE(Compound,BondIndex);

    typedef std::map<AtomIndex, Vec3> AtomTargetLocations;

    class SingleAtom;




    Compound();

    explicit Compound(const Name& name 
        );

    int getNumAtoms() const;

    size_t getNumBondCenters() const;

    size_t getNumBondCenters(Compound::AtomIndex atomIndex) const;

    size_t getNumBonds() const;

    // int getNumSubcompounds() const;

    bool hasAtom(const AtomPathName& name 
        ) const;

    bool hasBondCenter(const BondCenterPathName& bondCenter 
        ) const;

    // bool hasSubcompound(const Name& name ///< name of the subcompound relative to this compound.
    //     ) const;

    Compound::AtomIndex getBondAtomIndex(
        Compound::BondIndex bondIndex, 
        int which                      
        ) const;


    Compound& setBaseAtom(
        const Compound::AtomName& name,   
        const Element& element,           
        const Transform& location = Vec3(0)  
        );
    Compound& setBaseAtom(
        const Compound::AtomName& name,   
        const Biotype& biotype,           
        const Transform& location = Vec3(0)    
        );  

    Compound& setBaseAtom(
        const Compound::SingleAtom& c, 
        const Transform& = Transform() 
        ); 

    Compound& setBaseCompound(
        const Compound::Name& n,         
        const Compound& c,               
        const Transform& location = Transform()
        );

    Compound& bondAtom(
        const Compound::SingleAtom& c, 
        const BondCenterPathName& parentBondName, 
        mdunits::Length distance, 
        Angle dihedral = 0, 
        BondMobility::Mobility = BondMobility::Default 
        );

    Compound& bondAtom(
        const Compound::SingleAtom& c, 
        const BondCenterPathName& parentBondName 
        );

    Compound& bondCompound(
        const Compound::Name& n, 
        const Compound& c, 
        const BondCenterPathName& parentBondName, 
        mdunits::Length distance, 
        Angle dihedral = 180*Deg2Rad, 
        BondMobility::Mobility mobility = BondMobility::Default 
        );

    Compound& bondCompound(
        const Compound::Name& n, 
        const Compound& c,  
        const BondCenterPathName& parentBondName 
        );

    Compound& bondCompound(
        const Compound::Name& n, 
        const Compound& c, 
        const BondCenterPathName& parentBondName,  
        BondMobility::Mobility mobility 
        );

    Compound& setInboardBondCenter(
        const Compound::BondCenterName& centerName 
        );

    Compound& convertInboardBondCenterToOutboard();

    Compound& addFirstBondCenter(
        const Compound::BondCenterName& centerName, 
        const Compound::AtomPathName& atomName 
        );

    Compound& addSecondBondCenter(
        const Compound::BondCenterName& centerName, 
        const Compound::AtomName& atomName, 
        Angle bondAngle1 
        );

    Compound& addPlanarBondCenter(
        const Compound::BondCenterName& centerName, 
        const Compound::AtomName& atomName, 
        Angle bondAngle1, 
        Angle bondAngle2 
        );

    Compound& addRightHandedBondCenter(
        const Compound::BondCenterName& centerName, 
        const Compound::AtomName& atomName, 
        Angle bondAngle1, 
        Angle bondAngle2 
        );

    Compound& addLeftHandedBondCenter(
        const Compound::BondCenterName& centerName, 
        const Compound::AtomName& atomName, 
        Angle bondAngle1, 
        Angle bondAngle2 
        );

    Compound& addRingClosingBond(
        const Compound::BondCenterPathName& centerName1, 
        const Compound::BondCenterPathName& centerName2, 
        mdunits::Length bondLength, 
        Angle dihedral = 180*Deg2Rad, 
        BondMobility::Mobility mobility = BondMobility::Default 
        );

    Compound& addRingClosingBond(
        const Compound::BondCenterPathName& centerName1, 
        const Compound::BondCenterPathName& centerName2 
        );

    //const Compound& getSubcompound(Compound::Index ///< integer index of subcompound
    //    ) const;

    //Compound& updSubcompound(Compound::Index ///< integer index of subcompound
    //    );

    const Compound& getSubcompound(const Compound::Name& subcompoundName 
        ) const;

    Compound& updSubcompound(const Compound::Name& subcompoundName 
        );

    // end topology methods



    Vec3 calcDefaultAtomLocationInGroundFrame(const AtomPathName& atomName 
        ) const;

    Vec3 calcDefaultAtomLocationInCompoundFrame(const AtomPathName& atomName 
        ) const;

    Compound& setDefaultInboardBondLength(mdunits::Length 
        );

    Compound& setDefaultInboardDihedralAngle(Angle 
        );

    Compound& setDefaultBondAngle(
        Angle angle, 
        const AtomPathName& atom1, 
        const AtomPathName& atom2, 
        const AtomPathName& atom3 
        );

    Compound& setDefaultBondLength(
        mdunits::Length length, 
        const AtomPathName& atom1, 
        const AtomPathName& atom2 
        );

    Compound& setDefaultDihedralAngle(
        const DihedralName& dihedralName, 
        Angle angleInRadians 
        );

    Compound& setDefaultDihedralAngle(
                Angle angle, 
                Compound::AtomIndex atom1,
                Compound::AtomIndex atom2,
                Compound::AtomIndex atom3,
                Compound::AtomIndex atom4
                );
    
    Compound& setDefaultDihedralAngle(
                Angle angle, 
                const Compound::AtomName& atom1,
                const Compound::AtomName& atom2,
                const Compound::AtomName& atom3,
                const Compound::AtomName& atom4
                );
    
    Angle calcDefaultDihedralAngle(
        const DihedralName& dihedralName 
        ) const;

    Compound& setDihedralAngle(
        State& state, 
        const DihedralName& dihedralName, 
        Angle 
        );

    Angle calcDihedralAngle(
        const State& state, 
        const DihedralName& dihedralName 
        ) const;

    Transform calcDefaultAtomFrameInCompoundFrame(Compound::AtomIndex 
        ) const;

    Vec3 calcAtomLocationInGroundFrame(
        const State& state, 
        Compound::AtomIndex atomId 
        ) const;

    Vec3 calcAtomLocationInCompoundFrame(
        const State& state, 
        Compound::AtomIndex atomId 
        ) const;

    Vec3 calcAtomVelocityInGroundFrame(
        const State& state, 
        Compound::AtomIndex atomId  
        ) const;

    Vec3 calcAtomAccelerationInGroundFrame(
        const State& state, 
        Compound::AtomIndex atomId   
        ) const;

    Transform getSubcompoundFrameInParentFrame(
        const Compound::Name& subcompoundName 
        ) const;

    virtual AtomTargetLocations createAtomTargets(const class PdbStructure& targetStructure) const;
    virtual AtomTargetLocations createAtomTargets(const class PdbChain& targetChain) const;
    virtual AtomTargetLocations createAtomTargets(const class PdbResidue& targetResidue) const;

    Compound& matchDefaultAtomChirality(
            const AtomTargetLocations& atomTargets, 
            Angle planarityTolerance = 90.0 * Deg2Rad, 
            bool flipAll = true 
            );
    Compound& matchDefaultBondLengths(const AtomTargetLocations& atomTargets);
    Compound& matchDefaultBondAngles(const AtomTargetLocations& atomTargets);

    enum PlanarBondMatchingPolicy {KeepPlanarBonds, DistortPlanarBonds, FlipPlanarBonds};
    Compound& matchDefaultDihedralAngles(
        const AtomTargetLocations& atomTargets,
        PlanarBondMatchingPolicy policy = FlipPlanarBonds 
        );

    Compound& matchDefaultTopLevelTransform(const AtomTargetLocations& atomTargets);
    TransformAndResidual getTransformAndResidual(const Compound::AtomTargetLocations& atomTargets) const;
    
    Compound& matchDefaultConfiguration(
            const AtomTargetLocations& atomTargets, 
            MatchStratagem matchStratagem = Match_Exact) 
    {
        if (matchStratagem == Compound::Match_TopologyOnly)
            // do nothing, topology is already there
            return *this;
        
        else if (matchStratagem == Compound::Match_Exact) 
        {
            // low tolerance breaks planarity just about everywhere
            matchDefaultAtomChirality(atomTargets, 0.01, false);

            matchDefaultBondLengths(atomTargets);
            matchDefaultBondAngles(atomTargets);
            
            // Set dihedral angles even when bonded atoms are planar
            matchDefaultDihedralAngles(atomTargets, Compound::DistortPlanarBonds);

            matchDefaultTopLevelTransform(atomTargets);
            
            // No further optimization should be needed
            
            return *this;
        }
        
        else if (matchStratagem == Compound::Match_Idealized) 
        {
            // Break planarity constraint only when input atoms are 90 degrees out of plane
            matchDefaultAtomChirality(atomTargets, 90*Deg2Rad, true);
            
            // Avoid setting non zero/180 dihedral angles on bonded planar atoms
            matchDefaultDihedralAngles(atomTargets, Compound::FlipPlanarBonds);
            
            matchDefaultTopLevelTransform(atomTargets);

            // At this point the structure match is approximate and may have
            // many bad contacts.  Optimization is needed
            // TODO - ObservedPointFitter on internal coordinates
            fitDefaultConfiguration(atomTargets, 0.005);
            
            return *this;
        }
            
        else assert(false); // unknown match stratagem

        return *this;
    }

    Compound& fitDefaultConfiguration(
            const AtomTargetLocations& atomTargets,
            SimTK::Real targetRms);
    
    const Compound& populateDefaultPdbChain(
        class PdbChain&, 
        int& defaultNextResidueNumber,
        const Transform& transform = Transform()) const;

    const Compound& populatePdbChain(
        const State& state, 
        class PdbChain&, 
        int& defaultNextResidueNumber,
        const Transform& transform = Transform()) const;

    std::ostream& writeDefaultPdb(
        std::ostream& os, 
        const Transform& transform = Transform()  
        ) const;

    void writeDefaultPdb(
        const char* outFileName, 
        const Transform& transform
        ) const;

    std::ostream& writeDefaultPdb(
        std::ostream& os, 
        int& nextAtomSerialNumber, 
        const Transform& transform = Transform() 
        ) const;

    std::ostream& writePdb(
        const State& state, 
        std::ostream& os,  
        const Transform& transform = Transform() 
        ) const;


    std::ostream& writePdb(
        const State& state, 
        std::ostream& os,  
        int& nextAtomSerialNumber,  
        const Transform& transform = Transform() 
        ) const;

    // end configuration section



    Compound& setCompoundName(const Name& 
        );


    const Name& getCompoundName() const;

    Compound& addCompoundSynonym(const Name& 
        );

    const AtomPathName getAtomName(Compound::AtomIndex 
        ) const;

    const Element& getAtomElement(Compound::AtomIndex 
        ) const;

    const Element& getAtomElement(const Compound::AtomName& 
        ) const;

    Compound& nameAtom(
        const Compound::AtomName& newName, 
        const AtomPathName& oldName 
        );

    Compound& nameAtom(
        const Compound::AtomName& newName,  
        const AtomPathName& oldName,  
        BiotypeIndex biotype 
        );

    Compound& defineDihedralAngle(
        const Compound::DihedralName& angleName, 
        const Compound::AtomPathName& atom1, 
        const Compound::AtomPathName& atom2, 
        const Compound::AtomPathName& atom3, 
        const Compound::AtomPathName& atom4, 
        Angle offset = 0*Deg2Rad 
        );

    Compound& defineDihedralAngle(
        const Compound::DihedralName& angleName, 
        const Compound::BondCenterPathName& bond1, 
        const Compound::BondCenterPathName& bond2, 
        Angle offset = 0*Deg2Rad 
        );

    Compound& setPdbResidueNumber(int resNum);

    int getPdbResidueNumber() const;

    Compound& setPdbResidueName(const String& resName);

    const String& getPdbResidueName() const;

    Compound& setPdbChainId(char chainId 
        );

    char getPdbChainId() const;

    Compound& nameBondCenter(
        const Compound::BondCenterName& newName, 
        const BondCenterPathName& oldName 
        );

    Compound& inheritAtomNames(const Compound::Name& 
        );
    
    Compound& inheritCompoundSynonyms(const Compound& otherCompound);
    
    Compound& inheritBondCenterNames(const Compound::Name& 
        );

    Compound::AtomIndex getAtomIndex(const Compound::AtomPathName& 
        ) const;

    // end nomenclature methods



    //void setCompoundSystem(CompoundSystem& system, ///< The CompoundSystem that will take ownership of this Compound.
    //                       Compound::Index);       ///< The index within the new owner CompoundSystem that this Compound will have.

    void setMultibodySystem(MultibodySystem& system);

    Compound& setBondMobility(
        BondMobility::Mobility mobility, 
        const AtomPathName& atom1, 
        const AtomPathName& atom2 
        );

    Compound& setBondMobility(
        BondMobility::Mobility mobility,  
        Compound::BondIndex bondIndex 
        );

    Compound  & setCompoundBondMobility(BondMobility::Mobility mobility);

    MobilizedBodyIndex getAtomMobilizedBodyIndex(Compound::AtomIndex 
        ) const;

    Vec3 getAtomLocationInMobilizedBodyFrame(Compound::AtomIndex 
        ) const;

    Compound& setBiotypeIndex(
        const Compound::AtomPathName& atomName, 
        BiotypeIndex biotype 
        );
    
    Compound& setAtomBiotype(
            const Compound::AtomPathName& atomName,
            const String& biotypeResidueName,
            const String& biotypeAtomName,
            SimTK::Ordinality::Residue ordinality = SimTK::Ordinality::Any
            )
    {
        Compound::AtomIndex atomIndex = getAtomIndex(atomName);
        
        if ( ! Biotype::exists(biotypeResidueName, biotypeAtomName, ordinality) )
            Biotype::defineBiotype(
                getAtomElement(atomIndex),
                getNumBondCenters(atomIndex), 
                biotypeResidueName, 
                biotypeAtomName, 
                ordinality
            );
        
        const Biotype& biotype = Biotype::get(biotypeResidueName, biotypeAtomName, ordinality);
        
        assert( biotype.getElement() == getAtomElement(atomIndex) );
        assert( biotype.getValence() == getNumBondCenters(atomIndex) );
        
        BiotypeIndex biotypeIndex = biotype.getIndex();
        setBiotypeIndex(atomName, biotypeIndex);

        return *this;
    }

    BiotypeIndex getAtomBiotypeIndex(Compound::AtomIndex 
        ) const;

    // RUNTIME INTERFACE

    DuMM::AtomIndex getDuMMAtomIndex(Compound::AtomIndex) const;

    // end simulation methods

    Compound& setTopLevelTransform(const Transform& transform);

    const Transform& getTopLevelTransform() const;



protected:

    void setDuMMAtomIndex(
        Compound::AtomIndex, 
        DuMM::AtomIndex 
        );

    explicit Compound(CompoundRep* ip);
    friend class CompoundSystem;

private:
    // OBSOLETE: use getNumAtoms()
    int getNAtoms() const {return getNumAtoms();}
    // OBSOLETE: use getNumBondCenters()
    size_t getNBondCenters() const {return getNumBondCenters();}
    // OBSOLETE: use getNumBondCenters()
    size_t getNBondCenters(Compound::AtomIndex atomIndex) const {return getNumBondCenters(atomIndex);}
    // OBSOLETE: use getNumBonds()
    size_t getNBonds() const {return getNumBonds();}
};


SimTK_MOLMODEL_EXPORT std::ostream& operator<<(std::ostream& o, const Compound& c);

class  Compound::SingleAtom : public Compound {
public:
    SingleAtom(
        const Compound::AtomName& atomName, 
        const Element& element 
        ) 
    {
        setBaseAtom( atomName, element );

        setCompoundName("SingleAtom"); // should be overridden by derived class constructors
    }
};


class  UnivalentAtom : public Compound::SingleAtom {
public:
    UnivalentAtom(
        const Compound::AtomName& atomName, 
        const Element& element 
        ) 
        : Compound::SingleAtom(atomName, element)
    {
        addFirstBondCenter( "bond", atomName);
        setInboardBondCenter("bond");

        setCompoundName("UnivalentAtom"); // should be overridden by derived class constructors
    }
};


class  BivalentAtom : public Compound::SingleAtom {
public:
    BivalentAtom(
        const Compound::AtomName& atomName,  
        const Element& element,  
        Angle angle = 180*Deg2Rad 
        ) 
        : Compound::SingleAtom(atomName, element)
    {        
        // BondCenter1 dihedral will be relative to BondCenter2
        addFirstBondCenter( "bond1", atomName);
        // conversely, bond center 2 dihedral is relative to bond center 1
        addSecondBondCenter( "bond2", atomName, angle);
        setInboardBondCenter("bond1"); // without loss of generality

        setCompoundName("BivalentAtom"); // should be overridden by derived class constructors
    }
};


class  TrivalentAtom : public Compound::SingleAtom {
public:
    TrivalentAtom(
        const Compound::AtomName& atomName,   
        const Element& element,   
        Angle angle1 = 120*Deg2Rad, 
        Angle angle2 = 120*Deg2Rad 
        ) 
        : Compound::SingleAtom(atomName, element)
    {
        // BondCenter1 dihedral will be relative to BondCenter2
        addFirstBondCenter( "bond1", atomName );
        // bond centers 2 and 3 dihedrals relative to bond center 1
        addSecondBondCenter( "bond2", atomName,  angle1);
        addPlanarBondCenter( "bond3", atomName, angle2, 360*Deg2Rad - angle1 - angle2);

        // Choice of inboard bond may differ from bond priority - user may change this
        setInboardBondCenter("bond1");

        setCompoundName("TrivalentAtom"); // should be overridden by derived class constructors
    }
};

class  QuadrivalentAtom : public Compound::SingleAtom {
public:
    QuadrivalentAtom(
        const Compound::AtomName& atomName, 
        const Element& element 
        ) 
        : Compound::SingleAtom(atomName, element)
    {
        static const Angle TetrahedralAngle = 109.47 * Deg2Rad;

        // BondCenter1 dihedral will be relative to BondCenter2
        // Using Rotation() constructor that takes x axis and approximate y axis
        addFirstBondCenter( "bond1", atomName );
        // bond centers 2, 3, and 4 dihedrals will be relative to bond1
        addSecondBondCenter( "bond2", atomName, TetrahedralAngle );
        addLeftHandedBondCenter( "bond3", atomName, TetrahedralAngle, TetrahedralAngle );
        addRightHandedBondCenter( "bond4", atomName, TetrahedralAngle, TetrahedralAngle );

        // Choice of inboard bond may differ from bond priority - user may change this
        setInboardBondCenter("bond1");

        // default dihedral should be left blank, so it could be
        // defined from the other side.
        // the ultimate default will be 180 degrees anyway.
        // setDefaultInboardDihedralAngle(180*Deg2Rad);

        setCompoundName("QuadrivalentAtom"); // should be overridden by derived class constructors
    }
        
    QuadrivalentAtom(
        const Compound::AtomName& atomName, 
        const Element& element, 
        Angle bond12Angle, 
        Angle bond13Angle, 
        Angle bond14Angle, 
        Angle dihedral3, 
        Angle dihedral4 
        ) 
        : Compound::SingleAtom(atomName, element)
    {
        // BondCenter1 dihedral will be relative to BondCenter2
        // Using Rotation() constructor that takes x axis and approximate y axis
        addFirstBondCenter( "bond1", atomName );
        // bond centers 2, 3, and 4 dihedrals will be relative to bond1
        addSecondBondCenter( "bond2", atomName, bond12Angle );
        
        // Compute bond23Angle and bond24Angle
        SimTK::Vec3 v1(0,0,1);
        SimTK::Vec3 v2 = SimTK::Rotation(bond12Angle, ZAxis) * SimTK::Vec3(0,0,1);
        SimTK::Vec3 v3 = SimTK::Rotation(dihedral3, YAxis) * SimTK::Rotation(bond12Angle, ZAxis) * SimTK::Vec3(0,0,1);
        SimTK::Vec3 v4 = SimTK::Rotation(dihedral4, YAxis) * SimTK::Rotation(bond12Angle, ZAxis) * SimTK::Vec3(0,0,1);
        
        Angle bond23Angle = std::acos(SimTK::dot(v2, v3));
        Angle bond24Angle = std::acos(SimTK::dot(v2, v4));
        
        // restrict angle range to (-Pi,Pi)
        while (dihedral3 >   Pi) dihedral3 -= 2.0*Pi;
        while (dihedral3 <= -Pi) dihedral3 += 2.0*Pi;
        while (dihedral4 >   Pi) dihedral4 -= 2.0*Pi;
        while (dihedral4 <= -Pi) dihedral4 += 2.0*Pi;
        
        if (dihedral3 > 0)
            addLeftHandedBondCenter( "bond3", atomName, bond13Angle, bond23Angle );
        else
            addRightHandedBondCenter( "bond3", atomName, bond13Angle, bond23Angle );
        
        if (dihedral4 > 0)
            addLeftHandedBondCenter( "bond4", atomName, bond14Angle, bond24Angle );
        else 
            addRightHandedBondCenter( "bond4", atomName, bond14Angle, bond24Angle );

        // Choice of inboard bond may differ from bond priority - user may change this
        setInboardBondCenter("bond1");

        // default dihedral should be left blank, so it could be
        // defined from the other side.
        // the ultimate default will be 180 degrees anyway.
        // setDefaultInboardDihedralAngle(180*Deg2Rad);

        setCompoundName("QuadrivalentAtom"); // should be overridden by derived class constructors
    }
};


class  AliphaticHydrogen : public UnivalentAtom {
public:
    explicit AliphaticHydrogen(const AtomName& atomName = "H" 
        ) 
        : UnivalentAtom(atomName, Element::Hydrogen())
    {
        setDefaultInboardBondLength(0.1112); // for bonding to aliphatic carbon

        setCompoundName("AliphaticHydrogenAtom");
    }
};


class  AliphaticCarbon : public QuadrivalentAtom {
public:
    explicit AliphaticCarbon(const AtomName& atomName = "C" 
        ) 
        : QuadrivalentAtom(atomName, Element::Carbon()) 
    {
        // In case this bonds to another aliphatic carbon
        setDefaultInboardBondLength(0.15620); // for bonding to another aliphatic carbon

        setCompoundName("AliphaticCarbon");
    }
};


class  MethyleneGroup : public Compound {
public:
    MethyleneGroup() {
        static const mdunits::Length C_Hdistance = 0.1112; // nanometers

        setBaseAtom(AliphaticCarbon("C"));
        bondAtom(AliphaticHydrogen("H1"), "C/bond3");
        bondAtom(AliphaticHydrogen("H2"), "C/bond4");
        nameBondCenter("bond1", "C/bond1");
        nameBondCenter("bond2", "C/bond2");

        setDefaultInboardBondLength(0.15620); // for bonding to another aliphatic carbon

        setCompoundName("MethyleneGroup");
    }
};


class  MethylGroup : public Compound {
public:
    MethylGroup() {
        setBaseAtom(AliphaticCarbon("C"));
        bondAtom(AliphaticHydrogen("H1"), "C/bond2");
        bondAtom(AliphaticHydrogen("H2"), "C/bond3");
        bondAtom(AliphaticHydrogen("H3"), "C/bond4");
        nameBondCenter("bond", "C/bond1");

        setDefaultInboardBondLength(0.15620); // for bonding to another aliphatic carbon
        setCompoundName("MethylGroup");
    }
};

class  AromaticSixMemberedCHGroup : public Compound {
public:
    AromaticSixMemberedCHGroup() {
        static const mdunits::Length C_Hdistance = 0.1080; // in nanometers, from tinker amber99.dat
        static const mdunits::Length C_Cdistance = 0.1400; // in nanometers, from tinker amber99.dat

        setBaseAtom( TrivalentAtom("C", Element::Carbon(), 120*Deg2Rad, 120*Deg2Rad) );
        bondAtom( UnivalentAtom("H", Element::Hydrogen()), "C/bond3", C_Hdistance);

        // remember to change default length for bonding other other things
        setDefaultInboardBondLength(C_Cdistance); // for bonding to other aromatic CH

        setCompoundName("AromaticSixMemberedCHGroup");
    }
};


class  AlcoholOHGroup : public Compound {
public:
    AlcoholOHGroup() {
        static const mdunits::Length O_Hdistance = 0.0960; // nanometers

        setBaseAtom( BivalentAtom("O", Element::Oxygen(), 108.50 * Deg2Rad) );

        bondAtom( UnivalentAtom("H", Element::Hydrogen()), "O/bond2", O_Hdistance );
        nameBondCenter("bond", "O/bond1");

        // In case of bonding to aliphatic carbon
        setDefaultInboardBondLength(0.1410); // for bonding to aliphatic carbon

        setCompoundName("AlcoholOHGroup");
    }
};


class  PrimaryAmineGroup : public Compound {
public:
    PrimaryAmineGroup() {
        static const mdunits::Length H_Ndistance = 0.1010; // nanometers
        static const mdunits::Length N_Cdistance = 0.1471; // in nanometers, for bonding to aliphatic carbon

        setBaseAtom( QuadrivalentAtom("N", Element::Nitrogen()) );
        bondAtom( UnivalentAtom("H1", Element::Hydrogen()), "N/bond2", H_Ndistance);
        bondAtom( UnivalentAtom("H2", Element::Hydrogen()), "N/bond3", H_Ndistance);
        bondAtom( UnivalentAtom("H3", Element::Hydrogen()), "N/bond4", H_Ndistance);

        setDefaultInboardBondLength(N_Cdistance); // for bonding to aliphatic carbon
    }
};


class  CarboxylateGroup : public Compound {
public:
    CarboxylateGroup() {
        // parameters from Tinker amber99.param
        static const mdunits::Length O_Cdistance = 0.1250; // nanometers
        static const mdunits::Length C_CtDistance = 0.15220; // nanometers
        static const Angle O_C_Oangle = 126*Deg2Rad;

        // we actually use the inboard-C-O angle during construction
        static const Angle O_C_CtAngle = 180*Deg2Rad - 0.5*O_C_Oangle;

        setBaseAtom( TrivalentAtom("C", Element::Carbon(), O_C_CtAngle, O_C_CtAngle) );
        bondAtom( UnivalentAtom("O1", Element::Oxygen()), "C/bond2", O_Cdistance);
        bondAtom( UnivalentAtom("O2", Element::Oxygen()), "C/bond3", O_Cdistance);

        setDefaultInboardBondLength(C_CtDistance);

        nameBondCenter("bond", "C/bond1");
    }
};

class  SimTK_MOLMODEL_EXPORT Molecule : public Compound {
public:
    Molecule();
    SimTK_INSERT_DERIVED_HANDLE_DECLARATIONS(Molecule, CompoundRep, Compound);
protected:
    explicit Molecule(CompoundRep* rep);
};


class  Argon : public Molecule {
public:
    Argon() {
        setPdbResidueName("AR ");

        setBaseAtom( "Ar", Biotype::Argon() );

        setCompoundName("Argon");
    }
};


class  Methane : public Molecule {
public:
    Methane() 
    {
        setBaseCompound("methyl", MethylGroup());
        inheritAtomNames("methyl");

        // Ordinarily, methyl group bonds to aliphatic carbon,
        // and has a default bond length to match.
        // Here we turn off the methyl inboard bond, so the
        // AliphaticHydrogen will, as desired, dictate the bond length
        convertInboardBondCenterToOutboard();

        bondAtom(AliphaticHydrogen("H4"), "methyl/bond", 0.1112);

        setBiotypeIndex( "C", Biotype::MethaneC().getIndex() );
        setBiotypeIndex( "H1", Biotype::MethaneH().getIndex() );
        setBiotypeIndex( "H2", Biotype::MethaneH().getIndex() );
        setBiotypeIndex( "H3", Biotype::MethaneH().getIndex() );
        setBiotypeIndex( "H4", Biotype::MethaneH().getIndex() );

        setCompoundName("Methane");


    }
};

class  Ethane : public Molecule {
public:
    Ethane() 
    {
        setPdbResidueName("EHN");

        setBaseCompound("methyl1", MethylGroup());
        nameAtom("C1", "methyl1/C", Biotype::EthaneC().getIndex() );
        nameAtom("H1", "methyl1/H1", Biotype::EthaneH().getIndex() );
        nameAtom("H2", "methyl1/H2", Biotype::EthaneH().getIndex() );
        nameAtom("H3", "methyl1/H3", Biotype::EthaneH().getIndex() );

        // This first methyl is a base, not a decoration
        convertInboardBondCenterToOutboard();

        bondCompound("methyl2", MethylGroup(), "methyl1/bond");

        nameAtom("C2", "methyl2/C", Biotype::EthaneC().getIndex() );
        nameAtom("H4", "methyl2/H1", Biotype::EthaneH().getIndex() );
        nameAtom("H5", "methyl2/H2", Biotype::EthaneH().getIndex() );
        nameAtom("H6", "methyl2/H3", Biotype::EthaneH().getIndex() );

        defineDihedralAngle("torsion", "H1", "C1", "C2", "H4");

        setDefaultTorsionAngle(180*Deg2Rad);
        // setBondRotatable(false, "C1", "C2");

        setCompoundName("Ethane");

    }

    Ethane& setDefaultTorsionAngle(Angle angle 
        ) {
        setDefaultDihedralAngle("torsion", angle);

        return *this;
    }

    Angle calcDefaultTorsionAngle() const {
        return calcDefaultDihedralAngle("torsion");
    }
};

class BiopolymerResidueRep;
// class SimTK_MOLMODEL_EXPORT BiopolymerResidue : public Compound {
class SimTK_MOLMODEL_EXPORT BiopolymerResidue : public Compound {
public:   
    BiopolymerResidue(
        const Compound::Name& residueTypeName, 
        const String& threeLetterCode, 
        char oneLetterCode 
        );

    BiopolymerResidue& setOneLetterCode(char olc 
        );

    BiopolymerResidue& setThreeLetterCode(const String& tlc 
        );

    BiopolymerResidue& setResidueTypeName(const String& name 
        );
    
    char getOneLetterCode() const;
    const String& getThreeLetterCode() const;
    const String& getResidueTypeName() const;

    bool assignBiotypes(
        Ordinality::Residue ordinality = Ordinality::Any 
        );

    Compound::AtomIndex getParentCompoundAtomIndex(AtomIndex residueAtomIndex) const;

    SimTK_INSERT_DERIVED_HANDLE_DECLARATIONS(BiopolymerResidue,BiopolymerResidueRep,Compound);
};

class BiopolymerRep;

// This class to point to residue atoms in a Biopolymer
class ResidueInfo {
public:
    SimTK_DEFINE_UNIQUE_LOCAL_INDEX_TYPE(ResidueInfo,Index);
    SimTK_DEFINE_UNIQUE_LOCAL_INDEX_TYPE(ResidueInfo,AtomIndex);

    class AtomInfo {
    public:
        friend class ResidueInfo;

        AtomInfo(Compound::AtomIndex index, const Compound::AtomName& name) 
            : pdbAtomName(name), biopolymerAtomIndex(index)
        {
            synonyms.insert(name);
        }

        const std::set<Compound::AtomName>& getNames() const {
            return synonyms;
        }

    private:
        Compound::AtomIndex biopolymerAtomIndex;
        String pdbAtomName;
        std::set<Compound::AtomName> synonyms;
    };


    ResidueInfo(
        ResidueInfo::Index ix, 
        const Compound::Name& name, 
        const BiopolymerResidue& res,
        Compound::AtomIndex atomOffset,
        char insertionCode = ' ');

    AtomIndex addAtom(Compound::AtomIndex index, const Compound::AtomName& name) {
        ResidueInfo::AtomIndex answer(atoms.size());
        atoms.push_back(AtomInfo(index, name));
        atomIdsByName[name] = answer;
        return answer;
    }

    size_t getNumAtoms() const {
        return atoms.size();
    }

    const AtomInfo& getAtomInfo(AtomIndex a) const {
        return atoms[a];
    }
    AtomInfo& updAtomInfo(AtomIndex a) {
        return atoms[a];
    }

    const std::set<Compound::Name>& getNames() const {return synonyms;}

    const Compound::Name& getPdbResidueName() const {
        return pdbResidueName;
    }

    const int getPdbResidueNumber() const {
        return pdbResidueNumber;
    }

    const char getPdbInsertionCode() const {
        return pdbInsertionCode;
    }

    Compound::AtomIndex getAtomIndex(ResidueInfo::AtomIndex a) const {
        return atoms[a].biopolymerAtomIndex;
    }

    Compound::AtomIndex getAtomIndex(Compound::AtomName a) const {
        return atoms[atomIdsByName.find(a)->second].biopolymerAtomIndex;
    }

    const String& getAtomName(ResidueInfo::AtomIndex a) const {
        return atoms[a].pdbAtomName;
    }

    const std::set<Compound::AtomName>& getAtomSynonyms(ResidueInfo::AtomIndex a) const {
        return atoms[a].synonyms;
    }

    const Compound::Name& getName() const {return nameInCompound;}

    ResidueInfo& setPdbResidueNumber(int num) {
        pdbResidueNumber = num;
        return *this;
    }

    ResidueInfo::Index getIndex() const {return index;}

private:
    ResidueInfo::Index index;
    Compound::Name nameInCompound;
    Compound::Name pdbResidueName;
    std::set<Compound::Name> synonyms;
    int pdbResidueNumber;
    char pdbInsertionCode;
    std::vector<AtomInfo> atoms;
    std::map<const Compound::Name, AtomIndex> atomIdsByName;
    // TODO - put private data in ResidueInfoRep
};

class SimTK_MOLMODEL_EXPORT Biopolymer : public Molecule
{
public:
    typedef String Sequence;
    
    Biopolymer();

    int getNumResidues() const;

    const ResidueInfo& getResidue(
        ResidueInfo::Index residueIndex 
        ) const;

    ResidueInfo& updResidue(
        ResidueInfo::Index residueIndex 
        ); 

    const ResidueInfo& getResidue(
        Compound::Name residueName 
        ) const;

    ResidueInfo& updResidue(
        Compound::Name residueName 
        );

    const Compound::Name& getResidueName(
        ResidueInfo::Index residueIndex 
        ) const;

    // BiopolymerResidue createResidueCompound(ResidueInfo::Index r) const;

    Biopolymer& renumberPdbResidues(int firstPdbResidueNumber); 
    
    bool assignResidueBiotypes(ResidueInfo::Index, Ordinality::Residue);

    void assignBiotypes();

    ResidueInfo::Index appendResidue(
        const Compound::Name& resName, 
        const BiopolymerResidue& residue 
        );

    ResidueInfo::Index appendResidue(
        const Compound::Name& resName, 
        const BiopolymerResidue& residue, 
        BondMobility::Mobility mobility 
        );

    Biopolymer& setResidueBondMobility(ResidueInfo::Index, BondMobility::Mobility);
    MobilizedBodyIndex getResidueAtomMobilizedBodyIndex(ResidueInfo::Index res, ResidueInfo::AtomIndex a) const {
        return getAtomMobilizedBodyIndex(getResidue(res).getAtomIndex(a));
    }
    Vec3 getResidueAtomLocationInMobilizedBodyFrame(ResidueInfo::Index res, ResidueInfo::AtomIndex a) const {
        return getAtomLocationInMobilizedBodyFrame(getResidue(res).getAtomIndex(a));
    }

    virtual AtomTargetLocations createAtomTargets(const class PdbStructure& targetStructure) const;
    virtual AtomTargetLocations createAtomTargets(const class PdbChain& targetChain) const;

    SimTK_INSERT_DERIVED_HANDLE_DECLARATIONS(Biopolymer, BiopolymerRep, Molecule);

    // OBSOLETE; TODO: remove in SimTK 2.0
    int getNResidues() const {return getNumResidues();}

private:
};


} // namespace SimTK

#endif // SimTK_MOLMODEL_COMPOUND_H_