import sys
from simbios.simtk import *
from  simbios.simtk import Mobility as BondMobility

def main():
    system = CompoundSystem()
    matter = SimbodyMatterSubsystem(system)
    decorations = DecorationSubsystem(system)
    forceField = DuMMForceFieldSubsystem(system)

    forceField.loadAmber99Parameters()

    protein = Protein("SIMTK")
    protein.assignBiotypes()
    system.adoptCompound(protein)

    for b in range(protein.getNumBonds()):
        bondIx = Compound.BondIndex(b) 
        # set all bonds rigid
        protein.setBondMobility(
            BondMobility.Rigid, 
            bondIx)

    system.updDefaultSubsystem().addEventReporter(VTKEventReporter(system,
        0.020))

    # finalize multibody system
    system.modelCompounds() 

    # Maintain a constant temperature
    system.updDefaultSubsystem().addEventHandler(VelocityRescalingThermostat(
           system,  293.15, 0.1))

    # Instantiate simbody model
    system.realizeTopology()
    state = system.updDefaultState()

    # Relax the structure before dynamics run
    LocalEnergyMinimizer.minimizeEnergy(system, state, 15.0)

    # Simulate it.
    
    integ = VerletIntegrator(system)
    ts = TimeStepper(system, integ)
    ts.initialize(state)
    ts.stepTo(20.0)

if __name__ == '__main__':
    main()