Psi4

Psi4: run through a small driver script and decode the JSON it writes.

The codec writes job.json and a driver script (input.py, from _psi4_driver.py) and runs psi4 input.py -k -o stdout: Psi4 executes the script with its own interpreter, so the caller's Python environment never imports Psi4. The script calls the Psi4 driver for the calculation type and writes results.json; stdout carries Psi4's normal output.

Inputs map to Psi4 as follows:

• model.method and model.basis form the method string, e.g. "wb97x-d/def2-svp" for psi4.energy("wb97x-d/def2-svp"). Leave basis empty for composite methods or to set basis in keywords.

• keywords go to psi4.set_options unchanged, so every Psi4 option and its module-scoped form ("scf__maxiter") works as in Psi4's documentation. Two kinds of keyword are not Psi4 options and are routed instead:

• Driver arguments go to psi4.energy()/gradient()/frequency()/ optimize(): external_potentials (point charges for QM/MM electrostatic embedding, [[charge, [x, y, z]], ...] in bohr), dertype (0 or "energy" for derivatives by finite differences of energies, 1 or "gradient" for Hessians from finite differences of gradients), and project_trans/project_rot (whether frequency analysis removes translations and rotations; Psi4 decides by default).

• symmetry is the molecule block's point group, e.g. "c1" to turn symmetry off. By default Psi4 detects whatever symmetry the input frame already has.

• cmdline_args are appended to the psi4 command, e.g. ["-n", "4"] for four threads or ["--memory", "8GB"].

The molecule keeps the input's frame (no_com, no_reorient): Psi4 never moves or rotates it, so gradients, Hessians, normal modes, and trajectories are in the input's coordinates, and symmetry detection only finds symmetry the frame already has. No input field carries atom indices.

Calculation types: energy, gradient, hessian (psi4.frequency: also frequencies, normal modes, and the gradient), optimization and transition_state (psi4.optimize with optking; every step carries its energy and gradient). A transition-state search sets opt_type ts and, unless given, full_hess_every 0 (an exact Hessian at the first step). Every single point also records Mulliken charges in atomic_charges when the method supports them, the gradient on external point charges in point_charge_gradient, and all scalar Psi4 variables in extras["psi4_variables"].

Psi4 1.11 handles point charges correctly only for energies and gradients. Its Hessians with external_potentials either fail or come back wrong (they disagree with finite differences of its own gradients by about as much as the charges change them), and its optimizations reset the molecule to the starting geometry at every step, so they never move. The codec rejects those combinations; for QM/MM geometry optimizations, run geomeTRIC with Psi4 gradients as its child.

Examples

Each script runs as written. Copy it, delete the options you don't need, and run it with python script.py. They live in the repository's examples/ folder.

Single point: energy, gradient, Hessian

examples/psi4/single_point.py

"""Psi4 single point: energy, gradient, or Hessian.

Requires: ``pip install qccompute`` (it installs qccodec) and Psi4 on PATH, e.g.
``conda install -c conda-forge psi4`` (add ``pyddx`` for ddX implicit solvent).

How the codec talks to Psi4 (https://psicode.org/psi4manual/master/):

- ``model.method`` and ``model.basis`` become Psi4's method string, e.g. "wb97x-d/def2-svp".
  Leave ``basis`` out for composite methods such as "mp2/cc-pv[dt]z".
- ``keywords`` go to ``psi4.set_options`` unchanged, so use Psi4's option names.
  Module-scoped options are written "module__option", e.g. "scf__maxiter".
- A few keywords are not Psi4 options and are routed where Psi4 expects them:
  ``external_potentials`` (point charges), ``dertype``, ``project_trans``,
  ``project_rot`` (driver arguments) and ``symmetry`` (the molecule block).
- The molecule stays in your coordinates: Psi4 is told not to move or reorient it,
  so gradients and Hessians line up with ``structure``.
- Psi4 errors (unknown option, SCF not converged, ...) raise with Psi4's message.
"""

