PCM analytical IR dipole disagrees with finite-difference PCM dipoles

[corinwagen]

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.

Frequency (cm⁻¹)	Analytical	Finite difference	Relative difference
2068.5596	0.787536	0.594638	+32.44%
3974.7373	9.690058	11.701083	−17.19%
4109.9406	5.679714	8.989973	−36.82%

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")

[puzhichen]

Thanks for the detailed report and the minimal reproducer! I've looked into this and can confirm the root cause.

The analytical dipole derivative in gpu4pyscf/properties/ir.py currently only computes the gas-phase contributions. It does not include any PCM solvent-response contributions . So the analytical IR intensities are effectively gas-phase even when a PCM solvent is attached to the mean-field object.

This is less a bug and more a missing feature — PCM analytical IR intensities haven't been implemented yet. I'm going to relabel this issue accordingly and work on adding PCM support. In the meantime, finite-difference of mf.dip_moment() (as in your reference script) is the correct workaround.

[corinwagen]

Ah, makes sense... thanks!

[michalkrompiec]

Could the code raise an error for now, until the PCM part is implemented?

[puzhichen]

Thanks for the suggestion! Yes, we'll make sure the code raises a clear error, until the PCM part is fully implemented.