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

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

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

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

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

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

   system.adoptCompound(rna)
   system.modelCompounds()

   # Write pdb file to use for initial topology in VMD
   pdbOut = open("ThreeAdenylates.pdb", "w")
   rna.writeDefaultPdb(pdbOut)

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

   # Send coordinates to VMD
   system.updDefaultSubsystem().addEventReporter(PeriodicVmdReporter(system, 0.015, 3000, True))

   # Prepare for molecular dynamics
   system.realizeTopology()
   state = system.updDefaultState()

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

   integrator = VerletIntegrator(system)
   timeStepper = TimeStepper(system, integrator)
   timeStepper.initialize(state)

   # Start simulation
   timeStepper.stepTo(10.0)

if __name__ == "__main__":
    main()