def METHOD_NAME(theta, psi0):
    def energy1(theta, psi0):
        w = variational_wfn(theta, psi0)
        wl = tc.backend.conj(w)
        wr = tc.backend.stop_gradient(w)
        wl = tc.backend.reshape(wl, [1, -1])
        wr = tc.backend.reshape(wr, [-1, 1])
        e = wl @ h @ wr  # <\partial psi0|H| psi0>
        return tc.backend.real(e)[0, 0]
    def energy2(theta, psi0):
        w = variational_wfn(theta, psi0)
        wr0 = tc.backend.stop_gradient(w)
        wr0 = tc.backend.reshape(wr0, [-1, 1])
        wl0 = tc.backend.stop_gradient(w)
        wl0 = tc.backend.conj(wl0)
        wl0 = tc.backend.reshape(wl0, [1, -1])
        e0 = wl0 @ h @ wr0  # <psi0| H | psi0>
        wl = tc.backend.conj(w)
        wl = tc.backend.reshape(wl, [1, -1])
        w0 = wl @ wr0  # <\partial psi0| psi0>
        return tc.backend.real((w0 * e0)[0, 0])
    eg1 = tc.backend.grad(energy1, argnums=0)
    eg2 = tc.backend.grad(energy2, argnums=0)
    rhs1 = eg1(theta, psi0)
    rhs1 = tc.backend.imag(rhs1)
    rhs2 = eg2(theta, psi0)
    rhs2 = tc.backend.imag(rhs2)  # should be a imaginary number
    rhs = rhs1 - rhs2
    return rhs