Update schema in Q-Chem calculator (#954)

Closes #942.

---------

Co-authored-by: deepsource-autofix[bot] <62050782+deepsource-autofix[bot]@users.noreply.github.com>

src/quacc/calculators/qchem.py

@@ -10,8 +10,10 @@
 import numpy as np
 from ase import Atoms, units
 from ase.calculators.calculator import FileIOCalculator
+from emmet.core.tasks import _parse_custodian
 from monty.io import zopen
 from pymatgen.io.ase import AseAtomsAdaptor
+from pymatgen.io.qchem.inputs import QCInput
 from pymatgen.io.qchem.outputs import QCOutput
 from pymatgen.io.qchem.sets import QChemDictSet
 
@@ -20,6 +22,8 @@
 logger = logging.getLogger(__name__)
 
 
+# TODO: Would be more consistent to convert the Hessian to ASE units if it's being
+# stored in the first-class results.
 class QChem(FileIOCalculator):
  """
 
@@ -43,10 +47,29 @@ class QChem(FileIOCalculator):
  Returns
  -------
  Atoms
- The ASE Atoms object with attached Q-Chem calculator.
+ The ASE Atoms object with attached Q-Chem calculator. In calc.results,
+ the following properties are set:
+
+ - energy: electronic energy in eV
+ - forces: forces in eV/A
+ - hessian: Hessian in native Q-Chem units
+ - entropy: total enthalpy in eV
+ - enthalpy: total entropy in eV/K
+ - qc_output: Output from `pymatgen.io.qchem.outputs.QCOutput.data`
+ - qc_input: Input from `pymatgen.io.qchem.inputs.QCInput.as_dict()`
+ - custodian: custodian.json file metadata
  """
 
- implemented_properties = ["energy", "forces", "hessian"] # noqa: RUF012
+ implemented_properties = [
+ "energy",
+ "forces",
+ "hessian",
+ "enthalpy",
+ "entropy",
+ "qc_output",
+ "qc_input",
+ "custodian",
+ ] # noqa: RUF012
 
  def __init__(
  self,
@@ -177,8 +200,14 @@ def write_input(self, atoms, properties=None, system_changes=None):
  qcin.write("mol.qin")
 
  def read_results(self):
- data = QCOutput("mol.qout").data
- self.results["energy"] = data["final_energy"] * units.Hartree
+ qc_input = QCInput.from_file("mol.qin").as_dict()
+ qc_output = QCOutput("mol.qout").data
+ self.results["qc_output"] = qc_output
+ self.results["qc_input"] = qc_input
+ self.results["custodian"] = _parse_custodian(Path.cwd())
+
+ self.results["energy"] = qc_output["final_energy"] * units.Hartree
+
  if self.job_type in ["force", "opt"]:
  tmp_grad_data = []
  # Read the gradient scratch file in 8 byte chunks
@@ -198,10 +227,10 @@ def read_results(self):
  ]
  # Ensure that the scratch values match the correct values from the
  # output file but with higher precision
- if data["pcm_gradients"] is not None:
- gradient = data["pcm_gradients"][0]
+ if qc_output["pcm_gradients"] is not None:
+ gradient = qc_output["pcm_gradients"][0]
  else:
- gradient = data["gradients"][0]
+ gradient = qc_output["gradients"][0]
  for ii, subgrad in enumerate(grad):
  for jj, val in enumerate(subgrad):
  if abs(gradient[ii, jj] - val) > 1e-6:
@@ -213,14 +242,16 @@ def read_results(self):
  self.results["forces"] = gradient * (-units.Hartree / units.Bohr)
  else:
  self.results["forces"] = None
- self.prev_orbital_coeffs = []
+
  # Read orbital coefficients scratch file in 8 byte chunks
+ self.prev_orbital_coeffs = []
  with zopen("53.0", mode="rb") as file:
  binary = file.read()
  self.prev_orbital_coeffs.extend(
  struct.unpack("d", binary[ii * 8 : (ii + 1) * 8])[0]
  for ii in range(int(len(binary) / 8))
  )
+
  if self.job_type == "freq":
  tmp_hess_data = []
  # Read Hessian scratch file in 8 byte chunks
@@ -232,12 +263,13 @@ def read_results(self):
  )
  self.results["hessian"] = np.reshape(
  np.array(tmp_hess_data),
- (len(data["species"]) * 3, len(data["species"]) * 3),
+ (len(qc_output["species"]) * 3, len(qc_output["species"]) * 3),
+ )
+ self.results["enthalpy"] = qc_output["total_enthalpy"] * (
+ units.kcal / units.mol
+ )
+ self.results["entropy"] = qc_output["total_entropy"] * (
+ 0.001 * units.kcal / units.mol
  )
