4075885785

Build Type

push

github

Committed by Shyue Ping Ong

Commit Message

Merge branch 'master' of github.com:materialsproject/pymatgen

Run Details

96 of 96 new or added lines in 27 files covered. (100.0%)

81013 of 102710 relevant lines covered (78.88%)

0.79 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

16.02

/pymatgen/command_line/bader_caller.py

# Copyright (c) Pymatgen Development Team.
# Distributed under the terms of the MIT License.

"""
This module implements an interface to the Henkelmann et al.'s excellent
Fortran code for calculating a Bader charge analysis.

This module depends on a compiled bader executable available in the path.
Please download the library at http://theory.cm.utexas.edu/henkelman/code/bader/
and follow the instructions to compile the executable.

If you use this module, please cite the following:

G. Henkelman, A. Arnaldsson, and H. Jonsson, "A fast and robust algorithm for
Bader decomposition of charge density", Comput. Mater. Sci. 36, 254-360 (2006).
"""

from __future__ import annotations

import glob
import os
import shutil
import subprocess
import warnings
from shutil import which

import numpy as np
from monty.dev import requires
from monty.io import zopen
from monty.tempfile import ScratchDir

from pymatgen.io.common import VolumetricData
from pymatgen.io.vasp.inputs import Potcar
from pymatgen.io.vasp.outputs import Chgcar

__author__ = "shyuepingong"
__version__ = "0.1"
__maintainer__ = "Shyue Ping Ong"
__email__ = "shyuep@gmail.com"
__status__ = "Beta"
__date__ = "4/5/13"

BADEREXE = which("bader") or which("bader.exe")


