import numpy as np


eVtoHr = 27.2114
Nsample = 2000000

def sto_1s(r, zeta):
    N = (zeta**3/np.pi)**0.5
    return N * np.exp(-zeta * r)

def onsite_energy_mc(zeta, Nsample=100000, L=10.0):
    """
    モンテカルロ法で水素1s軌道のオンサイトエネルギーを近似
    H = -1/r （原子単位系）
    """
    pts = np.random.uniform(-L, L, size=(Nsample, 3))
    r = np.linalg.norm(pts, axis=1)
    phi = sto_1s(r, zeta)
    density = phi**2
    V = -1.0 / r
    f = density * V
    volume = (2*L)**3
    return volume * np.mean(f)

# 例: ζ=1.0 の水素1s
E_mc = onsite_energy_mc(1.0, Nsample=Nsample)
print("Onsite energy (MC approx) =", E_mc * eVtoHr, " eV")