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
@@ -97,16 +97,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")
@@ -171,34 +169,32 @@ 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)
+
+            subprocess.run(CHARGEMOL_EXE, capture_output=True, check=True)
+
+            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.
@@ -208,8 +204,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)
@@ -228,21 +223,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
 
@@ -252,8 +247,8 @@ def _from_data_dir(self, chargemol_output_path=None):
         else:
             self.cm5_charges = None
 
-    def get_charge_transfer(self, atom_index, charge_type="ddec"):
-        """Get the charge transferred for a particular atom. A positive value means
+    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
         cationic character). This is the same thing as the negative of the partial atomic
@@ -338,6 +333,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",
@@ -346,6 +342,8 @@ def _write_jobscript_for_chargemol(
         """Write 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.
@@ -364,6 +362,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."
@@ -381,7 +382,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
@@ -403,8 +404,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,58 +1,106 @@
 from __future__ import annotations
 
+import os
+import subprocess
+import zipfile
+from unittest.mock import patch
+from urllib.request import urlretrieve
+
+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
 
-TEST_DIR = f"{TEST_FILES_DIR}/command_line/chargemol"
+ddec_url = "https://sourceforge.net/projects/ddec/files/latest/download"
+download_path = TEST_FILES_DIR / "command_line" / "chargemol" / "ddec-latest.gz"
+extract_path = TEST_FILES_DIR / "command_line" / "chargemol"
+
+if not os.path.exists(download_path):
+    urlretrieve(ddec_url, download_path)
+
+if not os.path.exists(extract_path / "chargemol_09_26_2017"):
+    with zipfile.ZipFile(download_path, "r") as zip_ref:
+        zip_ref.extractall(extract_path)
+
+exe_path = extract_path / "chargemol_09_26_2017/chargemol_FORTRAN_09_26_2017/compiled_binaries/linux"
+
+current_path = os.environ.get("PATH", "")
+if str(exe_path) not in current_path:
+    os.environ["PATH"] = str(exe_path) + os.pathsep + current_path
+
+subprocess.call(["chmod", "+x", str(exe_path / "Chargemol_09_26_2017_linux_serial")])
+subprocess.call(["chmod", "+x", str(exe_path / "Chargemol_09_26_2017_linux_parallel")])
 
 
 class TestChargemolAnalysis:
     def test_parse_chargemol(self):
-        test_dir = f"{TEST_DIR}/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] == {
+        test_dir = f"{TEST_FILES_DIR}/command_line/chargemol/spin_unpolarized"
+        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_DIR}/spin_polarized"
-        ca = ChargemolAnalysis(path=test_dir, run_chargemol=False)
-        assert ca.ddec_spin_moments == [0.201595, 0.399203, 0.399203]
-        assert ca.dipoles == [[-0.0, 0.0, -0.127251], [0.0, -0.027857, -0.010944], [0.0, 0.027857, -0.010944]]
-        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
+        test_dir = f"{TEST_FILES_DIR}/command_line/chargemol/spin_polarized"
+        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}/command_line/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)
+
+    def test_chargemol_custom_path(self):
+        chg_mol = ChargemolAnalysis(
+            path=TEST_FILES_DIR / "command_line" / "bader",
+            atomic_densities_path=extract_path / "chargemol_09_26_2017" / "atomic_densities",
+            run_chargemol=True,
+        )
+        print(chg_mol.ddec_charges)