DuMMForceFieldSubsystem.h

Go to the documentation of this file.

#ifndef SimTK_MOLMODEL_DUMM_FORCE_FIELD_SUBSYSTEM_H_
#define SimTK_MOLMODEL_DUMM_FORCE_FIELD_SUBSYSTEM_H_

/* -------------------------------------------------------------------------- *
 *                      SimTK Core: SimTK Simbody(tm)                         *
 * -------------------------------------------------------------------------- *
 * 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-9 Stanford University and the Authors.         *
 * Authors: Michael Sherman                                                   *
 * Contributors: Chris Bruns, Peter Eastman, Randy Radmer                     *
 *                                                                            *
 * 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 "SimTKcommon.h"
#include "simbody/internal/ForceSubsystem.h"

#include "molmodel/internal/common.h"
#include "molmodel/internal/Biotype.h"

#include <cassert>


namespace SimTK {

class MolecularMechanicsSystem;

namespace DuMM {
SimTK_DEFINE_UNIQUE_INDEX_TYPE(AtomIndex);
SimTK_DEFINE_UNIQUE_INDEX_TYPE(AtomClassIndex);
SimTK_DEFINE_UNIQUE_INDEX_TYPE(ChargedAtomTypeIndex);
SimTK_DEFINE_UNIQUE_INDEX_TYPE(BondIndex);
SimTK_DEFINE_UNIQUE_INDEX_TYPE(ClusterIndex);



class CustomBondStretch {
public:
    virtual ~CustomBondStretch() {}
    virtual Real calcEnergy(Real distance) const = 0;
    virtual Real calcForce(Real distance) const = 0;
};

class CustomBondBend {
public:
    virtual ~CustomBondBend() {}
    virtual Real calcEnergy(Real bendAngle) const = 0;
    virtual Real calcTorque(Real bendAngle) const = 0;
};

class CustomBondTorsion {
public:
    virtual ~CustomBondTorsion() {}
    virtual Real calcEnergy(Real dihedralAngle) const = 0;
    virtual Real calcTorque(Real dihedralAngle) const = 0;
};

static const Real Ang2Nm  = (Real)0.1L; 
static const Real Nm2Ang  = (Real)10.L; 
static const Real Deg2Rad = (Real)SimTK_DEGREE_TO_RADIAN; 
static const Real Rad2Deg = (Real)SimTK_RADIAN_TO_DEGREE; 
static const Real KJ2Kcal = (Real)SimTK_KJOULE_TO_KCAL;   
static const Real Kcal2KJ = (Real)SimTK_KCAL_TO_KJOULE;   

static const Real Sigma2Radius = (Real)std::pow(2.L,  1.L/6.L);
static const Real Radius2Sigma = (Real)std::pow(2.L, -1.L/6.L);

} // namespace DuMM


class SimTK_MOLMODEL_EXPORT DuMMForceFieldSubsystem : public ForceSubsystem {
public:
    enum VdwMixingRule {
        WaldmanHagler       = 1,    
        HalgrenHHG          = 2,    
        Jorgensen           = 3,    
        LorentzBerthelot    = 4,    
        Kong                = 5     
    };

    DuMMForceFieldSubsystem();
    explicit DuMMForceFieldSubsystem(MolecularMechanicsSystem&);


    const Vec3& getAtomLocationInGround(const State& state, DuMM::AtomIndex atomIx);
    const Vec3& getAtomVelocityInGround(const State& state, DuMM::AtomIndex atomIx);

    const Vec3& getAtomBodyStationInGround(const State& state, DuMM::AtomIndex atomIx);



    DuMM::AtomIndex addAtom(DuMM::ChargedAtomTypeIndex chargedAtomTypeIx);

    DuMM::BondIndex addBond(DuMM::AtomIndex atom1Ix, DuMM::AtomIndex atom2Ix);

    DuMM::AtomIndex  getBondAtom(DuMM::BondIndex bond, int which) const;

    int getNumAtoms() const;
    int getNumBonds() const;

    Real   getAtomMass(DuMM::AtomIndex atomIx) const;
    int    getAtomElement(DuMM::AtomIndex atomIx) const;
    Real   getAtomRadius(DuMM::AtomIndex atomIx) const;
    MobilizedBodyIndex getAtomBody(DuMM::AtomIndex atomIx) const;
    Vec3   getAtomStationOnBody(DuMM::AtomIndex atomIx) const;

    Vec3   getAtomStationInCluster(DuMM::AtomIndex atomIx, DuMM::ClusterIndex clusterIx) const;

    Vec3   getElementDefaultColor(int atomicNumber) const;
    Vec3   getAtomDefaultColor(DuMM::AtomIndex atomIx) const;




    DuMM::ClusterIndex createCluster(const char* clusterName);

    void placeAtomInCluster(DuMM::AtomIndex atomIx, DuMM::ClusterIndex clusterIx, const Vec3& station);

    void placeClusterInCluster(DuMM::ClusterIndex childClusterIndex, DuMM::ClusterIndex parentClusterIndex, 
                               const Transform& X_PC);

    MassProperties calcClusterMassProperties(DuMM::ClusterIndex clusterIx, const Transform& X_BC = Transform()) const;

    void attachClusterToBody(DuMM::ClusterIndex clusterIx, MobilizedBodyIndex body, const Transform& X_BC = Transform());

    void attachAtomToBody(DuMM::AtomIndex atomIx, MobilizedBodyIndex body, const Vec3& station = Vec3(0));

    MobilizedBodyIndex getClusterBody(DuMM::ClusterIndex clusterIx) const;

    Transform getClusterPlacementOnBody(DuMM::ClusterIndex clusterIx) const;

    Transform getClusterPlacementInCluster(DuMM::ClusterIndex childClusterIndex, DuMM::ClusterIndex parentClusterIndex) const;



        // DEFINE FORCE FIELD PARAMETERS
    
    void defineAtomClass(DuMM::AtomClassIndex atomClassIx, const char* atomClassName,
                         int elementNumber, int expectedValence,
                         Real vdwRadiusInNm, Real vdwWellDepthInKJ) {
         defineIncompleteAtomClass(atomClassIx, atomClassName, elementNumber, expectedValence);
         setAtomClassVdwParameters(atomClassIx, vdwRadiusInNm, vdwWellDepthInKJ);
    }

    void defineAtomClass_KA(DuMM::AtomClassIndex atomClassIx, const char* atomClassName,
                            int element, int valence,
                            Real vdwRadiusInAng, Real vdwWellDepthInKcal)
    {
        defineAtomClass(atomClassIx, atomClassName, element, valence,
            vdwRadiusInAng*DuMM::Ang2Nm, vdwWellDepthInKcal*DuMM::Kcal2KJ);
    }
    // For backwards compatibility and compactness of expression, permit integer indices
    void defineAtomClass_KA(int atomClassIx, const char* atomClassName,
                            int element, int valence,
                            Real vdwRadiusInAng, Real vdwWellDepthInKcal)
    {
        defineAtomClass_KA((DuMM::AtomClassIndex)atomClassIx, atomClassName, element, valence, vdwRadiusInAng, vdwWellDepthInKcal);
    }

    bool hasAtomClass(DuMM::AtomClassIndex) const; // TODO
    bool hasAtomClass(const String& atomClassName) const; // TODO
    DuMM::AtomClassIndex getAtomClassIndex(const String& atomClassName) const;
    DuMM::AtomClassIndex getNextUnusedAtomClassIndex() const; // TODO
    
    DuMM::AtomClassIndex getAtomClassIndex(DuMM::AtomIndex atomIx) const;
    Real getVdwRadius(DuMM::AtomClassIndex atomClassIx) const;
    Real getVdwWellDepth(DuMM::AtomClassIndex atomClassIx) const;

    void defineChargedAtomType(DuMM::ChargedAtomTypeIndex atomTypeIx, const char* atomTypeName,
        DuMM::AtomClassIndex atomClassIx, Real partialChargeInE) {
            defineIncompleteChargedAtomType(atomTypeIx, atomTypeName, atomClassIx);
            setChargedAtomTypeCharge(atomTypeIx, partialChargeInE);
    }

    void defineChargedAtomType_KA(DuMM::ChargedAtomTypeIndex atomTypeIx, const char* atomTypeName,
                                  DuMM::AtomClassIndex atomClassIx, Real partialChargeInE)
    {
        defineChargedAtomType(atomTypeIx, atomTypeName, atomClassIx, partialChargeInE); // easy!
    }
    // For backwards compatibility and compactness of expression, permit integer indices
    void defineChargedAtomType_KA(int atomTypeIx, const char* atomTypeName,
                                  int atomClassIx, Real partialChargeInE)
    {
        defineChargedAtomType_KA((DuMM::ChargedAtomTypeIndex) atomTypeIx, atomTypeName, (DuMM::AtomClassIndex)atomClassIx, partialChargeInE);
    }

    bool hasChargedAtomType(DuMM::ChargedAtomTypeIndex) const; // TODO
    bool hasChargedAtomType(const String& chargedTypeName) const; // TODO
    DuMM::ChargedAtomTypeIndex getChargedAtomTypeIndex(const String& chargedTypeName) const; // TODO
    DuMM::ChargedAtomTypeIndex getNextUnusedChargedAtomTypeIndex() const; // TODO



    const char* getVdwMixingRuleName(VdwMixingRule) const;

    void setVdwMixingRule(VdwMixingRule);
    VdwMixingRule getVdwMixingRule() const;

    void setVdw12ScaleFactor(Real); 
    void setVdw13ScaleFactor(Real); 
    void setVdw14ScaleFactor(Real); 
    void setVdw15ScaleFactor(Real); 

    void setCoulomb12ScaleFactor(Real); 
    void setCoulomb13ScaleFactor(Real); 
    void setCoulomb14ScaleFactor(Real); 
    void setCoulomb15ScaleFactor(Real); 



    void loadAmber99Parameters();

    // void dumpCForcefieldParameters(std::ostream& os, const String methodName) const;


    void populateFromTinkerParameterFile(
        std::istream& 
        );

    void setBiotypeChargedAtomType(
        DuMM::ChargedAtomTypeIndex chargedAtomTypeIndex, 
        BiotypeIndex biotypeIx 
        );

    DuMM::ChargedAtomTypeIndex getBiotypeChargedAtomType(BiotypeIndex biotypeIx) const;

    std::ostream& generateBiotypeChargedAtomTypeSelfCode(std::ostream& os) const;



    void defineBondStretch(DuMM::AtomClassIndex class1, DuMM::AtomClassIndex class2,
                           Real stiffnessInKJperNmSq, Real nominalLengthInNm);

    void defineBondStretch_KA(DuMM::AtomClassIndex class1, DuMM::AtomClassIndex class2,
                              Real stiffnessInKcalPerAngSq, Real nominalLengthInAng)
    {
        defineBondStretch(class1, class2, 
                          stiffnessInKcalPerAngSq * DuMM::Kcal2KJ/square(DuMM::Ang2Nm),
                          nominalLengthInAng      * DuMM::Ang2Nm);
    }
    void defineBondStretch_KA(int class1, int class2,
                              Real stiffnessInKcalPerAngSq, Real nominalLengthInAng)
    {   defineBondStretch_KA((DuMM::AtomClassIndex)class1, (DuMM::AtomClassIndex)class2,
                             stiffnessInKcalPerAngSq, nominalLengthInAng); }

    void defineCustomBondStretch(DuMM::AtomClassIndex       class1, 
                                 DuMM::AtomClassIndex       class2,
                                 DuMM::CustomBondStretch*   bondStretchTerm);



    void defineBondBend(DuMM::AtomClassIndex class1, DuMM::AtomClassIndex class2, DuMM::AtomClassIndex class3,
                        Real stiffnessInKJPerRadSq, Real nominalAngleInDeg);

    void defineBondBend_KA(DuMM::AtomClassIndex class1, DuMM::AtomClassIndex class2, DuMM::AtomClassIndex class3,
                           Real stiffnessInKcalPerRadSq, Real nominalAngleInDeg) 
    {
        defineBondBend(class1,class2,class3,
                       stiffnessInKcalPerRadSq * DuMM::Kcal2KJ,
                       nominalAngleInDeg);
    }
    void defineBondBend_KA(int class1, int class2, int class3,
                           Real stiffnessInKcalPerRadSq, Real nominalAngleInDeg)
    {   defineBondBend_KA((DuMM::AtomClassIndex)class1,(DuMM::AtomClassIndex)class2,(DuMM::AtomClassIndex)class3,
                          stiffnessInKcalPerRadSq,nominalAngleInDeg); }

    void defineCustomBondBend(DuMM::AtomClassIndex      class1, 
                              DuMM::AtomClassIndex      class2, 
                              DuMM::AtomClassIndex      class3,
                              DuMM::CustomBondBend*     bondBendTerm);



    void defineBondTorsion
       (DuMM::AtomClassIndex class1, DuMM::AtomClassIndex class2, DuMM::AtomClassIndex class3, DuMM::AtomClassIndex class4, 
        int periodicity, Real ampInKJ, Real phaseInDegrees);
    void defineBondTorsion
       (DuMM::AtomClassIndex class1, DuMM::AtomClassIndex class2, DuMM::AtomClassIndex class3, DuMM::AtomClassIndex class4,  
        int periodicity1, Real amp1InKJ, Real phase1InDegrees,
        int periodicity2, Real amp2InKJ, Real phase2InDegrees);
    void defineBondTorsion
       (DuMM::AtomClassIndex class1, DuMM::AtomClassIndex class2, DuMM::AtomClassIndex class3, DuMM::AtomClassIndex class4, 
        int periodicity1, Real amp1InKJ, Real phase1InDegrees,
        int periodicity2, Real amp2InKJ, Real phase2InDegrees,
        int periodicity3, Real amp3InKJ, Real phase3InDegrees);

    void defineBondTorsion_KA
       (DuMM::AtomClassIndex class1, DuMM::AtomClassIndex class2, DuMM::AtomClassIndex class3, DuMM::AtomClassIndex class4, 
        int periodicity1, Real amp1InKcal, Real phase1InDegrees)
    { 
        defineBondTorsion(class1,class2,class3,class4,
                          periodicity1, amp1InKcal * DuMM::Kcal2KJ, phase1InDegrees);
    }
    void defineBondTorsion_KA
       (DuMM::AtomClassIndex class1, DuMM::AtomClassIndex class2, DuMM::AtomClassIndex class3, DuMM::AtomClassIndex class4,  
        int periodicity1, Real amp1InKcal, Real phase1InDegrees,
        int periodicity2, Real amp2InKcal, Real phase2InDegrees)
    {
        defineBondTorsion(class1,class2,class3,class4,
                          periodicity1, amp1InKcal * DuMM::Kcal2KJ, phase1InDegrees,
                          periodicity2, amp2InKcal * DuMM::Kcal2KJ, phase2InDegrees);
    }
    void defineBondTorsion_KA
       (DuMM::AtomClassIndex class1, DuMM::AtomClassIndex class2, DuMM::AtomClassIndex class3, DuMM::AtomClassIndex class4, 
        int periodicity1, Real amp1InKcal, Real phase1InDegrees,
        int periodicity2, Real amp2InKcal, Real phase2InDegrees,
        int periodicity3, Real amp3InKcal, Real phase3InDegrees)
    {
        defineBondTorsion(class1,class2,class3,class4,
                          periodicity1, amp1InKcal * DuMM::Kcal2KJ, phase1InDegrees,
                          periodicity2, amp2InKcal * DuMM::Kcal2KJ, phase2InDegrees,
                          periodicity3, amp3InKcal * DuMM::Kcal2KJ, phase3InDegrees);
    }
    void defineBondTorsion_KA
       (int class1, int class2, int class3, int class4, 
        int periodicity1, Real amp1InKcal, Real phase1InDegrees)
    {
        defineBondTorsion_KA
               ((DuMM::AtomClassIndex)class1, (DuMM::AtomClassIndex)class2, (DuMM::AtomClassIndex)class3, (DuMM::AtomClassIndex)class4, 
                periodicity1, amp1InKcal, phase1InDegrees);
    }
    void defineBondTorsion_KA
       (int class1, int class2, int class3, int class4, 
        int periodicity1, Real amp1InKcal, Real phase1InDegrees,
        int periodicity2, Real amp2InKcal, Real phase2InDegrees)
    {
        defineBondTorsion_KA
               ((DuMM::AtomClassIndex)class1, (DuMM::AtomClassIndex)class2, (DuMM::AtomClassIndex)class3, (DuMM::AtomClassIndex)class4, 
                periodicity1, amp1InKcal, phase1InDegrees,
                periodicity2, amp2InKcal, phase2InDegrees);
    }
    void defineBondTorsion_KA
       (int class1, int class2, int class3, int class4, 
        int periodicity1, Real amp1InKcal, Real phase1InDegrees,
        int periodicity2, Real amp2InKcal, Real phase2InDegrees,
        int periodicity3, Real amp3InKcal, Real phase3InDegrees)
    {
        defineBondTorsion_KA
               ((DuMM::AtomClassIndex)class1, (DuMM::AtomClassIndex)class2, (DuMM::AtomClassIndex)class3, (DuMM::AtomClassIndex)class4, 
                periodicity1, amp1InKcal, phase1InDegrees,
                periodicity2, amp2InKcal, phase2InDegrees,
                periodicity3, amp3InKcal, phase3InDegrees);
    }


    void defineCustomBondTorsion(DuMM::AtomClassIndex       class1, 
                                 DuMM::AtomClassIndex       class2, 
                                 DuMM::AtomClassIndex       class3,
                                 DuMM::AtomClassIndex       class4,
                                 DuMM::CustomBondTorsion*   bondTorsionTerm);



    void defineAmberImproperTorsion
       (DuMM::AtomClassIndex class1, DuMM::AtomClassIndex class2, DuMM::AtomClassIndex class3, DuMM::AtomClassIndex class4,
        int periodicity, Real ampInKJ, Real phaseInDegrees);
    void defineAmberImproperTorsion
       (DuMM::AtomClassIndex class1, DuMM::AtomClassIndex class2, DuMM::AtomClassIndex class3, DuMM::AtomClassIndex class4,
        int periodicity1, Real amp1InKJ, Real phase1InDegrees,
        int periodicity2, Real amp2InKJ, Real phase2InDegrees);
    void defineAmberImproperTorsion
       (DuMM::AtomClassIndex class1, DuMM::AtomClassIndex class2, DuMM::AtomClassIndex class3, DuMM::AtomClassIndex class4,
        int periodicity1, Real amp1InKJ, Real phase1InDegrees,
        int periodicity2, Real amp2InKJ, Real phase2InDegrees,
        int periodicity3, Real amp3InKJ, Real phase3InDegrees);

    void defineAmberImproperTorsion_KA
       (DuMM::AtomClassIndex class1, DuMM::AtomClassIndex class2, DuMM::AtomClassIndex class3, DuMM::AtomClassIndex class4,
        int periodicity1, Real amp1InKcal, Real phase1InDegrees)
    {
        defineAmberImproperTorsion(class1,class2,class3,class4,
                          periodicity1, amp1InKcal * DuMM::Kcal2KJ, phase1InDegrees);
    }
    void defineAmberImproperTorsion_KA
       (DuMM::AtomClassIndex class1, DuMM::AtomClassIndex class2, DuMM::AtomClassIndex class3, DuMM::AtomClassIndex class4,
        int periodicity1, Real amp1InKcal, Real phase1InDegrees,
        int periodicity2, Real amp2InKcal, Real phase2InDegrees)
    {
        defineAmberImproperTorsion(class1,class2,class3,class4,
                          periodicity1, amp1InKcal * DuMM::Kcal2KJ, phase1InDegrees,
                          periodicity2, amp2InKcal * DuMM::Kcal2KJ, phase2InDegrees);
    }
    void defineAmberImproperTorsion_KA
       (DuMM::AtomClassIndex class1, DuMM::AtomClassIndex class2, DuMM::AtomClassIndex class3, DuMM::AtomClassIndex class4,
        int periodicity1, Real amp1InKcal, Real phase1InDegrees,
        int periodicity2, Real amp2InKcal, Real phase2InDegrees,
        int periodicity3, Real amp3InKcal, Real phase3InDegrees)
    {
        defineAmberImproperTorsion(class1,class2,class3,class4,
                          periodicity1, amp1InKcal * DuMM::Kcal2KJ, phase1InDegrees,
                          periodicity2, amp2InKcal * DuMM::Kcal2KJ, phase2InDegrees,
                          periodicity3, amp3InKcal * DuMM::Kcal2KJ, phase3InDegrees);
    }


    void setSolventDielectric(Real);    // typically 80 for water
    void setSoluteDielectric(Real);     // typically 1 or 2 for protein
    Real getSolventDielectric() const;
    Real getSoluteDielectric() const;

    void setGbsaIncludeAceApproximation(bool);
    void setGbsaIncludeAceApproximationOn()  {setGbsaIncludeAceApproximation(true );}
    void setGbsaIncludeAceApproximationOff() {setGbsaIncludeAceApproximation(false);}

    void setVdwGlobalScaleFactor(Real);                 
    void setCoulombGlobalScaleFactor(Real);             
    void setGbsaGlobalScaleFactor(Real);                
    void setBondStretchGlobalScaleFactor(Real);         
    void setBondBendGlobalScaleFactor(Real);            
    void setBondTorsionGlobalScaleFactor(Real);         
    void setAmberImproperTorsionGlobalScaleFactor(Real);
    void setCustomBondStretchGlobalScaleFactor(Real);   
    void setCustomBondBendGlobalScaleFactor(Real);      
    void setCustomBondTorsionGlobalScaleFactor(Real);   

    Real getVdwGlobalScaleFactor() const;                   
    Real getCoulombGlobalScaleFactor() const;               
    Real getGbsaGlobalScaleFactor() const;                  
    Real getBondStretchGlobalScaleFactor() const;           
    Real getBondBendGlobalScaleFactor() const;              
    Real getBondTorsionGlobalScaleFactor() const;           
    Real getAmberImproperTorsionGlobalScaleFactor() const;  
    Real getCustomBondStretchGlobalScaleFactor() const;     
    Real getCustomBondBendGlobalScaleFactor() const;        
    Real getCustomBondTorsionGlobalScaleFactor() const;     

    void setAllGlobalScaleFactors(Real s) {
        setVdwGlobalScaleFactor(s);
        setCoulombGlobalScaleFactor(s);
        setGbsaGlobalScaleFactor(s);
        setBondStretchGlobalScaleFactor(s);
        setBondBendGlobalScaleFactor(s);
        setBondTorsionGlobalScaleFactor(s);
        setAmberImproperTorsionGlobalScaleFactor(s);
        setCustomBondStretchGlobalScaleFactor(s);
        setCustomBondBendGlobalScaleFactor(s);
        setCustomBondTorsionGlobalScaleFactor(s);
    }


    void setUseMultithreadedComputation(bool);
    bool getUseMultithreadedComputation() const;

    void setUseOpenMMAcceleration(bool);
    bool getUseOpenMMAcceleration() const;

    void setTraceOpenMM(bool);

    void setAllowOpenMMReference(bool);
    bool getAllowOpenMMReference() const;

    bool isUsingOpenMM() const;
    std::string getOpenMMPlatformInUse() const;



    long long getForceEvaluationCount() const;

    void dump() const;

    void dumpCForceFieldParameters(std::ostream& os, const String& methodName = "loadParameters") const;

    void loadTestMoleculeParameters();

    int getNAtoms() const {return getNumAtoms();}
    int getNBonds() const {return getNumBonds();}


protected:
    SimTK_PIMPL_DOWNCAST(DuMMForceFieldSubsystem, ForceSubsystem);
    SimTK_PIMPL_DOWNCAST(DuMMForceFieldSubsystem, Subsystem);

    void defineIncompleteAtomClass(
        DuMM::AtomClassIndex classIx, 
        const char* name, 
        int elementNumber, 
        int valence);

    void defineIncompleteAtomClass_KA(
        DuMM::AtomClassIndex classIx, 
        const char* name, 
        int elementNumber, 
        int valence) 
    {
        defineIncompleteAtomClass(
            classIx, 
            name, 
            elementNumber, 
            valence
            );
    }

    void setAtomClassVdwParameters(DuMM::AtomClassIndex atomClassIx, Real vdwRadiusInNm, Real vdwWellDepthInKJPerMol);
    void setAtomClassVdwParameters_KA(DuMM::AtomClassIndex atomClassIx, Real radiusInAng, Real wellDepthInKcal) {
        setAtomClassVdwParameters(atomClassIx, radiusInAng*DuMM::Ang2Nm, wellDepthInKcal*DuMM::Kcal2KJ);
    }

    bool isValidAtomClass(DuMM::AtomClassIndex) const;
    
    void defineIncompleteChargedAtomType(
        DuMM::ChargedAtomTypeIndex typeIx, 
        const char* name,
        DuMM::AtomClassIndex classIx);

    void defineIncompleteChargedAtomType_KA(
        DuMM::ChargedAtomTypeIndex typeIx, 
        const char* name,
        DuMM::AtomClassIndex classIx) 
    {
        defineIncompleteChargedAtomType(typeIx, name, classIx);
    }

    void setChargedAtomTypeCharge(DuMM::ChargedAtomTypeIndex, Real charge);
    void setChargedAtomTypeCharge_KA(DuMM::ChargedAtomTypeIndex chargedAtomTypeIx, Real charge) {
        setChargedAtomTypeCharge(chargedAtomTypeIx, charge);
    }

private:
    class DuMMForceFieldSubsystemRep& updRep();
    const DuMMForceFieldSubsystemRep& getRep() const;

friend class MolecularMechanicsSystem;
};

class Amber99ForceSubsystem : public DuMMForceFieldSubsystem {
public:
    explicit Amber99ForceSubsystem(MolecularMechanicsSystem& system) 
        : DuMMForceFieldSubsystem(system)
    {
        loadAmber99Parameters();
    }
};


} // namespace SimTK



#endif // SimTK_MOLMODEL_DUMM_FORCE_FIELD_SUBSYSTEM_H_