Feat: add clamped Pockels calculation capability

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,

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

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

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 (