a-SiO2

Example: α-SiO₂

In this example, we perform a detailed simulation of the infrared (IR) and Raman spectra of α-SiO₂ (quartz) in the P3₂21 spacegroup.

We calculate the band intensities, assign the modes to irreducible representations, and also include first-principles linewidths obtained from a high-quality Phono3py calculation.

We check the convergence of the calculated IR and Raman intensities with respect to the plane-wave cutoff, k-point sampling and exchange-correlation functional, and we also check the convergence of the linewidths with respect to q-point sampling and Brillouin-zone integration method.

The Phono(3)py calculations used as the base for this example were prepared by A. Togo.

Calculations

1. Preparation

The input files generated from the Phono(3)py calculations are the structure (POSCAR), the second- and third-order force constants (fc2.hdf5, fc3.hdf5), and the Born effective-charge tensors (BORN). For post-processing with phonopy, fc2.hdf5 needs to be copied/renamed to force_constants.hdf5.

1. Generate a mesh.yaml or mesh.hdf5 file containing the Γ-point phonon frequencies and eigenvectors: phonopy --dim="6 6 3" --readfc --hdf5 --fc-symmetry --mesh="1 1 1" --eigenvectors; for older versions of Phonopy, you may need to use --fc_symmetry=1 or fc-symmetry=1 in place of --fc-symmetry

2. Generate an irreps.yaml file containing the grouping of the modes by Mulliken symbol: phonopy --dim="6 6 3" --readfc --hdf5 --fc-symmetry --irreps="0 0 0"

3. Generate the phonon linewidths for the Γ-point modes: phono3py --dim="2 2 2" --dim_fc2="6 6 3" --fc2 --fc3 -v --br --thm --mesh="48 48 48" --write_gamma --gp=0

2. IR spectrum

Using the BORN file provided with the calculaton, generate a simulated IR spectrum and peak table with the room-temperature (300 K) linewidths:

phonopy-ir --ir-reps --linewidth-hdf5="kappa-m484848-g0.hdf5" --linewidth-temperature=300

Sample VASP input files for calculating the Born charges - used for the convergence tests - can be found in the VASP-Files folder.

3. Raman spectrum

From the irreps.yaml file, and with reference to the character table for the D₃ point group, the following modes may show Raman activity: 4-9, 11-15, 17, 18, 20-23, and 25-27.

1. Generate displaced structures for the Raman-intensity calculations: phonopy-raman -d --bands="4 5 6 7 8 9 11 12 13 14 15 17 18 20 21 22 23 25 26 27"

2. Run dielectric-constant calculations on the displaced structures (sample VASP input files in the VASP-Files folder)

3. Collect the dielectric constants calculated for the displaced structures: phonopy-raman -r OUTCAR.*

4. Post-process the results to generate a peak table and simulated Raman spectrum: phonopy-raman -p --ir-reps --linewidth-hdf5="kappa-m484848-g0.hdf5" --linewidth-temperature=300

Results

1. IR spectrum

The simulated spectrum was compared to ATR-IR data from the RRUFF database (RRUFF IDs: R040031-1, R050125-1, R050125-1).

The calculations predict much sharper lines than in the three measurements. However, given the much better reproduction of the linewidths of the Raman-active modes (see below), we believe this to be an experimental artefact.

There is also a slight shift in the peak positions, possibly due to the use of the LDA functional and/or the exclusion of finite-temperature effects from the calculations, but this is much less of an issue.

Nonetheless, the correspondance between the simulated and measured spectra is at least sufficient to assign the major bands.

2. Raman spectrum

The simulated spectrum was compared to depolarised 514 nm Raman spectra from the RRUFF database (RRUFF IDs: R040031-3, R050125-3, R060604-3).

Both the relative band intensities and peak widths are in excellent agreement with the experimental measurements, with the main discrepancies being in the two low-frequency bands below 250 cm^-1.

As in the IR spectra, there are some noticeable shifts in band frequencies, but the correspondance is sufficient to assign the spectrum.

3. Peak table

The irreducible representation of the Γ-point modes is obtained as 4 A₁ + 4 A₂ + 8 E. From the character table, we expect the A₁ modes to be Raman active, the A₂ modes to be IR active, and the E modes to be both Raman and IR active.

