Merge pull request #869 from materialsproject/add-run-ddec6-option-to…

…-calc-from-vasp-files

Add `run_ddec6` option to `Calculation.from_vasp_files`

emmet-core/emmet/core/vasp/calculation.py

@@ -3,13 +3,14 @@
 # mypy: ignore-errors
 
 import logging
+from datetime import datetime
 from pathlib import Path
 from typing import Any, Dict, List, Optional, Tuple, Union
-from datetime import datetime
 
 import numpy as np
 from pydantic import BaseModel, Extra, Field
 from pymatgen.command_line.bader_caller import bader_analysis_from_path
+from pymatgen.command_line.chargemol_caller import ChargemolAnalysis
 from pymatgen.core.lattice import Lattice
 from pymatgen.core.structure import Structure
 from pymatgen.core.trajectory import Trajectory
@@ -18,17 +19,17 @@
 from pymatgen.electronic_structure.dos import CompleteDos, Dos
 from pymatgen.io.vasp import (
     BSVasprun,
+    Kpoints,
     Locpot,
     Outcar,
     Poscar,
     Potcar,
     PotcarSingle,
     Vasprun,
     VolumetricData,
-    Kpoints,
 )
 
-from emmet.core.math import Matrix3D, Vector3D, ListMatrix3D
+from emmet.core.math import ListMatrix3D, Matrix3D, Vector3D
 from emmet.core.utils import ValueEnum
 from emmet.core.vasp.calc_types import (
     CalcType,
@@ -134,7 +135,9 @@ class CalculationInput(BaseModel):
     structure: Optional[Structure] = Field(
         None, description="Input structure for the calculation"
     )
-    is_hubbard: bool = Field(False, description="Is this a Hubbard +U calculation")
+    is_hubbard: bool = Field(
+        default=False, description="Is this a Hubbard +U calculation"
+    )
     hubbards: Optional[dict] = Field(None, description="The hubbard parameters used")
 
     @classmethod
@@ -585,7 +588,8 @@ class Calculation(BaseModel):
         description="Paths (relative to dir_name) of the VASP output files "
         "associated with this calculation",
     )
-    bader: Optional[dict] = Field(None, description="Output from the bader software")
+    bader: Optional[dict] = Field(None, description="Output from bader charge analysis")
+    ddec6: Optional[dict] = Field(None, description="Output from DDEC6 charge analysis")
     run_type: Optional[RunType] = Field(
         None, description="Calculation run type (e.g., HF, HSE06, PBE)"
     )
@@ -610,6 +614,7 @@ def from_vasp_files(
         parse_bandstructure: Union[str, bool] = False,
         average_locpot: bool = True,
         run_bader: bool = False,
+        run_ddec6: bool | str = False,
         strip_bandstructure_projections: bool = False,
         strip_dos_projections: bool = False,
         store_volumetric_data: Optional[Tuple[str]] = None,
@@ -658,8 +663,13 @@ def from_vasp_files(
 
         average_locpot
             Whether to store the average of the LOCPOT along the crystal axes.
-        run_bader
+        run_bader : bool = False
             Whether to run bader on the charge density.
+        run_ddec6 : bool | str = False
+            Whether to run DDEC6 on the charge density. If a string, it's interpreted
+            as the path to the atomic densities directory. Can also be set via the
+            DDEC6_ATOMIC_DENSITIES_DIR environment variable. The files are available at
+            https://sourceforge.net/projects/ddec/files.
         strip_dos_projections
             Whether to strip the element and site projections from the density of
             states. This can help reduce the size of DOS objects in systems with many
@@ -715,6 +725,11 @@ def from_vasp_files(
             suffix = "" if task_name == "standard" else f".{task_name}"
             bader = bader_analysis_from_path(dir_name, suffix=suffix)
 
+        ddec6 = None
+        if run_ddec6 and VaspObject.CHGCAR in output_file_paths:
+            densities_path = run_ddec6 if isinstance(run_ddec6, (str, Path)) else None
+            ddec6 = ChargemolAnalysis(path=dir_name, atomic_densities_path=densities_path).summary
+
         locpot = None
         if average_locpot:
             if VaspObject.LOCPOT in vasp_objects:
@@ -768,6 +783,7 @@ def from_vasp_files(
                     k.name.lower(): v for k, v in output_file_paths.items()
                 },
                 bader=bader,
+                ddec6=ddec6,
                 run_type=run_type(input_doc.parameters),
                 task_type=task_type(input_doc.model_dump()),
                 calc_type=calc_type(input_doc.model_dump(), input_doc.parameters),

emmet-core/tests/test_calculation.py

@@ -1,4 +1,5 @@
 import pytest
+
 from tests.conftest import assert_schemas_equal, get_test_object
 
 
@@ -151,6 +152,10 @@ def test_calculation(test_dir, object_name, task_name):
     assert_schemas_equal(test_doc, valid_doc)
     assert set(objects.keys()) == set(test_object.objects[task_name])
 
+    # check bader and ddec6 keys exist
+    assert test_doc.bader is None
+    assert test_doc.ddec6 is None
+
     # test document can be jsanitized
     d = jsanitize(test_doc, strict=True, enum_values=True, allow_bson=True)