# Copyright (C) 2024 The authors of Py-ChemShell
#
# This file is part of Py-ChemShell.
#
# Py-ChemShell is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as
# published by the Free Software Foundation, either version 3 of the
# License, or (at your option) any later version.
#
# Py-ChemShell is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public
# License along with Py-ChemShell. If not, see
# <http://www.gnu.org/licenses/>.
__author__ = "Joseph Thacker <joseph.thacker@stfc.ac.uk>"
import os
import subprocess
import shutil
from typing import Tuple
import numpy as np
import numpy.typing
from ..qm import _QM
from ...base import errors as chemsh_errors
from ...utils import testutils as chemsh_testutils
from ...objects.fragment import Fragment
MINIMUM_PSI4_VERSION = "1.9.1"
[docs]
class Psi4(_QM):
"""Interface to Psi4
Attributes
----------
frag : Fragment, default: Fragment()
The geometry and physical properties of the system to study
method : str, default: "rhf"
Defines the level of theory must be one of:
`["HF", "DFT", "UDFT", "RHF", "ROHF", "UHF", "CUHF", "RKS", "UKS"]`
basis : str, default: "3-21G"
Defines the basis set
charge : int, default: 0
Charge of the system
mult : int, default: 1
Multiplicity of the system
input_file_name: str, default: "_psi4.inp"
Input file name that ChemShell will write and overwrite
output_file_name: str, default: "_psi4.out"
Output file name that ChemShell will read from and will be
overwritten by Psi4
convergence: float, default: 1.0e-8
Defines the energy convergence parameter for the SCF cycle
functional: str, default: ""
Available only when using Density Functional Theory
maxiter: int, default: 30
Maximum number of SCF cycles
use_ri: bool, default: False
Whether to use Resolution of the Identity (RI)
nprocs: int, default: 1
Number of threads that PSi4 will use
syscmd: str, default: "psi4"
Defines the Psi4 executable (must be in path)
fcidump: bool, default: False
Should Psi4 output an FCIDUMP file?
set_options: dict, default: dict()
Additional options to be set in Psi4. They must take the form:
`set key value`
Notes
-----
This interface uses the Psi4 input files and executable,
not the Python API. The minimum required version of Psi4
is 1.9.1.
Examples
--------
.. highlight:: python
.. code-block:: python
#
# Psi4 HF single point energy
#
from chemsh import Fragment, Psi4, SP
# Define a fragment using a json file
water = Fragment(coords="water.cjson")
# set up a qm task
qm_task = Psi4(frag=water, method="hf", basis="3-21g")
# define a single point calculation
sp = SP(theory=qm_task, gradients=True)
# run the single point calculation
sp.run()
# print the energy & gradients
print(sp.result.energy)
print(sp.result.gradients)
.. highlight:: python
.. code-block:: python
#
# Psi4/Gulp QM/MM Water Optimisation, with electrostatic
# embedding
# see "chemsh-py/tests/qmmm/psi4-gulp/h2o_dimer_elec_opt.py"
#
from chemsh import Fragment, Psi4, GULP, QMMM, Opt
# Define the water dimer from json
water = Fragment(coords='h2o_dimer.cjson')
# Define the GULP FF parameters
ff = '''#
# molmec removes Coulomb terms between bonded atoms
# and between atoms bonded to a common third atom
keyword molmec
# Flexible TIP3P model
# last number in LJ term is a cutoff
lennard 12 6 kcal x13
O O 582000.0 595.0 99.9
# intramolecular force-field parameters for MM
# from QM calculation on water monomer
# Note GULP k is twice that of ChemShell format
harmonic bond kcal
H O 1059.162 0.96
# Note GULP k is twice that of ChemShell format
# note different ordering of atoms
three bond kcal
O H H 68.087 104.5
'''
# Define the QM region calculation
qm = Psi4(method="hf", basis="3-21g")
# Define the MM region calculation
mm = GULP(ff=ff, molecule=True)
# Define the QM/MM calculation
qmmm = QMMM(frag=water, qm_region=range(3), qm=qm, mm=mm,
embedding='electrostatic', coupling='covalent')
# Define and run the optimisation
optimisation = Opt(theory=qmmm)
optimisation.run()
"""
[docs]
def __init__(
self,
frag: Fragment = Fragment(),
method: str = "rhf",
basis: str = "3-21G",
charge: int = 0,
mult: int = 1,
input_file_name: str = "_psi4.inp",
output_file_name: str = "_psi4.out",
convergence: float = 1.0e-8,
functional: str = "",
maxiter: int = 30,
use_ri: bool = False,
nprocs: int = 1,
syscmd: str = "psi4",
fcidump: bool = False,
set_options: dict = {},
):
""" """
# START: Hack to allow for a proper initialiser
# Lets not use kwargs and instead use "locals" to grab all the input
# (local arguments at this point)
_syscmd = syscmd
kwargs = locals().copy()
kwargs["input"] = input_file_name
# Remove following from the dict to prevent
# ChemShell raising an error about kwargs
kwargs.pop("self")
kwargs.pop("__class__")
kwargs.pop("input_file_name")
kwargs.pop("output_file_name")
kwargs.pop("convergence")
kwargs.pop("mult")
kwargs.pop("use_ri")
kwargs.pop("nprocs")
kwargs.pop("syscmd")
kwargs.pop("set_options")
kwargs.pop("fcidump")
# END: Hack to allow for a proper initialiser
# Init Base Class
super().__init__(**kwargs)
# Convert method from ChemShell to Psi4 format
method_conversions = {"HF": "RHF", "DFT": "RKS", "UDFT": "UKS"}
if method.upper() in method_conversions.keys():
method = method_conversions[method.upper()]
# Perform checks on input args
allowed_methods = ["RHF", "ROHF", "UHF", "CUHF", "RKS", "UKS"]
if method.upper() not in allowed_methods:
raise RuntimeError(
f"User input 'method=\"{method.upper()}\"' is invalid.\n"
"Must be one of the following:\n"
+ ("\n".join(allowed_methods))
)
# If doing DFT, we set the default to "b-lyp" unless other functional
# is provided
if "KS" in method.upper():
functional = functional if functional != "" else "b-lyp"
else:
if functional != "":
raise RuntimeError(
"User asked for DFT functional when not performing DFT.\n"
f"Method: '{method}'\n"
f"Functional: '{functional}'"
)
# Public Member Variables
self.frag = frag
self.method = method
self.basis = basis
self.charge = charge
self.mult = mult
self.input = input_file_name
self.output = output_file_name
self.convergence = convergence
self.functional = functional
self.maxiter = maxiter
self.use_ri = use_ri
self.nprocs = nprocs
self.set_options = set_options
self.fcidump = fcidump
# Private Member Variables
self._syscmd = syscmd
self._file_string = ""
self._nuc_grads_fname = self.output + ".nuc_grads.npy"
self._extern_grads_fname = self.output + ".extern_grads.npy"
self._fcidump_fname = self.output + ".fcidump"
# Check Psi4 executable exists and is executable
psi4_available = shutil.which(self._syscmd) is not None
if not psi4_available:
chemsh_testutils.setUnsupported()
exit(chemsh_errors.UNSUPPORTED)
# Check for a minimum version of Psi4
self.version = self._get_psi4_version()
if self.version < MINIMUM_PSI4_VERSION:
raise RuntimeError(
f"Psi4 version ({self.version}) is less than"
f"the minimum supported version in ChemShell: {MINIMUM_PSI4_VERSION}"
)
def _get_psi4_version(self) -> bool:
"""
checks that the current Psi4 version is higher than the minimum
allowed version
"""
psi4_subprocess = subprocess.run(
[self._syscmd, "--version"], capture_output=True
)
version = psi4_subprocess.stdout.decode("utf-8")
return version
def _write_banner(self, inp_str: str) -> str:
"""
Returns the ChemShell Banner for the input file
"""
return (
inp_str + "#! This file has been generated automatically "
"by PyChemShell (https://chemshell.org/)\n\n"
)
def _write_molecule(self, inp_str: str) -> str:
"""Writes a string with format:
``
molecule name {
0 1
units au
symmetry c1
noreorient
element x y z
. ...
. ...
. ...
}
``
"""
inp_str += "molecule chemsh_mol {\n"
inp_str += f"{self.charge} {self.mult}\n"
inp_str += "units au\n"
inp_str += "symmetry c1\n"
inp_str += "noreorient\n"
coordstr = self.frag.coords2str(charges=None, unit="a.u.")
inp_str += coordstr
inp_str += "}\n"
return inp_str
def _write_set_options(self, inp_str: str, set_options: dict) -> str:
"""Loops through the dictionary and writes each option as:
``
set key value
``
"""
for key, val in set_options.items():
inp_str += f"set {key} {val}\n"
return inp_str
def _write_external_point_charges(
self, inp_str: str, extra_options: dict
) -> str:
"""
Returns the input required to set external point charges
in Psi4
"""
add_point_charges = bool(self.frag.nbqs)
if add_point_charges:
inp_str += "external_point_charges = np.array([\n"
for i in range(self.frag.nbqs):
inp_str += (
f"[{self.frag.bqs.charges[i]},"
f" {self.frag.bqs.coords[i, 0]},"
f" {self.frag.bqs.coords[i, 1]},"
f" {self.frag.bqs.coords[i, 2]}],\n"
)
inp_str += "]);\n"
extra_options["external_potentials"] = "external_point_charges"
return inp_str, extra_options
def _write_calc_type(
self, inp_str: str, gradients: bool, extra_options: dict
) -> str:
"""
Returns the input required to set the calculation type.
"""
starting_str = (
"_, wfn = energy('scf', return_wfn=True"
if gradients or self.fcidump
else "energy('scf'"
)
inp_str += starting_str
for k, v in extra_options.items():
inp_str += f", {k} = {v}"
inp_str += ")\n"
return inp_str
def _write_gradients(
self, inp_str: str, has_external_charges: bool
) -> str:
"""
Returns the input required to compute gradients for the nuclear
positions and the point charges.
"""
inp_str += "# Compute the gradients of the nuclear positions\n"
inp_str += (
"nuclear_gradients = np.asarray(gradient('scf'," " ref_wfn=wfn))\n"
)
inp_str += "# Write out the nuclear gradients to a numpy file\n"
inp_str += f"nuclear_gradients.tofile('{self._nuc_grads_fname}')\n"
if has_external_charges:
inp_str += "# Get the point charge gradients\n"
inp_str += (
"external_charge_gradients = "
"np.asarray(wfn.external_pot().gradient_on_charges())\n"
)
inp_str += (
"external_charge_gradients"
f".tofile('{self._extern_grads_fname}')\n"
)
return inp_str
def _write_fcidump(self, inp_str: str) -> str:
"""
Returns in input required to write an fcidump file
"""
if self.fcidump:
inp_str += f"fcidump(wfn, '{self._fcidump_fname}')"
return inp_str
def _read_energy_from_output(self) -> float:
"""
Scans the output file for the total energy
Raises
------
RuntimeError if the total energy is not found in the output file
"""
with open(self.output, "r") as output_file:
for line in output_file:
if "Total Energy =" in line:
return float(line.split()[-1])
raise RuntimeError(f"Energy not found in output file: {self.output}")
def _read_gradients_from_output(
self, has_external_charges: bool
) -> numpy.typing.NDArray:
"""
Scans the output file for the total gradients and the point charge
gradients
"""
natoms = self.frag.natoms
gradients = np.zeros((natoms, 3))
gradients = np.fromfile(self._nuc_grads_fname)
self._result.gradients = np.copy(gradients)
if has_external_charges:
gradients = np.fromfile(self._extern_grads_fname)
self.frag.bqs.gradients = np.copy(gradients)
def run(self, **kwargs):
"""
TODO: we should look to not use kwargs
"""
# Get the gradients flag, if it is set in the kwargs
if "_gradients" in kwargs.keys():
self._gradients = kwargs["_gradients"]
compute_grad = self._gradients
add_point_charges = bool(self.frag.nbqs)
# Define a dictionary for additional options for the calculation type
extra_options = {}
if "KS" in self.method.upper():
extra_options["dft_functional"] = "'" + self.functional + "'"
# These are multiple options written as: "set property value"
# Initialise with the any extra options provided by the user
set_options = self.set_options.copy()
# Now set any standard ChemShell options
set_options["scf_type"] = "DF" if self.use_ri else "PK"
set_options["maxiter"] = str(self.maxiter)
set_options["basis"] = self.basis
set_options["reference"] = self.method
set_options["e_convergence"] = self.convergence
# Write the input file string
inp_str = ""
inp_str = self._write_banner(inp_str)
inp_str = self._write_molecule(inp_str)
inp_str = self._write_set_options(inp_str, set_options)
inp_str, extra_options = self._write_external_point_charges(
inp_str, extra_options
)
inp_str = self._write_calc_type(inp_str, compute_grad, extra_options)
if compute_grad:
inp_str = self._write_gradients(inp_str, add_point_charges)
# Write string to the input file
with open(self.input, "w") as input_file:
input_file.write(inp_str)
# Run the Psi4 job
psi4_subprocess = subprocess.run(
[
self._syscmd,
"-i",
self.input,
"-o",
self.output,
"-n",
str(self.nprocs),
],
capture_output=True,
)
# Check to see if Psi4 has finished correctly
# if not then raise an error
if psi4_subprocess.returncode != 0:
raise RuntimeError(
f"Psi4 Failed with return code: {psi4_subprocess.returncode}\n"
"StdErr:" + psi4_subprocess.stderr.decode("utf-8") + "\n"
"StdOut:" + psi4_subprocess.stdout.decode("utf-8")
)
# Read in the energy from the output
self._result.energy = self._read_energy_from_output()
# Read the gradients from output, if requested
if compute_grad:
self._read_gradients_from_output(add_point_charges)