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()

    pdbFile = open("1AKGtrim.pdb", "r")
    protein = Protein(pdbFile)

    protein.assignBiotypes()
    system.adoptCompound(protein)
    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()
    LocalEnergyMinimizer.minimizeEnergy(system, state, 15.0)
    
    # Simulate it.
    integ = VerletIntegrator(system)
    integ.setAccuracy(1e-2)
    ts = TimeStepper(system, integ)
    ts.initialize(state)
    ts.stepTo(1.0)

if __name__ == "__main__":
    main()