6086389871

Committed 05 Sep 2023 03:02PM UTC coverage: 47.839% (+0.07%) from 47.768%

Build # 6086389871

Build Type

Pull #12

github-actions

Committed by

web-flow

Commit Message

Merge 84a3c20f6 into 0d9ba7ade

Pull Request Pull Request #12: Add Phonopy interface

Run Details

122 of 122 new or added lines in 2 files covered. (100.0%)

487 of 1018 relevant lines covered (47.84%)

0.48 hits per line

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

67.16

/atomistics/phonons/calculator.py

from typing import Optional
import posixpath

import numpy as np
import scipy.constants
from phonopy import Phonopy
from phonopy.units import VaspToTHz
from phonopy.file_IO import write_FORCE_CONSTANTS
import structuretoolkit

from atomistics.shared.calculator import Calculator
from atomistics.phonons.helper import (
    get_supercell_matrix,
    get_hesse_matrix,
    get_band_structure,
    plot_band_structure,
    plot_dos,
)


class PhonopyCalculator(Calculator):
    """
    Phonopy wrapper for the calculation of free energy in the framework of quasi harmonic approximation.

    Get output via `get_thermal_properties()`.

    Note:

    - This class does not consider the thermal expansion. For this, use `QuasiHarmonicJob` (find more in its
        docstring)
    - Depending on the value given in `job.input['interaction_range']`, this class automatically changes the number of
        atoms. The input parameters of the reference job might have to be set appropriately (e.g. use `k_mesh_density`
        for DFT instead of setting k-points directly).
    - The structure used in the reference job should be a relaxed structure.
    - Theory behind it: https://en.wikipedia.org/wiki/Quasi-harmonic_approximation
    """

    def __init__(
        self,
        structure,
        interaction_range=10,
        factor=VaspToTHz,
        displacement=0.01,
        dos_mesh=20,
        primitive_matrix=None,
        number_of_snapshots=None,
    ):
        self._interaction_range = interaction_range
        self._displacement = displacement
        self._dos_mesh = dos_mesh
        self._number_of_snapshots = number_of_snapshots
        self.structure = structure
        self._phonopy_unit_cell = structuretoolkit.common.atoms_to_phonopy(
            self.structure
        )
        self.phonopy = Phonopy(
            unitcell=self._phonopy_unit_cell,
            supercell_matrix=get_supercell_matrix(
                interaction_range=self._interaction_range,
                cell=self._phonopy_unit_cell.get_cell(),
            ),
            primitive_matrix=primitive_matrix,
            factor=factor,
        )
        self.phonopy.generate_displacements(
            distance=self._displacement,
            number_of_snapshots=self._number_of_snapshots,
        )
        self._dos_dict = {}
        self._mesh_dict = {}

    def generate_structures(self):
        return {
            ind: self._restore_magmoms(structuretoolkit.common.phonopy_to_atoms(sc))
            for ind, sc in enumerate(self.phonopy.supercells_with_displacements)
        }

    def _restore_magmoms(self, structure):
        """
        Args:
            structure (pyiron.atomistics.structure.atoms): input structure

        Returns:
            structure (pyiron_atomistics.atomistics.structure.atoms): output structure with magnetic moments
        """
        if self.structure.has("initial_magmoms"):
            magmoms = self.structure.get_initial_magnetic_moments()
            magmoms = np.tile(
                magmoms,
                np.prod(
                    np.diagonal(
                        get_supercell_matrix(
                            interaction_range=self._interaction_range,
                            cell=self._phonopy_unit_cell.get_cell(),
                        )
                    )
                ).astype(int),
            )
            structure.set_initial_magnetic_moments(magmoms)
        return structure

    def analyse_structures(self, output_dict):
        """

        Returns:

        """
        forces_lst = [output_dict[k] for k in sorted(output_dict.keys())]
        self.phonopy.forces = forces_lst
        self.phonopy.produce_force_constants(
            fc_calculator=None if self._number_of_snapshots is None else "alm"
        )
        self.phonopy.run_mesh(mesh=[self._dos_mesh] * 3)
        self._mesh_dict = self.phonopy.get_mesh_dict()
        self.phonopy.run_total_dos()
        self._dos_dict = self.phonopy.get_total_dos_dict()
        return self._mesh_dict, self._dos_dict

    def get_thermal_properties(self, t_min=1, t_max=1500, t_step=50, temperatures=None):
        """
        Returns thermal properties at constant volume in the given temperature range.  Can only be called after job
        successfully ran.

        Args:
            t_min (float): minimum sample temperature
            t_max (float): maximum sample temperature
            t_step (int):  tempeature sample interval
            temperatures (array_like, float):  custom array of temperature samples, if given t_min, t_max, t_step are
                                               ignored.

        Returns:
            :class:`Thermal`: thermal properties as returned by Phonopy
        """
        self.phonopy.run_thermal_properties(
            t_step=t_step, t_max=t_max, t_min=t_min, temperatures=temperatures
        )
        tp_dict = self.phonopy.get_thermal_properties_dict()
        kJ_mol_to_eV = 1000 / scipy.constants.Avogadro / scipy.constants.electron_volt
        tp_dict["free_energy"] *= kJ_mol_to_eV
        return tp_dict

    def get_dynamical_matrix(self):
        """

        Returns:

        """
        return np.real_if_close(self.phonopy.dynamical_matrix.dynamical_matrix)

    def dynamical_matrix_at_q(self, q):
        """

        Args:
            q:

        Returns:

        """
        return np.real_if_close(self.phonopy.get_dynamical_matrix_at_q(q))

    def write_phonopy_force_constants(self, file_name="FORCE_CONSTANTS", cwd=None):
        """

        Args:
            file_name:
            cwd:

        Returns:

        """
        if cwd is not None:
            file_name = posixpath.join(cwd, file_name)
        write_FORCE_CONSTANTS(
            force_constants=self.phonopy.force_constants, filename=file_name
        )

    def get_hesse_matrix(self):
        return get_hesse_matrix(force_constants=self.phonopy.force_constants)

    def get_band_structure(
        self, npoints=101, with_eigenvectors=False, with_group_velocities=False
    ):
        return get_band_structure(
            phonopy=self.phonopy,
            npoints=npoints,
            with_eigenvectors=with_eigenvectors,
            with_group_velocities=with_group_velocities,
        )

    def plot_band_structure(
        self, axis=None, *args, label: Optional[str] = None, **kwargs
    ):
        try:
            results = self.phonopy.get_band_structure_dict()
        except RuntimeError:
            results = self.get_band_structure()

        # HACK: strictly speaking this breaks phonopy API and could bite us
        path_connections = self.phonopy._band_structure.path_connections
        labels = self.phonopy._band_structure.labels
        return plot_band_structure(
            results=results,
            path_connections=path_connections,
            labels=labels,
            axis=axis,
            *args,
            label=label,
            **kwargs
        )

    def plot_dos(self, *args, axis=None, **kwargs):
        return plot_dos(
            dos_energies=self._dos_dict["frequency_points"],
            dos_total=self._dos_dict["total_dos"],
            *args,
            axis=axis,
            **kwargs
        )

1	from typing import Optional	1✔
2	import posixpath	1✔
3
4	import numpy as np	1✔
5	import scipy.constants	1✔
6	from phonopy import Phonopy	1✔
7	from phonopy.units import VaspToTHz	1✔
8	from phonopy.file_IO import write_FORCE_CONSTANTS	1✔
9	import structuretoolkit	1✔
10
11	from atomistics.shared.calculator import Calculator	1✔
12	from atomistics.phonons.helper import (	1✔
13	get_supercell_matrix,
14	get_hesse_matrix,
15	get_band_structure,
16	plot_band_structure,
17	plot_dos,
18	)
19
20
21	class PhonopyCalculator(Calculator):	1✔
22	"""
23	Phonopy wrapper for the calculation of free energy in the framework of quasi harmonic approximation.
24
25	Get output via `get_thermal_properties()`.
26
27	Note:
28
29	- This class does not consider the thermal expansion. For this, use `QuasiHarmonicJob` (find more in its
30	docstring)
31	- Depending on the value given in `job.input['interaction_range']`, this class automatically changes the number of
32	atoms. The input parameters of the reference job might have to be set appropriately (e.g. use `k_mesh_density`
33	for DFT instead of setting k-points directly).
34	- The structure used in the reference job should be a relaxed structure.
35	- Theory behind it: https://en.wikipedia.org/wiki/Quasi-harmonic_approximation
36	"""
37
38	def __init__(	1✔
39	self,
40	structure,
41	interaction_range=10,
42	factor=VaspToTHz,
43	displacement=0.01,
44	dos_mesh=20,
45	primitive_matrix=None,
46	number_of_snapshots=None,
47	):
48	self._interaction_range = interaction_range	1✔
49	self._displacement = displacement	1✔
50	self._dos_mesh = dos_mesh	1✔
51	self._number_of_snapshots = number_of_snapshots	1✔
52	self.structure = structure	1✔
53	self._phonopy_unit_cell = structuretoolkit.common.atoms_to_phonopy(	1✔
54	self.structure
55	)
56	self.phonopy = Phonopy(	1✔
57	unitcell=self._phonopy_unit_cell,
58	supercell_matrix=get_supercell_matrix(
59	interaction_range=self._interaction_range,
60	cell=self._phonopy_unit_cell.get_cell(),
61	),
62	primitive_matrix=primitive_matrix,
63	factor=factor,
64	)
65	self.phonopy.generate_displacements(	1✔
66	distance=self._displacement,
67	number_of_snapshots=self._number_of_snapshots,
68	)
69	self._dos_dict = {}	1✔
70	self._mesh_dict = {}	1✔
71
72	def generate_structures(self):	1✔
73	return {	1✔
74	ind: self._restore_magmoms(structuretoolkit.common.phonopy_to_atoms(sc))
75	for ind, sc in enumerate(self.phonopy.supercells_with_displacements)
76	}
77
78	def _restore_magmoms(self, structure):	1✔
79	"""
80	Args:
81	structure (pyiron.atomistics.structure.atoms): input structure
82
83	Returns:
84	structure (pyiron_atomistics.atomistics.structure.atoms): output structure with magnetic moments
85	"""
86	if self.structure.has("initial_magmoms"):	1✔
87	magmoms = self.structure.get_initial_magnetic_moments()	×
88	magmoms = np.tile(	×
89	magmoms,
90	np.prod(
91	np.diagonal(
92	get_supercell_matrix(
93	interaction_range=self._interaction_range,
94	cell=self._phonopy_unit_cell.get_cell(),
95	)
96	)
97	).astype(int),
98	)
99	structure.set_initial_magnetic_moments(magmoms)	×
100	return structure	1✔
101
102	def analyse_structures(self, output_dict):	1✔
103	"""
104
105	Returns:
106
107	"""
108	forces_lst = [output_dict[k] for k in sorted(output_dict.keys())]	1✔
109	self.phonopy.forces = forces_lst	1✔
110	self.phonopy.produce_force_constants(	1✔
111	fc_calculator=None if self._number_of_snapshots is None else "alm"
112	)
113	self.phonopy.run_mesh(mesh=[self._dos_mesh] * 3)	1✔
114	self._mesh_dict = self.phonopy.get_mesh_dict()	1✔
115	self.phonopy.run_total_dos()	1✔
116	self._dos_dict = self.phonopy.get_total_dos_dict()	1✔
117	return self._mesh_dict, self._dos_dict	1✔
118
119	def get_thermal_properties(self, t_min=1, t_max=1500, t_step=50, temperatures=None):	1✔
120	"""
121	Returns thermal properties at constant volume in the given temperature range. Can only be called after job
122	successfully ran.
123
124	Args:
125	t_min (float): minimum sample temperature
126	t_max (float): maximum sample temperature
127	t_step (int): tempeature sample interval
128	temperatures (array_like, float): custom array of temperature samples, if given t_min, t_max, t_step are
129	ignored.
130
131	Returns:
132	:class:`Thermal`: thermal properties as returned by Phonopy
133	"""
134	self.phonopy.run_thermal_properties(	×
135	t_step=t_step, t_max=t_max, t_min=t_min, temperatures=temperatures
136	)
137	tp_dict = self.phonopy.get_thermal_properties_dict()	×
138	kJ_mol_to_eV = 1000 / scipy.constants.Avogadro / scipy.constants.electron_volt	×
139	tp_dict["free_energy"] *= kJ_mol_to_eV	×
140	return tp_dict	×
141
142	def get_dynamical_matrix(self):	1✔
143	"""
144
145	Returns:
146
147	"""
148	return np.real_if_close(self.phonopy.dynamical_matrix.dynamical_matrix)	×
149
150	def dynamical_matrix_at_q(self, q):	1✔
151	"""
152
153	Args:
154	q:
155
156	Returns:
157
158	"""
159	return np.real_if_close(self.phonopy.get_dynamical_matrix_at_q(q))	×
160
161	def write_phonopy_force_constants(self, file_name="FORCE_CONSTANTS", cwd=None):	1✔
162	"""
163
164	Args:
165	file_name:
166	cwd:
167
168	Returns:
169
170	"""
171	if cwd is not None:	×
172	file_name = posixpath.join(cwd, file_name)	×
173	write_FORCE_CONSTANTS(	×
174	force_constants=self.phonopy.force_constants, filename=file_name
175	)
176
177	def get_hesse_matrix(self):	1✔
178	return get_hesse_matrix(force_constants=self.phonopy.force_constants)	1✔
179
180	def get_band_structure(	1✔
181	self, npoints=101, with_eigenvectors=False, with_group_velocities=False
182	):
183	return get_band_structure(	×
184	phonopy=self.phonopy,
185	npoints=npoints,
186	with_eigenvectors=with_eigenvectors,
187	with_group_velocities=with_group_velocities,
188	)
189
190	def plot_band_structure(	1✔
191	self, axis=None, args, label: Optional[str] = None, *kwargs
192	):
193	try:	×
194	results = self.phonopy.get_band_structure_dict()	×
195	except RuntimeError:	×
196	results = self.get_band_structure()	×
197
198	# HACK: strictly speaking this breaks phonopy API and could bite us
199	path_connections = self.phonopy._band_structure.path_connections	×
200	labels = self.phonopy._band_structure.labels	×
201	return plot_band_structure(	×
202	results=results,
203	path_connections=path_connections,
204	labels=labels,
205	axis=axis,
206	*args,
207	label=label,
208	**kwargs
209	)
210
211	def plot_dos(self, args, axis=None, *kwargs):	1✔
212	return plot_dos(	×
213	dos_energies=self._dos_dict["frequency_points"],
214	dos_total=self._dos_dict["total_dos"],
215	*args,
216	axis=axis,
217	**kwargs
218	)

pyiron / atomistics / 6086389871

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