diff --git a/atomistics/workflows/evcurve/workflow.py b/atomistics/workflows/evcurve/workflow.py index 0cdeb58b..38ae8568 100644 --- a/atomistics/workflows/evcurve/workflow.py +++ b/atomistics/workflows/evcurve/workflow.py @@ -1,17 +1,14 @@ import numpy as np from ase.atoms import Atoms from collections import OrderedDict -from typing import Optional, List from atomistics.shared.output import OutputEnergyVolumeCurve -from atomistics.workflows.evcurve.fit import EnergyVolumeFit from atomistics.workflows.interface import Workflow from atomistics.workflows.evcurve.debye import ( get_thermal_properties, OutputThermodynamic, ) from atomistics.workflows.evcurve.helper import ( - get_strains, get_volume_lst, generate_structures_helper, analyse_structures_helper, @@ -72,9 +69,6 @@ def analyse_structures( ) return self.fit_dict - def get_volume_lst(self) -> np.ndarray: - return get_volume_lst(structure_dict=self._structure_dict) - def get_thermal_properties( self, t_min: float = 1.0, @@ -94,10 +88,3 @@ def get_thermal_properties( constant_volume=constant_volume, output_keys=output_keys, ) - - def _get_strains(self): - return get_strains( - vol_range=self.vol_range, - num_points=self.num_points, - strain_lst=self.strains, - ) diff --git a/atomistics/workflows/phonons/helper.py b/atomistics/workflows/phonons/helper.py index 8b0215d3..5918e99a 100644 --- a/atomistics/workflows/phonons/helper.py +++ b/atomistics/workflows/phonons/helper.py @@ -1,7 +1,283 @@ from typing import Optional + +from ase.atoms import Atoms import numpy as np -import phonopy +from phonopy import Phonopy import scipy.constants +import structuretoolkit + +from atomistics.shared.output import OutputThermodynamic, OutputPhonons +from atomistics.workflows.phonons.units import VaspToTHz, kJ_mol_to_eV + + +class PhonopyProperties(object): + def __init__( + self, + phonopy_instance: Phonopy, + dos_mesh: np.ndarray, + shift=None, + is_time_reversal: bool = True, + is_mesh_symmetry: bool = True, + with_eigenvectors: bool = False, + with_group_velocities: bool = False, + is_gamma_center: bool = False, + number_of_snapshots: int = None, + sigma: float = None, + freq_min: float = None, + freq_max: float = None, + freq_pitch: float = None, + use_tetrahedron_method: bool = True, + npoints: int = 101, + ): + self._phonopy = phonopy_instance + self._sigma = sigma + self._freq_min = freq_min + self._freq_max = freq_max + self._freq_pitch = freq_pitch + self._use_tetrahedron_method = use_tetrahedron_method + self._npoints = npoints + self._with_eigenvectors = with_eigenvectors + self._with_group_velocities = with_group_velocities + self._dos_mesh = dos_mesh + self._shift = shift + self._is_time_reversal = is_time_reversal + self._is_mesh_symmetry = is_mesh_symmetry + self._with_eigenvectors = with_eigenvectors + self._with_group_velocities = with_group_velocities + self._is_gamma_center = is_gamma_center + self._number_of_snapshots = number_of_snapshots + self._total_dos = None + self._band_structure_dict = None + self._mesh_dict = None + self._force_constants = None + + def _calc_band_structure(self): + self._phonopy.auto_band_structure( + npoints=self._npoints, + with_eigenvectors=self._with_eigenvectors, + with_group_velocities=self._with_group_velocities, + ) + self._band_structure_dict = self._phonopy.get_band_structure_dict() + + def _calc_force_constants(self): + self._phonopy.produce_force_constants( + fc_calculator=None if self._number_of_snapshots is None else "alm" + ) + self._force_constants = self._phonopy.force_constants + + def mesh_dict(self) -> dict: + if self._force_constants is None: + self._calc_force_constants() + if self._mesh_dict is None: + self._phonopy.run_mesh( + mesh=[self._dos_mesh] * 3, + shift=self._shift, + is_time_reversal=self._is_time_reversal, + is_mesh_symmetry=self._is_mesh_symmetry, + with_eigenvectors=self._with_eigenvectors, + with_group_velocities=self._with_group_velocities, + is_gamma_center=self._is_gamma_center, + ) + self._mesh_dict = self._phonopy.get_mesh_dict() + return self._mesh_dict + + def band_structure_dict(self) -> dict: + if self._band_structure_dict is None: + self._calc_band_structure() + return self._band_structure_dict + + def total_dos_dict(self) -> dict: + if self._total_dos is None: + self._phonopy.run_total_dos( + sigma=self._sigma, + freq_min=self._freq_min, + freq_max=self._freq_max, + freq_pitch=self._freq_pitch, + use_tetrahedron_method=self._use_tetrahedron_method, + ) + self._total_dos = self._phonopy.get_total_dos_dict() + return self._total_dos + + def dynamical_matrix(self) -> np.ndarray: + if self._band_structure_dict is None: + self._calc_band_structure() + return self._phonopy.dynamical_matrix.dynamical_matrix + + def force_constants(self) -> np.ndarray: + if self._force_constants is None: + self._calc_force_constants() + return self._force_constants + + +class PhonopyThermalProperties(object): + def __init__(self, phonopy_instance: Phonopy): + self._phonopy = phonopy_instance + self._thermal_properties = phonopy_instance.get_thermal_properties_dict() + + def free_energy(self) -> np.ndarray: + return self._thermal_properties["free_energy"] * kJ_mol_to_eV + + def temperatures(self) -> np.ndarray: + return self._thermal_properties["temperatures"] + + def entropy(self) -> np.ndarray: + return self._thermal_properties["entropy"] + + def heat_capacity(self) -> np.ndarray: + return self._thermal_properties["heat_capacity"] + + def volumes(self) -> np.ndarray: + return np.array( + [self._phonopy.unitcell.get_volume()] + * len(self._thermal_properties["temperatures"]) + ) + + +def restore_magmoms( + structure_with_magmoms: Atoms, + structure: Atoms, + interaction_range: float, + cell: np.ndarray, +) -> Atoms: + """ + Args: + structure_with_magmoms (ase.atoms.Atoms): input structure with magnetic moments + structure (ase.atoms.Atoms): input structure without magnetic moments + interaction_range (float): + cell (np.ndarray): + + Returns: + structure (ase.atoms.Atoms): output structure with magnetic moments + """ + if structure_with_magmoms.has("initial_magmoms"): + magmoms = structure_with_magmoms.get_initial_magnetic_moments() + magmoms = np.tile( + magmoms, + np.prod( + np.diagonal( + get_supercell_matrix( + interaction_range=interaction_range, + cell=cell, + ) + ) + ).astype(int), + ) + structure.set_initial_magnetic_moments(magmoms) + return structure + + +def generate_structures_helper( + structure: Atoms, + primitive_matrix: np.ndarray = None, + displacement: float = 0.01, + number_of_snapshots: int = None, + interaction_range: float = 10.0, + factor: float = VaspToTHz, +): + unitcell = structuretoolkit.common.atoms_to_phonopy(structure) + phonopy_obj = Phonopy( + unitcell=unitcell, + supercell_matrix=get_supercell_matrix( + interaction_range=interaction_range, + cell=unitcell.get_cell(), + ), + primitive_matrix=primitive_matrix, + factor=factor, + ) + phonopy_obj.generate_displacements( + distance=displacement, + number_of_snapshots=number_of_snapshots, + ) + structure_dict = { + ind: restore_magmoms( + structure_with_magmoms=structure, + structure=structuretoolkit.common.phonopy_to_atoms(sc), + interaction_range=interaction_range, + cell=unitcell.get_cell(), + ) + for ind, sc in enumerate(phonopy_obj.supercells_with_displacements) + } + return phonopy_obj, structure_dict + + +def analyse_structures_helper( + phonopy: Phonopy, + output_dict: dict, + dos_mesh: int = 20, + number_of_snapshots: int = None, + output_keys: tuple[str] = OutputPhonons.keys(), +) -> dict: + """ + + Returns: + + """ + if "forces" in output_dict.keys(): + output_dict = output_dict["forces"] + forces_lst = [output_dict[k] for k in sorted(output_dict.keys())] + phonopy.forces = forces_lst + phono = PhonopyProperties( + phonopy_instance=phonopy, + dos_mesh=dos_mesh, + shift=None, + is_time_reversal=True, + is_mesh_symmetry=True, + with_eigenvectors=False, + with_group_velocities=False, + is_gamma_center=False, + number_of_snapshots=number_of_snapshots, + sigma=None, + freq_min=None, + freq_max=None, + freq_pitch=None, + use_tetrahedron_method=True, + npoints=101, + ) + return OutputPhonons(**{k: getattr(phono, k) for k in OutputPhonons.keys()}).get( + output_keys=output_keys + ) + + +def get_thermal_properties( + phonopy: Phonopy, + t_min: float = 1.0, + t_max: float = 1500.0, + t_step: float = 50.0, + temperatures: np.ndarray = None, + cutoff_frequency: float = None, + pretend_real: bool = False, + band_indices: np.ndarray = None, + is_projection: bool = False, + output_keys: tuple[str] = OutputThermodynamic.keys(), +) -> dict: + """ + 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 + """ + phonopy.run_thermal_properties( + t_step=t_step, + t_max=t_max, + t_min=t_min, + temperatures=temperatures, + cutoff_frequency=cutoff_frequency, + pretend_real=pretend_real, + band_indices=band_indices, + is_projection=is_projection, + ) + phono = PhonopyThermalProperties(phonopy_instance=phonopy) + return OutputThermodynamic( + **{k: getattr(phono, k) for k in OutputThermodynamic.keys()} + ).get(output_keys=output_keys) def get_supercell_matrix(interaction_range: float, cell: np.ndarray) -> np.ndarray: @@ -62,7 +338,7 @@ def plot_dos( def get_band_structure( - phonopy: phonopy.Phonopy, + phonopy: Phonopy, npoints: int = 101, with_eigenvectors: bool = False, with_group_velocities: bool = False, diff --git a/atomistics/workflows/phonons/workflow.py b/atomistics/workflows/phonons/workflow.py index 7ac8f1ac..7033d68b 100644 --- a/atomistics/workflows/phonons/workflow.py +++ b/atomistics/workflows/phonons/workflow.py @@ -3,144 +3,20 @@ from ase.atoms import Atoms import numpy as np -import phonopy -from phonopy import Phonopy from phonopy.file_IO import write_FORCE_CONSTANTS -import structuretoolkit from atomistics.shared.output import OutputThermodynamic, OutputPhonons from atomistics.workflows.interface import Workflow from atomistics.workflows.phonons.helper import ( - get_supercell_matrix, + analyse_structures_helper, + generate_structures_helper, + get_thermal_properties, get_hesse_matrix, get_band_structure, plot_band_structure, plot_dos, ) -from atomistics.workflows.phonons.units import VaspToTHz, kJ_mol_to_eV - - -class PhonopyProperties(object): - def __init__( - self, - phonopy_instance: phonopy.Phonopy, - dos_mesh: np.ndarray, - shift=None, - is_time_reversal: bool = True, - is_mesh_symmetry: bool = True, - with_eigenvectors: bool = False, - with_group_velocities: bool = False, - is_gamma_center: bool = False, - number_of_snapshots: int = None, - sigma: float = None, - freq_min: float = None, - freq_max: float = None, - freq_pitch: float = None, - use_tetrahedron_method: bool = True, - npoints: int = 101, - ): - self._phonopy = phonopy_instance - self._sigma = sigma - self._freq_min = freq_min - self._freq_max = freq_max - self._freq_pitch = freq_pitch - self._use_tetrahedron_method = use_tetrahedron_method - self._npoints = npoints - self._with_eigenvectors = with_eigenvectors - self._with_group_velocities = with_group_velocities - self._dos_mesh = dos_mesh - self._shift = shift - self._is_time_reversal = is_time_reversal - self._is_mesh_symmetry = is_mesh_symmetry - self._with_eigenvectors = with_eigenvectors - self._with_group_velocities = with_group_velocities - self._is_gamma_center = is_gamma_center - self._number_of_snapshots = number_of_snapshots - self._total_dos = None - self._band_structure_dict = None - self._mesh_dict = None - self._force_constants = None - - def _calc_band_structure(self): - self._phonopy.auto_band_structure( - npoints=self._npoints, - with_eigenvectors=self._with_eigenvectors, - with_group_velocities=self._with_group_velocities, - ) - self._band_structure_dict = self._phonopy.get_band_structure_dict() - - def _calc_force_constants(self): - self._phonopy.produce_force_constants( - fc_calculator=None if self._number_of_snapshots is None else "alm" - ) - self._force_constants = self._phonopy.force_constants - - def mesh_dict(self) -> dict: - if self._force_constants is None: - self._calc_force_constants() - if self._mesh_dict is None: - self._phonopy.run_mesh( - mesh=[self._dos_mesh] * 3, - shift=self._shift, - is_time_reversal=self._is_time_reversal, - is_mesh_symmetry=self._is_mesh_symmetry, - with_eigenvectors=self._with_eigenvectors, - with_group_velocities=self._with_group_velocities, - is_gamma_center=self._is_gamma_center, - ) - self._mesh_dict = self._phonopy.get_mesh_dict() - return self._mesh_dict - - def band_structure_dict(self) -> dict: - if self._band_structure_dict is None: - self._calc_band_structure() - return self._band_structure_dict - - def total_dos_dict(self) -> dict: - if self._total_dos is None: - self._phonopy.run_total_dos( - sigma=self._sigma, - freq_min=self._freq_min, - freq_max=self._freq_max, - freq_pitch=self._freq_pitch, - use_tetrahedron_method=self._use_tetrahedron_method, - ) - self._total_dos = self._phonopy.get_total_dos_dict() - return self._total_dos - - def dynamical_matrix(self) -> np.ndarray: - if self._band_structure_dict is None: - self._calc_band_structure() - return self._phonopy.dynamical_matrix.dynamical_matrix - - def force_constants(self) -> np.ndarray: - if self._force_constants is None: - self._calc_force_constants() - return self._force_constants - - -class PhonopyThermalProperties(object): - def __init__(self, phonopy_instance: phonopy.Phonopy): - self._phonopy = phonopy_instance - self._thermal_properties = phonopy_instance.get_thermal_properties_dict() - - def free_energy(self) -> np.ndarray: - return self._thermal_properties["free_energy"] * kJ_mol_to_eV - - def temperatures(self) -> np.ndarray: - return self._thermal_properties["temperatures"] - - def entropy(self) -> np.ndarray: - return self._thermal_properties["entropy"] - - def heat_capacity(self) -> np.ndarray: - return self._thermal_properties["heat_capacity"] - - def volumes(self) -> np.ndarray: - return np.array( - [self._phonopy.unitcell.get_volume()] - * len(self._thermal_properties["temperatures"]) - ) +from atomistics.workflows.phonons.units import VaspToTHz class PhonopyWorkflow(Workflow): @@ -175,55 +51,21 @@ def __init__( 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._primitive_matrix = primitive_matrix + self._factor = factor + self.phonopy = None self._phonopy_dict = {} def generate_structures(self) -> dict: - self.phonopy.generate_displacements( - distance=self._displacement, + self.phonopy, structure_dict = generate_structures_helper( + structure=self.structure, + primitive_matrix=self._primitive_matrix, + displacement=self._displacement, number_of_snapshots=self._number_of_snapshots, + interaction_range=self._interaction_range, + factor=self._factor, ) - return { - "calc_forces": { - 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: Atoms) -> Atoms: - """ - Args: - structure (ase.atoms.Atoms): input structure - - Returns: - structure (ase.atoms.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 + return {"calc_forces": structure_dict} def analyse_structures( self, output_dict: dict, output_keys: tuple[str] = OutputPhonons.keys() @@ -233,30 +75,13 @@ def analyse_structures( Returns: """ - if "forces" in output_dict.keys(): - output_dict = output_dict["forces"] - forces_lst = [output_dict[k] for k in sorted(output_dict.keys())] - self.phonopy.forces = forces_lst - phono = PhonopyProperties( - phonopy_instance=self.phonopy, + self._phonopy_dict = analyse_structures_helper( + phonopy=self.phonopy, + output_dict=output_dict, dos_mesh=self._dos_mesh, - shift=None, - is_time_reversal=True, - is_mesh_symmetry=True, - with_eigenvectors=False, - with_group_velocities=False, - is_gamma_center=False, number_of_snapshots=self._number_of_snapshots, - sigma=None, - freq_min=None, - freq_max=None, - freq_pitch=None, - use_tetrahedron_method=True, - npoints=101, + output_keys=output_keys, ) - self._phonopy_dict = OutputPhonons( - **{k: getattr(phono, k) for k in OutputPhonons.keys()} - ).get(output_keys=output_keys) return self._phonopy_dict def get_thermal_properties( @@ -285,20 +110,18 @@ def get_thermal_properties( Returns: :class:`Thermal`: thermal properties as returned by Phonopy """ - self.phonopy.run_thermal_properties( - t_step=t_step, - t_max=t_max, + return get_thermal_properties( + phonopy=self.phonopy, t_min=t_min, + t_max=t_max, + t_step=t_step, temperatures=temperatures, cutoff_frequency=cutoff_frequency, pretend_real=pretend_real, band_indices=band_indices, is_projection=is_projection, + output_keys=output_keys, ) - phono = PhonopyThermalProperties(phonopy_instance=self.phonopy) - return OutputThermodynamic( - **{k: getattr(phono, k) for k in OutputThermodynamic.keys()} - ).get(output_keys=output_keys) def get_dynamical_matrix(self, npoints: int = 101) -> np.ndarray: """ diff --git a/atomistics/workflows/quasiharmonic.py b/atomistics/workflows/quasiharmonic.py index cc982555..856bd39a 100644 --- a/atomistics/workflows/quasiharmonic.py +++ b/atomistics/workflows/quasiharmonic.py @@ -1,14 +1,22 @@ +from typing import Optional, List + from ase.atoms import Atoms import numpy as np from atomistics.shared.output import OutputThermodynamic from atomistics.workflows.evcurve.workflow import EnergyVolumeCurveWorkflow from atomistics.workflows.evcurve.helper import ( + get_strains, + get_volume_lst, fit_ev_curve, _strain_axes, ) -from atomistics.workflows.phonons.workflow import PhonopyWorkflow -from atomistics.workflows.phonons.helper import get_supercell_matrix +from atomistics.workflows.phonons.helper import ( + get_supercell_matrix, + analyse_structures_helper as analyse_structures_phonopy_helper, + get_thermal_properties as get_thermal_properties_phonopy, + generate_structures_helper as generate_structures_phonopy_helper, +) from atomistics.workflows.phonons.units import ( VaspToTHz, kJ_mol_to_eV, @@ -27,7 +35,8 @@ def get_free_energy_classical( def get_thermal_properties( eng_internal_dict: dict, phonopy_dict: dict, - volume_lst: np.ndarray, + structure_dict: dict, + repeat_vector: np.ndarray, fit_type: str, fit_order: int, t_min: float = 1.0, @@ -55,6 +64,9 @@ def get_thermal_properties( Returns: :class:`Thermal`: thermal properties as returned by Phonopy """ + volume_lst = np.array(get_volume_lst(structure_dict=structure_dict)) / np.prod( + repeat_vector + ) if quantum_mechanical: tp_collect_dict = _get_thermal_properties_quantum_mechanical( phonopy_dict=phonopy_dict, @@ -152,7 +164,8 @@ def _get_thermal_properties_quantum_mechanical( :class:`Thermal`: thermal properties as returned by Phonopy """ return { - strain: phono.get_thermal_properties( + strain: get_thermal_properties_phonopy( + phonopy=phono, t_step=t_step, t_max=t_max, t_min=t_min, @@ -200,9 +213,7 @@ def _get_thermal_properties_classical( t_property_lst = [] for t in temperatures: t_property = 0.0 - for freqs, w in zip( - phono.phonopy.mesh.frequencies, phono.phonopy.mesh.weights - ): + for freqs, w in zip(phono.mesh.frequencies, phono.mesh.weights): freqs = np.array(freqs) * THzToEv cond = freqs > cutoff_frequency t_property += ( @@ -211,9 +222,7 @@ def _get_thermal_properties_classical( ) * w ) - t_property_lst.append( - t_property / np.sum(phono.phonopy.mesh.weights) * EvTokJmol - ) + t_property_lst.append(t_property / np.sum(phono.mesh.weights) * EvTokJmol) tp_collect_dict[strain] = { "temperatures": temperatures, "free_energy": np.array(t_property_lst) * kJ_mol_to_eV, @@ -268,6 +277,79 @@ def volumes(self) -> np.ndarray: return self._volumes_selected_lst +def generate_structures_helper( + structure: Atoms, + vol_range: Optional[float] = None, + num_points: Optional[int] = None, + strain_lst: Optional[List[float]] = None, + displacement: float = 0.01, + number_of_snapshots: int = None, + interaction_range: float = 10.0, + factor: float = VaspToTHz, +): + # Phonopy internally repeats structures that are "too small" + # Here we manually guarantee that all structures passed are big enough + # This provides some computational efficiency for classical calculations + # And for quantum calculations _ensures_ that force matrices and energy/atom + # get treated with the same kmesh + repeat_vector = np.array( + np.diag( + get_supercell_matrix( + interaction_range=interaction_range, + cell=structure.cell.array, + ) + ), + dtype=int, + ) + structure_energy_dict, structure_forces_dict, phonopy_dict = {}, {}, {} + for strain in get_strains( + vol_range=vol_range, + num_points=num_points, + strain_lst=strain_lst, + ): + strain_ind = 1 + np.round(strain, 7) + basis = _strain_axes( + structure=structure, axes=("x", "y", "z"), volume_strain=strain + ) + structure_energy_dict[strain_ind] = basis.repeat(repeat_vector) + phonopy_obj, structure_task_dict = generate_structures_phonopy_helper( + structure=basis, + displacement=displacement, + number_of_snapshots=number_of_snapshots, + interaction_range=interaction_range, + factor=factor, + ) + phonopy_dict[strain_ind] = phonopy_obj + structure_forces_dict.update( + { + (strain_ind, key): structure_phono + for key, structure_phono in structure_task_dict.items() + } + ) + return phonopy_dict, repeat_vector, structure_energy_dict, structure_forces_dict + + +def analyse_structures_helper( + phonopy_dict: dict, + output_dict: dict, + dos_mesh: int = 20, + number_of_snapshots: int = None, + output_keys: tuple[str] = ("force_constants", "mesh_dict"), +): + eng_internal_dict = output_dict["energy"] + phonopy_collect_dict = { + strain: analyse_structures_phonopy_helper( + phonopy=phono, + output_dict={k: v for k, v in output_dict["forces"].items() if strain in k}, + dos_mesh=dos_mesh, + number_of_snapshots=number_of_snapshots, + output_keys=output_keys, + ) + for strain, phono in phonopy_dict.items() + } + return eng_internal_dict, phonopy_collect_dict + + class QuasiHarmonicWorkflow(EnergyVolumeCurveWorkflow): def __init__( self, @@ -300,66 +382,42 @@ def __init__( self._primitive_matrix = primitive_matrix self._phonopy_dict = {} self._eng_internal_dict = None - # Phonopy internally repeats structures that are "too small" - # Here we manually guarantee that all structures passed are big enough - # This provides some computational efficiency for classical calculations - # And for quantum calculations _ensures_ that force matrices and energy/atom - # get treated with the same kmesh - self._repeat_vector = np.array( - np.diag( - get_supercell_matrix( - interaction_range=interaction_range, - cell=structure.cell.array, - ) - ), - dtype=int, - ) + self._repeat_vector = None def generate_structures(self) -> dict: - task_dict = {"calc_forces": {}} - for strain in self._get_strains(): - strain_ind = 1 + np.round(strain, 7) - basis = _strain_axes( - structure=self.structure, axes=self.axes, volume_strain=strain - ) - structure_ev = basis.repeat(self._repeat_vector) - self._structure_dict[strain_ind] = structure_ev - self._phonopy_dict[strain_ind] = PhonopyWorkflow( - structure=basis, - interaction_range=self._interaction_range, - factor=self._factor, - displacement=self._displacement, - dos_mesh=self._dos_mesh, - primitive_matrix=self._primitive_matrix, - number_of_snapshots=self._number_of_snapshots, - ) - structure_task_dict = self._phonopy_dict[strain_ind].generate_structures() - task_dict["calc_forces"].update( - { - (strain_ind, key): structure_phono - for key, structure_phono in structure_task_dict[ - "calc_forces" - ].items() - } - ) - task_dict["calc_energy"] = self._structure_dict - return task_dict + ( + self._phonopy_dict, + self._repeat_vector, + structure_energy_dict, + structure_forces_dict, + ) = generate_structures_helper( + structure=self.structure, + vol_range=self.vol_range, + num_points=self.num_points, + strain_lst=self.strains, + displacement=self._displacement, + number_of_snapshots=self._number_of_snapshots, + interaction_range=self._interaction_range, + factor=self._factor, + ) + self._structure_dict = structure_energy_dict + return { + "calc_energy": structure_energy_dict, + "calc_forces": structure_forces_dict, + } def analyse_structures( self, output_dict: dict, output_keys: tuple[str] = ("force_constants", "mesh_dict"), ): - self._eng_internal_dict = output_dict["energy"] - phonopy_collect_dict = { - strain: phono.analyse_structures( - output_dict={ - k: v for k, v in output_dict["forces"].items() if strain in k - }, - output_keys=output_keys, - ) - for strain, phono in self._phonopy_dict.items() - } + self._eng_internal_dict, phonopy_collect_dict = analyse_structures_helper( + phonopy_dict=self._phonopy_dict, + output_dict=output_dict, + dos_mesh=self._dos_mesh, + number_of_snapshots=self._number_of_snapshots, + output_keys=output_keys, + ) return self._eng_internal_dict, phonopy_collect_dict def get_thermal_properties( @@ -396,7 +454,8 @@ def get_thermal_properties( return get_thermal_properties( eng_internal_dict=self._eng_internal_dict, phonopy_dict=self._phonopy_dict, - volume_lst=np.array(self.get_volume_lst()) / np.prod(self._repeat_vector), + structure_dict=self._structure_dict, + repeat_vector=self._repeat_vector, fit_type=self.fit_type, fit_order=self.fit_order, t_min=t_min, diff --git a/tests/test_phonons_lammps_functional.py b/tests/test_phonons_lammps_functional.py new file mode 100644 index 00000000..5bc89d1b --- /dev/null +++ b/tests/test_phonons_lammps_functional.py @@ -0,0 +1,114 @@ +import os + +from ase.build import bulk +from phonopy.units import VaspToTHz +import unittest + +from atomistics.workflows.phonons.helper import ( + get_hesse_matrix, + get_thermal_properties, + generate_structures_helper, + analyse_structures_helper, +) + +try: + from atomistics.calculators import ( + evaluate_with_lammps, get_potential_by_name + ) + + skip_lammps_test = False +except ImportError: + skip_lammps_test = True + + +@unittest.skipIf( + skip_lammps_test, "LAMMPS is not installed, so the LAMMPS tests are skipped." +) +class TestPhonons(unittest.TestCase): + def test_calc_phonons(self): + structure = bulk("Al", cubic=True) + df_pot_selected = get_potential_by_name( + potential_name='1999--Mishin-Y--Al--LAMMPS--ipr1', + resource_path=os.path.join(os.path.dirname(__file__), "static", "lammps"), + ) + result_dict = evaluate_with_lammps( + task_dict={"optimize_positions_and_volume": structure}, + potential_dataframe=df_pot_selected, + ) + phonopy_obj, structure_dict = generate_structures_helper( + structure=result_dict["structure_with_optimized_positions_and_volume"], + primitive_matrix=None, + number_of_snapshots=None, + displacement=0.01, + interaction_range=10.0, + factor=VaspToTHz, + ) + result_dict = evaluate_with_lammps( + task_dict={"calc_forces": structure_dict}, + potential_dataframe=df_pot_selected, + ) + phonopy_dict = analyse_structures_helper( + phonopy=phonopy_obj, + output_dict=result_dict, + dos_mesh=20, + number_of_snapshots=None, + ) + mesh_dict, dos_dict = phonopy_dict["mesh_dict"], phonopy_dict["total_dos_dict"] + self.assertEqual((324, 324), get_hesse_matrix(force_constants=phonopy_obj.force_constants).shape) + self.assertTrue('qpoints' in mesh_dict.keys()) + self.assertTrue('weights' in mesh_dict.keys()) + self.assertTrue('frequencies' in mesh_dict.keys()) + self.assertTrue('eigenvectors' in mesh_dict.keys()) + self.assertTrue('group_velocities' in mesh_dict.keys()) + self.assertTrue('frequency_points' in dos_dict.keys()) + self.assertTrue('total_dos' in dos_dict.keys()) + thermal_dict = get_thermal_properties( + phonopy=phonopy_obj, + t_min=1, + t_max=1500, + t_step=50, + temperatures=None, + cutoff_frequency=None, + pretend_real=False, + band_indices=None, + is_projection=False, + ) + for key in ["temperatures", "free_energy", "volumes", "entropy", "heat_capacity"]: + self.assertTrue(len(thermal_dict[key]), 31) + self.assertEqual(thermal_dict["temperatures"][0], 1.0) + self.assertEqual(thermal_dict["temperatures"][-1], 1501.0) + self.assertTrue(thermal_dict["free_energy"][0] < 0.2) + self.assertTrue(thermal_dict["free_energy"][0] > 0.1) + self.assertTrue(thermal_dict["free_energy"][-1] < -2.6) + self.assertTrue(thermal_dict["free_energy"][-1] > -2.7) + self.assertTrue(thermal_dict["entropy"][0] < 0.1) + self.assertTrue(thermal_dict["entropy"][0] > 0.0) + self.assertTrue(thermal_dict["entropy"][-1] < 271) + self.assertTrue(thermal_dict["entropy"][-1] > 270) + self.assertTrue(thermal_dict["heat_capacity"][0] < 0.1) + self.assertTrue(thermal_dict["heat_capacity"][0] > 0.0) + self.assertTrue(thermal_dict["heat_capacity"][-1] < 100) + self.assertTrue(thermal_dict["heat_capacity"][-1] > 99) + self.assertTrue(thermal_dict["volumes"][-1] < 66.5) + self.assertTrue(thermal_dict["volumes"][-1] > 66.4) + self.assertTrue(thermal_dict["volumes"][0] < 66.5) + self.assertTrue(thermal_dict["volumes"][0] > 66.4) + thermal_dict = get_thermal_properties( + phonopy=phonopy_obj, + t_min=1, + t_max=1500, + t_step=50, + temperatures=None, + cutoff_frequency=None, + pretend_real=False, + band_indices=None, + is_projection=False, + output_keys=["temperatures", "free_energy"] + ) + self.assertEqual(len(thermal_dict.keys()), 2) + self.assertEqual(thermal_dict["temperatures"][0], 1.0) + self.assertEqual(thermal_dict["temperatures"][-1], 1501.0) + self.assertTrue(thermal_dict["free_energy"][0] < 0.2) + self.assertTrue(thermal_dict["free_energy"][0] > 0.1) + self.assertTrue(thermal_dict["free_energy"][-1] < -2.6) + self.assertTrue(thermal_dict["free_energy"][-1] > -2.7) diff --git a/tests/test_quasiharmonic_lammps_functional.py b/tests/test_quasiharmonic_lammps_functional.py new file mode 100644 index 00000000..6220fa76 --- /dev/null +++ b/tests/test_quasiharmonic_lammps_functional.py @@ -0,0 +1,121 @@ +import os + +from ase.build import bulk +import numpy as np +from phonopy.units import VaspToTHz +import unittest + +from atomistics.workflows.quasiharmonic import ( + generate_structures_helper, + analyse_structures_helper, + get_thermal_properties, +) + +try: + from atomistics.calculators import ( + evaluate_with_lammps, get_potential_by_name + ) + + skip_lammps_test = False +except ImportError: + skip_lammps_test = True + + +@unittest.skipIf( + skip_lammps_test, "LAMMPS is not installed, so the LAMMPS tests are skipped." +) +class TestPhonons(unittest.TestCase): + def test_calc_phonons(self): + structure = bulk("Al", cubic=True) + df_pot_selected = get_potential_by_name( + potential_name='1999--Mishin-Y--Al--LAMMPS--ipr1', + resource_path=os.path.join(os.path.dirname(__file__), "static", "lammps"), + ) + result_dict = evaluate_with_lammps( + task_dict={"optimize_positions_and_volume": structure}, + potential_dataframe=df_pot_selected, + ) + phonopy_dict, repeat_vector, structure_energy_dict, structure_forces_dict = generate_structures_helper( + structure=result_dict["structure_with_optimized_positions_and_volume"], + vol_range=0.05, + num_points=11, + strain_lst=None, + displacement=0.01, + number_of_snapshots=None, + interaction_range=10, + factor=VaspToTHz, + ) + result_dict = evaluate_with_lammps( + task_dict={"calc_energy": structure_energy_dict, "calc_forces": structure_forces_dict}, + potential_dataframe=df_pot_selected, + ) + eng_internal_dict, phonopy_collect_dict = analyse_structures_helper( + phonopy_dict=phonopy_dict, + output_dict=result_dict, + dos_mesh=20, + number_of_snapshots=None, + ) + tp_collect_dict = get_thermal_properties( + eng_internal_dict=eng_internal_dict, + phonopy_dict=phonopy_dict, + structure_dict=structure_energy_dict, + repeat_vector=repeat_vector, + fit_type="polynomial", + fit_order=3, + t_min=1, + t_max=1500, + t_step=50, + temperatures=None, + ) + for key in ["temperatures", "free_energy", "volumes", "entropy", "heat_capacity"]: + self.assertTrue(len(tp_collect_dict[key]), 31) + self.assertEqual(tp_collect_dict["temperatures"][0], 1.0) + self.assertEqual(tp_collect_dict["temperatures"][-1], 1501.0) + self.assertTrue(tp_collect_dict["free_energy"][0] < 0.2) + self.assertTrue(tp_collect_dict["free_energy"][0] > 0.1) + self.assertTrue(tp_collect_dict["free_energy"][-1] < -2.6) + self.assertTrue(tp_collect_dict["free_energy"][-1] > -2.7) + self.assertTrue(tp_collect_dict["entropy"][0] < 0.1) + self.assertTrue(tp_collect_dict["entropy"][0] > 0.0) + self.assertTrue(tp_collect_dict["entropy"][-1] < 270) + self.assertTrue(tp_collect_dict["entropy"][-1] > 269) + self.assertTrue(tp_collect_dict["heat_capacity"][0] < 0.1) + self.assertTrue(tp_collect_dict["heat_capacity"][0] > 0.0) + self.assertTrue(tp_collect_dict["heat_capacity"][-1] < 100) + self.assertTrue(tp_collect_dict["heat_capacity"][-1] > 99) + self.assertTrue(tp_collect_dict["volumes"][-1] < 66.6) + self.assertTrue(tp_collect_dict["volumes"][-1] > 66.5) + self.assertTrue(tp_collect_dict["volumes"][0] < 66.5) + self.assertTrue(tp_collect_dict["volumes"][0] > 66.4) + thermal_properties_dict = get_thermal_properties( + eng_internal_dict=eng_internal_dict, + phonopy_dict=phonopy_dict, + structure_dict=structure_energy_dict, + repeat_vector=repeat_vector, + fit_type="polynomial", + fit_order=3, + temperatures=[100, 1000], + output_keys=["free_energy", "temperatures", "volumes"], + quantum_mechanical=True, + ) + temperatures_qh_qm, volumes_qh_qm = thermal_properties_dict["temperatures"], thermal_properties_dict["volumes"] + thermal_properties_dict = get_thermal_properties( + eng_internal_dict=eng_internal_dict, + phonopy_dict=phonopy_dict, + structure_dict=structure_energy_dict, + repeat_vector=repeat_vector, + fit_type="polynomial", + fit_order=3, + temperatures=[100, 1000], + output_keys=["temperatures", "volumes"], + quantum_mechanical=False, + ) + temperatures_qh_cl, volumes_qh_cl = thermal_properties_dict["temperatures"], thermal_properties_dict["volumes"] + self.assertEqual(len(eng_internal_dict.keys()), 11) + self.assertEqual(len(tp_collect_dict.keys()), 5) + self.assertEqual(len(temperatures_qh_qm), 2) + self.assertEqual(len(volumes_qh_qm), 2) + self.assertTrue(volumes_qh_qm[0] < volumes_qh_qm[-1]) + self.assertEqual(len(temperatures_qh_cl), 2) + self.assertEqual(len(volumes_qh_cl), 2) + self.assertTrue(volumes_qh_cl[0] < volumes_qh_cl[-1])