- data["enthalpy"] = data["total_enthalpy"] * (units.kcal / units.mol)
- data["entropy"] = data["total_entropy"] * (0.001 * units.kcal / units.mol)
- for k in ["total_enthalpy", "total_entropy"]:
- data.pop(k)
  else:
  self.results["hessian"] = None
- self.results["qc_output"] = data

tests/calculators/qchem/test_qchem.py

@@ -213,6 +213,4 @@ def test_qchem_read_results_freq(tmpdir):
  assert len(calc.results["qc_output"]["frequency_mode_vectors"]) == 36
  assert len(calc.results["qc_output"]["frequency_mode_vectors"][0]) == 14
  assert len(calc.results["qc_output"]["frequency_mode_vectors"][0][0]) == 3
- assert calc.results["qc_output"]["enthalpy"] == 2.647248450819514
- assert calc.results["qc_output"]["entropy"] == 0.003996739364205975
  tmpdir.chdir()

tests/recipes/qchem/qchem_examples/custodian.json

@@ -0,0 +1,27 @@
+[
+ {
+ "job": {
+ "@module": "custodian.qchem.jobs",
+ "@class": "QCJob",
+ "@version": "2023.3.10",
+ "qchem_command": ["qchem"],
+ "max_cores": 40,
+ "multimode": "openmp",
+ "input_file": "mol.qin",
+ "output_file": "mol.qout",
+ "qclog_file": "mol.qclog",
+ "suffix": "",
+ "calc_loc": "/tmp",
+ "nboexe": "/clusterfs/mp/software/nbo7/bin/nbo7.i4.exe",
+ "save_scratch": false,
+ "backup": true
+ },
+ "corrections": [],
+ "handler": null,
+ "validator": null,
+ "max_errors": false,
+ "max_errors_per_job": false,
+ "max_errors_per_handler": false,
+ "nonzero_return_code": false
+ }
+]

tests/recipes/qchem/test_qchem_recipes.py

@@ -70,18 +70,21 @@ def mock_execute1(_self, **kwargs):
  copy(os.path.join(QCHEM_DIR, "mol.qout.basic"), "mol.qout")
  copy(os.path.join(QCHEM_DIR, "131.0.basic"), "131.0")
  copy(os.path.join(QCHEM_DIR, "53.0.basic"), "53.0")
+ copy(os.path.join(QCHEM_DIR, "custodian.json"), "custodian.json")
 
 
 def mock_execute2(_self, **kwargs):
  copy(os.path.join(QCHEM_DIR, "mol.qout.intermediate"), "mol.qout")
  copy(os.path.join(QCHEM_DIR, "131.0.intermediate"), "131.0")
  copy(os.path.join(QCHEM_DIR, "53.0.intermediate"), "53.0")
+ copy(os.path.join(QCHEM_DIR, "custodian.json"), "custodian.json")
 
 
 def mock_execute3(_self, **kwargs):
  copy(os.path.join(QCHEM_DIR, "mol.qout.alternate"), "mol.qout")
  copy(os.path.join(QCHEM_DIR, "131.0.alternate"), "131.0")
  copy(os.path.join(QCHEM_DIR, "53.0.alternate"), "53.0")
+ copy(os.path.join(QCHEM_DIR, "custodian.json"), "custodian.json")
 
 
 def mock_execute4(self, **kwargs):
