Commit 596a8718 authored by Samuel Poncé's avatar Samuel Poncé

Modifiction to pp.py script and phonon self-energy

- Removal of the xml conversion from the pp.py script
- Creation of a separate pp-xml.py script
- Change of sign in the retarded phonon self-energy to match with the
correct Eq. 145 from RMP 89, 015003 (2017)
parent f94baa58
This diff is collapsed.
......@@ -4,191 +4,11 @@
# 14/07/2015 - Creation of the script - Samuel Ponce
# 14/03/2018 - Automatically reads the number of q-points - Michael Waters
# 14/03/2018 - Detect if SOC is included in the calculation - Samuel Ponce
# 13/11/2018 - Write dyn files in xml format for SOC case - Shunhong Zhang (USTC)
#
import numpy as np
import os
from xml.dom import minidom
# Convert the dyn files to the xml form, for SOC case - Shunhong Zhang (USTC)
def dyn2xml(prefix):
ndyn=int(os.popen('head -2 {0}.dyn0|tail -1'.format(prefix)).read())
for idyn in range(1,ndyn+1):
print '{0}.dyn{1} to {0}.dyn_q{1}.xml'.format(prefix,idyn)
dynmat=dyn(prefix,idyn)
dynmat._write_xml()
def get_geom_info():
if os.path.isfile('ph.out')==False:
print 'cannot extract geometry info from ph.out'
return 1
else:
volm=float(os.popen('grep -a volume ph.out 2>/dev/null|tail -1').readline().split()[-2])
get_at=os.popen('grep -a -A 3 "crystal axes" ph.out 2>/dev/null|tail -3').readlines()
at=np.array([[float(item) for item in line.split()[3:6]] for line in get_at])
get_bg=os.popen('grep -a -A 3 "reciprocal axes" ph.out 2>/dev/null|tail -3').readlines()
bg=np.array([[float(item) for item in line.split()[3:6]] for line in get_bg])
return volm,at,bg
class dyn(object):
def __init__(self,prefix,idyn):
self._prefix=prefix
self._idyn=idyn
fil='{0}.dyn{1}'.format(prefix,idyn)
f=open(fil)
self._comment=f.readline()
f.readline()
line=f.readline().split()
self._ntype=int(line[0])
self._natom=int(line[1])
self._ibrav=int(line[2])
self._nspin=1
self._cell_dim=np.array([float(ii) for ii in line[3:]])
self._volm=0
self._at=np.zeros((3,3),float)
self._bg=np.zeros((3,3),float)
try: self._volm,self._at,self._bg = get_geom_info()
except: print 'warning: lattice info not found'
self._species=[];
self._mass=[]
for i in range(self._ntype):
line=f.readline().split()
self._species.append(line[1].strip("'"))
self._mass.append(float(line[-1])/911.4442) # normalize to atomic mass
self._atom_type=np.zeros(self._natom,int)
self._pos=np.zeros((self._natom,3),float)
for i in range(self._natom):
line=f.readline().split()
self._atom_type[i]=int(line[1])
for j in range(3): self._pos[i,j]=float(line[j+2])
self._nqpt=int(os.popen('grep -c "Dynamical Matrix" {0}'.format(fil)).read().split()[0])
self._qpt=[]
self._dynmat=np.zeros((self._nqpt,self._natom,self._natom,3,3,2),float)
f.readline()
for iqpt in range(self._nqpt):
f.readline();
f.readline()
line=f.readline().split()
self._qpt.append(np.array([float(item) for item in line[3:6]]))
f.readline()
for i in range(self._natom):
for j in range(self._natom):
f.readline()
data=np.fromfile(f,sep=' ',count=18,dtype=float).reshape(3,3,2)
self._dynmat[iqpt,i,j]=data
self._qpt=np.array(self._qpt)
for i in range(5): f.readline()
self._freq=np.zeros((self._natom*3,2),float)
self._disp=np.zeros((self._natom*3,self._natom,3,2),float)
for i in range(self._natom*3):
line=f.readline().split()
self._freq[i,0]=float(line[4])
self._freq[i,1]=float(line[7])
for j in range(self._natom):
line=f.readline().split()[1:-1]
data=np.array([float(item) for item in line]).reshape(3,2)
self._disp[i,j]=data
def _write_xml(self):
doc=minidom.Document()
root = doc.createElement('Root')
doc.appendChild(root)
geom_info=doc.createElement('GEOMETRY_INFO')
tags=('NUMBER_OF_TYPES','NUMBER_OF_ATOMS','BRAVAIS_LATTICE_INDEX','SPIN_COMPONENTS')
numbers=(self._ntype,self._natom,self._ibrav,self._nspin)
for i,(tag,num) in enumerate(zip(tags,numbers)):
inode=doc.createElement(tag)
inode.setAttribute('type','integer')
inode.setAttribute('size','1')
inode.text=num
inode.appendChild(doc.createTextNode(str(num)))
geom_info.appendChild(inode)
cell_dim=doc.createElement('CELL_DIMENSIONS')
cell_dim.setAttribute('type','real')
cell_dim.setAttribute('size','6')
for i in range(6):
cell_dim.appendChild(doc.createTextNode('{0:16.10f}'.format(self._cell_dim[i])))
geom_info.appendChild(cell_dim)
tags=['AT','BG']
for tag,lat in zip(tags,(self._at,self._bg)):
inode=doc.createElement(tag)
inode.setAttribute('type','real')
inode.setAttribute('size','9')
inode.setAttribute('columns','3')
for i in range(3):
text=' '.join(['{0:16.10f}'.format(item) for item in lat[i]])
inode.appendChild(doc.createTextNode(text))
geom_info.appendChild(inode)
volm=doc.createElement('UNIT_CELL_VOLUME_AU')
volm.setAttribute('type','real')
volm.setAttribute('size','1')
volm.appendChild(doc.createTextNode('{0:16.10f}'.format(self._volm)))
geom_info.appendChild(volm)
for itype in range(self._ntype):
nt=doc.createElement('TYPE_NAME.{0}'.format(itype+1))
nt.setAttribute('type','character')
nt.setAttribute('size','1')
nt.setAttribute('len','3')
nt.appendChild(doc.createTextNode('{0}'.format(self._species[itype])))
na=doc.createElement('MASS.{0}'.format(itype+1))
na.setAttribute('type','real')
na.setAttribute('size','1')
na.appendChild(doc.createTextNode('{0:16.10f}'.format(self._mass[itype])))
geom_info.appendChild(nt)
geom_info.appendChild(na)
for iat in range(self._natom):
at=doc.createElement('ATOM.{0}'.format(iat+1))
at.setAttribute('SPECIES','{0}'.format(self._species[self._atom_type[iat]-1]))
at.setAttribute('INDEX',str(iat+1))
pos=' '.join(['{0:16.10f}'.format(item) for item in self._pos[iat]])
at.setAttribute('TAU',pos)
geom_info.appendChild(at)
nqpt=doc.createElement('NUMBER_OF_Q')
nqpt.setAttribute('type','integer')
nqpt.setAttribute('size','1')
nqpt.appendChild(doc.createTextNode(str(self._nqpt)))
geom_info.appendChild(nqpt)
root.appendChild(geom_info)
for iqpt in range(self._nqpt):
dynmat=doc.createElement('DYNAMICAL_MAT_.{0}'.format(iqpt+1))
qpt=doc.createElement('Q_POINT')
qpt.setAttribute('type','real')
qpt.setAttribute('size','3')
qpt.setAttribute('columns','3')
tnode=doc.createTextNode(' '.join(['{0:16.10f}'.format(item) for item in self._qpt[iqpt]]))
qpt.appendChild(tnode)
dynmat.appendChild(qpt)
for iat in range(self._natom):
for jat in range(self._natom):
ph=doc.createElement('PHI.{0}.{1}'.format(iat+1,jat+1))
ph.setAttribute('type','complex')
ph.setAttribute('size','9')
ph.setAttribute('columns','3')
for i in range(3):
for j in range(3):
text='{0:16.10f} {1:16.10f}'.format(self._dynmat[iqpt,iat,jat,i,j,0],self._dynmat[iqpt,iat,jat,i,j,1])
ph.appendChild(doc.createTextNode(text))
dynmat.appendChild(ph)
root.appendChild(dynmat)
mode=doc.createElement('FREQUENCIES_THZ_CMM1')
for iomega in range(self._natom*3):
inode=doc.createElement('OMEGA.{0}'.format(iomega+1))
inode.setAttribute('type','real')
inode.setAttribute('size','2')
inode.setAttribute('columns','2')
inode.appendChild(doc.createTextNode('{0:16.10f} {1:16.10f}'.format(self._freq[iomega,0],self._freq[iomega,1])))
idisp=doc.createElement('DISPLACEMENT.{0}'.format(iomega+1))
idisp.setAttribute('tpye','complex')
idisp.setAttribute('size','3')
for iat in range(self._natom):
for j in range(3):
tnode=doc.createTextNode('{0:16.10f} {1:16.10f}'.format(self._disp[iomega,iat,j,0],self._disp[iomega,iat,j,1]))
idisp.appendChild(tnode)
mode.appendChild(inode)
mode.appendChild(idisp)
root.appendChild(mode)
fp = open('{0}.dyn_q{1}.xml'.format(self._prefix,self._idyn), 'w')
doc.writexml(fp, addindent=' ', newl='\n')
# Return the number of q-points in the IBZ
def get_nqpt(prefix):
fname = '_ph0/' +prefix+'.phsave/control_ph.xml'
......@@ -233,21 +53,6 @@ prefix = str(user_input)
# Test if SOC
SOC = hasSOC(prefix)
# If SOC detected, but dyn is not in XML and we want to convert it
if SOC=='true':
user_input = raw_input('Calculation with SOC detected. Do you want to convert dyn in XML format [y/n]?\n')
if str(user_input) == 'y':
dyn2xml(prefix)
os.system('mv {0}.dyn*.xml save'.format(prefix))
# If no SOC detected, do you want to convert into XML format
if SOC=='false':
user_input = raw_input('Calculation without SOC detected. Do you want to convert to xml anyway [y/n]?\n')
if str(user_input) == 'y':
SOC = 'true'
dyn2xml(prefix)
os.system('mv {0}.dyn*.xml save'.format(prefix))
# Test if seq. or parallel run
SEQ = isSEQ(prefix)
......
......@@ -65,8 +65,6 @@
zi_all(:,:), &!
esigmar_all(:,:,:), &!
esigmai_all(:,:,:), &!
gammar_all(:,:,:), &! Real part of the Phonon self-energy (freq. dependent for spectral function)
gammai_all(:,:,:), &! Imaginary part of the Phonon self-energy (freq. dependent for spectral function)
jdos(:), &!
spectra(:,:,:,:,:,:), &! dipole absorption spectra, polarizations, nomega, nsmear, dme/vme, absorption/emission
sumr(:,:,:,:), &! to apply the ASR correction to every xq
......
......@@ -1472,6 +1472,13 @@
!
ENDDO
CLOSE(lambda_phself)
!
! SP - 03/2019
! \Gamma = 1/\tau = phonon lifetime
! \Gamma = - 2 * Im \Pi^R where \Pi^R is the retarted phonon self-energy.
! Im \Pi^R = pi*k-point weight*[f(E_k+q) - f(E_k)]*delta[E_k+q - E_k - w_q]
! Since gamma_all = pi*k-point weight*[f(E_k) - f(E_k+q)]*delta[E_k+q - E_k - w_q] we have
! \Gamma = 2 * gamma_all
OPEN(unit=linewidth_phself,file='linewidth.phself')
WRITE(linewidth_phself, '(a)') '# Phonon frequency and phonon lifetime in meV '
WRITE(linewidth_phself,'(a)') '# Q-point Mode Phonon freq (meV) Phonon linewidth (meV)'
......@@ -1479,7 +1486,7 @@
!
DO imode=1, nmodes
WRITE(linewidth_phself,'(i9,i6,E20.8,E22.10)') iqq,imode,&
ryd2mev*wf(imode,iqq),ryd2mev*REAL(gamma_all(imode,iqq,1))
ryd2mev*wf(imode,iqq), 2.0d0 * ryd2mev * REAL(gamma_all(imode,iqq,1))
ENDDO
!
ENDDO
......
......@@ -329,13 +329,13 @@
wgkq = wgauss( -ekq*inv_eptemp0, -99)
!
! = k-point weight * [f(E_k) - f(E_k+q)]/ [E_k+q - E_k -w_q + id]
! This is the imaginary part of the phonon self-energy, sans
! This is the imaginary part of minus the phonon self-energy, sans
! the matrix elements
!
!weight = wkf (ikk) * (wgkk - wgkq) * &
! aimag ( cone / ( ekq - ekk - wq - ci * degaussw0 ) )
!
! SP: The expression below (phonon self-energy) corresponds to
! SP: The expression below (minus phonon self-energy) corresponds to
! = pi*k-point weight*[f(E_k) - f(E_k+q)]*delta[E_k+q - E_k - w_q]
weight = pi * wkf (ikk) * (wgkk - wgkq)* &
w0gauss ( (ekq - ekk - wq) / degaussw0, 0) / degaussw0
......
This diff is collapsed.
......@@ -45,6 +45,8 @@ dos1=`grep "DOS =" $fname | awk '{print $3}'`
e2=`grep " E(" $fname | awk '{print $4}'`
rsig=`grep "Re\[Sigma\]=" $fname | awk '{print $7}'`
isig=`grep "Im\[Sigma\]=" $fname | awk '{print $10}'`
rpi=`grep "Re\[Pi\]=" $fname | awk '{print $7}'`
ipi=`grep "Im\[Pi\]=" $fname | awk '{print $10}'`
z1=`grep " Z=" $fname | awk '{print $13}'`
lam=`grep "lam= " $fname | awk '{print $15}'`
lambda=`grep " lambda___(" $fname | awk '{print $4}'`
......@@ -166,6 +168,16 @@ if test "$isig" != ""; then
for x in $isig; do echo $x; done
fi
if test "$rpi" != ""; then
echo rpi
for x in $rpi; do echo $x; done
fi
if test "$ipi" != ""; then
echo ipi
for x in $ipi; do echo $x; done
fi
if test "$z1" != ""; then
echo z1
for x in $z1; do echo $x; done
......
......@@ -70,6 +70,8 @@ tolerance = ( (1.0e-6, 5.0e-3, 'e1'),
(1.0e-3, 5.0e-3, 'e2'),
(1.0 , 2.0e-1, 'rsig'), # epw_base3 on desktop 0.054364
(1.5 , 5.0e-1, 'isig'),
(1.0e-2, 1.0e-2, 'rpi'), # epw_base3 on desktop 0.054364
(1.0e-2, 1.0e-2, 'ipi'),
(5.0e-1, 1.0e-1, 'z1'), # epw_soc on desktop
(1.0e-2, None , 'lam'), # epw_base3 on desktop 7e-06
(1.0e-5, 1.0e-5, 'lambda'),
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment