psi4.in

molecule {
O 0.000000000 0.000000000 0.369372944
H -0.783975899 0.000000000 -0.184666472
H 0.783975899 0.000000000 -0.184686472
}

set {
    basis sto-3g
    scf_type direct
    df_scf_guess false
}

set optking {
    print_opt_params true
}

optimize("hf")

pyberny_psi4_example.py

import logging
import subprocess as sp
from pathlib import Path
from typing import Any, Callable, Generator, Optional, Sequence, Tuple, Union

import psi4
import numpy as np
from berny import Berny, geomlib
from berny.coords import angstrom

Atoms = Sequence[Tuple[str, np.ndarray]]
GradRet = Tuple[float, np.ndarray]


def atoms_to_str(atoms: Atoms) -> str:
    """Convert the atoms received by a PyBerny Solver into a Psi4 geometry
    string.
    """
    return "\n".join(
        (f"{symbol} {coords[0]} {coords[1]} {coords[2]}" for symbol, coords in atoms)
    )


def Psi4GradientWrapper(
    atoms: Atoms,
    lattice: Optional[Any] = None,
    method: str = "hf",
    basis: str = "sto-3g",
) -> GradRet:
    """Calculate the energy and gradient using Psi4."""
    psi4.core.clean()
    psi4.set_memory("500 MB")
    psi4.core.set_output_file("output.dat", True)
    _ = psi4.geometry(atoms_to_str(atoms))
    psi4.set_options({"basis": basis, "scf_type": "direct", "df_scf_guess": False})
    ret = psi4.gradient(method, return_wfn=True)
    energy = ret[1].energy()
    gradient = ret[0].np
    return energy, gradient


def GenericGradientSolver(
    f: Callable[..., GradRet], *args, **kwargs
) -> Generator[GradRet, GradRet, None]:
    """Given a function that ..."""
    atoms, lattice = yield 0.0, np.asanyarray([])
    while True:
        energy, gradients = f(atoms, *args, **kwargs)
        atoms, lattice = yield energy, gradients / angstrom


def optimize(optimizer, solver, trajectory=None):
    """Optimize a geometry with respect to a solver.

    Args:
        optimizer (:class:`~collections.abc.Generator`): Optimizer object with
            the same generator interface as :class:`~berny.Berny`
        solver (:class:`~collections.abc.Generator`): unprimed generator that
            receives geometry as a 2-tuple of a list of 2-tuples of the atom
            symbol and coordinate (as a 3-tuple), and of a list of lattice
            vectors (or :data:`None` if molecule), and yields the energy and
            gradients (as a :math:`N`-by-3 matrix or :math:`(N+3)`-by-3 matrix
            in case of a crystal geometry).

            See :class:`~berny.solvers.MopacSolver` for an example.
        trajectory (str): filename for the XYZ trajectory

    Returns:
        The optimized geometry.

    The function is equivalent to::

        next(solver)
        for geom in optimizer:
            energy, gradients = solver.send((list(geom), geom.lattice))
            optimizer.send((energy, gradients))
    """
    if trajectory:
        trajectory = open(trajectory, "w")
    try:
        next(solver)
        for geom in optimizer:
            energy, gradients = solver.send((list(geom), geom.lattice))
            if trajectory:
                geom.dump(trajectory, "xyz")
            optimizer._log.debug(f"gradient:\n{gradients * angstrom}")
            debug = optimizer.send((energy, gradients))
    finally:
        if trajectory:
            trajectory.close()
    return geom


def optimize_xyz(path_to_xyz: Union[str, Path]) -> str:
    geom = geomlib.readfile(str(path_to_xyz))
    berny = Berny(geom, debug=True, trust=0.5)
    solver = GenericGradientSolver(Psi4GradientWrapper)
    final_geom = optimize(berny, solver)
    return final_geom.dumps("xyz")


if __name__ == "__main__":
    logging.basicConfig(level=logging.DEBUG)
    sp.call(["rm", "output.dat"])
    xyz = optimize_xyz("water_start.xyz")

water_start.xyz

3

O 0.000000000 0.000000000 0.369372944
H -0.783975899 0.000000000 -0.184666472
H 0.783975899 0.000000000 -0.184686472