@@ -97,6 +100,7 @@ def mock_execute5(_self, **kwargs):
  copy(os.path.join(QCHEM_DIR, "mol.qout.freq"), "mol.qout")
  copy(os.path.join(QCHEM_DIR, "132.0.freq"), "132.0")
  copy(os.path.join(QCHEM_DIR, "53.0.freq"), "53.0")
+ copy(os.path.join(QCHEM_DIR, "custodian.json"), "custodian.json")
 
 
 def mock_read(self, **kwargs):
@@ -106,7 +110,6 @@ def mock_read(self, **kwargs):
 
 def test_static_job_v1(monkeypatch, tmpdir):
  tmpdir.chdir()
-
  monkeypatch.setattr(FileIOCalculator, "execute", mock_execute1)
  charge, spin_multiplicity = check_charge_and_spin(TEST_ATOMS)
  output = static_job(TEST_ATOMS, charge, spin_multiplicity)
@@ -119,10 +122,12 @@ def test_static_job_v1(monkeypatch, tmpdir):
  assert output["parameters"]["spin_multiplicity"] == 1
  assert output["results"]["energy"] == pytest.approx(-606.1616819641 * units.Hartree)
  assert output["results"]["forces"][0][0] == pytest.approx(-1.3826330655069403)
+ assert output["results"]["custodian"][0]["job"]["max_cores"] == 40
 
  qcin = QCInput.from_file("mol.qin.gz")
  ref_qcin = QCInput.from_file(os.path.join(QCHEM_DIR, "mol.qin.basic"))
  qcinput_nearly_equal(qcin, ref_qcin)
+ qcinput_nearly_equal(ref_qcin, QCInput.from_dict(output["results"]["qc_input"]))
 
 
 def test_static_job_v2(monkeypatch, tmpdir):
@@ -152,6 +157,7 @@ def test_static_job_v2(monkeypatch, tmpdir):
  qcin = QCInput.from_file("mol.qin.gz")
  ref_qcin = QCInput.from_file(os.path.join(QCHEM_DIR, "mol.qin.intermediate"))
  qcinput_nearly_equal(qcin, ref_qcin)
+ qcinput_nearly_equal(ref_qcin, QCInput.from_dict(output["results"]["qc_input"]))
 
 
 def test_static_job_v3(monkeypatch, tmpdir):
@@ -180,6 +186,7 @@ def test_static_job_v3(monkeypatch, tmpdir):
  qcin = QCInput.from_file("mol.qin.gz")
  ref_qcin = QCInput.from_file(os.path.join(QCHEM_DIR, "mol.qin.alternate"))
  qcinput_nearly_equal(qcin, ref_qcin)
+ qcinput_nearly_equal(ref_qcin, QCInput.from_dict(output["results"]["qc_input"]))
 
 
 def test_static_job_v4(monkeypatch, tmpdir):
@@ -229,6 +236,7 @@ def test_relax_job_v1(monkeypatch, tmpdir):
  os.path.join(QCHEM_DIR, "mol.qin.basic.sella_opt_iter1")
  )
  qcinput_nearly_equal(qcin, ref_qcin)
+ assert len(output["results"]["qc_input"]) > 1
 
 
 @pytest.mark.skipif(
@@ -327,8 +335,15 @@ def test_freq_job_v1(monkeypatch, tmpdir):
  assert output["parameters"]["spin_multiplicity"] == 2
  assert output["results"]["energy"] == pytest.approx(-605.6859554019 * units.Hartree)
  assert output["results"]["hessian"] is not None
- assert output["results"]["qc_output"]["enthalpy"] == pytest.approx(
- 61.047 * (units.kcal / units.mol)
+ assert output["results"]["enthalpy"] == pytest.approx(
+ output["results"]["qc_output"]["total_enthalpy"] * (units.kcal / units.mol)
+ )
+ assert (
+ output["results"]["entropy"]
+ == output["results"]["qc_output"]["total_entropy"]
+ * 0.001
+ * units.kcal
+ / units.mol
  )