Phonopy

Tutorial: Phonon Calculations using Phonopy and CRYSTAL23

1. Geometry Optimization (CRYSTAL23)

Before calculating phonons, you must start with a tightly optimized structure. Any residual forces in the unit cell will result in imaginary frequencies (negative phonons).

Create a standard CRYSTAL23 input file (e.g., crystal.d12) for your primitive or conventional unit cell.

Ensure you have high convergence criteria for the geometry optimization and SCF cycles (e.g., TOLDEE 8, TOLINTEG 8 8 8 8 16).

Run the optimization:

runcrystal crystal

Once completed, verify that the optimization converged. The output file should be changed to crystal.o and will be used as the structure input for Phonopy.

2. Generate displaced supercells

Now, we will use Phonopy to create supercells with small atomic displacements.

NOTE: Before running this step, make sure there is a TEMPLATE file in the directory containing your basis set and Hamiltonian parameters. An example of Hamiltonian parameters of the used in this example is shown in subsection 7. Configuration Files Reference. 

To generate the displacements in supercells, run:

phonopy -v --crystal --dim="8 8 1" -d

There is a formatting issue when Phonopy generates the .ext files for CRYSTAL: if there is a negative symbol, it sometimes omits the space, creating two badly joined columns.

To automatically fix this formatting error, run:

./parse_supercells.py

NOTE: Make sure parse_supercells.py is executable and in your working directory. This script will scan for all .ext files, find merged columns (like a digit followed immediately by a minus sign), and insert the missing space.

3. Calculate Forces (CRYSTAL23)

Next, you need to run a single-point energy and gradient calculation on every displaced supercell you just generated. Usually, having a first-guess wave function or charge density from a previous calculation helps convergence. 

We suggest to first run a single-point energy on the supercell.d12 ground structure, and then introduce the final wave function information to the further displacements. The keyword in CRYSTAL23 for this is GUESSP.

STEP 1: Run supercell.d12 without GRADCAL keyword in .d12 file. Once finished, the wave function information will be written in in file fort.9.

STEP 2: Run displacement calculations with GRADCAL keyword in .d12 file. The wave function information will be read from file fort.9. In your CRYSTAL23 slurm you will need to add this lines:

export GUESS_FILE=SnSe2-supercell.f9

cp $DIR/$JOB.d12 $scratch/$JOB/INPUT

cp $DIR/$JOB.ext $scratch/$JOB/fort.34

cp $DIR/$GUESS_FILE $scratch/$JOB/fort.20

After running step 2, a TEMPLATE file will be read. Then the following lines will be added when running phonopy:

TOLDEE

10

GRADCAL

END

Add your basis set and Hamiltonian parameters in the TEMPLATE file, erase GRADCAL keyword in .d12 file from the STEP 1. Add GUESSP keyword in .d12 file for further displacement calculations. Run CRYSTAL23 on all of them to generate the .o output files.

4. Create the FORCE_SETS File

Once CRYSTAL23 finishes calculating the forces for all supercells, we gather these forces into a format Phonopy can read.

Run the following command to parse the CRYSTAL outputs and create the FORCE_SETS file:

phonopy -v --crystal -f supercell-*o

5. Calculate FORCE CONSTANTS

With the forces extracted, you can now generate the force constants for your system.

phonopy -v --crystal --dim="8 8 1" --writefc

6. Calculate Phonon properties

Now that the force constants are written, you can calculate and plot various phonon properties using your configuration files.

Phonon Mesh / Thermodynamic Properties:

phonopy -p mesh.conf

Phonon Dispersion (Band Structure):

phonopy -p -s band.conf

Partial Density of States (PDOS):

phonopy -p -s pdos.conf

Phonon Animation:

To visualize the phonon modes, you can generate an animation file:

phonopy --readfc -v --dim="3 3 3" --anime="0.5 0.5 0.5"

Combined Band Structure and PDOS Plot:

If you want to plot the dispersion bands alongside the Density of States, use the phonopy-bandplot tool:

phonopy-bandplot --dos projected_dos.dat --fmin '-0.25' --fmax '9' -i '1, 2 3' --legacy -l -t '1T-SnSe_{2}' -o bands_pdos.pdf

Dynamical Matrix:

