PCM analytical IR dipole disagrees with finite-difference PCM dipoles

Using GPU4PySCF 1.8.0 with PySCF 2.13.1, the analytical dipole derivative used by gpu4pyscf.properties.ir.eval_ir_freq_intensity disagrees substantially with central finite differences of PCM dipole moments.

Minimal system:

H  0 0 0
O  0 0 1
H  0 1 1

Settings: B3LYP/STO-3G, IEF-PCM with ε = 2.3741.

The analytical and numerical Cartesian dipole derivatives were projected through exactly the same externally computed normal modes from the same analytical Hessian to avoid differences from mode ordering, normalization, or translation/rotation projection.

The finite-difference step was 0.005 Å. A gas-phase control agrees within 0.4%, which suggests to me that the discrepancy is specific to solvent-response contributions to the analytical dipole derivative (or similar).

Calculations run through Modal on an A10 GPU.

Here's some minimal scripts to show the setup, although these don't do the mode-projection work mentioned above so it's not quite a 1:1 comparison. (Disclaimer: Codex was used to write these scripts.)

"""Minimal analytical PCM IR reproducer for GPU4PySCF 1.8.0."""

import numpy as np
from gpu4pyscf.dft import rks
from gpu4pyscf.properties import ir
from pyscf import gto

ATOM = """
H  0.0  0.0  0.0
O  0.0  0.0  1.0
H  0.0  1.0  1.0
"""


def asnumpy(value):
    get = getattr(value, "get", None)
    return np.asarray(get() if callable(get) else value)


mol = gto.M(atom=ATOM, basis="sto-3g", unit="Angstrom", verbose=0)
mf = rks.RKS(mol, xc="B3LYP").density_fit().PCM()
mf.with_solvent.method = "IEF-PCM"
mf.with_solvent.eps = 2.3741  # toluene
mf.grids.atom_grid = (99, 590)
mf.level_shift = 0.1
mf.conv_tol = 1e-10
mf.conv_tol_cpscf = 1e-7
assert np.isfinite(mf.kernel()) and mf.converged

frequency, intensity = ir.eval_ir_freq_intensity(mf, mf.Hessian())
for freq, inten in zip(asnumpy(frequency), asnumpy(intensity), strict=True):
    print(f"{freq:12.6f} cm^-1  {inten:12.6f} km/mol")

"""Finite-difference PCM dipole reference for GPU4PySCF 1.8.0."""

import numpy as np
from gpu4pyscf.dft import rks
from pyscf import gto
from pyscf.hessian import thermo

SYMBOLS = ("H", "O", "H")
COORDS = np.array(
    [
        [0.0, 0.0, 0.0],
        [0.0, 0.0, 1.0],
        [0.0, 1.0, 1.0],
    ]
)
STEP_ANGSTROM = 0.005
IR_CONVERSION = 42.255  # (D / Angstrom / sqrt(amu))^2 -> km/mol


def asnumpy(value):
    get = getattr(value, "get", None)
    return np.asarray(get() if callable(get) else value)


def run_scf(coords):
    mol = gto.M(
        atom=list(zip(SYMBOLS, coords, strict=True)),
        basis="sto-3g",
        unit="Angstrom",
        verbose=0,
    )
    mf = rks.RKS(mol, xc="B3LYP").density_fit().PCM()
    mf.with_solvent.method = "IEF-PCM"
    mf.with_solvent.eps = 2.3741  # toluene
    mf.grids.atom_grid = (99, 590)
    mf.level_shift = 0.1
    mf.conv_tol = 1e-10
    mf.conv_tol_cpscf = 1e-7
    assert np.isfinite(mf.kernel()) and mf.converged
    return mf


reference = run_scf(COORDS)
raw_hessian = asnumpy(reference.Hessian().kernel())
analysis = thermo.harmonic_analysis(
    reference.mol,
    raw_hessian,
    imaginary_freq=False,
)

# d(mu_Debye) / d(x_Angstrom), axes (atom, coordinate, dipole component)
dipole_gradient = np.empty((len(SYMBOLS), 3, 3))
for atom in range(len(SYMBOLS)):
    for axis in range(3):
        plus = COORDS.copy()
        minus = COORDS.copy()
        plus[atom, axis] += STEP_ANGSTROM
        minus[atom, axis] -= STEP_ANGSTROM
        mu_plus = asnumpy(run_scf(plus).dip_moment(verbose=0))
        mu_minus = asnumpy(run_scf(minus).dip_moment(verbose=0))
        dipole_gradient[atom, axis] = (mu_plus - mu_minus) / (2 * STEP_ANGSTROM)

# PySCF norm_mode axes are (mode, atom, coordinate), mass-normalized.
dmu_dq = np.einsum("mac,acd->md", analysis["norm_mode"], dipole_gradient)
intensity = IR_CONVERSION * np.einsum("md,md->m", dmu_dq, dmu_dq)

for freq, inten in zip(analysis["freq_wavenumber"], intensity, strict=True):
    print(f"{freq:12.6f} cm^-1  {inten:12.6f} km/mol")