class BaderAnalysis:
    """
    Bader analysis for Cube files and VASP outputs.

    .. attribute: data

        Atomic data parsed from bader analysis. Essentially a list of dicts
        of the form::

        [
            {
                "atomic_vol": 8.769,
                "min_dist": 0.8753,
                "charge": 7.4168,
                "y": 1.1598,
                "x": 0.0079,
                "z": 0.8348
            },
            ...
        ]

    .. attribute: vacuum_volume

        Vacuum volume of the Bader analysis.

    .. attribute: vacuum_charge

        Vacuum charge of the Bader analysis.

    .. attribute: nelectrons

        Number of electrons of the Bader analysis.

    .. attribute: chgcar

        Chgcar object associated with input CHGCAR file.

    .. attribute: atomic_densities

        list of charge densities for each atom centered on the atom
        excess 0's are removed from the array to reduce the size of the array
        the charge densities are dicts with the charge density map,
        the shift vector applied to move the data to the center, and the original dimension of the charge density map
        charge:
            {
            "data": charge density array
            "shift": shift used to center the atomic charge density
            "dim": dimension of the original charge density map
            }
    """

    @requires(
        which("bader") or which("bader.exe"),
        "BaderAnalysis requires the executable bader to be in the path."
        " Please download the library at http://theory.cm.utexas"
        ".edu/vasp/bader/ and compile the executable.",
    )
    def __init__(
        self,
        chgcar_filename=None,
        potcar_filename=None,
        chgref_filename=None,
        parse_atomic_densities=False,
        cube_filename=None,
    ):
        """
        Initializes the Bader caller.

        Args:
            chgcar_filename (str): The filename of the CHGCAR.
            potcar_filename (str): The filename of the POTCAR.
            chgref_filename (str): The filename of the reference charge density.
            parse_atomic_densities (bool): Optional. turns on atomic partition of the charge density
                charge densities are atom centered
            cube_filename (str): Optional. The filename of the cube file.
        """
        if not BADEREXE:
            raise RuntimeError(
                "BaderAnalysis requires the executable bader to be in the path."
                " Please download the library at http://theory.cm.utexas"
                ".edu/vasp/bader/ and compile the executable."
            )

        if not (cube_filename or chgcar_filename):
            raise ValueError("You must provide either a cube file or a CHGCAR")
        if cube_filename and chgcar_filename:
            raise ValueError("You cannot parse a cube and a CHGCAR at the same time.")
        self.parse_atomic_densities = parse_atomic_densities

        if chgcar_filename:
            fpath = os.path.abspath(chgcar_filename)
            self.is_vasp = True
            self.chgcar = Chgcar.from_file(fpath)
            self.structure = self.chgcar.structure
            self.potcar = Potcar.from_file(potcar_filename) if potcar_filename is not None else None
            self.natoms = self.chgcar.poscar.natoms
            chgrefpath = os.path.abspath(chgref_filename) if chgref_filename else None
            self.reference_used = bool(chgref_filename)

            # List of nelects for each atom from potcar
            potcar_indices = []
            for i, v in enumerate(self.natoms):
                potcar_indices += [i] * v
            self.nelects = (
                [self.potcar[potcar_indices[i]].nelectrons for i in range(len(self.structure))] if self.potcar else []
            )

        else:
            fpath = os.path.abspath(cube_filename)
            self.is_vasp = False
            self.cube = VolumetricData.from_cube(fpath)
            self.structure = self.cube.structure
            self.nelects = None
            chgrefpath = os.path.abspath(chgref_filename) if chgref_filename else None
            self.reference_used = bool(chgref_filename)

        with ScratchDir("."):
            tmpfile = "CHGCAR" if chgcar_filename else "CUBE"
            with zopen(fpath, "rt") as f_in:
                with open(tmpfile, "w") as f_out:
                    shutil.copyfileobj(f_in, f_out)
            args = [BADEREXE, tmpfile]
            if chgref_filename:
                with zopen(chgrefpath, "rt") as f_in:
                    with open("CHGCAR_ref", "w") as f_out:
                        shutil.copyfileobj(f_in, f_out)
                args += ["-ref", "CHGCAR_ref"]
            if parse_atomic_densities:
                args += ["-p", "all_atom"]
            with subprocess.Popen(args, stdout=subprocess.PIPE, stdin=subprocess.PIPE, close_fds=True) as rs:
                stdout, _ = rs.communicate()
            if rs.returncode != 0:
                raise RuntimeError(
                    f"bader exited with return code {rs.returncode}. Please check your bader installation."
                )

            try:
                self.version = float(stdout.split()[5])
            except ValueError:
                self.version = -1  # Unknown
            if self.version < 1.0:
                warnings.warn(
                    "Your installed version of Bader is outdated, calculation of vacuum charge may be incorrect.",
                    UserWarning,
                )

            data = []
            with open("ACF.dat") as f:
                raw = f.readlines()
                headers = ("x", "y", "z", "charge", "min_dist", "atomic_vol")
                raw.pop(0)
                raw.pop(0)
                while True:
                    l = raw.pop(0).strip()
                    if l.startswith("-"):
                        break
                    vals = map(float, l.split()[1:])
                    data.append(dict(zip(headers, vals)))
                for l in raw:
                    toks = l.strip().split(":")
                    if toks[0] == "VACUUM CHARGE":
                        self.vacuum_charge = float(toks[1])
                    elif toks[0] == "VACUUM VOLUME":
                        self.vacuum_volume = float(toks[1])
                    elif toks[0] == "NUMBER OF ELECTRONS":
                        self.nelectrons = float(toks[1])
            self.data = data

            if self.parse_atomic_densities:
                # convert the charge density for each atom spit out by Bader into Chgcar objects for easy parsing
                atom_chgcars = [
                    Chgcar.from_file(f"BvAt{str(i).zfill(4)}.dat") for i in range(1, len(self.chgcar.structure) + 1)
                ]

                atomic_densities = []
                # For each atom in the structure
                for _, loc, chg in zip(
                    self.chgcar.structure,
                    self.chgcar.structure.frac_coords,
                    atom_chgcars,
                ):
                    # Find the index of the atom in the charge density atom
                    index = np.round(np.multiply(loc, chg.dim))

                    data = chg.data["total"]
                    # Find the shift vector in the array
                    shift = (np.divide(chg.dim, 2) - index).astype(int)

                    # Shift the data so that the atomic charge density to the center for easier manipulation
                    shifted_data = np.roll(data, shift, axis=(0, 1, 2))

                    # Slices a central window from the data array
                    def slice_from_center(data, xwidth, ywidth, zwidth):
                        x, y, z = data.shape
                        startx = x // 2 - (xwidth // 2)
                        starty = y // 2 - (ywidth // 2)
                        startz = z // 2 - (zwidth // 2)
                        return data[
                            startx : startx + xwidth,
                            starty : starty + ywidth,
                            startz : startz + zwidth,
                        ]

                    # Finds the central encompassing volume which holds all the data within a precision
                    def find_encompassing_vol(data):
                        total = np.sum(data)
                        for i in range(np.max(data.shape)):
                            sliced_data = slice_from_center(data, i, i, i)
                            if total - np.sum(sliced_data) < 0.1:
                                return sliced_data
                        return None

                    d = {
                        "data": find_encompassing_vol(shifted_data),
                        "shift": shift,
                        "dim": self.chgcar.dim,
                    }
                    atomic_densities.append(d)
                self.atomic_densities = atomic_densities

    def get_charge(self, atom_index):
        """
        Convenience method to get the charge on a particular atom. This is the "raw"
        charge generated by the Bader program, not a partial atomic charge. If the cube file
        is a spin-density file, then this will return the spin density per atom with
        positive being spin up and negative being spin down.

        Args:
            atom_index:
                Index of atom.

        Returns:
            Charge associated with atom from the Bader analysis.
        """
        return self.data[atom_index]["charge"]

    def get_charge_transfer(self, atom_index, nelect=None):
        """
        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). If the arg nelect is not supplied, then POTCAR must be
        supplied to determine nelect.

        Args:
            atom_index:
                Index of atom.
            nelect:
                number of electrons associated with an isolated atom at this index.
                For most DFT codes this corresponds to the number of valence electrons
                associated with the pseudopotential (e.g. ZVAL for VASP).

        Returns:
            Charge transfer associated with atom from the Bader analysis.
            Given by bader charge on atom - nelect for associated atom.
        """
        if not self.nelects and nelect is None:
            raise ValueError("No NELECT info! Need POTCAR for VASP or nelect argument for cube file")
        return self.data[atom_index]["charge"] - (nelect if nelect is not None else self.nelects[atom_index])

    def get_partial_charge(self, atom_index, nelect=None):
        """
        Convenience method to get the partial charge on a particular atom. This is
        simply the negative value of the charge transferred. A positive value indicates
        that the atom has cationic character, whereas a negative value indicates the
        site has anionic character.

        Args:
            atom_index:
                Index of atom.
            nelect:
                number of electrons associated with an isolated atom at this index.
                For most DFT codes this corresponds to the number of valence electrons
                associated with the pseudopotential (e.g. ZVAL for VASP).

        Returns:
            Charge associated with atom from the Bader analysis.
        """
        return -self.get_charge_transfer(atom_index, nelect)

    def get_charge_decorated_structure(self):
        """
        Returns a charge decorated structure

        Note, this assumes that the Bader analysis was correctly performed on a file
        with electron densities
        """
        charges = [-self.get_charge(i) for i in range(len(self.structure))]
        struct = self.structure.copy()
        struct.add_site_property("charge", charges)
        return struct

    def get_oxidation_state_decorated_structure(self, nelects=None):
        """
        Returns an oxidation state decorated structure based on bader analysis results.
        Each site is assigned a charge based on the computed partial atomic charge from bader.

        Note, this assumes that the Bader analysis was correctly performed on a file
        with electron densities
        """
        charges = [self.get_partial_charge(i, None if not nelects else nelects[i]) for i in range(len(self.structure))]
        struct = self.structure.copy()
        struct.add_oxidation_state_by_site(charges)
        return struct

    def get_decorated_structure(self, property_name, average=False):
        """
        Get a property-decorated structure from the Bader analysis.

        This is distinct from getting charge decorated structure, which assumes
        the "standard" Bader analysis of electron densities followed by converting
        electron count to charge. The expected way to use this is to call Bader on
        a non-charge density file such as a spin density file, electrostatic potential
        file, etc., while using the charge density file as the reference (chgref_filename)
        so that the partitioning is determined via the charge, but averaging or integrating
        is done for another property.

        User warning: Bader analysis cannot automatically determine what property is
        inside of the file. So if you want to use this for a non-conventional property
        like spin, you must ensure that you have the file is for the appropriate
        property and you have an appropriate reference file.

        Args:
            property_name: name of the property to assign to the structure, note that
                if name is "spin" this is handled as a special case, and the appropriate
                spin properties are set on the species in the structure
            average: whether or not to return the average of this property, rather
                than the total, by dividing by the atomic volume.

        Returns:
            structure with site properties assigned via Bader Analysis
        """
        vals = [self.get_charge(i) for i in range(len(self.structure))]
        struct = self.structure.copy()
        if average:
            vals = np.divide(vals, [d["atomic_vol"] for d in self.data])
        struct.add_site_property(property_name, vals)
        if property_name == "spin":
            struct.add_spin_by_site(vals)
        return struct

    @property
    def summary(self):
        """
        :return: Dict summary of key analysis, e.g., atomic volume, charge, etc.
        """
        summary = {
            "min_dist": [d["min_dist"] for d in self.data],
            "charge": [d["charge"] for d in self.data],
            "atomic_volume": [d["atomic_vol"] for d in self.data],
            "vacuum_charge": self.vacuum_charge,
            "vacuum_volume": self.vacuum_volume,
            "reference_used": self.reference_used,
            "bader_version": self.version,
        }

        if self.parse_atomic_densities:
            summary["charge_densities"] = self.atomic_densities

        if self.potcar:
            charge_transfer = [self.get_charge_transfer(i) for i in range(len(self.data))]
            summary["charge_transfer"] = charge_transfer

        return summary

    @classmethod
    def from_path(cls, path, suffix=""):
        """
        Convenient constructor that takes in the path name of VASP run
        to perform Bader analysis.

        Args:
            path (str): Name of directory where VASP output files are
                stored.
            suffix (str): specific suffix to look for (e.g. '.relax1'
                for 'CHGCAR.relax1.gz').
        """

        def _get_filepath(filename):
            name_pattern = filename + suffix + "*" if filename != "POTCAR" else filename + "*"
            paths = glob.glob(os.path.join(path, name_pattern))
            fpath = None
            if len(paths) >= 1:
                # using reverse=True because, if multiple files are present,
                # they likely have suffixes 'static', 'relax', 'relax2', etc.
                # and this would give 'static' over 'relax2' over 'relax'
                # however, better to use 'suffix' kwarg to avoid this!
                paths.sort(reverse=True)
                warning_msg = f"Multiple files detected, using {os.path.basename(paths[0])}" if len(paths) > 1 else None
                fpath = paths[0]
            else:
                warning_msg = f"Could not find {filename}"
                if filename in ["AECCAR0", "AECCAR2"]:
                    warning_msg += ", cannot calculate charge transfer."
                elif filename == "POTCAR":
                    warning_msg += ", interpret Bader results with caution."
            if warning_msg:
                warnings.warn(warning_msg)
            return fpath

        chgcar_filename = _get_filepath("CHGCAR")
        if chgcar_filename is None:
            raise FileNotFoundError("Could not find CHGCAR!")
        potcar_filename = _get_filepath("POTCAR")
        aeccar0 = _get_filepath("AECCAR0")
        aeccar2 = _get_filepath("AECCAR2")
        if aeccar0 and aeccar2:
            # `chgsum.pl AECCAR0 AECCAR2` equivalent to obtain chgref_file
            chgref = Chgcar.from_file(aeccar0) + Chgcar.from_file(aeccar2)
            chgref_filename = "CHGREF"
            chgref.write_file(chgref_filename)
        else:
            chgref_filename = None
        return cls(
            chgcar_filename=chgcar_filename,
            potcar_filename=potcar_filename,
            chgref_filename=chgref_filename,
        )