If you need to explicitly extract and look at the dynamical matrix, you can tell Phonopy to write it out for specific q-points. You do this using the --writedm and --qpoints flags.

For example, to get the dynamical matrix at the Gamma point (0, 0, 0), you would run a command like this:

phonopy --crystal --dim="8 8 1" --qpoints="0 0 0" --writedm

This will generate a file called qpoints.yaml, which contains the full dynamical matrix (both real and imaginary parts) for that specific q-point.

7. Configuration Files Reference

To run the property calculations in Step 6, Phonopy requires specific configuration files (.conf). These files define your system, the q-point grids, and the high-symmetry paths in the Brillouin zone.Below are the exact configuration files used for this $SnSe_{2}$ tutorial. Save these in your working directory.

A. Band Structure (band.conf)

This file defines the high-symmetry path for the phonon dispersion plot.

ATOM_NAME = Sn Se

DIM = 8 8 1

BAND = 0.0 0.0 0.0  0.5 0.0 0.0  0.333333 0.333333 0.0  0.0 0.0 0.0

BAND_LABELS = $\Gamma$ M K $\Gamma$

EIGENVECTORS=.TRUE.

BAND_POINTS = 202

DIM: Matches our supercell dimensions of 8x8x1.

BAND / BAND_LABELS: Traces the path Γ - M - K - Γ.

BAND_POINTS: Calculates 202 points along this path for a smooth curve.

EIGENVECTORS: Set to .TRUE. to save the vibrational modes for each point.

B. Phonon Mesh (mesh.conf)

This file defines a uniform q-point mesh across the Brillouin zone, typically used for thermodynamic properties or total Density of States.

ATOM_NAME = Sn Se

DIM = 8 8 1

MP = 91 91 1

MESH_SYMMETRY = .FALSE.

EIGENVECTORS = .FALSE.

MP: Specifies a highly dense 91x91x1 Monkhorst-Pack mesh suitable for 2D/layered materials.

MESH_SYMMETRY: Set to .FALSE. to calculate the full mesh without reducing it by symmetry.

C. Partial Density of States (pdos.conf)

This file defines a uniform q-point mesh across the Brillouin zone, typically used for thermodynamic properties or total Density of States.

ATOM_NAME = Sn Se

DIM = 8 8 1

MP = 30 30 1

PDOS = 1, 2 3

MP: Uses a slightly lighter 30x30x1 grid for the PDOS calculation.

PDOS: The syntax 1, 2 3 tells Phonopy to calculate the projected DOS for atom 1 (Sn) separately from the combined projected DOS of atoms 2 and 3 (Se).

D. TEMPLATE example

TEMPLATE file contains your basis set and Hamiltonian parameters.

[BASIS SETS]

99 0

END

DFT

SPIN

r2SCANh-D3

XXLGRID

END

TOLINTEG

9 9 9 11 38

SHRINK

3 6

SCFDIR

BIPOSIZE

110000000

EXCHSIZE

110000000

MAXCYCLE

800

FMIXING

30

DIIS

PPAN

8. Interpreting the Results

Once you run the phonopy-bandplot plotting script with combined band structure and PDOS, it will generate a combined PDF file (e.g., bands_pdos.pdf). The graph is divided in two halves:

The Phonon Dispersion (Left Side)

This plots the phonon frequencies against the momentum (q-points) along the high-symmetry path you defined (Γ - M - K - Γ).

Acoustic Modes: You will see three bands starting exactly at zero frequency at the Γ point. These correspond to the acoustic vibrations (sound waves) where the atoms move together in the same direction.

Optical Modes: The higher-frequency bands are the optical modes, where atoms in the unit cell move out of phase with one another.

Structural Stability: If you see any bands dipping below the zero-frequency line (these are called "imaginary frequencies"), it means your structure is dynamically unstable. If all bands are at or above zero, you have a stable local minimum.

The Partial Density of States - PDOS (Right Side)

This plot shows the number of available phonon states at each frequency, broken down by atomic contribution.

Because you separated the atoms in your pdos.conf file (1 for Sn, and 2 3 for Se), you will see two distinct curves.

Mass differences: Generally, heavier atoms (like Tin) will dominate the peaks at lower frequencies, while lighter atoms (like Selenium) will dominate the peaks at higher frequencies.