from qccompute import compute
from qcdata import CalcType, Model, ProgramInput, Structure

calctype = CalcType.gradient  # or CalcType.energy, CalcType.hessian

# Coordinates are in bohr (Structure.from_xyz / Structure.open read angstrom files).
structure = Structure(
    symbols=["O", "H", "H"],
    geometry=[[0.0, 0.0, 0.0], [0.524, 1.687, 0.48], [1.146, -0.45, -1.354]],
    charge=0,
    multiplicity=1,
)

prog_input = ProgramInput(
    program="psi4",
    calctype=calctype,
    structure=structure,
    model=Model(method="wb97x-d", basis="def2-svp"),  # e.g. "b3lyp-d3bj"/"def2-tzvp"
    keywords={
        # Any Psi4 option, exactly as in Psi4's documentation:
        # "scf_type": "pk",  # exact integrals instead of density fitting
        "e_convergence": 1e-8,
        "d_convergence": 1e-8,
        "scf__maxiter": 100,  # module-scoped option
        # "reference": "uhf",  # open shells: also set structure multiplicity
        #
        # Implicit solvent (ddX; needs the pyddx package next to Psi4):
        "ddx": True,
        "ddx_model": "pcm",  # or "cosmo", "lpb"
        "ddx_solvent": "water",
        #
        # Point charges for QM/MM electrostatic embedding, [charge, [x, y, z]] in
        # bohr. Psi4 1.11 handles them only for energies and gradients; the
        # gradient on the charges is returned as point_charge_gradient.
        # "external_potentials": [[-0.8, [0.0, 0.0, 5.0]], [0.4, [1.0, 1.0, 6.0]]],
        #
        # Symmetry: Psi4 detects any symmetry your frame already has without
        # moving atoms. Name a point group to override, or "c1" to turn it off.
        # "symmetry": "c1",
        #
        # Hessians: frequency analysis projects out translations and rotations
        # only at stationary points unless told otherwise.
        # "project_rot": True,
        # "dertype": 1,  # Hessian by finite differences of analytic gradients
    },
    # Arguments for the psi4 command itself: threads and memory.
    cmdline_args=["-n", "1", "--memory", "2GB"],
)

output = compute(
    prog_input,
    # timeout=3600,  # seconds before Psi4 is killed
    # collect_files=True,  # keep every output file in output.results.files
    # executables={"psi4": "/path/to/psi4"},  # or set QCCOMPUTE_PSI4_EXECUTABLE
)

results = output.results
print("Energy (Hartree):", results.energy)
if results.gradient is not None:
    print("Gradient (Hartree/bohr):\n", results.gradient)
if results.hessian is not None:
    print("Frequencies (cm^-1):", results.freqs_wavenumber)
print("Mulliken charges:", results.atomic_charges)
print("Dipole (e*bohr):", results.scf_dipole_moment)
print("Psi4 variables:", sorted(results.extras["psi4_variables"])[:5], "...")

Optimization

examples/psi4/optimization.py

"""Psi4 geometry optimization with its built-in optimizer, optking.

Requires: ``pip install qccompute`` (it installs qccodec) and Psi4 on PATH, e.g.
``conda install -c conda-forge psi4``.

Notes for Psi4 optimizations:

- optking options are ordinary Psi4 options: ``g_convergence``, ``geom_maxiter``,
  ``step_type``, ``opt_coordinates``, ``full_hess_every``, ... (Psi4 manual, "Geometry
  Optimization").
- Every step in ``results.trajectory`` carries its structure, energy, and gradient;
  the last step also carries the final charges, orbitals, and Psi4 variables.
- The first step is exactly your input structure: Psi4 does not move or reorient it.
- Point charges (``external_potentials``) do not work in Psi4 1.11 optimizations
  (it never leaves the starting geometry), so the codec rejects them. For QM/MM
  optimizations use geomeTRIC with Psi4 gradients; see examples/geometric/.
- For transition states use ``CalcType.transition_state`` (transition_state.py).
"""

