bugfix: chargemol runtime directory and os.path.exists fix

pymatgen/command_line/chargemol_caller.py

@@ -47,7 +47,7 @@
 import warnings
 from glob import glob
 from shutil import which
-from typing import TYPE_CHECKING
+from typing import TYPE_CHECKING, Literal
 
 import numpy as np
 from monty.tempfile import ScratchDir
@@ -96,16 +96,14 @@ def __init__(
             run_chargemol (bool): Whether to run the Chargemol analysis. If False,
                 the existing Chargemol output files will be read from path. Default: True.
         """
-        path = path or os.getcwd()
+        path = os.path.abspath(path) if path else os.getcwd()
         if run_chargemol and not CHARGEMOL_EXE:
             raise OSError(
                 "ChargemolAnalysis requires the Chargemol executable to be in PATH."
                 " Please download the library at https://sourceforge.net/projects/ddec/files"
                 "and follow the instructions."
             )
-        if atomic_densities_path == "":
-            atomic_densities_path = os.getcwd()
-        self._atomic_densities_path = atomic_densities_path
+        self._atomic_densities_path = atomic_densities_path or os.getcwd()
 
         self._chgcar_path = self._get_filepath(path, "CHGCAR")
         self._potcar_path = self._get_filepath(path, "POTCAR")
@@ -170,34 +168,34 @@ def _execute_chargemol(self, **job_control_kwargs):
         """Internal function to run Chargemol.
 
         Args:
+            path (str): Path to the directory to run Chargemol in.
+                Default: None (current working directory).
             atomic_densities_path (str): Path to the atomic densities directory
             required by Chargemol. If None, Pymatgen assumes that this is
             defined in a "DDEC6_ATOMIC_DENSITIES_DIR" environment variable.
                 Default: None.
             job_control_kwargs: Keyword arguments for _write_jobscript_for_chargemol.
         """
-        with ScratchDir("."):
-            try:
-                os.symlink(self._chgcar_path, "./CHGCAR")
-                os.symlink(self._potcar_path, "./POTCAR")
-                os.symlink(self._aeccar0_path, "./AECCAR0")
-                os.symlink(self._aeccar2_path, "./AECCAR2")
-            except OSError as exc:
-                print(f"Error creating symbolic link: {exc}")
-
-            # write job_script file:
-            self._write_jobscript_for_chargemol(**job_control_kwargs)
-
-            # Run Chargemol
-            with subprocess.Popen(CHARGEMOL_EXE, stdout=subprocess.PIPE, stdin=subprocess.PIPE, close_fds=True) as rs:
-                _stdout, stderr = rs.communicate()
-            if rs.returncode != 0:
-                raise RuntimeError(
-                    f"{CHARGEMOL_EXE} exit code: {rs.returncode}, error message: {stderr!s}. "
-                    "Please check your Chargemol installation."
-                )
-
-            self._from_data_dir()
+        with ScratchDir(".") as temp_dir:
+            # with tempfile.TemporaryDirectory() as temp_dir:
+
+            os.symlink(self._chgcar_path, "CHGCAR")
+            os.symlink(self._potcar_path, "POTCAR")
+            os.symlink(self._aeccar0_path, "AECCAR0")
+            os.symlink(self._aeccar2_path, "AECCAR2")
+
+            lines = self._write_jobscript_for_chargemol(
+                chgcar_path=os.path.join(temp_dir, "CHGCAR"),
+                **job_control_kwargs,
+            )
+            with open(os.path.join(temp_dir, "job_control.txt"), mode="w") as file:
+                file.write(lines)
+
+            response = subprocess.run(CHARGEMOL_EXE, capture_output=True, check=False)
+            _stdout = response.stdout.decode("utf-8")
+            response.stderr.decode("utf-8")
+
+            self._from_data_dir(chargemol_output_path=temp_dir)
 
     def _from_data_dir(self, chargemol_output_path=None):
         """Internal command to parse Chargemol files from a directory.
@@ -207,8 +205,7 @@ def _from_data_dir(self, chargemol_output_path=None):
             Chargemol output files.
                 Default: None (current working directory).
         """
-        if chargemol_output_path is None:
-            chargemol_output_path = "."
+        chargemol_output_path = chargemol_output_path or "."
 
         charge_path = f"{chargemol_output_path}/DDEC6_even_tempered_net_atomic_charges.xyz"
         self.ddec_charges = self._get_data_from_xyz(charge_path)
@@ -227,21 +224,21 @@ def _from_data_dir(self, chargemol_output_path=None):
         else:
             self.ddec_spin_moments = None
 
-        rsquared_path = f"{chargemol_output_path}/DDEC_atomic_Rsquared_moments.xyz"
-        if os.path.isfile(rsquared_path):
-            self.ddec_rsquared_moments = self._get_data_from_xyz(rsquared_path)
+        r2_path = f"{chargemol_output_path}/DDEC_atomic_Rsquared_moments.xyz"
+        if os.path.isfile(r2_path):
+            self.ddec_rsquared_moments = self._get_data_from_xyz(r2_path)
         else:
             self.ddec_rsquared_moments = None
 
-        rcubed_path = f"{chargemol_output_path}/DDEC_atomic_Rcubed_moments.xyz"
-        if os.path.isfile(rcubed_path):
-            self.ddec_rcubed_moments = self._get_data_from_xyz(rcubed_path)
+        r3_path = f"{chargemol_output_path}/DDEC_atomic_Rcubed_moments.xyz"
+        if os.path.isfile(r3_path):
+            self.ddec_rcubed_moments = self._get_data_from_xyz(r3_path)
         else:
             self.ddec_rcubed_moments = None
 
-        rfourth_path = f"{chargemol_output_path}/DDEC_atomic_Rfourth_moments.xyz"
-        if os.path.isfile(rfourth_path):
-            self.ddec_rfourth_moments = self._get_data_from_xyz(rfourth_path)
+        r4_path = f"{chargemol_output_path}/DDEC_atomic_Rfourth_moments.xyz"
+        if os.path.isfile(r4_path):
+            self.ddec_rfourth_moments = self._get_data_from_xyz(r4_path)
         else:
             self.ddec_rfourth_moments = None
 
@@ -251,7 +248,7 @@ def _from_data_dir(self, chargemol_output_path=None):
         else:
             self.cm5_charges = None
 
-    def get_charge_transfer(self, atom_index, charge_type="ddec"):
+    def get_charge_transfer(self, atom_index: int, charge_type: Literal["cm5", "ddec"] = "ddec") -> float:
         """Returns the charge transferred for a particular atom. A positive value means
         that the site has gained electron density (i.e. exhibits anionic character)
         whereas a negative value means the site has lost electron density (i.e. exhibits
@@ -265,13 +262,11 @@ def get_charge_transfer(self, atom_index, charge_type="ddec"):
         Returns:
             float: charge transferred at atom_index
         """
-        if charge_type.lower() not in ["ddec", "cm5"]:
-            raise ValueError(f"Invalid {charge_type=}")
         if charge_type.lower() == "ddec":
-            charge_transfer = -self.ddec_charges[atom_index]
-        elif charge_type.lower() == "cm5":
-            charge_transfer = -self.cm5_charges[atom_index]
-        return charge_transfer
+            return -self.ddec_charges[atom_index]
+        if charge_type.lower() == "cm5":
+            return -self.cm5_charges[atom_index]
+        raise ValueError(f"Invalid {charge_type=}")
 
     def get_charge(self, atom_index, nelect=None, charge_type="ddec"):
         """Convenience method to get the charge on a particular atom using the same
@@ -337,6 +332,7 @@ def get_bond_order(self, index_from, index_to):
 
     def _write_jobscript_for_chargemol(
         self,
+        chgcar_path=None,
         net_charge=0.0,
         periodicity=(True, True, True),
         method="ddec6",
@@ -345,6 +341,8 @@ def _write_jobscript_for_chargemol(
         """Writes job_script.txt for Chargemol execution.
 
         Args:
+            chgcar_path (str): Path to the CHGCAR file. If None, CHGCAR is assumed to be
+                in the directory where the job_script.txt is written.
             net_charge (float): Net charge of the system.
                 Defaults to 0.0.
             periodicity (tuple[bool]): Periodicity of the system.
@@ -363,6 +361,9 @@ def _write_jobscript_for_chargemol(
         if net_charge:
             lines += f"<net charge>\n{net_charge}\n</net charge>\n"
 
+        if chgcar_path:
+            lines += f"<input filename>\n{chgcar_path}\n</input filename>\n"
+
         # Periodicity
         if periodicity:
             per_a = ".true." if periodicity[0] else ".false."
@@ -380,7 +381,7 @@ def _write_jobscript_for_chargemol(
                 "The DDEC6_ATOMIC_DENSITIES_DIR environment variable must be set or the atomic_densities_path must"
                 " be specified"
             )
-        if not os.path.isfile(atomic_densities_path):
+        if not os.path.exists(atomic_densities_path):
             raise FileNotFoundError(f"{atomic_densities_path=} does not exist")
 
         # This is to fix a Chargemol filepath nuance
@@ -402,8 +403,7 @@ def _write_jobscript_for_chargemol(
             bo = ".true." if compute_bond_orders else ".false."
             lines += f"\n<compute BOs>\n{bo}\n</compute BOs>\n"
 
-        with open("job_control.txt", mode="w") as file:
-            file.write(lines)
+        return lines
 
     @staticmethod
     def _get_dipole_info(filepath):

tests/command_line/test_chargemol_caller.py

@@ -1,55 +1,74 @@
 from __future__ import annotations
 
+from unittest.mock import patch
+
+import pytest
+
 from pymatgen.command_line.chargemol_caller import ChargemolAnalysis
-from pymatgen.core import Element
+from pymatgen.core import Element, Structure
 from pymatgen.util.testing import TEST_FILES_DIR
 
 
 class TestChargemolAnalysis:
     def test_parse_chargemol(self):
         test_dir = f"{TEST_FILES_DIR}/chargemol/spin_unpolarized"
-        ca = ChargemolAnalysis(path=test_dir, run_chargemol=False)
-        assert ca.ddec_charges == [0.8432, -0.8432]
-        assert ca.get_partial_charge(0) == 0.8432
-        assert ca.get_partial_charge(0, charge_type="cm5") == 0.420172
-        assert ca.get_charge_transfer(0) == -0.8432
-        assert ca.get_charge_transfer(0, charge_type="cm5") == -0.420172
-        assert ca.get_charge(0, nelect=1) == 1 - 0.8432
-        assert ca.get_charge(0, nelect=1, charge_type="cm5") == 1 - 0.420172
-        assert ca.dipoles == [[0.0, 0.0, 0.0], [0.0, 0.0, 0.0]]
-        assert ca.ddec_spin_moments is None
-        assert ca.bond_order_sums == [0.53992, 0.901058]
-        assert ca.ddec_rsquared_moments == [8.261378, 34.237274]
-        assert ca.ddec_rcubed_moments == [14.496002, 88.169236]
-        assert ca.ddec_rfourth_moments == [37.648248, 277.371929]
-        assert ca.cm5_charges == [0.420172, -0.420172]
-        assert ca.summary["ddec"]["partial_charges"] == ca.ddec_charges
-        assert ca.summary["ddec"]["dipoles"] == ca.dipoles
-        assert ca.summary["ddec"]["bond_order_sums"] == ca.bond_order_sums
-        assert ca.summary["ddec"]["rsquared_moments"] == ca.ddec_rsquared_moments
-        assert ca.summary["ddec"]["rcubed_moments"] == ca.ddec_rcubed_moments
-        assert ca.summary["ddec"]["rfourth_moments"] == ca.ddec_rfourth_moments
-        assert ca.summary["cm5"]["partial_charges"] == ca.cm5_charges
-        assert ca.summary["ddec"]["bond_order_dict"] == ca.bond_order_dict
-        assert ca.summary["ddec"].get("spin_moments") is None
-        assert ca.natoms == [1, 1]
-        assert ca.structure is not None
-        assert len(ca.bond_order_dict) == 2
-        assert ca.bond_order_dict[0]["bonded_to"][0] == {
+        chg_mol = ChargemolAnalysis(path=test_dir, run_chargemol=False)
+        assert chg_mol.ddec_charges == [0.8432, -0.8432]
+        assert chg_mol.get_partial_charge(0) == 0.8432
+        assert chg_mol.get_partial_charge(0, charge_type="cm5") == 0.420172
+        assert chg_mol.get_charge_transfer(0) == -0.8432
+        assert chg_mol.get_charge_transfer(0, charge_type="cm5") == -0.420172
+        assert chg_mol.get_charge(0, nelect=1) == 1 - 0.8432
+        assert chg_mol.get_charge(0, nelect=1, charge_type="cm5") == 1 - 0.420172
+        assert chg_mol.dipoles == [[0.0, 0.0, 0.0], [0.0, 0.0, 0.0]]
+        assert chg_mol.ddec_spin_moments is None
+        assert chg_mol.bond_order_sums == [0.53992, 0.901058]
+        assert chg_mol.ddec_rsquared_moments == [8.261378, 34.237274]
+        assert chg_mol.ddec_rcubed_moments == [14.496002, 88.169236]
+        assert chg_mol.ddec_rfourth_moments == [37.648248, 277.371929]
+        assert chg_mol.cm5_charges == [0.420172, -0.420172]
+        assert chg_mol.summary["ddec"]["partial_charges"] == chg_mol.ddec_charges
+        assert chg_mol.summary["ddec"]["dipoles"] == chg_mol.dipoles
+        assert chg_mol.summary["ddec"]["bond_order_sums"] == chg_mol.bond_order_sums
+        assert chg_mol.summary["ddec"]["rsquared_moments"] == chg_mol.ddec_rsquared_moments
+        assert chg_mol.summary["ddec"]["rcubed_moments"] == chg_mol.ddec_rcubed_moments
+        assert chg_mol.summary["ddec"]["rfourth_moments"] == chg_mol.ddec_rfourth_moments
+        assert chg_mol.summary["cm5"]["partial_charges"] == chg_mol.cm5_charges
+        assert chg_mol.summary["ddec"]["bond_order_dict"] == chg_mol.bond_order_dict
+        assert chg_mol.summary["ddec"].get("spin_moments") is None
+        assert chg_mol.natoms == [1, 1]
+        assert isinstance(chg_mol.structure, Structure)
+        assert len(chg_mol.structure) == 2
+        assert chg_mol.structure.formula == "Na1 Cl1"
+        assert len(chg_mol.bond_order_dict) == 2
+        assert chg_mol.bond_order_dict[0]["bonded_to"][0] == {
             "spin_polarization": 0.0,
             "index": 1,
             "direction": (-1, -1, 0),
             "bond_order": 0.0882,
             "element": Element("Cl"),
         }
-        assert ca.bond_order_dict[1]["bonded_to"][-1]["direction"] == (-1, 0, 0)
-        assert ca.get_property_decorated_structure().site_properties["partial_charge_ddec6"] == ca.ddec_charges
+        assert chg_mol.bond_order_dict[1]["bonded_to"][-1]["direction"] == (-1, 0, 0)
+        # check that partial charges are written to structure site properties
+        charge_decorated_struct = chg_mol.get_property_decorated_structure()
+        struct_partial_charge = charge_decorated_struct.site_properties["partial_charge_ddec6"]
+        assert struct_partial_charge == chg_mol.ddec_charges
 
     def test_parse_chargemol2(self):
         test_dir = f"{TEST_FILES_DIR}/chargemol/spin_polarized"
-        ca = ChargemolAnalysis(path=test_dir, run_chargemol=False)
-        assert ca.ddec_spin_moments == [0.201595, 0.399203, 0.399203]
-        assert ca.summary["ddec"]["bond_order_dict"][0]["bonded_to"][0]["spin_polarization"] == 0.0490
-        assert ca.summary["ddec"]["spin_moments"] == ca.ddec_spin_moments
-        assert ca.natoms is None
-        assert ca.structure is None
+        chg_mol = ChargemolAnalysis(path=test_dir, run_chargemol=False)
+        assert chg_mol.ddec_spin_moments == [0.201595, 0.399203, 0.399203]
+        assert chg_mol.summary["ddec"]["bond_order_dict"][0]["bonded_to"][0]["spin_polarization"] == 0.0490
+        assert chg_mol.summary["ddec"]["spin_moments"] == chg_mol.ddec_spin_moments
+        assert chg_mol.natoms is None
+        assert chg_mol.structure is None
+
+    def test_missing_exe_error(self):
+        # monkeypatch CHARGEMOL_EXE to raise in ChargemolAnalysis.__init__
+        patch.dict("os.environ", {"CHARGEMOL_EXE": "non_existent"})
+
+        test_dir = f"{TEST_FILES_DIR}/chargemol/spin_unpolarized"
+        with pytest.raises(
+            OSError, match="ChargemolAnalysis requires the Chargemol executable to be in PATH. Please download"
+        ):
+            ChargemolAnalysis(path=test_dir, run_chargemol=True)