Calculate the free energy and related quantities from the phonon density of states

fqha.py

#!/usr/bin/env python
# Calculate the free energy and related quantities from the phonon density
# of states. Optionally, renormalize "imaginary" frequencies according to
# a simple scheme.
# Davide Ceresoli <dceresoli@gmail.com>

import numpy as np
import sys, argparse

def omega_to_Ry(hbar_omega):   # omega in cm^-1
    return hbar_omega * 9.1126705e-06

kB = 6.3336231e-06
def kbT_to_Ry(T):              # T in K
    return T * kB

# internal energy
def U(omega, T):
    omega_over_kbT = omega_to_Ry(omega+1e-6) / kbT_to_Ry(T)
    return omega_to_Ry(omega)/2.0 + omega_to_Ry(omega)/(np.exp(omega_over_kbT)-1.0)

# entropy
def S(omega, T):
    omega_over_kbT = omega_to_Ry(omega+1e-6) / kbT_to_Ry(T)
    return -kB*np.log(1- np.exp(-omega_over_kbT)) + (omega_to_Ry(omega)/T) / (np.exp(omega_over_kbT)-1.0)

# F = U - T*S
def F(omega, T):
    return U(omega,T) - T*S(omega,T)

# C_V
def C_V(omega, T):
    omega_over_kbT = omega_to_Ry(omega+1e-6) / kbT_to_Ry(T)
    return kB * omega_over_kbT**2 * np.exp(-omega_over_kbT) / (np.exp(-omega_over_kbT) - 1.0)**2


#======================================================================
# main program
#======================================================================
parser = argparse.ArgumentParser(description='Calculate free energy from phonon DOS')
parser.add_argument('-f', '--fildos', dest='fildos', action='store', required=True,
                    help='file containing the phonon DOS')
parser.add_argument('--Tmin', dest='Tmin', action='store', default=0,
                    help='initial temperature')
parser.add_argument('--Tmax', dest='Tmax', action='store', default=1000,
                    help='final temperature')
parser.add_argument('--deltaT', dest='deltaT', action='store', default=100,
                    help='temperature increment')
parser.add_argument('-u', '--units', dest='units', default='cminv', action='store',
                    help='frequency units: cminv or THz')
parser.add_argument('-r', '--renormalize', dest='renorm', default=False, action='store_true',
                    help='renormalize imaginary frequencies')
args = parser.parse_args()

# read Phonon DOS
phdos = np.loadtxt(args.fildos)

# omega and dos
omega = phdos[:,0]
dos = phdos[:,1]

# convert omega to cminv?
if args.units == 'THz':
    omega *= 33.35641
    dos /= 33.35641
elif args.units == 'cminv':
    pass
else:
    sys.stderr.write('unknown units: THz or cminv\n')
    sys.exit(1)
    
# keep original omega, just for the integration
omega0 = omega.copy()   

# renormalize
if np.any(omega < 0.0):
    if args.renorm:
        count = np.count_nonzero(omega<0.0)
        print('# WARNING: {0} imaginary freqs, renormalized by -sqrt(2)'.format(count))
        omega[omega<0.0] = -np.sqrt(2)*omega[omega<0.0]
    else:
        discard = len(omega[omega<0.0])
        print('# WARNING: discarding {0} imaginary freqs'.format(discard))
        omega = omega[discard:]
        omega0 = omega0[discard:]
        dos = dos[discard:]

# number of modes
n = np.trapz(dos, omega0)
print('# number of phonon modes: {0} ({1:.3f} atoms)'.format(n, n/3))

# zero point energy
zpe = 0.5 * omega_to_Ry(np.trapz(dos*omega, omega0))
print('# zero point energy:', zpe)

print('#    T (K)     F (Ry)  S (mRy/K) Cv (mRy/K)     E (Ry)')
T = float(args.Tmin)
while T <= float(args.Tmax):
    if abs(T) <= 1e-2:
       F_t = zpe
       U_t = zpe
       S_t = 0.0
       CV_t = 0.0
    else:
       F_t = np.trapz(F(omega,T)*dos, omega0)
       U_t = np.trapz(U(omega,T)*dos, omega0)
       S_t = np.trapz(S(omega,T)*dos, omega0) * 1000
       CV_t = np.trapz(C_V(omega,T)*dos, omega0) * 1000
       
    print('{0:10.4f} {1:10.4f} {2:10.4f} {3:10.4f} {4:10.4f}'.format(T, F_t, S_t, CV_t, U_t))

    T += float(args.deltaT)