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

def main():
   # Initialize simbody objects
   system = simtk.CompoundSystem()
   matter = simtk.SimbodyMatterSubsystem(system)
   decoration = simtk.DecorationSubsystem(system)
   forces = simtk.DuMMForceFieldSubsystem(system)
   forces.loadAmber99Parameters()

   # Initialize molmodel objects
   rna = simtk.RNA("AAA")

   # Set first residue to Euclidean mobilities
   rna.setResidueBondMobility(simtk.ResidueInfo.Index(0), BondMobility.Free)

   # Leave second residue at default combination of Torsion and Rigid mobilities

   # Set third residue to Rigid
   rna.setResidueBondMobility(simtk.ResidueInfo.Index(2), BondMobility.Rigid)

   # Finalize the multibody system
   system.adoptCompound(rna)
   system.modelCompounds()

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

   # Show me a movie
   system.updDefaultSubsystem().addEventReporter(simtk.VTKEventReporter(system, 0.015) )

   system.realizeTopology()
   state = system.updDefaultState()

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

   # Prepare for molecular dynamics
   integrator = simtk.VerletIntegrator(system)
   integrator.setAccuracy(0.001)
   timeStepper = simtk.TimeStepper(system, integrator)
   timeStepper.initialize(state)

   # Start simulation
   timeStepper.stepTo(10.0)

if __name__ == '__main__':
    main()