<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Tully model 1 from JCP 93, 1061 (1990). Details of this usage of the MMST matches that of Ananth et al. (JCP 127, 084114 (2007)). 
The potentials are
Step1:  Running a trajectory using just a single surface
Step2:  Running a trajectory using MMST with no coupling
This should give the same results as the single surface.
Step3:  Running a trajectory with the real MMST
<ASSISTANT_TASK:>
Python Code:
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import simtk.openmm as mm
import simtk.openmm.app as app
import simtk.unit as unit
sys11 = mm.openmm.System()
sys12 = mm.openmm.System()
sys22 = mm.openmm.System()
sys00 = mm.openmm.System()
for sys in [sys11, sys12, sys22, sys00]:
    mass = 1980.0 * unit.amu
    sys.addParticle(mass)
V11 = mm.openmm.CustomExternalForce("V0*tanh(a*x)")
V11.addGlobalParameter("V0", 0.01)
V11.addGlobalParameter("a", 1.6)
V11.addParticle(0, ())
V22 = mm.openmm.CustomExternalForce("-V0*tanh(a*x)")
V22.addGlobalParameter("V0", 0.01)
V22.addGlobalParameter("a", 1.6)
V22.addParticle(0, ())
V12 = mm.openmm.CustomExternalForce("C*exp(-D*(x+f))")
V12.addGlobalParameter("C", 0.005)
V12.addGlobalParameter("D", 1.0)
V12.addGlobalParameter("f", 0.0)
V12.addParticle(0, ())
V00 = mm.openmm.CustomExternalForce("0.0*x")
V00.addParticle(0, ())
sys00.addForce(V00)
sys11.addForce(V11)
sys12.addForce(V12)
sys22.addForce(V22)
topology = app.Topology()
dt = 5*46.0 * unit.femtoseconds 
init_pos = np.array([[-5.0, 0.0, 0.0]]) #* unit.nanometer
#init_vel = np.array([[19.9/1980.0, 0.0, 0.0]]) #* unit.nanometer / unit.picosecond
init_vel = np.array([[0.0022, 0.0, 0.0]])
Explanation: Tully model 1 from JCP 93, 1061 (1990). Details of this usage of the MMST matches that of Ananth et al. (JCP 127, 084114 (2007)). 
The potentials are:
$V_{11}(R) = V_0 \tanh(a R) \
V_{22}(R) = -V_{11}(R) \
V_{12}(R) = V_{21}(R) = C e^{-D(R+f)^2}$ 
with $V_0 = 0.01$, $a=1.6$, $C=0.005$, $D=1$, and $f=0$.
End of explanation
integ = mm.VerletIntegrator(dt)
simulation = app.Simulation(topology, sys11, integ)
simulation.context.setPositions(init_pos)
simulation.context.setVelocities(init_vel)
traj = []
forces = []
energies = []
for i in range(40000):
    state = simulation.context.getState(getPositions=True,
        getForces=True, getEnergy=True
    )
    pos = state.getPositions(asNumpy=True)
    force = state.getForces(asNumpy=True)
    energy = state.getPotentialEnergy()
    forces.append(force[0][0] / force.unit)
    energies.append(energy / energy.unit)
    traj.append(pos[0][0] / pos.unit)
    simulation.step(1)
plt.plot(traj, forces)
plt.plot(traj, energies)
plt.plot(traj)
Explanation: Running a trajectory using just a single surface
End of explanation
import openmm_mmst as mmst
Hmat = [[sys11, sys00], [sys00, sys22]]
#mmst_integ = NonadiabaticIntegrator(Hmat)
Explanation: Running a trajectory using MMST with no coupling
This should give the same results as the single surface.
End of explanation
Hmat = [[sys11, sys12], [sys12, sys22]]
#mmst_integ = NonadiabaticIntegrator(Hmat)
f = simulation.system.getForce(0)
Explanation: Running a trajectory with the real MMST
End of explanation
<END_TASK>