Vibrational spectra of Si#
In this tutorial you will learn how to calculate the infrared and Raman spectra of silicon using aiida-vibroscopy.
In this tutorial we will make use of the silicon structure to give you an overall understanding of the usage of the package. If you are interested in other features, please have a look at the other tutorials.
Let’s get started!
Finite electric fields#
In the Placzek approximations (good for insulators), all the needed quantities for Raman can be computed as finite difference using the electric fields, as demonstrated in the previous tutorial for the dielectric tensor. Now we a more challenging task, as a third order derivative is needed: the Raman tensor, defined as:
(1)#\[\begin{equation}
\frac{\partial \chi_{ij}}{\partial \tau_{K,k}}
=
\frac{1}{\Omega}
\frac{\partial^2 F_{K,k}}{\partial \mathcal{E}_i \partial \mathcal{E}_j}
\end{equation}\]
The same theory we saw in previous tutorials applies here. Note that for computing second order derivatives of forces we need a slightly different formula:
(2)#\[\begin{equation}
\left . \frac{\partial^2 f(x)}{\partial x^2} \right|_{x=0}
=
\frac{1}{h^2}
\left [
f(h) -2f(0) +f(-h)
\right ]
+ \mathcal{O}(h^2)
\end{equation}\]
Have a look at the finite difference coefficients for coefficients of any accuracy order.
Show code cell content
from local_module import load_temp_profile
# If you download this file, you can run it with your own profile.
# Put these lines instead:
# from aiida import load_profile
# load_profile()
data = load_temp_profile(
name="raman-tutorial",
add_computer=True,
add_pw_code=True,
add_sssp=True,
)
The IRamanSpectraWorkChain#
For computing the spectra we need both phonons and Raman tensors. For understanding better the ingredients, refer to the formulation and to the references
Note
Also the Raman tensors can be computed using perturbative approaches, such as the so called 2n+1 theorem and the second order response of the density matrix. Quantum ESPRESSO can carry out such simulation using the PHonon code (ph.x routine). Nevertheless, only LDA and norm-conserving pseudopotential can be employed.
Let’s import the WorkChain and run it! We use the get_builder_from_protocol to get a prefilled builder with all inputs.
Note
These inputs should be considered not as converged parameters, but as a good starting point. You may also need to tweak some parameters, e.g. add magnetization etc., depending on your case.
Show code cell content
from aiida.plugins import DbImporterFactory
CodDbImporter = DbImporterFactory('cod')
cod = CodDbImporter()
results = cod.query(id='1526655') # Si 1526655
structure = results[0].get_aiida_structure() # it has 8 atoms
from aiida.plugins import WorkflowFactory
from aiida.engine import run_get_node
IRamanSpectraWorkChain = WorkflowFactory("vibroscopy.spectra.iraman")
builder = IRamanSpectraWorkChain.get_builder_from_protocol(
code=data.pw_code,
structure=structure,
protocol="fast",
)
results, calc = run_get_node(builder)
Show code cell output
05/24/2023 11:02:11 AM <4148882> aiida.orm.nodes.process.workflow.workchain.WorkChainNode: [REPORT] [106|IRamanSpectraWorkChain|run_spectra]: submitting `HarmonicWorkChain` <PK=108>
05/24/2023 11:02:11 AM <4148882> aiida.engine.processes.functions: [WARNING] function `generate_preprocess_data` has invalid type hints: unsupported operand type(s) for |: 'AbstractNodeMeta' and 'NoneType'
05/24/2023 11:02:11 AM <4148882> aiida.engine.processes.functions: [WARNING] function `generate_phonopy_data` has invalid type hints: unsupported operand type(s) for |: 'AbstractNodeMeta' and 'NoneType'
05/24/2023 11:02:12 AM <4148882> aiida.orm.nodes.process.workflow.workchain.WorkChainNode: [REPORT] [108|HarmonicWorkChain|run_phonon]: submitting `PhononWorkChain` <PK=114>
05/24/2023 11:02:12 AM <4148882> aiida.orm.nodes.process.workflow.workchain.WorkChainNode: [REPORT] [108|HarmonicWorkChain|run_dielectric]: submitting `DielectricWorkChain` <PK=118>
05/24/2023 11:02:14 AM <4148882> aiida.orm.nodes.process.workflow.workchain.WorkChainNode: [REPORT] [118|DielectricWorkChain|run_base_scf]: launching base scf PwBaseWorkChain<131>
05/24/2023 11:02:15 AM <4148882> aiida.orm.nodes.process.workflow.workchain.WorkChainNode: [REPORT] [131|PwBaseWorkChain|run_process]: launching PwCalculation<134> iteration #1
05/24/2023 11:02:16 AM <4148882> aiida.orm.nodes.process.workflow.workchain.WorkChainNode: [REPORT] [114|PhononWorkChain|run_base_supercell]: launching base supercell scf PwBaseWorkChain<140>
05/24/2023 11:02:16 AM <4148882> aiida.orm.nodes.process.workflow.workchain.WorkChainNode: [REPORT] [140|PwBaseWorkChain|run_process]: launching PwCalculation<143> iteration #1
05/24/2023 11:02:35 AM <4148882> aiida.orm.nodes.process.workflow.workchain.WorkChainNode: [REPORT] [131|PwBaseWorkChain|sanity_check_insufficient_bands]: PwCalculation<134> run with smearing and highest band is occupied
05/24/2023 11:02:35 AM <4148882> aiida.orm.nodes.process.workflow.workchain.WorkChainNode: [REPORT] [131|PwBaseWorkChain|sanity_check_insufficient_bands]: BandsData<147> has invalid occupations: Occupation of 0.009307178051580323 at last band lkn<0,0,20>
05/24/2023 11:02:35 AM <4148882> aiida.orm.nodes.process.workflow.workchain.WorkChainNode: [REPORT] [131|PwBaseWorkChain|sanity_check_insufficient_bands]: PwCalculation<134> had insufficient bands
05/24/2023 11:02:35 AM <4148882> aiida.orm.nodes.process.workflow.workchain.WorkChainNode: [REPORT] [131|PwBaseWorkChain|sanity_check_insufficient_bands]: Action taken: increased number of bands to 24 and restarting from the previous charge density.
05/24/2023 11:02:35 AM <4148882> aiida.orm.nodes.process.workflow.workchain.WorkChainNode: [REPORT] [131|PwBaseWorkChain|inspect_process]: PwCalculation<134> finished successfully but a handler was triggered, restarting
05/24/2023 11:02:35 AM <4148882> aiida.orm.nodes.process.workflow.workchain.WorkChainNode: [REPORT] [131|PwBaseWorkChain|run_process]: launching PwCalculation<152> iteration #2
05/24/2023 11:02:36 AM <4148882> aiida.orm.nodes.process.workflow.workchain.WorkChainNode: [REPORT] [140|PwBaseWorkChain|sanity_check_insufficient_bands]: PwCalculation<143> run with smearing and highest band is occupied
05/24/2023 11:02:36 AM <4148882> aiida.orm.nodes.process.workflow.workchain.WorkChainNode: [REPORT] [140|PwBaseWorkChain|sanity_check_insufficient_bands]: BandsData<154> has invalid occupations: Occupation of 0.009307178051003007 at last band lkn<0,0,20>
05/24/2023 11:02:36 AM <4148882> aiida.orm.nodes.process.workflow.workchain.WorkChainNode: [REPORT] [140|PwBaseWorkChain|sanity_check_insufficient_bands]: PwCalculation<143> had insufficient bands
05/24/2023 11:02:36 AM <4148882> aiida.orm.nodes.process.workflow.workchain.WorkChainNode: [REPORT] [140|PwBaseWorkChain|sanity_check_insufficient_bands]: Action taken: increased number of bands to 24 and restarting from the previous charge density.
05/24/2023 11:02:36 AM <4148882> aiida.orm.nodes.process.workflow.workchain.WorkChainNode: [REPORT] [140|PwBaseWorkChain|inspect_process]: PwCalculation<143> finished successfully but a handler was triggered, restarting
05/24/2023 11:02:37 AM <4148882> aiida.orm.nodes.process.workflow.workchain.WorkChainNode: [REPORT] [140|PwBaseWorkChain|run_process]: launching PwCalculation<160> iteration #2
05/24/2023 11:02:47 AM <4148882> aiida.orm.nodes.process.workflow.workchain.WorkChainNode: [REPORT] [131|PwBaseWorkChain|results]: work chain completed after 2 iterations
05/24/2023 11:02:47 AM <4148882> aiida.orm.nodes.process.workflow.workchain.WorkChainNode: [REPORT] [131|PwBaseWorkChain|on_terminated]: remote folders will not be cleaned
05/24/2023 11:02:48 AM <4148882> aiida.orm.nodes.process.workflow.workchain.WorkChainNode: [REPORT] [118|DielectricWorkChain|run_nscf]: launching base scf PwBaseWorkChain<168>
05/24/2023 11:02:49 AM <4148882> aiida.orm.nodes.process.workflow.workchain.WorkChainNode: [REPORT] [168|PwBaseWorkChain|run_process]: launching PwCalculation<171> iteration #1
05/24/2023 11:02:52 AM <4148882> aiida.orm.nodes.process.workflow.workchain.WorkChainNode: [REPORT] [140|PwBaseWorkChain|results]: work chain completed after 2 iterations
05/24/2023 11:02:52 AM <4148882> aiida.orm.nodes.process.workflow.workchain.WorkChainNode: [REPORT] [140|PwBaseWorkChain|on_terminated]: remote folders will not be cleaned
05/24/2023 11:02:54 AM <4148882> aiida.orm.nodes.process.workflow.workchain.WorkChainNode: [REPORT] [114|PhononWorkChain|run_forces]: submitting `PwBaseWorkChain` <PK=180> with supercell n.o 1
05/24/2023 11:02:56 AM <4148882> aiida.orm.nodes.process.workflow.workchain.WorkChainNode: [REPORT] [180|PwBaseWorkChain|run_process]: launching PwCalculation<183> iteration #1
05/24/2023 11:03:07 AM <4148882> aiida.orm.nodes.process.workflow.workchain.WorkChainNode: [REPORT] [168|PwBaseWorkChain|results]: work chain completed after 1 iterations
05/24/2023 11:03:07 AM <4148882> aiida.orm.nodes.process.workflow.workchain.WorkChainNode: [REPORT] [168|PwBaseWorkChain|on_terminated]: remote folders will not be cleaned
05/24/2023 11:03:12 AM <4148882> aiida.orm.nodes.process.workflow.workchain.WorkChainNode: [REPORT] [118|DielectricWorkChain|run_null_field_scfs]: launching PwBaseWorkChain<197> with null electric field 0
05/24/2023 11:03:12 AM <4148882> aiida.orm.nodes.process.workflow.workchain.WorkChainNode: [REPORT] [118|DielectricWorkChain|run_null_field_scfs]: launching PwBaseWorkChain<198> with null electric field 1
05/24/2023 11:03:14 AM <4148882> aiida.orm.nodes.process.workflow.workchain.WorkChainNode: [REPORT] [197|PwBaseWorkChain|run_process]: launching PwCalculation<201> iteration #1
05/24/2023 11:03:15 AM <4148882> aiida.orm.nodes.process.workflow.workchain.WorkChainNode: [REPORT] [198|PwBaseWorkChain|run_process]: launching PwCalculation<204> iteration #1
05/24/2023 11:04:04 AM <4148882> aiida.orm.nodes.process.workflow.workchain.WorkChainNode: [REPORT] [180|PwBaseWorkChain|sanity_check_insufficient_bands]: PwCalculation<183> run with smearing and highest band is occupied
05/24/2023 11:04:04 AM <4148882> aiida.orm.nodes.process.workflow.workchain.WorkChainNode: [REPORT] [180|PwBaseWorkChain|sanity_check_insufficient_bands]: BandsData<208> has invalid occupations: Occupation of 0.008059168787033277 at last band lkn<0,0,20>
05/24/2023 11:04:04 AM <4148882> aiida.orm.nodes.process.workflow.workchain.WorkChainNode: [REPORT] [180|PwBaseWorkChain|sanity_check_insufficient_bands]: PwCalculation<183> had insufficient bands
05/24/2023 11:04:04 AM <4148882> aiida.orm.nodes.process.workflow.workchain.WorkChainNode: [REPORT] [180|PwBaseWorkChain|sanity_check_insufficient_bands]: Action taken: increased number of bands to 24 and restarting from the previous charge density.
05/24/2023 11:04:04 AM <4148882> aiida.orm.nodes.process.workflow.workchain.WorkChainNode: [REPORT] [180|PwBaseWorkChain|inspect_process]: PwCalculation<183> finished successfully but a handler was triggered, restarting
05/24/2023 11:04:05 AM <4148882> aiida.orm.nodes.process.workflow.workchain.WorkChainNode: [REPORT] [180|PwBaseWorkChain|run_process]: launching PwCalculation<213> iteration #2
05/24/2023 11:04:53 AM <4148882> aiida.orm.nodes.process.workflow.workchain.WorkChainNode: [REPORT] [180|PwBaseWorkChain|results]: work chain completed after 2 iterations
05/24/2023 11:04:53 AM <4148882> aiida.orm.nodes.process.workflow.workchain.WorkChainNode: [REPORT] [180|PwBaseWorkChain|on_terminated]: cleaned remote folders of calculations: 183 213
05/24/2023 11:04:57 AM <4148882> aiida.orm.nodes.process.workflow.workchain.WorkChainNode: [REPORT] [114|PhononWorkChain|on_terminated]: cleaned remote folders of calculations: 143 160 183 213
05/24/2023 11:05:02 AM <4148882> aiida.orm.nodes.process.workflow.workchain.WorkChainNode: [REPORT] [197|PwBaseWorkChain|results]: work chain completed after 1 iterations
05/24/2023 11:05:03 AM <4148882> aiida.orm.nodes.process.workflow.workchain.WorkChainNode: [REPORT] [197|PwBaseWorkChain|on_terminated]: remote folders will not be cleaned
05/24/2023 11:05:10 AM <4148882> aiida.orm.nodes.process.workflow.workchain.WorkChainNode: [REPORT] [198|PwBaseWorkChain|results]: work chain completed after 1 iterations
05/24/2023 11:05:10 AM <4148882> aiida.orm.nodes.process.workflow.workchain.WorkChainNode: [REPORT] [198|PwBaseWorkChain|on_terminated]: remote folders will not be cleaned
05/24/2023 11:05:11 AM <4148882> aiida.orm.nodes.process.workflow.workchain.WorkChainNode: [REPORT] [118|DielectricWorkChain|run_electric_field_scfs]: launching PwBaseWorkChain<232> with electric field index 2 and sign 1.0 iteration #0
05/24/2023 11:05:13 AM <4148882> aiida.orm.nodes.process.workflow.workchain.WorkChainNode: [REPORT] [118|DielectricWorkChain|run_electric_field_scfs]: launching PwBaseWorkChain<235> with electric field index 3 and sign 1.0 iteration #0
05/24/2023 11:05:15 AM <4148882> aiida.orm.nodes.process.workflow.workchain.WorkChainNode: [REPORT] [232|PwBaseWorkChain|run_process]: launching PwCalculation<238> iteration #1
05/24/2023 11:05:16 AM <4148882> aiida.orm.nodes.process.workflow.workchain.WorkChainNode: [REPORT] [235|PwBaseWorkChain|run_process]: launching PwCalculation<241> iteration #1
05/24/2023 11:15:20 AM <4148882> aiida.orm.nodes.process.workflow.workchain.WorkChainNode: [REPORT] [232|PwBaseWorkChain|results]: work chain completed after 1 iterations
05/24/2023 11:15:20 AM <4148882> aiida.orm.nodes.process.workflow.workchain.WorkChainNode: [REPORT] [232|PwBaseWorkChain|on_terminated]: remote folders will not be cleaned
05/24/2023 11:15:34 AM <4148882> aiida.orm.nodes.process.workflow.workchain.WorkChainNode: [REPORT] [235|PwBaseWorkChain|results]: work chain completed after 1 iterations
05/24/2023 11:15:34 AM <4148882> aiida.orm.nodes.process.workflow.workchain.WorkChainNode: [REPORT] [235|PwBaseWorkChain|on_terminated]: remote folders will not be cleaned
05/24/2023 11:15:35 AM <4148882> aiida.orm.nodes.process.workflow.workchain.WorkChainNode: [REPORT] [118|DielectricWorkChain|run_electric_field_scfs]: launching PwBaseWorkChain<254> with electric field index 2 and sign 1.0 iteration #1
05/24/2023 11:15:36 AM <4148882> aiida.orm.nodes.process.workflow.workchain.WorkChainNode: [REPORT] [118|DielectricWorkChain|run_electric_field_scfs]: launching PwBaseWorkChain<257> with electric field index 3 and sign 1.0 iteration #1
05/24/2023 11:15:38 AM <4148882> aiida.orm.nodes.process.workflow.workchain.WorkChainNode: [REPORT] [254|PwBaseWorkChain|run_process]: launching PwCalculation<260> iteration #1
05/24/2023 11:15:39 AM <4148882> aiida.orm.nodes.process.workflow.workchain.WorkChainNode: [REPORT] [257|PwBaseWorkChain|run_process]: launching PwCalculation<263> iteration #1
05/24/2023 11:18:40 AM <4148882> aiida.orm.nodes.process.workflow.workchain.WorkChainNode: [REPORT] [254|PwBaseWorkChain|results]: work chain completed after 1 iterations
05/24/2023 11:18:40 AM <4148882> aiida.orm.nodes.process.workflow.workchain.WorkChainNode: [REPORT] [254|PwBaseWorkChain|on_terminated]: remote folders will not be cleaned
05/24/2023 11:20:25 AM <4148882> aiida.orm.nodes.process.workflow.workchain.WorkChainNode: [REPORT] [257|PwBaseWorkChain|results]: work chain completed after 1 iterations
05/24/2023 11:20:25 AM <4148882> aiida.orm.nodes.process.workflow.workchain.WorkChainNode: [REPORT] [257|PwBaseWorkChain|on_terminated]: remote folders will not be cleaned
05/24/2023 11:20:27 AM <4148882> aiida.orm.nodes.process.workflow.workchain.WorkChainNode: [REPORT] [118|DielectricWorkChain|run_numerical_derivatives]: launching NumericalDerivativesWorkChain<281> for computing numerical derivatives.
/aiida-plugins/aiida-vibroscopy/src/aiida_vibroscopy/calculations/symmetry.py:162: UserWarning: Symmetry of dChi/dr tensors is largely broken.
The max difference is:
21.971663701634725
warnings.warn('\n'.join(lines))
/aiida-plugins/aiida-vibroscopy/src/aiida_vibroscopy/calculations/symmetry.py:162: UserWarning: Symmetry of dChi/dr tensors is largely broken.
The max difference is:
94.52835279718761
warnings.warn('\n'.join(lines))
/aiida-plugins/aiida-vibroscopy/src/aiida_vibroscopy/calculations/symmetry.py:162: UserWarning: Symmetry of dChi/dr tensors is largely broken.
The max difference is:
118.71391582903925
warnings.warn('\n'.join(lines))
05/24/2023 11:20:32 AM <4148882> aiida.orm.nodes.process.workflow.workchain.WorkChainNode: [REPORT] [118|DielectricWorkChain|on_terminated]: cleaned remote folders of calculations: 134 152 171 201 204 238 241 260 263
05/24/2023 11:20:36 AM <4148882> aiida.orm.nodes.process.workflow.workchain.WorkChainNode: [REPORT] [108|HarmonicWorkChain|on_terminated]: cleaned remote folders of calculations: 143 160 183 213 134 152 171 201 204 238 241 260 263
These are the results:
results
{'output_dielectric': AttributeDict({'fields_data': AttributeDict({'field_index_2': AttributeDict({'0': <TrajectoryData: uuid: d443ca1f-b30b-4b75-a209-e680a1f732bb (pk: 246)>, '1': <TrajectoryData: uuid: 91c958f3-48da-4f89-8fb9-8515069cdf03 (pk: 268)>}), 'field_index_3': AttributeDict({'0': <TrajectoryData: uuid: d54a99d5-5790-43f5-8461-1241d1e1577e (pk: 250)>, '1': <TrajectoryData: uuid: 9d1234d1-e77e-442f-b01e-58e92518a8bc (pk: 272)>})}), 'tensors': AttributeDict({'numerical_accuracy_2_step_2': <ArrayData: uuid: 5bac9170-8bd0-4bd7-b122-5a54c70e7e9f (pk: 294)>, 'numerical_accuracy_2_step_1': <ArrayData: uuid: 01953989-f02d-414a-83ac-88efb148add0 (pk: 296)>, 'numerical_accuracy_4': <ArrayData: uuid: 9b1701ed-4157-421a-b7dc-2ab709ed323f (pk: 298)>}), 'critical_electric_field': <Float: uuid: cfebef1d-8916-43b7-b6a3-c76297cd305d (pk: 190) value: 0.00091230365735765>, 'electric_field_step': <Float: uuid: c9355f6c-6283-4053-af16-07383a8e565f (pk: 194) value: 0.00045615182867883>, 'units': <Dict: uuid: 4ffca53f-8eac-4ae6-9341-658e47444649 (pk: 292)>}),
'output_phonon': AttributeDict({'supercells': AttributeDict({'supercell_1': <StructureData: uuid: 64812fa5-1e5b-4f87-9517-e9d978fab7f0 (pk: 178)>}), 'supercells_forces': AttributeDict({'forces_0': <TrajectoryData: uuid: 8da64d92-77a3-45b1-9be9-ac796b97d39b (pk: 174)>, 'forces_1': <TrajectoryData: uuid: 74a664bb-7bad-4b4c-b2a8-c565368ea7b5 (pk: 217)>}), 'phonopy_data': <PhonopyData: uuid: caa40178-e9ff-46af-9d77-1bd8a68e71fd (pk: 221)>}),
'vibrational_data': AttributeDict({'numerical_accuracy_2_step_2': <VibrationalData: uuid: 08aab549-cc64-4c71-a153-e61dd7f8e0db (pk: 302)>, 'numerical_accuracy_2_step_1': <VibrationalData: uuid: e2b78311-d5e9-461c-929d-ec45c21dc342 (pk: 306)>, 'numerical_accuracy_4': <VibrationalData: uuid: b382f41d-2a41-457b-a510-b7cb7d471411 (pk: 310)>})}
The VibrationalData#
In aiida-vibroscopy we design a data type able to store all the phonon and dielectric properties, to ease the share of the data and make it as costistent as possible.
In this data, you can access to multiple information:
Structure (cell, sites, …)
Symmetry
Dielectric tensor (\(\epsilon^{\infty}\))
Born effective charges (\(Z^*\))
Raman tensors (\(\partial \chi / \partial \tau\))
Non-linear optical susceptibility (\(\chi^{(2)}\))
Again, Born effective charges and non-linear optical susceptibility are null in silicon, due to symmetry. In the next tutorial we are going to compute such quantities for AlAs, which are non vanishing.
vibro = calc.outputs.vibrational_data.numerical_accuracy_4
print("The VibrationalData: ", vibro, "\n")
print("The dielectric tensor: ", "\n", vibro.dielectric.round(5), "\n")
print("The Raman tensors (in 1/Angstrom): ", "\n", vibro.raman_tensors.round(5), "\n")
The VibrationalData: uuid: b382f41d-2a41-457b-a510-b7cb7d471411 (pk: 310)
The dielectric tensor:
[[11.08524 0. -0. ]
[ 0. 11.08524 -0. ]
[-0. -0. 11.08524]]
The Raman tensors (in 1/Angstrom):
[[[[-0. 0. -0. ]
[ 0. -0. -0.08896]
[-0. -0.08896 -0. ]]
[[ 0. -0. -0.08896]
[ 0. 0. -0. ]
[-0.08896 -0. -0. ]]
[[-0. -0.08896 -0. ]
[-0.08896 0. -0. ]
[-0. -0. -0. ]]]
[[[ 0. 0. 0. ]
[-0. 0. 0.08896]
[ 0. 0.08896 0. ]]
[[-0. 0. 0.08896]
[ 0. -0. -0. ]
[ 0.08896 0. 0. ]]
[[ 0. 0.08896 0. ]
[ 0.08896 0. 0. ]
[ 0. 0. 0. ]]]]
And we can also check how much the Raman tensors, which are third order derivatives, are sensible to the numerical accuracy chosen. Usually, the Raman tensor is expressed in Å\(^2\). To get it, one can use the convention in the code:
(3)#\[\begin{equation}
\left [
\frac{\partial \chi}{\partial \tau}
\right] (\text{Å}^2)
=
\Omega_{unitcell}
\left [
\frac{\partial \chi}{\partial \tau}
\right] (\text{Å}^{-1})
\end{equation}\]
vol = calc.outputs.vibrational_data.numerical_accuracy_4.get_unitcell().get_cell_volume()
print(
calc.outputs.vibrational_data.numerical_accuracy_4.raman_tensors[1,0,1,2]*vol,
calc.outputs.vibrational_data.numerical_accuracy_2_step_1.raman_tensors[1,0,1,2]*vol,
calc.outputs.vibrational_data.numerical_accuracy_2_step_2.raman_tensors[1,0,1,2]*vol,
)
13.86057079085759 13.895477711258131 14.00019847245988
We can see that the values are converged within 0.1 Å\(^2\) with all three steps.
Plotting the powder Raman spectra#
Now, we can get directly from the data the Raman intensities.
polarized_intensities, depolarized_intensities, frequencies, labels = vibro.run_powder_raman_intensities(frequency_laser=532, temperature=300)
print(frequencies, "\n", labels)
[531.91425155 531.91425155 531.91425155]
['T2g', 'T2g', 'T2g']
We have 3 degenerate modes. Quite expected for silicon, which is very symmetric. To plot the actual spectra:
from aiida_vibroscopy.utils.plotting import get_spectra_plot
total_intensities = polarized_intensities + depolarized_intensities
plt = get_spectra_plot(frequencies, total_intensities)
plt.show()
Analysing the workflow#
Look how the complexity of computing a Raman spectra has been handled fully automatically.
%verdi process status {calc.pk}
IRamanSpectraWorkChain<106> Finished [0] [2:if_(should_run_average)]
└── HarmonicWorkChain<108> Finished [0] [4:if_(should_run_phonopy)]
├── generate_preprocess_data<109> Finished [0]
├── PhononWorkChain<114> Finished [0] [7:if_(should_run_phonopy)]
│ ├── generate_preprocess_data<119> Finished [0]
│ ├── get_supercell<135> Finished [0]
│ ├── create_kpoints_from_distance<137> Finished [0]
│ ├── PwBaseWorkChain<140> Finished [0] [3:results]
│ │ ├── PwCalculation<143> Finished [0]
│ │ └── PwCalculation<160> Finished [0]
│ ├── get_supercells_with_displacements<177> Finished [0]
│ ├── PwBaseWorkChain<180> Finished [0] [3:results]
│ │ ├── PwCalculation<183> Finished [0]
│ │ └── PwCalculation<213> Finished [0]
│ └── generate_phonopy_data<220> Finished [0]
├── DielectricWorkChain<118> Finished [0] [11:results]
│ ├── create_kpoints_from_distance<120> Finished [0]
│ ├── create_directional_kpoints<124> Finished [0]
│ ├── create_directional_kpoints<127> Finished [0]
│ ├── PwBaseWorkChain<131> Finished [0] [3:results]
│ │ ├── PwCalculation<134> Finished [0]
│ │ └── PwCalculation<152> Finished [0]
│ ├── PwBaseWorkChain<168> Finished [0] [3:results]
│ │ └── PwCalculation<171> Finished [0]
│ ├── compute_critical_electric_field<189> Finished [0]
│ ├── get_accuracy_from_critical_field<191> Finished [0]
│ ├── get_electric_field_step<193> Finished [0]
│ ├── PwBaseWorkChain<197> Finished [0] [3:results]
│ │ └── PwCalculation<201> Finished [0]
│ ├── PwBaseWorkChain<198> Finished [0] [3:results]
│ │ └── PwCalculation<204> Finished [0]
│ ├── PwBaseWorkChain<232> Finished [0] [3:results]
│ │ └── PwCalculation<238> Finished [0]
│ ├── PwBaseWorkChain<235> Finished [0] [3:results]
│ │ └── PwCalculation<241> Finished [0]
│ ├── PwBaseWorkChain<254> Finished [0] [3:results]
│ │ └── PwCalculation<260> Finished [0]
│ ├── PwBaseWorkChain<257> Finished [0] [3:results]
│ │ └── PwCalculation<263> Finished [0]
│ ├── subtract_residual_forces<276> Finished [0]
│ └── NumericalDerivativesWorkChain<281> Finished [0] [None]
│ ├── generate_preprocess_data<282> Finished [0]
│ ├── compute_nac_parameters<284> Finished [0]
│ ├── compute_susceptibility_derivatives<288> Finished [0]
│ ├── join_tensors<293> Finished [0]
│ ├── join_tensors<295> Finished [0]
│ └── join_tensors<297> Finished [0]
├── elaborate_tensors<299> Finished [0]
├── generate_vibrational_data_from_phonopy<301> Finished [0]
├── elaborate_tensors<303> Finished [0]
├── generate_vibrational_data_from_phonopy<305> Finished [0]
├── elaborate_tensors<307> Finished [0]
└── generate_vibrational_data_from_phonopy<309> Finished [0]
Convergence with k points#
As a final remark, we show here how the convergence with k points is speeded up by sampling in this directional manner, compared to the most ‘traditional’ uniform mesh.
Show code cell source
import numpy as np
import matplotlib.pyplot as plt
# Some options to make the plot nice
plt.rcParams['font.size'] = 18
plt.rcParams['font.family'] = "serif"
plt.rcParams['text.usetex'] = False
plt.rcParams['xtick.major.size'] = 7.0
plt.rcParams['xtick.minor.size'] = 4.
plt.rcParams['ytick.major.size'] = 7.0
plt.rcParams['ytick.minor.size'] = 4.
plt.rcParams['axes.linewidth'] = 1.2
plt.rcParams['legend.fontsize'] = 15
# Directional sampling
raman_new = np.array([18.445, 19.943, 20.442, 20.442])
kpoints_new = np.array([3*6*6, 3*12*12, 3*24*24, 3*58*58])
# Uniform sampling
raman = np.array([15.610, 18.762, 20.436])
kpoints = np.array([3*3*3, 6*6*6, 12*12*12])
# Plotting
fig, ax = plt.subplots(figsize=(6.5,5)) # dichiaro il canvas e i subplots
ax.plot(kpoints_new, raman_new, color='red', marker='o', fillstyle='none', ms=10, linewidth=1.5, linestyle='-', label='Directional sampling')
ax.plot(kpoints, raman, color='black', marker='o', fillstyle='none', ms=10, linewidth=1.5, linestyle='-', label='Uniform sampling')
ax.set_xlabel('# k points')
ax.set_ylabel(r'Raman tensor ($\AA^2$)')
ax.legend(bbox_to_anchor=(0., 1.02, 1., .102), loc='lower left',
ncol=1, mode="expand", borderaxespad=0., frameon=False)
# plt.show()
<matplotlib.legend.Legend at 0x7f7e01cbeac0>
Final considerations#
We managed to compute the Raman spectra of silicon fully ab initio! :tada:
While for silicon the spectra is trivial, have a look to the tutorial on different functionals, where you will see the power of using extended Hubbard functionals.
Learn more and in details
To learn the full sets of inputs, to use proficiently the get_builder_from_protocol and more, have a look at the following sections: