diff --git a/src/aiida_vibroscopy/calculations/spectra_utils.py b/src/aiida_vibroscopy/calculations/spectra_utils.py index 0cef018..be11f00 100644 --- a/src/aiida_vibroscopy/calculations/spectra_utils.py +++ b/src/aiida_vibroscopy/calculations/spectra_utils.py @@ -10,7 +10,7 @@ from __future__ import annotations from copy import deepcopy -from typing import Union +from typing import Tuple, Union from aiida import orm from aiida.engine import calcfunction @@ -29,6 +29,7 @@ 'compute_raman_susceptibility_tensors', 'compute_polarization_vectors', 'compute_complex_dielectric', + 'compute_clamped_pockels_tensor', 'get_supercells_for_hubbard', 'elaborate_susceptibility_derivatives', 'generate_vibrational_data_from_forces', @@ -420,6 +421,82 @@ def compute_complex_dielectric( return diel_infinity + prefactor * (complex_diel) / phonopy_instance.unitcell.volume +def compute_clamped_pockels_tensor( + phonopy_instance: Phonopy, + raman_tensors: np.ndarray, + nlo_susceptibility: np.ndarray, + nac_direction: None | list[float, float, float] = None, + imaginary_thr: float = -5.0 / UNITS.thz_to_cm, + skip_frequencies: int = 3, +) -> Tuple[np.ndarray, np.ndarray, np.ndarray]: + """Compute the clamped Pockels tensor in Cartesian coordinates. + + .. note:: Units are in pm/V + + :param nac_direction: non-analytical direction in Cartesian coordinates; + (3,) shape :class:`list` or :class:`numpy.ndarray` + :param degeneracy_tolerance: degeneracy tolerance for irreducible representation + :param imaginary_thr: threshold for activating warnings on negative frequencies (in Hz) + :param skip_frequencies: number of frequencies to not include (i.e. the acoustic modes) + + :return: tuple of (r_ion + r_el, r_el, r_ion), each having (3, 3, 3) shape array + """ + borns = phonopy_instance.nac_params['born'] * UNITS.elementary_charge_si # convert to Coulomb + dielectric = phonopy_instance.nac_params['dielectric'] + + dielectric_inv = np.linalg.inv(dielectric) + raman = raman_tensors * 1.0e10 # convert to 1/m + + q_reduced = None + # Have a look in compute_active_modes for the notation + if nac_direction is not None: + q_reduced = np.dot(phonopy_instance.primitive.cell, + nac_direction) / (2. * np.pi) # in reduced/crystal (PRIMITIVE) coordinates + + phonopy_instance.run_qpoints(q_points=[0, 0, 0], nac_q_direction=q_reduced, with_eigenvectors=True) + frequencies = phonopy_instance.qpoints.frequencies[0] * 1.0e+12 # THz -> Hz + eigvectors = phonopy_instance.qpoints.eigenvectors.T.real + + if frequencies.min() < imaginary_thr: + raise ValueError('Negative frequencies detected.') + + masses = phonopy_instance.masses * UNITS.atomic_mass_si + sqrt_masses = np.array([[np.sqrt(mass)] for mass in masses]) + + shape = (len(frequencies), len(masses), 3) # (modes, atoms, 3) + eigvectors = eigvectors.reshape(shape) + norm_eigvectors = np.array([eigv / sqrt_masses for eigv in eigvectors]) + + # norm_eigvectors shape|indices = (modes, atoms, 3) | (m, a, p) + # raman shape|indices = (atoms, 3, 3, 3) | (a, p, i, j) + # The contraction is performed over a and k, resulting in (m, i, j) raman susceptibility tensors. + alpha = np.tensordot(norm_eigvectors, raman, axes=([1, 2], [0, 1])) + + # borns charges shape|indices = (atoms, 3, 3) | (a, k, p) + # norm_eigvectors shape|indices = (modes, atoms, 3) | (m, a, p) + # polarization vectors shape | indices = (3, modes) | (k, m) + polvec = np.tensordot(borns, norm_eigvectors, axes=([0, 2], [1, 2])) # (k, m) + + ir_contribution = polvec / frequencies**2 # (k, m) + + # sets the first `skip_frequencies` IR contributions to zero (i.e. the acoustic modes) + for i in range(skip_frequencies): + ir_contribution[i] = 0 + r_ion_inner = np.tensordot(alpha, ir_contribution, axes=([0], [1])) # (i, j, k) + + r_ion_left = np.dot(dielectric_inv, np.transpose(r_ion_inner, axes=[1, 0, 2])) + r_ion_transposed = np.dot(np.transpose(r_ion_left, axes=[0, 2, 1]), dielectric_inv) + r_ion = -np.transpose(r_ion_transposed, axes=[0, 2, 1]) * 1.0e+12 # pm/V + + r_el_left = np.dot(dielectric_inv, np.transpose(nlo_susceptibility, axes=[1, 0, 2])) + r_el_transposed = -2 * np.dot(np.transpose(r_el_left, axes=[0, 2, 1]), dielectric_inv) + r_el = np.transpose(r_el_transposed, axes=[0, 2, 1]) + + r = r_el + r_ion + + return r, r_el, r_ion + + @calcfunction def get_supercells_for_hubbard( preprocess_data: PreProcessData, diff --git a/src/aiida_vibroscopy/common/constants.py b/src/aiida_vibroscopy/common/constants.py index ce1058a..6eab6c2 100644 --- a/src/aiida_vibroscopy/common/constants.py +++ b/src/aiida_vibroscopy/common/constants.py @@ -43,8 +43,12 @@ # The intensity must be computed using frequencies in cm^-1 and normalized eigenvectors # by atomic masses expressed in atomic mass unit (Dalton). # IMPORTANT: still misss the units from the Dirac delta - raman_xsection=1.0e24 * 1.054571817e-34 / (2.0 * units.SpeedOfLight**4 * units.AMU * units.THzToCm**3 - ) # removed 1/4pi due to convention on Chi2 for the correction + raman_xsection=1.0e24 * 1.054571817e-34 / + (2.0 * units.SpeedOfLight**4 * units.AMU * + units.THzToCm**3), # removed 1/4pi due to convention on Chi2 for the correction + elementary_charge_si=1.602176634e-19, # elementary charge in Coulomb + electron_mass_si=units.Me, # electron mass in kg + atomic_mass_si=units.AMU, # atomic mass unit in kg # to be defined: # * kelvin to eV # * nm to eV diff --git a/src/aiida_vibroscopy/data/vibro_mixin.py b/src/aiida_vibroscopy/data/vibro_mixin.py index 544412d..af7df7d 100644 --- a/src/aiida_vibroscopy/data/vibro_mixin.py +++ b/src/aiida_vibroscopy/data/vibro_mixin.py @@ -19,6 +19,7 @@ compute_raman_space_average, compute_raman_susceptibility_tensors, ) +from aiida_vibroscopy.common import UNITS from aiida_vibroscopy.utils.integration.lebedev import LebedevScheme from aiida_vibroscopy.utils.spectra import raman_prefactor @@ -662,3 +663,60 @@ def get_available_quadrature_order_schemes(): """Return the available orders for quadrature integration on the nac direction unitary sphere.""" from aiida_vibroscopy.utils.integration.lebedev import get_available_quadrature_order_schemes get_available_quadrature_order_schemes() + + def run_clamped_pockels_tensor( + self, + nac_direction: tuple[float, float, float] = lambda: [0, 0, 0], + imaginary_thr: float = -5.0 / UNITS.thz_to_cm, + skip_frequencies: int = 3, + asr_sum_rules: bool = False, + symmetrize_fc: bool = False, + **kwargs, + ) -> np.ndarray: + """Compute the clamped Pockels tensor in Cartesian coordinates. + + .. note:: Units are in pm/V + + :param nac_direction: non-analytical direction in Cartesian coordinates; + (3,) shape :class:`list` or :class:`numpy.ndarray` + :param degeneracy_tolerance: degeneracy tolerance for irreducible representation + :param imaginary_thr: threshold for activating warnings on negative frequencies (in Hz) + :param skip_frequencies: number of frequencies to not include (i.e. the acoustic modes) + :param asr_sum_rules: whether to apply acoustic sum rules to the force constants + :param symmetrize_fc: whether to symmetrize the force constants using space group + :param kwargs: see also the :func:`~aiida_phonopy.data.phonopy.get_phonopy_instance` method + + * subtract_residual_forces: + whether or not subract residual forces (if set); + bool, defaults to False + * symmetrize_nac: + whether or not to symmetrize the nac parameters + using point group symmetry; bool, defaults to self.is_symmetry + + :return: tuple of (r_ion + r_el, r_el, r_ion), each having (3, 3, 3) shape array + """ + from aiida_vibroscopy.calculations.spectra_utils import compute_clamped_pockels_tensor + + if not isinstance(symmetrize_fc, bool) or not isinstance(asr_sum_rules, bool): + raise TypeError('the input is not of the correct type') + + phonopy_instance = self.get_phonopy_instance(**kwargs) + + if phonopy_instance.force_constants is None: + phonopy_instance.produce_force_constants() + + if asr_sum_rules: + phonopy_instance.symmetrize_force_constants() + if symmetrize_fc: + phonopy_instance.symmetrize_force_constants_by_space_group() + + results = compute_clamped_pockels_tensor( + phonopy_instance=phonopy_instance, + raman_tensors=self.raman_tensors, + nlo_susceptibility=self.nlo_susceptibility, + nac_direction=nac_direction, + imaginary_thr=imaginary_thr, + skip_frequencies=skip_frequencies, + ) + + return results diff --git a/tests/calculations/test_spectra.py b/tests/calculations/test_spectra.py index 3868aba..7c3e05e 100644 --- a/tests/calculations/test_spectra.py +++ b/tests/calculations/test_spectra.py @@ -24,13 +24,13 @@ def generate_phonopy_instance(): - symmetry info """ - def _generate_phonopy_instance(): + def _generate_phonopy_instance(name='AlAs'): """Return AlAs Phonopy instance.""" import os import phonopy - filename = 'phonopy_AlAs.yaml' + filename = f'phonopy_{name}.yaml' basepath = os.path.dirname(os.path.abspath(__file__)) phyaml = os.path.join(basepath, filename) @@ -151,6 +151,30 @@ def test_compute_raman_susceptibility_tensors(generate_phonopy_instance, generat assert np.abs(abs(alpha_comp) - abs(alpha_theo)) < 1e-3 +def test_compute_clamped_pockels_tensor(generate_phonopy_instance, generate_third_rank_tensors): + """Test the `compute_clamped_pockels_tensor` function.""" + from aiida_vibroscopy.calculations.spectra_utils import compute_clamped_pockels_tensor + + ph = generate_phonopy_instance() + raman, chi2 = generate_third_rank_tensors() + + pockels, pockels_el, pockels_ion = compute_clamped_pockels_tensor( + phonopy_instance=ph, + raman_tensors=raman, + nlo_susceptibility=chi2, + ) + + if DEBUG: + print('\n', '================================', '\n') + print('\t', 'DEBUG') + print(pockels) + print('\t', 'DEBUG EL') + print(pockels_el) + print('\t', 'DEBUG ION') + print(pockels_ion) + print('\n', '================================', '\n') + + def test_compute_methods(generate_phonopy_instance, generate_third_rank_tensors, ndarrays_regression): """Test the post-processing methods with data regression techniques.""" from aiida_vibroscopy.calculations.spectra_utils import (