subroutine CalcEnergy(Pos, M, L, rc, PEnergy, Dim, NAtom) implicit none integer, intent(in) :: M, Dim, NAtom real(8), intent(in), dimension(0:NAtom-1, 0:Dim-1) :: Pos real(8), intent(in) :: L, rc real(8), intent(out) :: PEnergy real(8), parameter :: k = 3000., r0 = 1. !======== YOUR CODE HERE ======== end subroutine subroutine CalcEnergyForces(Pos, M, L, rc, PEnergy, Forces, Dim, NAtom) implicit none integer, intent(in) :: M, Dim, NAtom real(8), intent(in), dimension(0:NAtom-1, 0:Dim-1) :: Pos real(8), intent(in) :: L, rc real(8), intent(out) :: PEnergy real(8), intent(inout), dimension(0:NAtom-1, 0:Dim-1) :: Forces !f2py intent(in, out) :: Forces real(8), parameter :: k = 3000., r0 = 1. !======== YOUR CODE HERE ======== end subroutine subroutine VVIntegrate(Pos, Vel, Accel, M, L, rc, dt, KEnergy, PEnergy, Dim, NAtom) implicit none integer, intent(in) :: M, Dim, NAtom real(8), intent(in) :: L, rc, dt real(8), intent(inout), dimension(0:NAtom-1, 0:Dim-1) :: Pos, Vel, Accel !f2py intent(in,out) :: Pos, Vel, Accel real(8), intent(out) :: KEnergy, PEnergy external :: CalcEnergyForces Pos = Pos + dt * Vel + 0.5 * dt*dt * Accel Vel = Vel + 0.5 * dt * Accel call CalcEnergyForces(Pos, M, L, rc, PEnergy, Accel, Dim, NAtom) Vel = Vel + 0.5 * dt * Accel KEnergy = 0.5 * sum(Vel*Vel) end subroutine