diff --git a/pymatgen/command_line/chargemol_caller.py b/pymatgen/command_line/chargemol_caller.py
index d9b17a27404..0fea13f5b90 100644
--- a/pymatgen/command_line/chargemol_caller.py
+++ b/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"\n{net_charge}\n\n"
+ if chgcar_path:
+ lines += f"\n{chgcar_path}\n\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\n{bo}\n\n"
- with open("job_control.txt", mode="w") as file:
- file.write(lines)
+ return lines
@staticmethod
def _get_dipole_info(filepath):
diff --git a/tests/command_line/test_chargemol_caller.py b/tests/command_line/test_chargemol_caller.py
index 2c641f7c862..a5a21ba3567 100644
--- a/tests/command_line/test_chargemol_caller.py
+++ b/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)