from qccompute import compute
from qcdata import CalcType, Model, ProgramInput, Structure

# Coordinates are in bohr.
structure = Structure(
    symbols=["O", "H", "H"],
    geometry=[[0.0, 0.0, 0.0], [0.524, 1.687, 0.48], [1.146, -0.45, -1.354]],
    charge=0,
    multiplicity=1,
)

prog_input = ProgramInput(
    program="psi4",
    calctype=CalcType.optimization,
    structure=structure,
    model=Model(method="wb97x-d", basis="def2-svp"),
    keywords={
        # optking settings:
        "g_convergence": "gau",  # or "qchem", "gau_tight", "gau_verytight", ...
        "geom_maxiter": 50,
        # "opt_coordinates": "both",  # redundant internals plus Cartesians
        # "full_hess_every": 0,  # exact Hessian at the first step
        #
        # Implicit solvent (ddX; needs pyddx next to Psi4):
        # "ddx": True,
        # "ddx_model": "pcm",
        # "ddx_solvent": "water",
        # "symmetry": "c1",  # keep an asymmetric path open from a symmetric start
    },
    cmdline_args=["-n", "1", "--memory", "2GB"],
)

output = compute(prog_input)  # timeout=..., collect_files=..., executables=...

results = output.results
print("Steps:", len(results.trajectory))
print("Energies (Hartree):", results.energies)
print("Final structure (bohr):\n", results.final_structure.geometry)
# Save the trajectory, e.g. to view in a molecule viewer:
# results.save("optimization.json"); results.final_structure.save("final.xyz")

Transition state

examples/psi4/transition_state.py

"""Psi4 transition-state search with optking.

Requires: ``pip install qccompute`` (it installs qccodec) and Psi4 on PATH, e.g.
``conda install -c conda-forge psi4``.

Notes for Psi4 transition-state searches:

- The codec sets ``opt_type ts`` and, unless you set it, ``full_hess_every 0``: an
  exact Hessian at the first step, which optking needs to find a saddle point.
  Psi4 has analytic Hessians only for Hartree-Fock, so for DFT it takes finite
  differences of gradients, which it runs by calling ``psi4`` again: keep ``psi4``
  on PATH (not only in ``executables``).
  Setting ``opt_type`` to anything else is an error for this calculation type.
- Start close to the saddle point. This example starts from nearly planar ammonia;
  its inversion transition state is planar, with one imaginary frequency.
- Check the result with a Hessian (single_point.py with CalcType.hessian): a
  transition state has exactly one imaginary (negative) frequency.
- Point charges are not supported here in Psi4 1.11; use geomeTRIC with Psi4
  gradients for QM/MM transition states.
"""

from qccompute import compute
from qcdata import CalcType, Model, ProgramInput, Structure

# Nearly planar ammonia, coordinates in bohr.
structure = Structure(
    symbols=["N", "H", "H", "H"],
    geometry=[
        [0.0, 0.0, 0.25],
        [1.77, 0.0, -0.10],
        [-0.885, 1.533, -0.10],
        [-0.885, -1.533, -0.10],
    ],
)

prog_input = ProgramInput(
    program="psi4",
    calctype=CalcType.transition_state,
    structure=structure,
    model=Model(method="wb97x-d", basis="def2-svp"),
    keywords={
        "geom_maxiter": 50,
        # "full_hess_every": 3,  # recompute the exact Hessian every 3 steps
        # "step_type": "rs_i_rfo",  # optking's default for transition states
    },
    cmdline_args=["-n", "1", "--memory", "2GB"],
)

output = compute(prog_input)

results = output.results
print("Steps:", len(results.trajectory))
print("Transition-state energy (Hartree):", results.energies[-1])
print("Transition-state structure (bohr):\n", results.final_structure.geometry)