diff --git a/Modules/DynamicalLanczos.py b/Modules/DynamicalLanczos.py index 51c8a322..d079e57b 100644 --- a/Modules/DynamicalLanczos.py +++ b/Modules/DynamicalLanczos.py @@ -986,6 +986,68 @@ def load_from_input_files(self, root_name = "tdscha", directory="."): + def _build_raman_vector(self, pol_vec_in = np.array([1,0,0]), + pol_vec_out = np.array([1,0,0]), mixed = False, + pol_in_2 = None, pol_out_2 = None, + unpolarized = None, normalized = True): + """Build a unit-cell Raman perturbation vector. + + ``unpolarized=None`` returns the contraction with the supplied light + polarizations. When ``mixed`` is true, the second contraction is + added coherently. + + Unpolarized indices 0--6 return, in order, the combinations + ``trace``, ``xx-yy``, ``xx-zz``, ``yy-zz``, ``xy``, ``xz``, and + ``yz``. With ``normalized=True`` these are multiplied by ``1/3``, + ``1/sqrt(2)`` (indices 1--3), and ``sqrt(3)`` (indices 4--6), as used + by :meth:`prepare_raman`. ``normalized=False`` retains the legacy + raw-component convention of :meth:`prepare_unpolarized_raman`. + """ + assert self.dyn.raman_tensor is not None, ( + "Error, no Raman tensor found. Cannot initialize the Raman response") + + if unpolarized is None: + raman_v = np.array( + self.dyn.GetRamanVector(pol_vec_in, pol_vec_out), copy=True) + if mixed: + raman_v += self.dyn.GetRamanVector(pol_in_2, pol_out_2) + return raman_v + + if unpolarized not in range(7): + raise ValueError( + "Error, unpolarized must be between [0, ..., 6] got invalid {}." + .format(unpolarized)) + + px = np.array([1, 0, 0]) + py = np.array([0, 1, 0]) + pz = np.array([0, 0, 1]) + + if unpolarized == 0: + raman_v = (self.dyn.GetRamanVector(px, px) + + self.dyn.GetRamanVector(py, py) + + self.dyn.GetRamanVector(pz, pz)) + elif unpolarized == 1: + raman_v = (self.dyn.GetRamanVector(px, px) + - self.dyn.GetRamanVector(py, py)) + elif unpolarized == 2: + raman_v = (self.dyn.GetRamanVector(px, px) + - self.dyn.GetRamanVector(pz, pz)) + elif unpolarized == 3: + raman_v = (self.dyn.GetRamanVector(py, py) + - self.dyn.GetRamanVector(pz, pz)) + elif unpolarized == 4: + raman_v = self.dyn.GetRamanVector(px, py) + elif unpolarized == 5: + raman_v = self.dyn.GetRamanVector(px, pz) + else: + raman_v = self.dyn.GetRamanVector(py, pz) + + if normalized: + scales = [1 / 3] + [1 / np.sqrt(2)] * 3 + [np.sqrt(3)] * 3 + raman_v = raman_v * scales[unpolarized] + + return np.array(raman_v, copy=True) + def prepare_raman(self, pol_vec_in = np.array([1,0,0]), pol_vec_out = np.array([1,0,0]), mixed = False, pol_in_2 = None, pol_out_2 = None, unpolarized: int = None): """ PREPARE LANCZOS FOR RAMAN SPECTRUM @@ -1001,6 +1063,11 @@ def prepare_raman(self, pol_vec_in = np.array([1,0,0]), pol_vec_out = np.array([ The polarization vector of the incoming light pol_vec_out : ndarray (size = 3) The polarization vector for the outcoming light + mixed : bool + If True, coherently add the contraction specified by + pol_in_2 and pol_out_2. + pol_in_2, pol_out_2 : ndarray (size = 3) or None + The second pair of light polarizations when mixed=True. unpolarized : int or None The perturbation for unpolarized raman (if different from None, overrides the behaviour of pol_vec_in and pol_vec_out). Indices goes from 0 to 6 (included). @@ -1016,97 +1083,18 @@ def prepare_raman(self, pol_vec_in = np.array([1,0,0]), pol_vec_out = np.array([ The total unpolarized raman intensity is 45 alpha^2 + 7 beta^2 """ - - # Check if the raman tensor is present - assert not self.dyn.raman_tensor is None, "Error, no Raman tensor found. Cannot initialize the Raman responce" - - # Get the raman vector (apply the ASR and contract the raman tensor with the polarization vectors) - raman_v = self.dyn.GetRamanVector(pol_vec_in, pol_vec_out) - if mixed: print('Prepare Raman') print('Adding other component of the Raman tensor') - raman_v += self.dyn.GetRamanVector(pol_in_2, pol_out_2) - # Get the raman vector in the supercelld - n_supercell = np.prod(self.dyn.GetSupercell()) + raman_v = self._build_raman_vector( + pol_vec_in=pol_vec_in, pol_vec_out=pol_vec_out, mixed=mixed, + pol_in_2=pol_in_2, pol_out_2=pol_out_2, + unpolarized=unpolarized, normalized=True) - if unpolarized is None: - # Get the raman vector - raman_v = self.dyn.GetRamanVector(pol_vec_in, pol_vec_out) - - # Get the raman vector in the supercelld - new_raman_v = np.tile(raman_v.ravel(), n_supercell) - - # Convert in the polarization basis and store the intensity - self.prepare_perturbation(new_raman_v, masses_exp=-1) - else: - px = np.array([1,0,0]) - py = np.array([0,1,0]) - pz = np.array([0,0,1]) - - if unpolarized == 0: - # Alpha - raman_v = self.dyn.GetRamanVector(px, px) - new_raman_v = np.tile(raman_v.ravel(), n_supercell) / 3 - self.prepare_perturbation(new_raman_v, masses_exp=-1) - - raman_v = self.dyn.GetRamanVector(py, py) - new_raman_v = np.tile(raman_v.ravel(), n_supercell) / 3 - self.prepare_perturbation(new_raman_v, masses_exp=-1, add = True) - - raman_v = self.dyn.GetRamanVector(pz, pz) - new_raman_v = np.tile(raman_v.ravel(), n_supercell) / 3 - self.prepare_perturbation(new_raman_v, masses_exp=-1, add = True) - elif unpolarized == 1: - # (xx -yy)^2 / 2 - raman_v = self.dyn.GetRamanVector(px, px) - new_raman_v = np.tile(raman_v.ravel(), n_supercell) / np.sqrt(2) - self.prepare_perturbation(new_raman_v, masses_exp=-1) - - raman_v = self.dyn.GetRamanVector(py, py) - new_raman_v = - np.tile(raman_v.ravel(), n_supercell) / np.sqrt(2) - self.prepare_perturbation(new_raman_v, masses_exp=-1, add = True) - elif unpolarized == 2: - # beta_2 = (xx -zz)^2 / 2 - raman_v = self.dyn.GetRamanVector(px, px) - new_raman_v = np.tile(raman_v.ravel(), n_supercell) / np.sqrt(2) - self.prepare_perturbation(new_raman_v, masses_exp=-1) - - raman_v = self.dyn.GetRamanVector(pz, pz) - new_raman_v = - np.tile(raman_v.ravel(), n_supercell) / np.sqrt(2) - self.prepare_perturbation(new_raman_v, masses_exp=-1, add = True) - elif unpolarized == 3: - # beta_2 = (yy -zz)^2 / 2 - raman_v = self.dyn.GetRamanVector(py, py) - new_raman_v = np.tile(raman_v.ravel(), n_supercell) / np.sqrt(2) - self.prepare_perturbation(new_raman_v, masses_exp=-1) - - raman_v = self.dyn.GetRamanVector(pz, pz) - new_raman_v = - np.tile(raman_v.ravel(), n_supercell) / np.sqrt(2) - self.prepare_perturbation(new_raman_v, masses_exp=-1, add = True) - elif unpolarized == 4: - # beta_2 = 3 xy^2 - raman_v = self.dyn.GetRamanVector(px, py) - new_raman_v = np.tile(raman_v.ravel(), n_supercell) * np.sqrt(3) - self.prepare_perturbation(new_raman_v, masses_exp=-1) - elif unpolarized == 5: - # beta_2 = 3 yz^2 - raman_v = self.dyn.GetRamanVector(py, pz) - new_raman_v = np.tile(raman_v.ravel(), n_supercell) * np.sqrt(3) - self.prepare_perturbation(new_raman_v, masses_exp=-1) - elif unpolarized == 6: - # beta_2 = 3 xz^2 - raman_v = self.dyn.GetRamanVector(px, pz) - new_raman_v = np.tile(raman_v.ravel(), n_supercell) * np.sqrt(3) - self.prepare_perturbation(new_raman_v, masses_exp=-1) - else: - raise ValueError("Error, unpolarized must be between [0, ... ,6] got invalid {}.".format(unpolarized)) - - - - - # Convert in the polarization basis and store the intensity + # Tile the unit-cell perturbation over the supercell. + n_supercell = np.prod(self.dyn.GetSupercell()) + new_raman_v = np.tile(raman_v.ravel(), n_supercell) self.prepare_perturbation(new_raman_v, masses_exp=-1) def get_prefactors_unpolarized_raman(self, index): @@ -1150,43 +1138,8 @@ def prepare_unpolarized_raman(self, index = 0, debug = False): + 7/2 [(xx-yy)^2 + (xx-zz)^2 + (yy-zz)^2] + 7 * 3 [(xy)^2 + (yz)^2 + (xz)^2] """ - # Check if the raman tensor is present - assert not self.dyn.raman_tensor is None, "Error, no Raman tensor found. Cannot initialize the Raman responce" - - labels = [i for i in range(7)] - if not(index in labels): - raise ValueError('{} should be in {}'.format(index, labels)) - - epols = {'x' : np.array([1,0,0]),\ - 'y' : np.array([0,1,0]),\ - 'z' : np.array([0,0,1])} - - # (xx + yy + zz)^2 - if index == 0: - raman_v = self.dyn.GetRamanVector(epols['x'], epols['x']) - raman_v += self.dyn.GetRamanVector(epols['y'], epols['y']) - raman_v += self.dyn.GetRamanVector(epols['z'], epols['z']) - # (xx - yy)^2 - elif index == 1: - raman_v = self.dyn.GetRamanVector(epols['x'], epols['x']) - raman_v -= self.dyn.GetRamanVector(epols['y'], epols['y']) - # (xx - zz)^2 - elif index == 2: - raman_v = self.dyn.GetRamanVector(epols['x'], epols['x']) - raman_v -= self.dyn.GetRamanVector(epols['z'], epols['z']) - # (yy - zz)^2 - elif index == 3: - raman_v = self.dyn.GetRamanVector(epols['y'], epols['y']) - raman_v -= self.dyn.GetRamanVector(epols['z'], epols['z']) - # (xy)^2 - elif index == 4: - raman_v = self.dyn.GetRamanVector(epols['x'], epols['y']) - # (xz)^2 - elif index == 5: - raman_v = self.dyn.GetRamanVector(epols['x'], epols['z']) - # (yz)^2 - elif index == 6: - raman_v = self.dyn.GetRamanVector(epols['y'], epols['z']) + raman_v = self._build_raman_vector( + unpolarized=index, normalized=False) if debug: np.save('raman_v_{}'.format(index), raman_v) @@ -1995,11 +1948,14 @@ def prepare_perturbation(self, vector, masses_exp = 1, add = False): # THIS IS OK IN THE WIGNER REPRESENTATION BECAUSE # THE PERTUBATION ENTERS ONLY IN THE R SECTOR - self.perturbation_modulus = new_v.dot(new_v) - if self.symmetrize: self.symmetrize_psi() + # Account for the full accumulated one-phonon perturbation. In + # particular, add=True must include cross terms with earlier vectors. + perturbation = self.psi[:self.n_modes] + self.perturbation_modulus = perturbation.dot(perturbation) + def prepare_mode(self, index): @@ -7360,6 +7316,3 @@ def min_stdes(func, args, x0, step = 1e-2, n_iters = 100): print("F: {} | G: {}".format( f, np.sqrt(np.sum(grad**2)))) return x - - - diff --git a/Modules/QSpaceLanczos.py b/Modules/QSpaceLanczos.py index 1584b057..eb6612d3 100644 --- a/Modules/QSpaceLanczos.py +++ b/Modules/QSpaceLanczos.py @@ -1247,87 +1247,20 @@ def prepare_raman(self, pol_vec_in=np.array([1.0, 0.0, 0.0]), pol_vec_out=np.arr The total unpolarized raman intensity is 45 alpha^2 + 7 beta^2 """ - # Check if the raman tensor is present - assert not self.dyn.raman_tensor is None, "Error, no Raman tensor found. Cannot initialize Raman response" + if mixed: + print('Prepare Raman') + print('Adding other component of the Raman tensor') + raman_v = self._build_raman_vector( + pol_vec_in=pol_vec_in, pol_vec_out=pol_vec_out, mixed=mixed, + pol_in_2=pol_in_2, pol_out_2=pol_out_2, + unpolarized=unpolarized, normalized=True) + + # A constant real-space perturbation is a Gamma vector whose + # unit-cell amplitude is multiplied by sqrt(number of cells). n_cell = np.prod(self.dyn.GetSupercell()) - - if unpolarized is None: - # Get the raman vector (apply the ASR and contract the raman tensor with the polarization vectors) - raman_v = self.dyn.GetRamanVector(pol_vec_in, pol_vec_out) - - if mixed: - print('Prepare Raman') - print('Adding other component of the Raman tensor') - raman_v += self.dyn.GetRamanVector(pol_in_2, pol_out_2) - - # Scale for Γ-point constant perturbation - new_raman_v = raman_v.ravel() * np.sqrt(n_cell) - - # Convert in the polarization basis - self.prepare_perturbation_q(0, new_raman_v) - else: - px = np.array([1, 0, 0]) - py = np.array([0, 1, 0]) - pz = np.array([0, 0, 1]) - - if unpolarized == 0: - # Alpha = (xx + yy + zz)^2/9 - raman_v = self.dyn.GetRamanVector(px, px) - new_raman_v = raman_v.ravel() * np.sqrt(n_cell) / 3 - self.prepare_perturbation_q(0, new_raman_v) - - raman_v = self.dyn.GetRamanVector(py, py) - new_raman_v = raman_v.ravel() * np.sqrt(n_cell) / 3 - self.prepare_perturbation_q(0, new_raman_v, add=True) - - raman_v = self.dyn.GetRamanVector(pz, pz) - new_raman_v = raman_v.ravel() * np.sqrt(n_cell) / 3 - self.prepare_perturbation_q(0, new_raman_v, add=True) - elif unpolarized == 1: - # (xx - yy)^2 / 2 - raman_v = self.dyn.GetRamanVector(px, px) - new_raman_v = raman_v.ravel() * np.sqrt(n_cell) / np.sqrt(2) - self.prepare_perturbation_q(0, new_raman_v) - - raman_v = self.dyn.GetRamanVector(py, py) - new_raman_v = -raman_v.ravel() * np.sqrt(n_cell) / np.sqrt(2) - self.prepare_perturbation_q(0, new_raman_v, add=True) - elif unpolarized == 2: - # (xx - zz)^2 / 2 - raman_v = self.dyn.GetRamanVector(px, px) - new_raman_v = raman_v.ravel() * np.sqrt(n_cell) / np.sqrt(2) - self.prepare_perturbation_q(0, new_raman_v) - - raman_v = self.dyn.GetRamanVector(pz, pz) - new_raman_v = -raman_v.ravel() * np.sqrt(n_cell) / np.sqrt(2) - self.prepare_perturbation_q(0, new_raman_v, add=True) - elif unpolarized == 3: - # (yy - zz)^2 / 2 - raman_v = self.dyn.GetRamanVector(py, py) - new_raman_v = raman_v.ravel() * np.sqrt(n_cell) / np.sqrt(2) - self.prepare_perturbation_q(0, new_raman_v) - - raman_v = self.dyn.GetRamanVector(pz, pz) - new_raman_v = -raman_v.ravel() * np.sqrt(n_cell) / np.sqrt(2) - self.prepare_perturbation_q(0, new_raman_v, add=True) - elif unpolarized == 4: - # 3 xy^2 - raman_v = self.dyn.GetRamanVector(px, py) - new_raman_v = raman_v.ravel() * np.sqrt(n_cell) * np.sqrt(3) - self.prepare_perturbation_q(0, new_raman_v) - elif unpolarized == 5: - # 3 yz^2 - raman_v = self.dyn.GetRamanVector(py, pz) - new_raman_v = raman_v.ravel() * np.sqrt(n_cell) * np.sqrt(3) - self.prepare_perturbation_q(0, new_raman_v) - elif unpolarized == 6: - # 3 xz^2 - raman_v = self.dyn.GetRamanVector(px, pz) - new_raman_v = raman_v.ravel() * np.sqrt(n_cell) * np.sqrt(3) - self.prepare_perturbation_q(0, new_raman_v) - else: - raise ValueError(f"Error, unpolarized must be between [0, ..., 6] got invalid {unpolarized}.") + new_raman_v = raman_v.ravel() * np.sqrt(n_cell) + self.prepare_perturbation_q(0, new_raman_v) def prepare_unpolarized_raman(self, index=0, debug=False): """ @@ -1347,45 +1280,9 @@ def prepare_unpolarized_raman(self, index=0, debug=False): Note: This method prepares the raw components WITHOUT prefactors. Use get_prefactors_unpolarized_raman() to get the correct prefactors. """ - # Check if the raman tensor is present - assert not self.dyn.raman_tensor is None, "Error, no Raman tensor found. Cannot initialize the Raman response" - - labels = [i for i in range(7)] - if index not in labels: - raise ValueError(f'{index} should be in {labels}') - - epols = {'x': np.array([1, 0, 0]), - 'y': np.array([0, 1, 0]), - 'z': np.array([0, 0, 1])} - + raman_v = self._build_raman_vector( + unpolarized=index, normalized=False) n_cell = np.prod(self.dyn.GetSupercell()) - - # (xx + yy + zz)^2 - if index == 0: - raman_v = self.dyn.GetRamanVector(epols['x'], epols['x']) - raman_v += self.dyn.GetRamanVector(epols['y'], epols['y']) - raman_v += self.dyn.GetRamanVector(epols['z'], epols['z']) - # (xx - yy)^2 - elif index == 1: - raman_v = self.dyn.GetRamanVector(epols['x'], epols['x']) - raman_v -= self.dyn.GetRamanVector(epols['y'], epols['y']) - # (xx - zz)^2 - elif index == 2: - raman_v = self.dyn.GetRamanVector(epols['x'], epols['x']) - raman_v -= self.dyn.GetRamanVector(epols['z'], epols['z']) - # (yy - zz)^2 - elif index == 3: - raman_v = self.dyn.GetRamanVector(epols['y'], epols['y']) - raman_v -= self.dyn.GetRamanVector(epols['z'], epols['z']) - # (xy)^2 - elif index == 4: - raman_v = self.dyn.GetRamanVector(epols['x'], epols['y']) - # (xz)^2 - elif index == 5: - raman_v = self.dyn.GetRamanVector(epols['x'], epols['z']) - # (yz)^2 - elif index == 6: - raman_v = self.dyn.GetRamanVector(epols['y'], epols['z']) if debug: np.save(f'raman_v_{index}', raman_v) @@ -1423,7 +1320,9 @@ def prepare_perturbation_q(self, iq, vector, add=False): v_scaled = vector / np.sqrt(m) R1 = np.conj(self.pols_q[:, :, iq]).T @ v_scaled # (n_bands,) complex self.psi[:self.n_bands] += R1 - self.perturbation_modulus = np.real(np.conj(R1) @ R1) + perturbation = self.psi[:self.n_bands] + self.perturbation_modulus = np.real( + np.conj(perturbation) @ perturbation) def reset_q(self): """Reset the Lanczos state for q-space.""" diff --git a/docs/examples.md b/docs/examples.md index 3095a14c..9b480f55 100644 --- a/docs/examples.md +++ b/docs/examples.md @@ -159,6 +159,12 @@ np.savetxt("raman_unpolarized_total.dat", print("Total unpolarized Raman intensity computed") ``` +This example uses the raw-component +`prepare_unpolarized_raman(index=i)` convention. Equivalently, call +`prepare_raman(unpolarized=i)` and replace `prefactors` with +`[45, 7, 7, 7, 7, 7, 7]`; the resulting weighted channel intensities are +identical. + ## Example 4: StaticHessian Calculation Compute free energy Hessian for stability analysis: @@ -650,4 +656,4 @@ See the `Examples/Comparison/` directory for: And `Examples/example_IR_Raman_2p/` for: - `IR_UNPOL/` - Unpolarized IR with 1ph/2ph processes - `RAMAN_UNPOL/` - Unpolarized Raman calculations -- Complete README files with instructions \ No newline at end of file +- Complete README files with instructions diff --git a/docs/usage.md b/docs/usage.md index 0315fac9..5eaad395 100644 --- a/docs/usage.md +++ b/docs/usage.md @@ -71,6 +71,11 @@ for i in range(7): lanczos.save_status(f"raman_unpolarized_{i}.npz") ``` +!!! warning "Raman data produced before the 1.7 hotfix" + + Saved `prepare_raman(unpolarized=i)` calculations for channels 0–3 used + incomplete diagonal combinations and must be recomputed. Channels 4–6 + were unaffected and can be reused. Then you can plot the unpolarized Raman spectrum by summing the contributions of the 7 components. This is done in the following way: @@ -86,8 +91,9 @@ w_ry = w/CC.Units.RY_TO_CM # Convert in Ry (the internal unit of tdscha) smearing = 2/CC.Units.RY_TO_CM # Smearing in cm⁻¹ raman_signal = np.zeros_like(w) +weights = [45, 7, 7, 7, 7, 7, 7] -# Load the 7 unpolarized Raman components and sum them. +# Load and combine the seven normalized invariant components. for i in range(7): lanczos = DL.Lanczos() lanczos.load_status(f"raman_unpolarized_{i}.npz") @@ -97,7 +103,7 @@ for i in range(7): # The response is proportional to the imaginary part of the Green's function. # The '-' sign selects the retarded response, which is the one relevant for Raman scattering. - raman_signal += -np.imag(gf) + raman_signal += weights[i] * -np.imag(gf) # Then, we can just plot the data @@ -107,6 +113,14 @@ plt.ylabel("Unpolarized Raman Intensity (arb. units)") plt.show() ``` +`prepare_raman(unpolarized=i)` prepares normalized invariants and therefore +uses weights `[45, 7, 7, 7, 7, 7, 7]`. The legacy +`prepare_unpolarized_raman(index=i)` method prepares the corresponding raw +Cartesian combinations. When using that API, multiply each spectrum by +`lanczos.get_prefactors_unpolarized_raman(i)`, currently +`[5, 7/2, 7/2, 7/2, 21, 21, 21]`. The two conventions give the same total +unpolarized intensity. + ## Parallel Execution Modes TD-SCHA supports four computation modes: diff --git a/tests/test_qspace/test_raman_ir_invariants.py b/tests/test_qspace/test_raman_ir_invariants.py new file mode 100644 index 00000000..f54c8d99 --- /dev/null +++ b/tests/test_qspace/test_raman_ir_invariants.py @@ -0,0 +1,257 @@ +"""Fast regression tests for Raman invariants and IR backend scaling.""" + +import importlib + +import numpy as np +import pytest + +import tdscha.DynamicalLanczos as DL +import tdscha.QSpaceLanczos as QL + + +N_CELL = 4 +N_ATOMS = 2 +N_CART = 3 * N_ATOMS +MASSES = np.array([1.0, 4.0]) + + +class _FakeDyn: + def __init__(self): + self.raman_tensor = np.zeros((3, 3, N_CART)) + components = { + (0, 0): [1.0, -2.0, 3.0, 4.0, -5.0, 6.0], + (1, 1): [-3.0, 5.0, 2.0, -1.0, 7.0, 4.0], + (2, 2): [8.0, 1.0, -4.0, 3.0, 2.0, -6.0], + (0, 1): [2.0, 3.0, -1.0, 5.0, -2.0, 7.0], + (0, 2): [-5.0, 4.0, 6.0, 2.0, 1.0, -3.0], + (1, 2): [7.0, -6.0, 5.0, -4.0, 3.0, -2.0], + } + for (i, j), values in components.items(): + self.raman_tensor[i, j] = values + self.raman_tensor[j, i] = values + + self.effective_charges = np.array([ + [[1.0, 2.0, -1.0], [3.0, -2.0, 4.0], [5.0, 1.0, 2.0]], + [[-2.0, 1.0, 3.0], [4.0, 2.0, -3.0], [1.0, -5.0, 6.0]], + ]) + + def GetRamanVector(self, pol_in, pol_out): + return np.einsum( + "i,j,ijk->k", pol_in, pol_out, self.raman_tensor) + + @staticmethod + def GetSupercell(): + return np.array([2, 1, 2]) + + +class _FakeStructure: + @staticmethod + def get_masses_array(): + return MASSES.copy() + + +class _RealHarness(DL.Lanczos): + def __init__(self, dyn): + self.dyn = dyn + masses_uc = np.repeat(MASSES, 3) + self.m = np.tile(masses_uc, N_CELL) + self.n_modes = N_CART + gamma = np.ones(N_CELL) / np.sqrt(N_CELL) + self.pols = np.kron(gamma[:, None], np.eye(N_CART)) + self.psi = np.zeros(self.n_modes) + self.symmetrize = False + self.ignore_small_w = True + + def reset(self): + pass + + +class _QSpaceHarness(QL.QSpaceLanczos): + def __init__(self, dyn): + self.dyn = dyn + self.uci_structure = _FakeStructure() + self.n_bands = N_CART + self.pols_q = np.eye(N_CART, dtype=np.complex128)[:, :, None] + self.psi = np.zeros(self.n_bands, dtype=np.complex128) + + def build_q_pair_map(self, iq): + assert iq == 0 + + def reset_q(self): + self.psi = np.zeros(self.n_bands, dtype=np.complex128) + + +@pytest.fixture +def dyn(): + return _FakeDyn() + + +def _raw_invariants(tensor): + return [ + tensor[0, 0] + tensor[1, 1] + tensor[2, 2], + tensor[0, 0] - tensor[1, 1], + tensor[0, 0] - tensor[2, 2], + tensor[1, 1] - tensor[2, 2], + tensor[0, 1], + tensor[0, 2], + tensor[1, 2], + ] + + +def _projected_unit_cell(vector): + masses = np.repeat(MASSES, 3) + return np.sqrt(N_CELL) * vector / np.sqrt(masses) + + +def test_shared_builder_returns_all_seven_invariants(dyn): + lanc = _RealHarness(dyn) + raw = _raw_invariants(dyn.raman_tensor) + scales = [1 / 3] + [1 / np.sqrt(2)] * 3 + [np.sqrt(3)] * 3 + + for index, expected_raw in enumerate(raw): + np.testing.assert_allclose( + lanc._build_raman_vector( + unpolarized=index, normalized=False), + expected_raw) + np.testing.assert_allclose( + lanc._build_raman_vector( + unpolarized=index, normalized=True), + expected_raw * scales[index]) + + +@pytest.mark.parametrize("backend", [_RealHarness, _QSpaceHarness]) +def test_polarized_and_coherently_mixed_raman(backend, dyn): + lanc = backend(dyn) + pol_in = np.array([0.5, -1.0, 2.0]) + pol_out = np.array([1.5, 0.25, -0.75]) + pol_in_2 = np.array([-0.5, 2.0, 1.0]) + pol_out_2 = np.array([1.0, -1.5, 0.5]) + + direct = dyn.GetRamanVector(pol_in, pol_out) + lanc.prepare_raman(pol_vec_in=pol_in, pol_vec_out=pol_out) + np.testing.assert_allclose(lanc.psi[:N_CART], _projected_unit_cell(direct)) + np.testing.assert_allclose( + lanc.perturbation_modulus, + np.vdot(_projected_unit_cell(direct), + _projected_unit_cell(direct)).real) + + direct_mixed = direct + dyn.GetRamanVector(pol_in_2, pol_out_2) + lanc.prepare_raman( + pol_vec_in=pol_in, pol_vec_out=pol_out, mixed=True, + pol_in_2=pol_in_2, pol_out_2=pol_out_2) + np.testing.assert_allclose( + lanc.psi[:N_CART], _projected_unit_cell(direct_mixed)) + np.testing.assert_allclose( + lanc.perturbation_modulus, + np.vdot(_projected_unit_cell(direct_mixed), + _projected_unit_cell(direct_mixed)).real) + + +@pytest.mark.parametrize("backend", [_RealHarness, _QSpaceHarness]) +def test_both_unpolarized_apis_have_equivalent_weighted_intensities( + backend, dyn): + normalized_weights = [45] + [7] * 6 + lanc = backend(dyn) + + for index, normalized_weight in enumerate(normalized_weights): + lanc.prepare_raman(unpolarized=index) + normalized_intensity = ( + normalized_weight * lanc.perturbation_modulus) + + lanc.prepare_unpolarized_raman(index=index) + raw_intensity = ( + lanc.get_prefactors_unpolarized_raman(index) + * lanc.perturbation_modulus) + + np.testing.assert_allclose(normalized_intensity, raw_intensity) + + +def test_real_and_qspace_raman_vectors_and_moduli_match(dyn): + real = _RealHarness(dyn) + qspace = _QSpaceHarness(dyn) + cases = [ + {"pol_vec_in": np.array([0.5, -1.0, 2.0]), + "pol_vec_out": np.array([1.5, 0.25, -0.75])}, + {"pol_vec_in": np.array([1.0, 0.0, 0.0]), + "pol_vec_out": np.array([0.0, 1.0, 0.0]), + "mixed": True, + "pol_in_2": np.array([0.0, 0.0, 1.0]), + "pol_out_2": np.array([0.0, 1.0, 0.0])}, + ] + cases.extend({"unpolarized": index} for index in range(7)) + + for kwargs in cases: + real.prepare_raman(**kwargs) + qspace.prepare_raman(**kwargs) + np.testing.assert_allclose( + qspace.psi[:N_CART], real.psi[:N_CART]) + np.testing.assert_allclose( + qspace.perturbation_modulus, real.perturbation_modulus) + + +def test_generic_add_uses_complete_accumulated_perturbation(dyn): + vector_1 = np.arange(1, N_CART + 1, dtype=float) + vector_2 = np.array([-2.0, 1.0, 4.0, -3.0, 5.0, 2.0]) + + real = _RealHarness(dyn) + vector_1_sc = np.tile(vector_1, N_CELL) + vector_2_sc = np.tile(vector_2, N_CELL) + real.prepare_perturbation(vector_1_sc, masses_exp=-1) + real.prepare_perturbation(vector_2_sc, masses_exp=-1, add=True) + expected = _projected_unit_cell(vector_1 + vector_2) + np.testing.assert_allclose(real.psi[:N_CART], expected) + np.testing.assert_allclose(real.perturbation_modulus, expected @ expected) + + qspace = _QSpaceHarness(dyn) + qspace.prepare_perturbation_q(0, vector_1 * np.sqrt(N_CELL)) + qspace.prepare_perturbation_q( + 0, vector_2 * np.sqrt(N_CELL), add=True) + np.testing.assert_allclose(qspace.psi[:N_CART], expected) + np.testing.assert_allclose( + qspace.perturbation_modulus, np.vdot(expected, expected).real) + + +def test_ir_real_qspace_parity_and_powder_average(dyn): + real = _RealHarness(dyn) + qspace = _QSpaceHarness(dyn) + moduli = [] + + for pol_vec in np.eye(3): + direct = np.einsum( + "abc,b->ac", dyn.effective_charges, pol_vec).ravel() + expected = _projected_unit_cell(direct) + + real.prepare_ir(pol_vec=pol_vec) + qspace.prepare_ir(pol_vec=pol_vec) + np.testing.assert_allclose(real.psi[:N_CART], expected) + np.testing.assert_allclose(qspace.psi[:N_CART], expected) + np.testing.assert_allclose( + qspace.perturbation_modulus, real.perturbation_modulus) + moduli.append(real.perturbation_modulus) + + powder_average = sum(moduli) / 3 + direct_powder_average = sum( + np.vdot( + _projected_unit_cell( + np.einsum( + "abc,b->ac", dyn.effective_charges, pol_vec).ravel()), + _projected_unit_cell( + np.einsum( + "abc,b->ac", dyn.effective_charges, pol_vec).ravel()) + ).real + for pol_vec in np.eye(3) + ) / 3 + np.testing.assert_allclose(powder_average, direct_powder_average) + + +def test_atom_fourier_lanczos_inherits_qspace_raman_implementation(): + try: + module = importlib.import_module("tdscha.QSpaceAtomFourier") + except ImportError: + pytest.skip("QSpaceAtomFourier is only present on the 1.8 branch") + + cls = module.QSpaceAtomFourierLanczos + assert issubclass(cls, QL.QSpaceLanczos) + assert "prepare_raman" not in cls.__dict__ + assert "prepare_unpolarized_raman" not in cls.__dict__ + assert "prepare_perturbation_q" not in cls.__dict__