from simbios.simtk import *

def main():
    # Load the PDB file and construct the system.
    system = CompoundSystem()
    matter = SimbodyMatterSubsystem(system)
    decorations = DecorationSubsystem(system)
    forceField = DuMMForceFieldSubsystem(system)
    forceField.loadAmber99Parameters()

    pdb = PDBReader("1AKG.pdb")
    pdb.createCompounds(system)
    system.modelCompounds()

    system.updDefaultSubsystem().addEventHandler(VelocityRescalingThermostat(
        system, 293.15, 0.1))
    system.updDefaultSubsystem().addEventReporter(VTKEventReporter(system,
        0.025))

    system.realizeTopology()
    
    # Create an initial state for the simulation.
    
    state = system.updDefaultState()
    pdb.createState(system, state)
    LocalEnergyMinimizer.minimizeEnergy(system, state, 15.0)
    
    # Simulate it.
    
    integ = VerletIntegrator(system)
    integ.setAccuracy(1e-2)
    ts = TimeStepper(system, integ)
    ts.initialize(state)
    ts.stepTo(2.0)

if __name__ == '__main__':
    main()