vibe-qc/Documentation

Density-dependent cavities and pressure with MOIST

Install MOIST’s Python/native API first. A cavity recipe and a list of model components determine the environment; the electronic method uses the same coupling regardless of that choice. This integration is experimental pending the matched scientific benchmark campaign.

import moist
import vibeqc as vq

recipe = vq.MoistOptions(
    cavity=moist.DROP(lsf=moist.Isodensity()),
    components=[moist.ModelComponentCPCM(78.39)],
    tol_density=1e-8,
    tol_potential=1e-8,
)
result = vq.run_job(mol, method="rks", functional="pbe0", basis="def2-svp",
    solvation=vq.SolvationOptions("moist", model=recipe), output="density_cavity")

Here mol is a molecular vibeqc.Molecule, for example the water geometry in the backend tutorial. MOIST receives the actual Gaussian basis and spin-summed density. At every density update it rebuilds the isodensity cavity. The returned density weights contribute to the Fock matrix, alongside the potential adjoints. Omitting those weights would differentiate a frozen cavity and solve a different electronic problem.

For an open-shell molecule, keep the multiplicity and reference explicit:

oh = vq.Molecule([vq.Atom(8, (0., 0., 0.)),
                  vq.Atom(1, (0., 0., 1.83))], 0, 2)
radical = vq.run_job(oh, method="uhf", basis="sto-3g",
    solvation=vq.SolvationOptions("moist", "water"), output="oh_water")

Both spin blocks receive the environment operator built from their summed density. method="uks", functional="pbe" selects the analogous DFT route. ROHF and excited-state solvent response are separate capabilities and are not supplied by this interface. For nearly degenerate radicals, compare electronic roots and converge the cavity angular grid: a converged SCF solution need not be the lowest root.

For pressure, supply pressure explicitly in Hartree/bohr cubed:

gpa_to_au = 1e9 * (5.29177210903e-11)**3 / 4.3597447222071e-18
pressure = 1.0 * gpa_to_au
recipe = vq.MoistOptions(
    solvent="vacuum",
    cavity=moist.DROP(lsf=moist.SvdW()),
    components=[moist.ModelComponentGOSTSHYP(pressure)],
)

ModelComponentPV(pressure) instead adds the pressure-volume component. On a geometry-only cavity PV has no density derivative; on a density-following cavity it does. GOSTSHYP uses Gaussian density moments and their amplitude response. These are different pressure models and should be studied separately before composing them with a PCM component.

The complete nuclear derivative currently uses energy finite differences. Converge the electronic and cavity tolerances more tightly than the desired gradient accuracy and compare at two displacement sizes. Pressure studies also need cavity and basis convergence, reference geometries and explicit pressure units. The executable example is examples/molecular/solvation/input-moist-pressure.py.