1a_calc_vel_ops.py — Compute Velocity Operators

Computes the velocity-operator matrix elements at each irreducible q-point using phono3py/WTE and saves them to velocity_operators/save_{iq}.hdf5. These files are required by all subsequent steps.

The velocity operator is built from the dynamical matrix and its momentum-space derivatives evaluated on a Γ-centred BZ mesh. One HDF5 file is written per irreducible q-point.

Prerequisites: phono3py environment (source activate_phono3py.sh).

cd tutorials/diffusivity
python 1a_calc_vel_ops.py

  1#!/usr/bin/env python
  2
  3import sys, os
  4import time
  5import numpy as np
  6import h5py
  7from tqdm import tqdm
  8
  9import phonopy
 10from phonopy import Phonopy
 11from phonopy import file_IO
 12from phonopy.interface.calculator import read_crystal_structure
 13VaspToTHz = phonopy.physical_units.get_physical_units().DefaultToTHz
 14
 15from wte.velocity_operator import VelocityOperator
 16
 17from ase.io import read
 18from ase.units import J, _hplanck  # conversion factors
 19from ase.dft.kpoints import bandpath
 20
 21
 22global start_time
 23
 24# if len(sys.argv) > 1:
 25#     iq = int(sys.argv[1])
 26# else:
 27#     print('error: iq not given')
 28#     exit()
 29
 30withBORN = False
 31
 32
 33def write__to_hdf5(filename, frequencies, qpoints, weights, velocity_operators):
 34    compression = "gzip"
 35    with h5py.File(filename, "w") as w:
 36        w.create_dataset("frequency", data=frequencies, compression=compression)
 37        w.create_dataset("qpoint", data=qpoints)
 38        w.create_dataset("weight", data=weights)
 39        w.create_dataset("velocity_operator", data=velocity_operators, compression=compression)
 40
 41
 42
 43start_time = time.time()
 44eV2Hz = 1 / (J * _hplanck)
 45
 46# input
 47primitive_filename = "POSCAR"
 48supercell = np.diag((2, 2, 2))
 49mesh = [5, 5, 5]
 50interpolation_mesh_size = mesh[0]
 51
 52atoms, string = read_crystal_structure(filename=primitive_filename)
 53at_positions = atoms.scaled_positions
 54nat = len(at_positions)
 55nat3 = 3 * nat
 56
 57print(nat3)
 58
 59
 60# setup FC2
 61phonons = Phonopy(atoms, supercell)
 62
 63force_constants = file_IO.read_force_constants_hdf5(filename="fc2.hdf5")
 64primitive = phonons.primitive
 65
 66if withBORN:
 67    nac_params = file_IO.parse_BORN(primitive, filename="BORN")
 68    phonons.nac_params = nac_params
 69
 70phonons.force_constants = force_constants
 71phonons.symmetrize_force_constants()
 72
 73
 74
 75# diagonalize the dynamical matrix
 76phonons.run_mesh(mesh, is_gamma_center=True, with_group_velocities=False)
 77mesh_dict = phonons.get_mesh_dict()
 78qpoints = mesh_dict['qpoints']
 79weights = mesh_dict['weights']
 80frequencies = np.asarray(mesh_dict['frequencies'])
 81n_q_pt = len(qpoints)
 82time_absolute_1 = time.time()
 83partial_time = time.time() - start_time
 84
 85print('time 1=', partial_time)
 86
 87
 88
 89# calculation of vel op elements
 90os.makedirs('velocity_operators', exist_ok=True)
 91for iq in tqdm(range(n_q_pt), desc='q-points'):
 92    q_pt = qpoints[iq, :]
 93    vel_op = VelocityOperator(
 94        dynamical_matrix=phonons.dynamical_matrix,
 95        q_length=5e-6,
 96        symmetry=phonons.primitive_symmetry,
 97        frequency_factor_to_THz=VaspToTHz,
 98    )
 99    vel_op.run([q_pt])
100    write__to_hdf5('velocity_operators/save_%d.hdf5' % (iq), frequencies[iq, :], qpoints[iq, :], weights[iq], vel_op.velocity_operators[0, :, :, :])
101
102
103print('max_f [cm^-1]=', np.amax(frequencies) * 33.356)
104print('min_f [cm^-1]=', np.amin(frequencies) * 33.356)
105
106time_absolute_2 = time.time()
107total_time = time_absolute_2 - start_time
108
109print('diagonalization time=', time_absolute_2 - time_absolute_1)
110print('total time=', total_time)