def bader_analysis_from_path(path, suffix=""):
    """
    Convenience method to run Bader analysis on a folder containing
    typical VASP output files.

    This method will:

    1. Look for files CHGCAR, AECCAR0, AECCAR2, POTCAR or their gzipped
    counterparts.
    2. If AECCAR* files are present, constructs a temporary reference
    file as AECCAR0 + AECCAR2
    3. Runs Bader analysis twice: once for charge, and a second time
    for the charge difference (magnetization density).

    :param path: path to folder to search in
    :param suffix: specific suffix to look for (e.g. '.relax1' for 'CHGCAR.relax1.gz'
    :return: summary dict
    """

    def _get_filepath(filename, warning, path=path, suffix=suffix):
        paths = glob.glob(os.path.join(path, filename + suffix + "*"))
        if not paths:
            warnings.warn(warning)
            return None
        if len(paths) > 1:
            # using reverse=True because, if multiple files are present,
            # they likely have suffixes 'static', 'relax', 'relax2', etc.
            # and this would give 'static' over 'relax2' over 'relax'
            # however, better to use 'suffix' kwarg to avoid this!
            paths.sort(reverse=True)
            warnings.warn(f"Multiple files detected, using {os.path.basename(path)}")
        path = paths[0]
        return path

    chgcar_path = _get_filepath("CHGCAR", "Could not find CHGCAR!")
    chgcar = Chgcar.from_file(chgcar_path)

    aeccar0_path = _get_filepath("AECCAR0", "Could not find AECCAR0, interpret Bader results with caution.")
    aeccar0 = Chgcar.from_file(aeccar0_path) if aeccar0_path else None

    aeccar2_path = _get_filepath("AECCAR2", "Could not find AECCAR2, interpret Bader results with caution.")
    aeccar2 = Chgcar.from_file(aeccar2_path) if aeccar2_path else None

    potcar_path = _get_filepath("POTCAR", "Could not find POTCAR, cannot calculate charge transfer.")
    potcar = Potcar.from_file(potcar_path) if potcar_path else None

    return bader_analysis_from_objects(chgcar, potcar, aeccar0, aeccar2)


