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 Molmoldel(tm)                            *
 * -------------------------------------------------------------------------- *
 * This is part of the SimTK 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/home/molmodel. *
 *                                                                            *
 * Portions copyright (c) 2007-11 Stanford University and the Authors.        *
 * Authors: Michael Sherman                                                   *
 * Contributors: Chris Bruns, Peter Eastman, Randy Radmer, Samuel Flores      *
 *                                                                            *
 * 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(IncludedAtomIndex);
SimTK_DEFINE_UNIQUE_INDEX_TYPE(NonbondAtomIndex);
SimTK_DEFINE_UNIQUE_INDEX_TYPE(BondIndex);
SimTK_DEFINE_UNIQUE_INDEX_TYPE(ClusterIndex);
SimTK_DEFINE_UNIQUE_INDEX_TYPE(AtomClassIndex);
SimTK_DEFINE_UNIQUE_INDEX_TYPE(ChargedAtomTypeIndex);



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&);






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;
void clearIncludedNonbondAtomList();

void clearIncludedBondList();

void resetIncludedNonbondAtomListToDefault();

void resetIncludedBondListToDefault();

void includeNonbondAtom(DuMM::AtomIndex atom);

void includeAllNonbondAtomsForOneBody(MobilizedBodyIndex mobod);

void includeAllInterbodyBondsForOneAtom(DuMM::AtomIndex atom);

void includeAllInterbodyBondsWithBothAtoms(DuMM::AtomIndex atom1, 
                                           DuMM::AtomIndex atom2);

void includeAllInterbodyBondsWithBothAtoms(DuMM::BondIndex bond);

void includeAllInterbodyBondsForOneBody(MobilizedBodyIndex mobod);

void includeAllInterbodyBondsBetweenTwoBodies
   (MobilizedBodyIndex mobod1, MobilizedBodyIndex mobod2);

int getNumIncludedAtoms() const;
DuMM::AtomIndex getAtomIndexOfIncludedAtom
   (DuMM::IncludedAtomIndex incAtomIx) const;

int getNumNonbondAtoms() const;
DuMM::IncludedAtomIndex getIncludedAtomIndexOfNonbondAtom
   (DuMM::NonbondAtomIndex nonbondAtomIx) const;
DuMM::AtomIndex getAtomIndexOfNonbondAtom
   (DuMM::NonbondAtomIndex nonbondAtomIx) const
{   return getAtomIndexOfIncludedAtom
            (getIncludedAtomIndexOfNonbondAtom(nonbondAtomIx)); }
    // DEFINE FORCE FIELD PARAMETERS
    
void defineAtomClass(DuMM::AtomClassIndex     atomClassIx, 
                     const char*              atomClassName,
                     int  atomicNumber,  int  expectedValence,
                     Real vdwRadiusInNm, Real vdwWellDepthInKJ) 
{   defineIncompleteAtomClass(atomClassIx, atomClassName, atomicNumber, 
                              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);
}

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;
bool hasAtomClass(const String& atomClassName) const;
DuMM::AtomClassIndex getAtomClassIndex(const String& atomClassName) const;
DuMM::AtomClassIndex getNextUnusedAtomClassIndex() const;
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!
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;
bool hasChargedAtomType(const String& chargedTypeName) const;
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 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 setTracing(bool);

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

void setNumThreadsRequested(int);
int getNumThreadsRequested() const;

bool isUsingMultithreadedComputation() const;

int getNumThreadsInUse() const;

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


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();


void setTraceOpenMM(bool shouldTrace) {setTracing(shouldTrace);}
protected: // don't let doxygen see these
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

 // End of MolecularMechanics module

#endif // SimTK_MOLMODEL_DUMM_FORCE_FIELD_SUBSYSTEM_H_