New utilities, nac and polarization directions in Cartesian, and load…

…s of tests (#51)

The following improvements are applied:
- Add post-processing of the complex dielectric funciton in the infrared
- Add post-processing of the normal reflectivity in the infrared
- Change the convention for non-analytical and polarization directions in Cartesian coordinates
- Add loads of tests, and where possible use data regression

docs/source/howto/overrides.md

@@ -1,6 +1,6 @@
 (howto-overrides)=
 
-# How-to use protocols and overrides
+# Protocols and overrides
 
 :::{important}
 The following how-to assume that you are familiar with submitting the workflows of the package.
@@ -17,7 +17,7 @@ as you probably noticed in the [tutorials](tutorials).
 Not necesseraly all the WorkChains have the `get_builder_from_protocol`. E.g. The {class}`~aiida_vibroscopy.workflows.dielectric.numerical_derivatives`.
 :::
 
-## How to use in general
+## General usage
 
 Depending on the WorkChain, you will need to specify some minimal, _but required_, inputs. Here is an example:
 
@@ -48,7 +48,7 @@ Out[1]:
 
 ```
 
-## How to use overrides (beginner)
+## Overrides (beginner)
 
 As stated at the beginning, you might need to tweak your builder inputs before submitting. Instead of modifying the inputs via the builder, you can do it via `overrides`. The way the overrides must be secified __is WorkChain dependent__. The structure should be the same as the input specification of the specific WorkChain, so always refer to the inputs of the particular workflow, which you can find [here](topics-workflows).
 
@@ -169,7 +169,7 @@ builder = DielectricWorkChain.get_builder_from_protocol(code=code, structure=str
 
 You got the gist, so it will be the same for the `phonon` inputs, which corresponds to the `PhononWorkChain`.
 
-## How to overrides (advanced)
+## Overrides (advanced)
 
 Once you got used to using the overrides, you will notice your python script will become quite a mess.
 The idea is to write the overrides in a separate file, using for example the convenient `YAML` format,
@@ -199,7 +199,7 @@ phonon:
     scf:
         pw:
             kpoints_distance: 0.2
-            pw:
+            parameters:
                 ...
 dielectric:
     central_difference:
@@ -209,7 +209,7 @@ dielectric:
             ...
 ```
 
-## Get automated inputs for insulators, metals, and magnetism
+## Automated inputs for insulators, metals, and magnetism
 
 Inputs might change depending on the nature of the material: insulating, metallic, a form of magnetism.
 This can be the need of specifying a smearing or a magnetic configuration.

docs/source/howto/postprocess.md

@@ -1,6 +1,6 @@
 (howto-postprocess)=
 
-# How-to post-process data
+# Post-process data
 
 Here you will find some _how-tos_ on how to post-process the `VibrationalData` and `PhonopyData`.
 These can be generated via the `HarmonicWorkChain`, the `IRamanSpectraWorkChain` or the `PhononWorkChain`.
@@ -10,13 +10,12 @@ Please, have a look at the [tutorials](tutorials) for an overview, and at the sp
 You can always access to the detailed description of a function/method putting a `?` question mark in front
 of the function and press enter. This will print the description of the function, with detailed information
 on inputs, outputs and purpose.
-``` {note}
+
 A rendered version can also be found in the documentation referring to
 the dedicated [API section](reference).
-```
 :::
 
-## How to get powder spectra
+## Powder spectra
 
 You can get the powder infrared and/or Raman spectra from a computed {{ vibrational_data }}. For Raman:
 
@@ -28,8 +27,8 @@ vibro = load_node(IDENTIFIER) # a VibrationalData node
 polarized_intensities, depolarized_intensities, frequencies, labels = vibro.run_powder_raman_intensities(frequency_laser=532, temperature=300)
 ```
 
-The total powder intensity will be the some of the polarized (backscattering geometry)
-and depolarized (90º geometry) intensities:
+The total powder intensity will be the some of the polarized (backscattering geometry, HH/VV setup)
+and depolarized (90º geometry, HV setup) intensities:
 
 ```python
 total = polarized_intensities + depolarized_intensities
@@ -44,7 +43,7 @@ intensities, frequencies, labels = vibro.run_powder_ir_intensities()
 The `labels` output is referring to the irreducible representation labels of the modes.
 
 
-## How to get single crystal spectra
+## Single crystal polarized spectra
 
 ::: {important}
 It is extremely important that you match the experimental conditions, meaning the direction of
@@ -53,7 +52,7 @@ on the convention used, both in the code and in the experiments.
 :::
 
 ::: {note} Convention
-In the following methods, the direction of the light must be given in crystal/fractional coordinates.
+In the following methods, the direction of the light must be given in **Cartesian coordinates**.
 :::
 
 You can get the single crystal infrared and/or Raman spectra from a computed {{ vibrational_data }}. For Raman:
@@ -86,7 +85,7 @@ incoming = [0,0,1] # light polatization of the incident laser beam
 intensities, frequencies, labels = vibro.run_single_crystal_ir_intensities(pol_incoming=incoming)
 ```
 
-::: {admonition} Cartesian coordinates
+<!-- ::: {admonition} Cartesian coordinates
 :class: hint
 To get the direction in Cartesian coordinates expressed in crystal coordinates of the primitive cell, you can use the following
 snippet
@@ -99,14 +98,9 @@ incoming_cartesian = [0,0,1]
 inv_cell = np.linalg.inv(cell)
 incoming = np.dot(invcell.T, incoming_cartesian)
 ```
-For the q-direction instead you need to transform them into reciprocal space crystal coordinates, as follows
-```python
-q_nac_cartesian = [0,0,1]
-q_nac_crystal = np.dot(cell, q_nac_cartesian)
-```
-:::
+::: -->
 
-## How to get the IR/Raman active modes
+## IR/Raman active modes
 
 To get the active modes from a computed {{ vibrational_data }}:
 

docs/source/howto/understand.md

@@ -1,6 +1,6 @@
 (howto-understand)=
 
-# How-to understand the input/builder structure
+# Understand the input/builder structure
 
 In AiiDA the CalcJobs and WorkChains have usually nested inputs and different options on how to run the calculation
 and/or workflows. To understand the nested input structure of CalcJobs/Workflows, we made them available in a clickable

src/aiida_vibroscopy/calculations/numerical_derivatives_utils.py

@@ -34,6 +34,32 @@
     'compute_nac_parameters'
 )
 