def bader_analysis_from_objects(chgcar, potcar=None, aeccar0=None, aeccar2=None):
    """
    Convenience method to run Bader analysis from a set
    of pymatgen Chgcar and Potcar objects.

    This method will:

    1. If aeccar objects are present, constructs a temporary reference
    file as AECCAR0 + AECCAR2
    2. Runs Bader analysis twice: once for charge, and a second time
    for the charge difference (magnetization density).

    :param chgcar: Chgcar object
    :param potcar: (optional) Potcar object
    :param aeccar0: (optional) Chgcar object from aeccar0 file
    :param aeccar2: (optional) Chgcar object from aeccar2 file
    :return: summary dict
    """

    with ScratchDir(".") as temp_dir:
        if aeccar0 and aeccar2:
            # construct reference file
            chgref = aeccar0.linear_add(aeccar2)
            chgref_path = os.path.join(temp_dir, "CHGCAR_ref")
            chgref.write_file(chgref_path)
        else:
            chgref_path = None

        chgcar.write_file("CHGCAR")
        chgcar_path = os.path.join(temp_dir, "CHGCAR")

        if potcar:
            potcar.write_file("POTCAR")
            potcar_path = os.path.join(temp_dir, "POTCAR")
        else:
            potcar_path = None

        ba = BaderAnalysis(
            chgcar_filename=chgcar_path,
            potcar_filename=potcar_path,
            chgref_filename=chgref_path,
        )

        summary = {
            "min_dist": [d["min_dist"] for d in ba.data],
            "charge": [d["charge"] for d in ba.data],
            "atomic_volume": [d["atomic_vol"] for d in ba.data],
            "vacuum_charge": ba.vacuum_charge,
            "vacuum_volume": ba.vacuum_volume,
            "reference_used": bool(chgref_path),
            "bader_version": ba.version,
        }

        if potcar:
            charge_transfer = [ba.get_charge_transfer(i) for i in range(len(ba.data))]
            summary["charge_transfer"] = charge_transfer

        if chgcar.is_spin_polarized:
            # write a CHGCAR containing magnetization density only
            chgcar.data["total"] = chgcar.data["diff"]
            chgcar.is_spin_polarized = False
            chgcar.write_file("CHGCAR_mag")

            chgcar_mag_path = os.path.join(temp_dir, "CHGCAR_mag")
            ba = BaderAnalysis(
                chgcar_filename=chgcar_mag_path,
                potcar_filename=potcar_path,
                chgref_filename=chgref_path,
            )
            summary["magmom"] = [d["charge"] for d in ba.data]

        return summary

