Water RHF with native Gaussian CPCM¶
Coordinates and custom radii use bohr. The iSwiG default uses 110 Lebedev nodes per atom, CPCM radii and a switching cutoff of 1e-10. This example uses 26 nodes for the matched benchmark.
from vibeqc import Atom, Molecule, SolvationOptions, run_job
from vibeqc.solvation import ISwiGCavity, MoistNativeOptions, NativePCM
water = Molecule([Atom(8, (0, 0, 0)), Atom(1, (0, 1.498, -1.159)),
Atom(1, (0, -1.498, -1.159))], 0, 1)
model = MoistNativeOptions(cavity=ISwiGCavity(nleb=26),
components=(NativePCM("cpcm", 78.39),))
result = run_job(water, basis="sto-3g", method="rhf", density_fit=False,
solvation=SolvationOptions("moist-native", model=model))
print(result.energy, result.solvent_result.e_solv)
Run the matched example once in each environment to compare the two backends.
Capability |
Implementation |
Numerical qualification |
|---|---|---|
iSwiG and Gaussian CPCM |
Dense solve, fitted positive-weight grids through 974 nodes |
Water RHF/STO-3G, 26 nodes per atom, matched through the AO operator and final energy |
COSMO |
Same matrix with COSMO dielectric scaling |
Water RHF and the conductor limit on H2+ |
PV |
Switched divergence-theorem volume times pressure in atomic units |
Water RHF, H2 DROP, and signed pressure on H2+ |
IEF-PCM, SS(V)PE, D-PCM and its original row scaling |
Dense Gaussian boundary operators and transpose solves |
Water RHF on iSwiG and H2 RHF on both level-2 DROP surfaces |
LPB and nonlocal PCM |
Gaussian Yukawa layers and coupled dielectric equations |
Water RHF on iSwiG and H2 RHF on both level-2 DROP surfaces |
SLIC |
Dense Newton solve with relative residual and Armijo search |
Water RHF on iSwiG and H2 RHF on both level-2 DROP surfaces |
GOSTSHYP |
Normalized pressure Gaussians with signed amplitude and regularization derivatives |
Water RHF on iSwiG and H2 RHF on both level-2 DROP surfaces |
SvdW-DROP and isodensity-DROP |
Single-branch level 2 projection, mapped weights, normals and curvature |
H2 RHF with every listed component and water RHF with CPCM |
RHF, UHF, RKS, UKS and frozen-reference correlation |
Shared molecular host |
CPCM water RHF, PBE/RKS, MP2 and CCSD, and H2+ UHF and PBE/UKS. OH remains unresolved |
The surface follows MOIST’s cavity/iswig.f90, cavity/bornzeta.f90 and
data/radii.f90. Each unpruned area is the solid-angle weight times the
squared sphere radius. The pair switch uses complementary error functions.
A positive area cutoff replaces the switching cutoff. The conductor kernel
in pcm/kernel/gauss/conductor.f90 uses a switched self term and unit-charge
Gaussians. The native off-diagonal uses SciPy’s error function and its
coincident-point limit in place of upstream polynomial branches.
See Lange and Herbert, sections II and III,
doi:10.1063/1.3511297.
Cavity and component settings¶
Use DROPCavity(lsf=SvdW(), nleb=194, projection_level=2) for a smooth
geometrical cavity. Use Isodensity(rho_iso=0.001, scale=1000.) as the
level set to rebuild the surface at every density. Both use single-branch
closest-point projection with SLSQP and Newton refinement. The critical
gradient and focal switches enter the mapped quadrature weights. The
Gaussian widths use those weights before reference switching. Native DROP
uses level 2. Upstream defaults to conditional multistart at level 3.
Levels 3 through 9, ECP core filling and marching-tetrahedra grids are
outside this implementation.
The isodensity Fock combines the fixed-surface AO operator with central
differences of the cavity contribution for each independent density element.
The physical density stays fixed while the two perturbed surfaces move.
The configurable response_step defaults to 1e-5. This is slower than
the upstream analytic surface adjoint and requires stable point membership.
Nuclear gradients use reconverged central differences of the total energy.
Pass a tuple of additive components to MoistNativeOptions. Dielectrics
are dimensionless. NativePCM accepts cpcm, cosmo, iefpcm, ssvpe,
dpcm and dpcmoriginal. NativeLPB(epsilon, kappa) takes the inverse
Debye length in inverse bohr. NativeNonlocal(epsilon, epsilon_inf, length_bohr, epsilon_solute=1.) takes the solvent correlation length.
NativeSLIC(epsilon, alpha, beta, gamma) also accepts mu, phi_static,
newton_tol and newton_maxiter. Its net-charge term uses the molecular
charge. Concave seams of an iSwiG surface can give multiple SLIC roots.
Dense real Gaussian surface operators are implemented. Point kernels,
matrix-free and iterative linear solvers, complex dielectric response,
excited states and Hessians are outside the native scope.
NativePV(pressure) and NativeGOSTSHYP(pressure) use hartree per cubic
bohr. GOSTSHYP accepts the trace regularization interval and the upstream
negative-amplitude suppression option. Its Fock differentiates both the
Gaussian overlap and the amplitude set by the normal force.
Source and equation crosswalk¶
The port uses the pinned MOIST source.
The kernel filenames below are relative to
model/continuum/component/pcm/kernel/gauss there. The Python module names
are relative to vibeqc.solvation.
Capability |
Upstream source |
Governing equations and native code |
|---|---|---|
Radii, Gaussian width and iSwiG |
|
|
SvdW-DROP |
|
|
Isodensity-DROP |
|
|
CPCM and COSMO |
|
|
IEF-PCM and SS(V)PE |
|
|
D-PCM |
|
|
LPB |
|
|
SLIC |
|
|
Nonlocal PCM |
|
|
PV |
|
|
GOSTSHYP |
|
|
The DROP SI files and the full original D-PCM and Juffer LPB papers were unavailable in the library collection used for this port. The numerical conventions for those cases come from the pinned implementation. The DROP comparisons use explicit molecular inputs, independently of the missing SI geometry fixtures.
Both cavity recipes accept radius_model with the upstream cpcm, smd,
d3, cosmo, bondi, rahm and gauss tables. Alternatively supply
radii_bohr in molecular atom order. The tables retain upstream placeholder
values and element limits. Missing Bondi entries raise an error.
Matched numerical comparisons¶
At the geometry in the example, native and external RHF/STO-3G water gave -74.97093711261074 hartree. Of the initial 78 nodes, both retained 70. Their first difference was the Lebedev ordering at the twelfth retained node. Matching nodes by location gave maximum coordinate and switched-area errors of 1.34e-15 bohr and 5.33e-15 bohr squared. At the common AO density D = I, the maximum potential difference was 2.23e-15 hartree per unit charge, the largest charge difference was 8.33e-16 elementary charge, the environment energy differed by 2.23e-15 hartree and the AO reaction operator by at most 1.78e-15 hartree.
Across the eleven model variants on water iSwiG and both H2 DROP cavities, the largest absolute total-energy difference was 8.53e-14 hartree. The largest native directional density derivative error was 2.70e-11 hartree per unit displacement. These checks cover the explicitly listed small grids and inputs. Native nuclear gradients are numerical. On H2 at 1.4 bohr the native CPCM gradient differed by 5.38e-8 hartree/bohr from an independently reconverged bond displacement at a different step size.
The existing OH doublet qualification remains open. Its 26-node UHF runs differed by 3.90e-6 hartree between hosts and the tight UKS/PBE run exhausted 100 outer cycles. The open-shell comparisons in the table use H2+.