phonopy

import numpy as np from phonopy import Phonopy from phonopy.structure.atoms import PhonopyAtoms # Define a simple unit cell (e.g., Silicon diamond structure) a = 5.404 # Lattice constant in Angstrom unitcell = PhonopyAtoms( symbols=['Si'] * 8, cell=np.diag([a, a, a]), scaled_positions=[ [0, 0, 0], [0, 0.5, 0.5], [0.5, 0, 0.5], [0.5, 0.5, 0], [0.25, 0.25, 0.25], [0.25, 0.75, 0.75], [0.75, 0.25, 0.75], [0.75, 0.75, 0.25] ] ) # Define the supercell matrix (e.g., 2x2x2 supercell) supercell_matrix = [[2, 0, 0], [0, 2, 0], [0, 0, 2]] # Initialize Phonopy object # The primitive matrix is often set based on the crystal structure, e.g., for FCC/diamond-like cells # For Si, the primitive cell is usually R-centered hexagonal or face-centered orthorhombic. # A common choice for diamond is: [[0, 0.5, 0.5], [0.5, 0, 0.5], [0.5, 0.5, 0]] primitive_matrix = [[0, 0.5, 0.5], [0.5, 0, 0.5], [0.5, 0.5, 0]] phonon = Phonopy( unitcell, supercell_matrix=supercell_matrix, primitive_matrix=primitive_matrix ) # Generate displacements for force calculations # distance defines the displacement magnitude (default 0.01 Angstrom for VASP-like calculators) phonon.generate_displacements(distance=0.03) # Get the supercells with displacements. These would typically be used to generate input files for DFT calculations. supercells_with_displacements = phonon.supercells_with_displacements print(f"Number of displaced supercells: {len(supercells_with_displacements)}") print("First displaced supercell (symbols and scaled positions):") print(supercells_with_displacements[0].symbols) print(supercells_with_displacements[0].scaled_positions) # In a real workflow, you would then run DFT calculations on these supercells to get forces, # and then use phonon.set_forces() and further methods to compute phonon properties.