ν [cm-1]	Ir. Rep.	IIR [e2 amu-1]	IRaman [e2 amu-1]	ΓT=300K [cm-1]	ρ
127.37	E	0.000	0.834	2.63	0.75
223.52	A1	Inactive	6.930	9.52	0.00
255.08	E	0.021	0.404	1.96	0.75
337.02	A1	Inactive	1.057	1.69	0.42
341.54	A2	0.202	Inactive	1.44	-
374.19	E	0.268	0.624	1.58	0.75
434.99	E	0.685	0.718	2.13	0.75
454.96	A1	Inactive	35.085	4.83	0.00
482.03	A2	0.362	Inactive	2.75	-
691.77	E	0.074	1.359	5.15	0.75
769.87	A2	0.213	Inactive	3.63	-
792.15	E	0.298	1.554	3.43	0.75
1070.57	E	3.441	1.366	4.15	0.75
1080.66	A2	1.804	Inactive	4.17	-
1085.88	A1	Inactive	1.113	3.56	0.72
1148.71	E	0.078	2.811	5.85	0.75

Some convergence checks

1. IR intensities

A set of Born effective charges were calculated using LDA as part of the Phono(3)py calculation. We further checked the convergence of the calculated IR intensities with respect to the plane-wave cutoff and k-point sampling used to calculate the charges, and we also tested the effect of using DFT functionals at different levels of approximation.

a. Plane-wave cutoff and k-point sampling

Cutoff and k-point tests were performed from a baseline of a 520 eV cutoff and a 6×6×6 Γ-centred mesh, which are the settings used for this compound in the Materials Project database (mp-6930).[1] Sample VASP input files for each parameter set can be found in the VASP-Files folder.

#	Cutoff [eV]	k-points
-	Togo	Togo
1	520	6×6×6
2	600	6×6×6
3	700	6×6×6
4	800	6×6×6
5	600	8×8×8
6	600	10×10×10
7	600	12×12×12

In this case, we observed negligable quantitative variation in the calculated intensities on either increasing the plane-wave cutoff or using denser k-point sampling. This is most likely because a) the baseline calculation is well converged, and b) SiO₂ is an insulator, and as such does not require a dense sampling mesh.

ν [cm-1]	Ir. Rep.	IIR [e2 amu-1]
		-	1	2	3	4	5	6	7
127.37	E	< 10-3	< 10-3	< 10-3	< 10-3	< 10-3	< 10-3	< 10-3	< 10-3
255.08	E	0.021	0.021	0.021	0.021	0.021	0.021	0.021	0.021
341.54	A2	0.202	0.202	0.202	0.202	0.202	0.202	0.202	0.202
374.19	E	0.268	0.268	0.268	0.268	0.268	0.268	0.268	0.268
434.99	E	0.685	0.685	0.685	0.685	0.685	0.685	0.685	0.685
482.03	A2	0.362	0.362	0.362	0.362	0.362	0.362	0.362	0.362
691.77	E	0.074	0.074	0.074	0.074	0.074	0.074	0.074	0.074
769.87	A2	0.213	0.213	0.213	0.213	0.213	0.213	0.213	0.213
792.15	E	0.298	0.298	0.298	0.298	0.298	0.298	0.298	0.298
1070.57	E	3.441	3.442	3.441	3.440	3.440	3.441	3.441	3.441
1080.66	A2	1.804	1.804	1.804	1.803	1.803	1.804	1.804	1.804
1148.71	E	0.078	0.078	0.078	0.078	0.078	0.078	0.078	0.078

b. Effect of functional

To test the effect of using different exchange-correlation functionals to calculate the Born charges, we compared the intensities obtained from one of the LDA calculations (600 eV cutoff with a 6×6×6 k-point mesh) to those calculated using the PBE GGA and TPSS meta-GGA functionals with equivalent settings.

We also performed calculations using the HSE06 hybrid functional with a relaxed 4×4×4 k-point mesh and a reduced total-energy convergence criteria of 10^-6 eV (vs. the usual 10^-8 eV).

Sample VASP input files for these tests can be found in the VASP-Files directory.

