import numpy as np

def sto_1s(r, zeta):
    # 規格化定数（1s STO）
    N = (zeta**3/np.pi)**0.5
    return N * np.exp(-zeta * r)

def overlap_mc(zeta, R, Nsample=100000, L=5.0):
    """
    モンテカルロ法で異中心1s軌道のオーバーラップを近似
    核A=(0,0,0), 核B=(0,0,R)
    """
    # サンプル点を立方体 [-L,L]^3 に一様分布で生成
    pts = np.random.uniform(-L, L, size=(Nsample, 3))
    # 核Aからの距離
    rA = np.linalg.norm(pts, axis=1)
    # 核Bからの距離
    rB = np.linalg.norm(pts - np.array([0,0,R]), axis=1)
    # 軌道値
    phiA = sto_1s(rA, zeta)
    phiB = sto_1s(rB, zeta)
    f = phiA * phiB
    # 平均値 × 体積
    volume = (2*L)**3
    return volume * np.mean(f)

# 例: H2, ζ=1.0, R=1.4 Bohr
zeta = 1.0
R = 1.4
S_mc = overlap_mc(zeta, R, Nsample=200000)
S_ana = np.exp(-zeta*R)*(1+zeta*R+(zeta*R)**2/3)

print("Monte Carlo overlap =", S_mc)
print("Analytic overlap    =", S_ana)