def METHOD_NAME(self, epsilon, state):
    axpy = linalg.blas.get_blas_funcs("axpy", dtype=self._dtype)
    pot = self._potential
    q_new = state.q.data.copy()
    p_new = state.p.data.copy()
    v_new = np.empty_like(q_new)
    q_new_grad = np.empty_like(q_new)
    dt = 0.5 * epsilon
    # p is already stored in p_new
    # p_new = p + dt * q_grad
    axpy(state.q_grad, p_new, a=dt)
    pot.velocity(p_new, out=v_new)
    # q is already stored in q_new
    # q_new = q + epsilon * v_new
    axpy(v_new, q_new, a=epsilon)
    p_new = RaveledVars(p_new, state.p.point_map_info)
    q_new = RaveledVars(q_new, state.q.point_map_info)
    logp = self._logp_dlogp_func(q_new, grad_out=q_new_grad)
    # p_new = p_new + dt * q_new_grad
    axpy(q_new_grad, p_new.data, a=dt)
    kinetic = pot.velocity_energy(p_new.data, v_new)
    energy = kinetic - logp
    return State(
        q_new,
        p_new,
        v_new,
        q_new_grad,
        energy,
        logp,
        state.index_in_trajectory + int(np.sign(epsilon)),
    )