Again, we observed very little quantitative difference in the calculated intensities of the IR-active modes. We note, of course, that different exchange-correlation functionals would very likely have a more significant impact on other parts of the calculation (e.g. the calculation of the phonon frequencies and/or linewidths).

ν [cm-1]	Ir. Rep.	IIR [e2 amu-1]
		LDA	PBE	TPSS*	HSE06*,**
127.37	E	< 10-3	< 10-3	< 10-3	< 10-3
255.08	E	0.021	0.021	0.021	0.023
341.54	A2	0.202	0.200	0.199	0.208
374.19	E	0.268	0.265	0.266	0.269
434.99	E	0.685	0.678	0.680	0.705
482.03	A2	0.362	0.359	0.360	0.367
691.77	E	0.074	0.072	0.073	0.075
769.87	A2	0.213	0.209	0.210	0.215
792.15	E	0.298	0.293	0.293	0.309
1070.57	E	3.441	3.453	3.488	3.418
1080.66	A2	1.804	1.811	1.829	1.790
1148.71	E	0.078	0.077	0.078	0.079

* These calculations were performed using the finite-field method (LCALCEPS = .TRUE. in VASP) with the default applied field of 10^-2 eV Å^-1
** As noted above, for the HSE06 hybrid calculations, the k-point mesh was reduced to 4×4×4 subdivisions for computational efficiency

2. Raman intensities

As for the IR intensities, we checked the convergence of the calculated Raman intensities with respect to the plane-wave cutoff and k-point sampling by performing tests with a subset of the parameters employed in the previous section (the same VASP input deck is used for calculating ε_static and the Born effective-charge tensors; sample input files are in the VASP-Files folder).

#	Cutoff [eV]	k-points
1	520	6×6×6
2	600	6×6×6
5	600	8×8×8
6	600	10×10×10
7	600	12×12×12

We again observed no significant variation in the calculated intensities, and therefore opted to use the base parameters (520 eV cutoff with 6×6×6 k-point sampling) for the "production" simulated spectrum.

We stress again, however, that the weak dependence of the calculated intensities on the technical parameters, in particular the k-point sampling, may be due to the insulating nature of α-SiO₂.

ν [cm-1]	Ir. Rep.	IRaman [Å4 amu-1]
		1	2	5	6	7
127.37	E	0.834	0.834	0.834	0.834	0.834
223.52	A1	6.931	6.928	6.928	6.928	6.928
255.08	E	0.404	0.403	0.403	0.403	0.403
337.02	A1	1.057	1.056	1.055	1.056	1.056
374.19	E	0.624	0.623	0.624	0.623	0.623
434.99	E	0.718	0.718	0.717	0.717	0.717
454.96	A1	35.085	35.090	35.096	35.093	35.099
691.77	E	1.359	1.359	1.358	1.359	1.359
792.15	E	1.554	1.556	1.556	1.556	1.556
1070.57	E	1.366	1.365	1.365	1.366	1.366
1085.88	A1	1.113	1.113	1.112	1.113	1.113
1148.71	E	2.811	2.810	2.810	2.810	2.810

3. Mode linewidths

Finally, we also checked the convergence of the Γ-point linewidths with respect to the isotropic q-point grid and integration method used to evaluate the three-phonon interaction strengths during the Phono3py post processing.

We tested three different integration methods, viz. the tetrahedron method (THM) and Gaussian smearing with σ = 0.1 and 0.01 THz. The tetrahedron method and Gaussian smearing with σ = 0.1 THz give broadly similar results, whereas for most of the modes the convergence with the Gaussian smearing using σ = 0.01 THz is erratic, particularly at sparser q-point grids.

In this particular set of calculations, q-point meshes with 56×56×56 subdivisions (1.55 × 10⁴ triplets) led to an anomalous spike in the calculated linewidths, indicating that explicit convergence testing ought to be considered a necessity.

In the present case, the tetrahedron method and a q-point sampling of 48×48×48 subdivisions (9.84 × 10³ triplets) appears to give well-converged results, so we employed these settings to obtain the linewidths used to generate the simulated spectra.

a. A₁ modes

b. A₂ modes

c. E modes

References

1. https://materialsproject.org/materials/mp-6930/