+# def map_polarization(polarization: np.ndarray, cell: np.ndarray, sign: Literal[-1, 1]) -> np.ndarray:
+#     """Map the polarization within a quantum of polarization.
+
+#     It maps P(dE) in [0, Pq] and P(-dE) in [-Pq, 0].
+
+#     :param polarization: (3,) vector in Cartesian coordinates, Ry atomic units
+#     :param cell: (3, 3) matrix cell in Cartesian coordinates
+#         (rows are lattice vectors, i.e. cell[0] = v1, ...)
+#     :param sign: sign (1 or -1) to select the side of the branch
+#     :return: (3,) vector in Cartesian coordinates, Ry atomic units
+#     """
+#     inv_cell = np.linalg.inv(cell)
+#     lengths = np.sqrt(np.sum(cell**2, axis=1)) / CONSTANTS.bohr_to_ang  # in Bohr
+#     volume = float(abs(np.dot(np.cross(cell[0], cell[1]), cell[2]))) / CONSTANTS.bohr_to_ang**3  # in Bohr^3
+
+#     pol_quantum = np.sqrt(2) * lengths / volume
+#     pol_crys = lengths * np.dot(polarization, inv_cell)  # in Bohr
+
+#     pol_branch = pol_quantum * (pol_crys // pol_quantum)
+#     pol_crys -= pol_branch
+
+#     if sign < 0:
+#         pol_crys -= pol_quantum
+
+#     return np.dot(pol_crys / lengths, cell)
+
 
 def get_central_derivatives_coefficients(accuracy: int, order: int) -> list[int]:
     r"""Return an array with the central derivatives coefficients.
@@ -174,7 +200,7 @@ def compute_susceptibility_derivatives(
         * Units are pm/V for non-linear optical susceptibility
         * Raman tensors indecis: (atomic,  atomic displacement, electric field, electric field)
 
-    :return: dictionary with :class:`~aiida.orm.ArrayData` having keys:
+    :return: dictionaries of numerical accuracies with :class:`~aiida.orm.ArrayData` having keys:
         * `raman_tensors` containing (num_atoms, 3, 3, 3) arrays;
         * `nlo_susceptibility` containing (3, 3, 3) arrays;
         * `units` as :class:`~aiida.orm.Dict` containing the units of the tensors.
@@ -211,10 +237,11 @@ def compute_susceptibility_derivatives(
         data = get_trajectories_from_symmetries(
             preprocess_data=preprocess_data, data=raw_data, data_0=data_0, accuracy_order=accuracy_order.value
         )
+    else:
+        data = raw_data
 
     # Conversion factors
-    # dchi_factor = evang_to_rybohr * CONSTANTS.bohr_to_ang**2  # --> angstrom^2
-    dchi_factor = (evang_to_rybohr * CONSTANTS.bohr_to_ang**2) / volume  # --> angstrom^-1
+    dchi_factor = (evang_to_rybohr * CONSTANTS.bohr_to_ang**2) / volume  # --> 4*pi / angstrom
     chi2_factor = 0.5 * (4 * np.pi) * 100 / (volume_au_units * efield_au_to_si)  # --> pm/Volt
 
     # Variables
@@ -403,6 +430,8 @@ def compute_nac_parameters(
         data = get_trajectories_from_symmetries(
             preprocess_data=preprocess_data, data=raw_data, data_0=data_0, accuracy_order=accuracy_order.value
         )
+    else:
+        data = raw_data
 
     # Conversion factors
     bec_factor = evang_to_rybohr / np.sqrt(2)