import simbios.simtk as simtk
from simbios.simtk import Transform
import sys

def main():
    # molecule-specialized simbody System
    system = simtk.CompoundSystem()

    matter = simtk.SimbodyMatterSubsystem(system)
    decorations = simtk.DecorationSubsystem(system)

    # molecular force field
    forceField = simtk.DuMMForceFieldSubsystem(system)
    forceField.loadAmber99Parameters()

    protein = simtk.Protein("SIMTK")
    protein.writeDefaultPdb(sys.stdout)

    protein.assignBiotypes()
    system.adoptCompound(protein)

    # finalize multibody system
    system.modelCompounds() 

    # show me a movie
    system.updDefaultSubsystem().addEventReporter(
        simtk.VTKEventReporter(system, 0.020))

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

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

    system.realize(state, simtk.Stage.Position)
    protein.writePdb(state, sys.stdout)

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

    system.realize(state, simtk.Stage.Position)
    protein.writePdb(state, sys.stdout)

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

    return 0


if __name__ == '__main__':
    main()