1	# Copyright (c) Pymatgen Development Team.
2	# Distributed under the terms of the MIT License.
3
4	"""	1✔
5	This module implements an interface to the Henkelmann et al.'s excellent
6	Fortran code for calculating a Bader charge analysis.
7
8	This module depends on a compiled bader executable available in the path.
9	Please download the library at http://theory.cm.utexas.edu/henkelman/code/bader/
10	and follow the instructions to compile the executable.
11
12	If you use this module, please cite the following:
13
14	G. Henkelman, A. Arnaldsson, and H. Jonsson, "A fast and robust algorithm for
15	Bader decomposition of charge density", Comput. Mater. Sci. 36, 254-360 (2006).
16	"""
17
18	from __future__ import annotations	1✔
19
20	import glob	1✔
21	import os	1✔
22	import shutil	1✔
23	import subprocess	1✔
24	import warnings	1✔
25	from shutil import which	1✔
26
27	import numpy as np	1✔
28	from monty.dev import requires	1✔
29	from monty.io import zopen	1✔
30	from monty.tempfile import ScratchDir	1✔
31
32	from pymatgen.io.common import VolumetricData	1✔
33	from pymatgen.io.vasp.inputs import Potcar	1✔
34	from pymatgen.io.vasp.outputs import Chgcar	1✔
35
36	__author__ = "shyuepingong"	1✔
37	__version__ = "0.1"	1✔
38	__maintainer__ = "Shyue Ping Ong"	1✔
39	__email__ = "shyuep@gmail.com"	1✔
40	__status__ = "Beta"	1✔
41	__date__ = "4/5/13"	1✔
42
43	BADEREXE = which("bader") or which("bader.exe")	1✔
44
45
46	class BaderAnalysis:	1✔
47	"""
48	Bader analysis for Cube files and VASP outputs.
49
50	.. attribute: data
51
52	Atomic data parsed from bader analysis. Essentially a list of dicts
53	of the form::
54
55	[
56	{
57	"atomic_vol": 8.769,
58	"min_dist": 0.8753,
59	"charge": 7.4168,
60	"y": 1.1598,
61	"x": 0.0079,
62	"z": 0.8348
63	},
64	...
65	]
66
67	.. attribute: vacuum_volume
68
69	Vacuum volume of the Bader analysis.
70
71	.. attribute: vacuum_charge
72
73	Vacuum charge of the Bader analysis.
74
75	.. attribute: nelectrons
76
77	Number of electrons of the Bader analysis.
78
79	.. attribute: chgcar
80
81	Chgcar object associated with input CHGCAR file.
82
83	.. attribute: atomic_densities
84
85	list of charge densities for each atom centered on the atom
86	excess 0's are removed from the array to reduce the size of the array
87	the charge densities are dicts with the charge density map,
88	the shift vector applied to move the data to the center, and the original dimension of the charge density map
89	charge:
90	{
91	"data": charge density array
92	"shift": shift used to center the atomic charge density
93	"dim": dimension of the original charge density map
94	}
95	"""
96
97	@requires(	1✔
98	which("bader") or which("bader.exe"),
99	"BaderAnalysis requires the executable bader to be in the path."
100	" Please download the library at http://theory.cm.utexas"
101	".edu/vasp/bader/ and compile the executable.",
102	)
103	def __init__(	1✔
104	self,
105	chgcar_filename=None,
106	potcar_filename=None,
107	chgref_filename=None,
108	parse_atomic_densities=False,
109	cube_filename=None,
110	):
111	"""
112	Initializes the Bader caller.
113
114	Args:
115	chgcar_filename (str): The filename of the CHGCAR.
116	potcar_filename (str): The filename of the POTCAR.
117	chgref_filename (str): The filename of the reference charge density.
118	parse_atomic_densities (bool): Optional. turns on atomic partition of the charge density
119	charge densities are atom centered
120	cube_filename (str): Optional. The filename of the cube file.
121	"""
122	if not BADEREXE:	×
123	raise RuntimeError(	×
124	"BaderAnalysis requires the executable bader to be in the path."
125	" Please download the library at http://theory.cm.utexas"
126	".edu/vasp/bader/ and compile the executable."
127	)
128
129	if not (cube_filename or chgcar_filename):	×
130	raise ValueError("You must provide either a cube file or a CHGCAR")	×
131	if cube_filename and chgcar_filename:	×
132	raise ValueError("You cannot parse a cube and a CHGCAR at the same time.")	×
133	self.parse_atomic_densities = parse_atomic_densities	×
134
135	if chgcar_filename:	×
136	fpath = os.path.abspath(chgcar_filename)	×
137	self.is_vasp = True	×
138	self.chgcar = Chgcar.from_file(fpath)	×
139	self.structure = self.chgcar.structure	×
140	self.potcar = Potcar.from_file(potcar_filename) if potcar_filename is not None else None	×
141	self.natoms = self.chgcar.poscar.natoms	×
142	chgrefpath = os.path.abspath(chgref_filename) if chgref_filename else None	×
143	self.reference_used = bool(chgref_filename)	×
144
145	# List of nelects for each atom from potcar
146	potcar_indices = []	×
147	for i, v in enumerate(self.natoms):	×
148	potcar_indices += [i] * v	×
149	self.nelects = (	×
150	[self.potcar[potcar_indices[i]].nelectrons for i in range(len(self.structure))] if self.potcar else []
151	)
152
153	else:
154	fpath = os.path.abspath(cube_filename)	×
155	self.is_vasp = False	×
156	self.cube = VolumetricData.from_cube(fpath)	×
157	self.structure = self.cube.structure	×
158	self.nelects = None	×
159	chgrefpath = os.path.abspath(chgref_filename) if chgref_filename else None	×
160	self.reference_used = bool(chgref_filename)	×
161
162	with ScratchDir("."):	×
163	tmpfile = "CHGCAR" if chgcar_filename else "CUBE"	×
164	with zopen(fpath, "rt") as f_in:	×
165	with open(tmpfile, "w") as f_out:	×
166	shutil.copyfileobj(f_in, f_out)	×
167	args = [BADEREXE, tmpfile]	×
168	if chgref_filename:	×
169	with zopen(chgrefpath, "rt") as f_in:	×
170	with open("CHGCAR_ref", "w") as f_out:	×
171	shutil.copyfileobj(f_in, f_out)	×
172	args += ["-ref", "CHGCAR_ref"]	×
173	if parse_atomic_densities:	×
174	args += ["-p", "all_atom"]	×
175	with subprocess.Popen(args, stdout=subprocess.PIPE, stdin=subprocess.PIPE, close_fds=True) as rs:	×
176	stdout, _ = rs.communicate()	×
177	if rs.returncode != 0:	×
178	raise RuntimeError(	×
179	f"bader exited with return code {rs.returncode}. Please check your bader installation."
180	)
181
182	try:	×
183	self.version = float(stdout.split()[5])	×
184	except ValueError:	×
185	self.version = -1 # Unknown	×
186	if self.version < 1.0:	×
187	warnings.warn(	×
188	"Your installed version of Bader is outdated, calculation of vacuum charge may be incorrect.",
189	UserWarning,
190	)
191
192	data = []	×
193	with open("ACF.dat") as f:	×
194	raw = f.readlines()	×
195	headers = ("x", "y", "z", "charge", "min_dist", "atomic_vol")	×
196	raw.pop(0)	×
197	raw.pop(0)	×
198	while True:
199	l = raw.pop(0).strip()	×
200	if l.startswith("-"):	×
201	break	×
202	vals = map(float, l.split()[1:])	×
203	data.append(dict(zip(headers, vals)))	×
204	for l in raw:	×
205	toks = l.strip().split(":")	×
206	if toks[0] == "VACUUM CHARGE":	×
207	self.vacuum_charge = float(toks[1])	×
208	elif toks[0] == "VACUUM VOLUME":	×
209	self.vacuum_volume = float(toks[1])	×
210	elif toks[0] == "NUMBER OF ELECTRONS":	×
211	self.nelectrons = float(toks[1])	×
212	self.data = data	×
213
214	if self.parse_atomic_densities:	×
215	# convert the charge density for each atom spit out by Bader into Chgcar objects for easy parsing
216	atom_chgcars = [	×
217	Chgcar.from_file(f"BvAt{str(i).zfill(4)}.dat") for i in range(1, len(self.chgcar.structure) + 1)
218	]
219
220	atomic_densities = []	×
221	# For each atom in the structure
222	for _, loc, chg in zip(	×
223	self.chgcar.structure,
224	self.chgcar.structure.frac_coords,
225	atom_chgcars,
226	):
227	# Find the index of the atom in the charge density atom
228	index = np.round(np.multiply(loc, chg.dim))	×
229
230	data = chg.data["total"]	×
231	# Find the shift vector in the array
232	shift = (np.divide(chg.dim, 2) - index).astype(int)	×
233
234	# Shift the data so that the atomic charge density to the center for easier manipulation
235	shifted_data = np.roll(data, shift, axis=(0, 1, 2))	×
236
237	# Slices a central window from the data array
238	def slice_from_center(data, xwidth, ywidth, zwidth):	×
239	x, y, z = data.shape	×
240	startx = x // 2 - (xwidth // 2)	×
241	starty = y // 2 - (ywidth // 2)	×
242	startz = z // 2 - (zwidth // 2)	×
243	return data[	×
244	startx : startx + xwidth,
245	starty : starty + ywidth,
246	startz : startz + zwidth,
247	]
248
249	# Finds the central encompassing volume which holds all the data within a precision
250	def find_encompassing_vol(data):	×
251	total = np.sum(data)	×
252	for i in range(np.max(data.shape)):	×
253	sliced_data = slice_from_center(data, i, i, i)	×
254	if total - np.sum(sliced_data) < 0.1:	×
255	return sliced_data	×
256	return None	×
257
258	d = {	×
259	"data": find_encompassing_vol(shifted_data),
260	"shift": shift,
261	"dim": self.chgcar.dim,
262	}
263	atomic_densities.append(d)	×
264	self.atomic_densities = atomic_densities	×
265
266	def get_charge(self, atom_index):	1✔
267	"""
268	Convenience method to get the charge on a particular atom. This is the "raw"
269	charge generated by the Bader program, not a partial atomic charge. If the cube file
270	is a spin-density file, then this will return the spin density per atom with
271	positive being spin up and negative being spin down.
272
273	Args:
274	atom_index:
275	Index of atom.
276
277	Returns:
278	Charge associated with atom from the Bader analysis.
279	"""
280	return self.data[atom_index]["charge"]	×
281
282	def get_charge_transfer(self, atom_index, nelect=None):	1✔
283	"""
284	Returns the charge transferred for a particular atom. A positive value means
285	that the site has gained electron density (i.e. exhibits anionic character)
286	whereas a negative value means the site has lost electron density (i.e. exhibits
287	cationic character). If the arg nelect is not supplied, then POTCAR must be
288	supplied to determine nelect.
289
290	Args:
291	atom_index:
292	Index of atom.
293	nelect:
294	number of electrons associated with an isolated atom at this index.
295	For most DFT codes this corresponds to the number of valence electrons
296	associated with the pseudopotential (e.g. ZVAL for VASP).
297
298	Returns:
299	Charge transfer associated with atom from the Bader analysis.
300	Given by bader charge on atom - nelect for associated atom.
301	"""
302	if not self.nelects and nelect is None:	×
303	raise ValueError("No NELECT info! Need POTCAR for VASP or nelect argument for cube file")	×
304	return self.data[atom_index]["charge"] - (nelect if nelect is not None else self.nelects[atom_index])	×
305
306	def get_partial_charge(self, atom_index, nelect=None):	1✔
307	"""
308	Convenience method to get the partial charge on a particular atom. This is
309	simply the negative value of the charge transferred. A positive value indicates
310	that the atom has cationic character, whereas a negative value indicates the
311	site has anionic character.
312
313	Args:
314	atom_index:
315	Index of atom.
316	nelect:
317	number of electrons associated with an isolated atom at this index.
318	For most DFT codes this corresponds to the number of valence electrons
319	associated with the pseudopotential (e.g. ZVAL for VASP).
320
321	Returns:
322	Charge associated with atom from the Bader analysis.
323	"""
324	return -self.get_charge_transfer(atom_index, nelect)	×
325
326	def get_charge_decorated_structure(self):	1✔
327	"""
328	Returns a charge decorated structure
329
330	Note, this assumes that the Bader analysis was correctly performed on a file
331	with electron densities
332	"""
333	charges = [-self.get_charge(i) for i in range(len(self.structure))]	×
334	struct = self.structure.copy()	×
335	struct.add_site_property("charge", charges)	×
336	return struct	×
337
338	def get_oxidation_state_decorated_structure(self, nelects=None):	1✔
339	"""
340	Returns an oxidation state decorated structure based on bader analysis results.
341	Each site is assigned a charge based on the computed partial atomic charge from bader.
342
343	Note, this assumes that the Bader analysis was correctly performed on a file
344	with electron densities
345	"""
346	charges = [self.get_partial_charge(i, None if not nelects else nelects[i]) for i in range(len(self.structure))]	×
347	struct = self.structure.copy()	×
348	struct.add_oxidation_state_by_site(charges)	×
349	return struct	×
350
351	def get_decorated_structure(self, property_name, average=False):	1✔
352	"""
353	Get a property-decorated structure from the Bader analysis.
354
355	This is distinct from getting charge decorated structure, which assumes
356	the "standard" Bader analysis of electron densities followed by converting
357	electron count to charge. The expected way to use this is to call Bader on
358	a non-charge density file such as a spin density file, electrostatic potential
359	file, etc., while using the charge density file as the reference (chgref_filename)
360	so that the partitioning is determined via the charge, but averaging or integrating
361	is done for another property.
362
363	User warning: Bader analysis cannot automatically determine what property is
364	inside of the file. So if you want to use this for a non-conventional property
365	like spin, you must ensure that you have the file is for the appropriate
366	property and you have an appropriate reference file.
367
368	Args:
369	property_name: name of the property to assign to the structure, note that
370	if name is "spin" this is handled as a special case, and the appropriate
371	spin properties are set on the species in the structure
372	average: whether or not to return the average of this property, rather
373	than the total, by dividing by the atomic volume.
374
375	Returns:
376	structure with site properties assigned via Bader Analysis
377	"""
378	vals = [self.get_charge(i) for i in range(len(self.structure))]	×
379	struct = self.structure.copy()	×
380	if average:	×
381	vals = np.divide(vals, [d["atomic_vol"] for d in self.data])	×
382	struct.add_site_property(property_name, vals)	×
383	if property_name == "spin":	×
384	struct.add_spin_by_site(vals)	×
385	return struct	×
386
387	@property	1✔
388	def summary(self):	1✔
389	"""
390	:return: Dict summary of key analysis, e.g., atomic volume, charge, etc.
391	"""
392	summary = {	×
393	"min_dist": [d["min_dist"] for d in self.data],
394	"charge": [d["charge"] for d in self.data],
395	"atomic_volume": [d["atomic_vol"] for d in self.data],
396	"vacuum_charge": self.vacuum_charge,
397	"vacuum_volume": self.vacuum_volume,
398	"reference_used": self.reference_used,
399	"bader_version": self.version,
400	}
401
402	if self.parse_atomic_densities:	×
403	summary["charge_densities"] = self.atomic_densities	×
404
405	if self.potcar:	×
406	charge_transfer = [self.get_charge_transfer(i) for i in range(len(self.data))]	×
407	summary["charge_transfer"] = charge_transfer	×
408
409	return summary	×
410
411	@classmethod	1✔
412	def from_path(cls, path, suffix=""):	1✔
413	"""
414	Convenient constructor that takes in the path name of VASP run
415	to perform Bader analysis.
416
417	Args:
418	path (str): Name of directory where VASP output files are
419	stored.
420	suffix (str): specific suffix to look for (e.g. '.relax1'
421	for 'CHGCAR.relax1.gz').
422	"""
423
424	def _get_filepath(filename):	×
425	name_pattern = filename + suffix + "" if filename != "POTCAR" else filename + ""	×
426	paths = glob.glob(os.path.join(path, name_pattern))	×
427	fpath = None	×
428	if len(paths) >= 1:	×
429	# using reverse=True because, if multiple files are present,
430	# they likely have suffixes 'static', 'relax', 'relax2', etc.
431	# and this would give 'static' over 'relax2' over 'relax'
432	# however, better to use 'suffix' kwarg to avoid this!
433	paths.sort(reverse=True)	×
434	warning_msg = f"Multiple files detected, using {os.path.basename(paths[0])}" if len(paths) > 1 else None	×
435	fpath = paths[0]	×
436	else:
437	warning_msg = f"Could not find {filename}"	×
438	if filename in ["AECCAR0", "AECCAR2"]:	×
439	warning_msg += ", cannot calculate charge transfer."	×
440	elif filename == "POTCAR":	×
441	warning_msg += ", interpret Bader results with caution."	×
442	if warning_msg:	×
443	warnings.warn(warning_msg)	×
444	return fpath	×
445
446	chgcar_filename = _get_filepath("CHGCAR")	×
447	if chgcar_filename is None:	×
448	raise FileNotFoundError("Could not find CHGCAR!")	×
449	potcar_filename = _get_filepath("POTCAR")	×
450	aeccar0 = _get_filepath("AECCAR0")	×
451	aeccar2 = _get_filepath("AECCAR2")	×
452	if aeccar0 and aeccar2:	×
453	# `chgsum.pl AECCAR0 AECCAR2` equivalent to obtain chgref_file
454	chgref = Chgcar.from_file(aeccar0) + Chgcar.from_file(aeccar2)	×
455	chgref_filename = "CHGREF"	×
456	chgref.write_file(chgref_filename)	×
457	else:
458	chgref_filename = None	×
459	return cls(	×
460	chgcar_filename=chgcar_filename,
461	potcar_filename=potcar_filename,
462	chgref_filename=chgref_filename,
463	)
464
465
466	def bader_analysis_from_path(path, suffix=""):	1✔
467	"""
468	Convenience method to run Bader analysis on a folder containing
469	typical VASP output files.
470
471	This method will:
472
473	1. Look for files CHGCAR, AECCAR0, AECCAR2, POTCAR or their gzipped
474	counterparts.
475	2. If AECCAR* files are present, constructs a temporary reference
476	file as AECCAR0 + AECCAR2
477	3. Runs Bader analysis twice: once for charge, and a second time
478	for the charge difference (magnetization density).
479
480	:param path: path to folder to search in
481	:param suffix: specific suffix to look for (e.g. '.relax1' for 'CHGCAR.relax1.gz'
482	:return: summary dict
483	"""
484
485	def _get_filepath(filename, warning, path=path, suffix=suffix):	×
486	paths = glob.glob(os.path.join(path, filename + suffix + "*"))	×
487	if not paths:	×
488	warnings.warn(warning)	×
489	return None	×
490	if len(paths) > 1:	×
491	# using reverse=True because, if multiple files are present,
492	# they likely have suffixes 'static', 'relax', 'relax2', etc.
493	# and this would give 'static' over 'relax2' over 'relax'
494	# however, better to use 'suffix' kwarg to avoid this!
495	paths.sort(reverse=True)	×
496	warnings.warn(f"Multiple files detected, using {os.path.basename(path)}")	×
497	path = paths[0]	×
498	return path	×
499
500	chgcar_path = _get_filepath("CHGCAR", "Could not find CHGCAR!")	×
501	chgcar = Chgcar.from_file(chgcar_path)	×
502
503	aeccar0_path = _get_filepath("AECCAR0", "Could not find AECCAR0, interpret Bader results with caution.")	×
504	aeccar0 = Chgcar.from_file(aeccar0_path) if aeccar0_path else None	×
505
506	aeccar2_path = _get_filepath("AECCAR2", "Could not find AECCAR2, interpret Bader results with caution.")	×
507	aeccar2 = Chgcar.from_file(aeccar2_path) if aeccar2_path else None	×
508
509	potcar_path = _get_filepath("POTCAR", "Could not find POTCAR, cannot calculate charge transfer.")	×
510	potcar = Potcar.from_file(potcar_path) if potcar_path else None	×
511
512	return bader_analysis_from_objects(chgcar, potcar, aeccar0, aeccar2)	×
513
514
515	def bader_analysis_from_objects(chgcar, potcar=None, aeccar0=None, aeccar2=None):	1✔
516	"""
517	Convenience method to run Bader analysis from a set
518	of pymatgen Chgcar and Potcar objects.
519
520	This method will:
521
522	1. If aeccar objects are present, constructs a temporary reference
523	file as AECCAR0 + AECCAR2
524	2. Runs Bader analysis twice: once for charge, and a second time
525	for the charge difference (magnetization density).
526
527	:param chgcar: Chgcar object
528	:param potcar: (optional) Potcar object
529	:param aeccar0: (optional) Chgcar object from aeccar0 file
530	:param aeccar2: (optional) Chgcar object from aeccar2 file
531	:return: summary dict
532	"""
533
534	with ScratchDir(".") as temp_dir:	×
535	if aeccar0 and aeccar2:	×
536	# construct reference file
537	chgref = aeccar0.linear_add(aeccar2)	×
538	chgref_path = os.path.join(temp_dir, "CHGCAR_ref")	×
539	chgref.write_file(chgref_path)	×
540	else:
541	chgref_path = None	×
542
543	chgcar.write_file("CHGCAR")	×
544	chgcar_path = os.path.join(temp_dir, "CHGCAR")	×
545
546	if potcar:	×
547	potcar.write_file("POTCAR")	×
548	potcar_path = os.path.join(temp_dir, "POTCAR")	×
549	else:
550	potcar_path = None	×
551
552	ba = BaderAnalysis(	×
553	chgcar_filename=chgcar_path,
554	potcar_filename=potcar_path,
555	chgref_filename=chgref_path,
556	)
557
558	summary = {	×
559	"min_dist": [d["min_dist"] for d in ba.data],
560	"charge": [d["charge"] for d in ba.data],
561	"atomic_volume": [d["atomic_vol"] for d in ba.data],
562	"vacuum_charge": ba.vacuum_charge,
563	"vacuum_volume": ba.vacuum_volume,
564	"reference_used": bool(chgref_path),
565	"bader_version": ba.version,
566	}
567
568	if potcar:	×
569	charge_transfer = [ba.get_charge_transfer(i) for i in range(len(ba.data))]	×
570	summary["charge_transfer"] = charge_transfer	×
571
572	if chgcar.is_spin_polarized:	×
573	# write a CHGCAR containing magnetization density only
574	chgcar.data["total"] = chgcar.data["diff"]	×
575	chgcar.is_spin_polarized = False	×
576	chgcar.write_file("CHGCAR_mag")	×
577
578	chgcar_mag_path = os.path.join(temp_dir, "CHGCAR_mag")	×
579	ba = BaderAnalysis(	×
580	chgcar_filename=chgcar_mag_path,
581	potcar_filename=potcar_path,
582	chgref_filename=chgref_path,
583	)
584	summary["magmom"] = [d["charge"] for d in ba.data]	×
585
586	return summary	×

materialsproject / pymatgen / 4075885785

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous