Select q-point extension to random grids.

Also modification of the pp.py script to be able to convert dyn file from
text to XML format. Courtesy of Shunhong Zhang.
This commit is contained in:
Samuel Ponce 2018-11-17 16:45:36 +00:00
parent 50ed822831
commit 145fc74df1
14 changed files with 420 additions and 157 deletions

View File

@ -1,13 +1,194 @@
#!/usr/bin/python
#
# Post-processing script from of PH data in format used by EPW
# 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'
@ -52,6 +233,21 @@ 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)
@ -67,7 +263,7 @@ else:
except ValueError:
raise Exception('The value you enter is not an integer!')
os.system('mkdir save')
os.system('mkdir save 2>/dev/null')
for iqpt in np.arange(1,nqpt+1):
label = str(iqpt)
@ -118,3 +314,4 @@ for iqpt in np.arange(1,nqpt+1):
else:
os.system('cp _ph0/'+prefix+'.q_'+str(iqpt)+'/'+prefix+'.dvscf1 save/'+prefix+'.dvscf_q'+label)
os.system('rm _ph0/'+prefix+'.q_'+str(iqpt)+'/*wfc*' )

View File

@ -99,6 +99,8 @@
!! Check wheter this is the first cycle after a restart.
LOGICAL :: first_time
!! Check wheter this is the first timeafter a restart.
LOGICAL :: homogeneous
!! Check if the k and q grids are homogenous and commensurate.
!
CHARACTER (len=256) :: filint
!! Name of the file to write/read
@ -840,46 +842,45 @@
! Determines which q-points falls within the fsthick windows
! Store the result in the selecq.fmt file
! If the file exists, automatically restart from the file
! This is only done in the case of homogeneous grids.
! -----------------------------------------------------------------------
totq = 0
!
IF ( (nkf1 /= 0) .AND. (nkf2 /= 0) .AND. (nkf3 /= 0) .AND. (nqf1 /= 0) .AND. (nqf2 /= 0) .AND. (nqf3 /= 0) ) THEN
!
! Check if the file has been pre-computed
IF (mpime == ionode_id) THEN
INQUIRE(FILE='selecq.fmt',EXIST=exst)
! Check if the grids are homogeneous and commensurate
homogeneous = .FALSE.
IF ( (nkf1 /= 0) .AND. (nkf2 /= 0) .AND. (nkf3 /= 0) .AND. &
(nqf1 /= 0) .AND. (nqf2 /= 0) .AND. (nqf3 /= 0) .AND. &
(MOD(nkf1,nqf1) == 0) .AND. (MOD(nkf2,nqf2) == 0) .AND. (MOD(nkf3,nqf3) == 0) ) THEN
homogeneous = .TRUE.
ELSE
homogeneous = .FALSE.
ENDIF
!
totq = 0
! Check if the file has been pre-computed
IF (mpime == ionode_id) THEN
INQUIRE(FILE='selecq.fmt',EXIST=exst)
ENDIF
CALL mp_bcast(exst, ionode_id, world_comm)
!
IF (exst) THEN
IF (selecqread) THEN
WRITE(stdout,'(5x,a)')' '
WRITE(stdout,'(5x,a)')'Reading selecq.fmt file. '
CALL qwindow(exst, nrr_k, dims, totq, selecq, irvec_r, ndegen_k, cufkk, cufkq, homogeneous)
ELSE
WRITE(stdout,'(5x,a)')' '
WRITE(stdout,'(5x,a)')'A selecq.fmt file was found but re-created because selecqread == .false. '
CALL qwindow(.FALSE., nrr_k, dims, totq, selecq, irvec_r, ndegen_k, cufkk, cufkq, homogeneous)
ENDIF
CALL mp_bcast(exst, ionode_id, world_comm)
!
IF (exst) THEN
IF (selecqread) THEN
WRITE(stdout,'(5x,a)')' '
WRITE(stdout,'(5x,a)')'Reading selecq.fmt file. '
CALL qwindow(exst, nrr_k, dims, totq, selecq, irvec_r, ndegen_k, cufkk, cufkq)
ELSE
WRITE(stdout,'(5x,a)')' '
WRITE(stdout,'(5x,a)')'A selecq.fmt file was found but re-created because selecqread == .false. '
CALL qwindow(.FALSE., nrr_k, dims, totq, selecq, irvec_r, ndegen_k, cufkk, cufkq)
ENDIF
ELSE ! exst
IF (selecqread) THEN
CALL errore( 'ephwann_shuffle', 'Variable selecqread == .true. but file selecq.fmt not found.',1 )
ELSE
CALL qwindow(exst, nrr_k, dims, totq, selecq, irvec_r, ndegen_k, cufkk, cufkq)
ENDIF
ELSE ! exst
IF (selecqread) THEN
CALL errore( 'ephwann_shuffle', 'Variable selecqread == .true. but file selecq.fmt not found.',1 )
ELSE
CALL qwindow(exst, nrr_k, dims, totq, selecq, irvec_r, ndegen_k, cufkk, cufkq, homogeneous)
ENDIF
!
WRITE(stdout,'(5x,a,i8,a)')'We only need to compute ',totq, ' q-points'
WRITE(stdout,'(5x,a)')' '
ELSE
! If Random points or points read from files, then take all.
totq = nqf
ALLOCATE(selecq(totq))
DO iq = 1, totq
selecq(iq) = iq
ENDDO
ENDIF ! homogeneous grids
ENDIF
!
WRITE(stdout,'(5x,a,i8,a)')'We only need to compute ',totq, ' q-points'
WRITE(stdout,'(5x,a)')' '
!
! -----------------------------------------------------------------------
! Possible restart during step 1)

View File

@ -100,6 +100,8 @@
!! Check wheter this is the first cycle after a restart.
LOGICAL :: first_time
!! Check wheter this is the first timeafter a restart.
LOGICAL :: homogeneous
!! Check if the k and q grids are homogenous and commensurate.
!
CHARACTER (len=256) :: filint
!! Name of the file to write/read
@ -809,45 +811,46 @@
! Determines which q-points falls within the fsthick windows
! Store the result in the selecq.fmt file
! If the file exists, automatically restart from the file
! This is only done in the case of homogeneous grids.
! -----------------------------------------------------------------------
!
! Check if the grids are homogeneous and commensurate
homogeneous = .FALSE.
IF ( (nkf1 /= 0) .AND. (nkf2 /= 0) .AND. (nkf3 /= 0) .AND. &
(nqf1 /= 0) .AND. (nqf2 /= 0) .AND. (nqf3 /= 0) .AND. &
(MOD(nkf1,nqf1) == 0) .AND. (MOD(nkf2,nqf2) == 0) .AND. (MOD(nkf3,nqf3) == 0) ) THEN
homogeneous = .TRUE.
ELSE
homogeneous = .FALSE.
ENDIF
!
totq = 0
!
IF ( (nkf1 /= 0) .AND. (nkf2 /= 0) .AND. (nkf3 /= 0) .AND. (nqf1 /= 0) .AND. (nqf2 /= 0) .AND. (nqf3 /= 0) ) THEN
!
! Check if the file has been pre-computed
IF (mpime == ionode_id) THEN
INQUIRE(FILE='selecq.fmt',EXIST=exst)
! Check if the file has been pre-computed
IF (mpime == ionode_id) THEN
INQUIRE(FILE='selecq.fmt',EXIST=exst)
ENDIF
CALL mp_bcast(exst, ionode_id, world_comm)
!
IF (exst) THEN
IF (selecqread) THEN
WRITE(stdout,'(5x,a)')' '
WRITE(stdout,'(5x,a)')'Reading selecq.fmt file. '
CALL qwindow(exst, nrr_k, dims, totq, selecq, irvec_r, ndegen_k, cufkk, cufkq, homogeneous)
ELSE
WRITE(stdout,'(5x,a)')' '
WRITE(stdout,'(5x,a)')'A selecq.fmt file was found but re-created because selecqread == .false. '
CALL qwindow(.FALSE., nrr_k, dims, totq, selecq, irvec_r, ndegen_k, cufkk, cufkq, homogeneous)
ENDIF
CALL mp_bcast(exst, ionode_id, world_comm)
!
IF (exst) THEN
IF (selecqread) THEN
WRITE(stdout,'(5x,a)')' '
WRITE(stdout,'(5x,a)')'Reading selecq.fmt file. '
CALL qwindow(exst, nrr_k, dims, totq, selecq, irvec_r, ndegen_k, cufkk, cufkq)
ELSE
WRITE(stdout,'(5x,a)')' '
WRITE(stdout,'(5x,a)')'A selecq.fmt file was found but re-created because selecqread == .false. '
CALL qwindow(.FALSE., nrr_k, dims, totq, selecq, irvec_r, ndegen_k, cufkk, cufkq)
ENDIF
ELSE ! exst
IF (selecqread) THEN
CALL errore( 'ephwann_shuffle', 'Variable selecqread == .true. but file selecq.fmt not found.',1 )
ELSE
CALL qwindow(exst, nrr_k, dims, totq, selecq, irvec_r, ndegen_k, cufkk, cufkq)
ENDIF
ELSE ! exst
IF (selecqread) THEN
CALL errore( 'ephwann_shuffle', 'Variable selecqread == .true. but file selecq.fmt not found.',1 )
ELSE
CALL qwindow(exst, nrr_k, dims, totq, selecq, irvec_r, ndegen_k, cufkk, cufkq, homogeneous)
ENDIF
!
WRITE(stdout,'(5x,a,i8,a)')'We only need to compute ',totq, ' q-points'
WRITE(stdout,'(5x,a)')' '
ELSE
totq = nqf
ALLOCATE(selecq(totq))
DO iq = 1, totq
selecq(iq) = iq
ENDDO
ENDIF ! homogeneous grids
ENDIF
!
WRITE(stdout,'(5x,a,i8,a)')'We only need to compute ',totq, ' q-points'
WRITE(stdout,'(5x,a)')' '
!
! -----------------------------------------------------------------------
! Possible restart during step 1)

View File

@ -17,7 +17,8 @@
CONTAINS
!
!-----------------------------------------------------------------------
SUBROUTINE qwindow(exst, nrr_k, dims, totq, selecq, irvec_r, ndegen_k, cufkk, cufkq)
SUBROUTINE qwindow( exst, nrr_k, dims, totq, selecq, irvec_r, ndegen_k, &
cufkk, cufkq, homogeneous )
!-----------------------------------------------------------------------
!!
!! This subroutine pre-computes the q-points that falls within the fstichk
@ -44,6 +45,8 @@
!
LOGICAL, INTENT(in) :: exst
!! If the file exist
LOGICAL, INTENT(in) :: homogeneous
!! Check if the grids are homogeneous and commensurate
INTEGER, INTENT(IN) :: nrr_k
!! Number of WS points for electrons
INTEGER, INTENT(IN) :: dims
@ -105,6 +108,8 @@
!! $r\cdot k$
REAL(kind=DP) :: etf_loc(nbndsub, nkf)
!! Eigen-energies all full k-grid.
REAL(kind=DP) :: etf_locq(nbndsub, nkf)
!! Eigen-energies all full k-grid.
REAL(kind=DP) :: etf_all(nbndsub, nkqtotf/2)
!! Eigen-energies all full k-grid.
REAL(kind=DP) :: xkf_tmp (3, nkqtotf)
@ -138,104 +143,153 @@
CALL mp_bcast(selecq, ionode_id, world_comm )
IF (nqtot /= nqtotf) THEN
CALL errore( 'qwindow', 'Cannot read from selecq.fmt, the q-point grid or &
fsthick window are different from read one. Remove the selecq.fmt file and restart. ',1 )
& fsthick window are different from read one. Remove the selecq.fmt file and restart.',1 )
ENDIF
!
ELSE
ALLOCATE(selecq(nqf))
selecq(:) = 0
etf_loc(:,:) = zero
etf_loc(:,:) = zero
etf_locq(:,:) = zero
etf_all(:,:) = zero
!
! First store eigen energies on full grid.
DO ik = 1, nkf
ikk = 2 * ik - 1
xkk = xkf(:, ikk)
CALL dgemv('t', 3, nrr_k, twopi, irvec_r, 3, xkk, 1, 0.0_DP, rdotk, 1 )
IF (use_ws) THEN
DO iw=1, dims
DO iw2=1, dims
DO ir = 1, nrr_k
IF (ndegen_k(ir,iw2,iw) > 0) THEN
cfac(ir,iw2,iw) = exp( ci*rdotk(ir) ) / ndegen_k(ir,iw2,iw)
ENDIF
ENDDO
ENDDO
ENDDO
ELSE
cfac(:,1,1) = exp( ci*rdotk(:) ) / ndegen_k(:,1,1)
ENDIF
CALL hamwan2bloch ( nbndsub, nrr_k, cufkk, etf_loc(:, ik), chw, cfac, dims)
ENDDO
CALL poolgather ( nbndsub, nkqtotf/2, nkf, etf_loc, etf_all )
!
! In case of k-point symmetry
IF (mp_mesh_k) THEN
BZtoIBZ(:) = 0
s_BZtoIBZ(:,:,:) = 0
!
IF ( mpime == ionode_id ) THEN
!
CALL set_sym_bl( )
!
! What we get from this call is BZtoIBZ
CALL kpoint_grid_epw ( nrot, time_reversal, .false., s, t_rev, bg, nkf1*nkf2*nkf3, &
nkf1,nkf2,nkf3, nkqtotf_tmp, xkf_tmp, wkf_tmp,BZtoIBZ,s_BZtoIBZ)
!
IF (iterative_bte) THEN
BZtoIBZ_tmp(:) = 0
DO ikbz=1, nkf1*nkf2*nkf3
BZtoIBZ_tmp(ikbz) = map_rebal( BZtoIBZ( ikbz ) )
ENDDO
BZtoIBZ(:) = BZtoIBZ_tmp(:)
ENDIF
!
ENDIF ! mpime
CALL mp_bcast( BZtoIBZ, ionode_id, inter_pool_comm )
!
ENDIF ! mp_mesh_k
!
DO iq=1, nqf
xxq = xqf (:, iq)
!
found(:) = 0
DO ik = 1, nkf
IF (homogeneous) THEN
! First store eigen energies on full grid.
DO ik = 1, nkf
ikk = 2 * ik - 1
xkk = xkf(:, ikk)
xkq = xkk + xxq
!
CALL kpmq_map( xkk, (/0d0,0d0,0d0/), 1, ind1 )
CALL kpmq_map( xkk, xxq, 1, ind2 )
!
! Use k-point symmetry
IF (mp_mesh_k) THEN
IF ( (( minval ( abs(etf_all(:, BZtoIBZ(ind1)) - ef) ) < fsthick ) .and. &
( minval ( abs(etf_all(:, BZtoIBZ(ind2)) - ef) ) < fsthick )) ) THEN
found(my_pool_id+1) = 1
EXIT ! exit the loop
ENDIF
CALL dgemv('t', 3, nrr_k, twopi, irvec_r, 3, xkk, 1, 0.0_DP, rdotk, 1 )
IF (use_ws) THEN
DO iw=1, dims
DO iw2=1, dims
DO ir = 1, nrr_k
IF (ndegen_k(ir,iw2,iw) > 0) THEN
cfac(ir,iw2,iw) = exp( ci*rdotk(ir) ) / ndegen_k(ir,iw2,iw)
ENDIF
ENDDO
ENDDO
ENDDO
ELSE
IF ( (( minval ( abs(etf_all(:, ind1) - ef) ) < fsthick ) .and. &
( minval ( abs(etf_all(:, ind2) - ef) ) < fsthick )) ) THEN
cfac(:,1,1) = exp( ci*rdotk(:) ) / ndegen_k(:,1,1)
ENDIF
CALL hamwan2bloch ( nbndsub, nrr_k, cufkk, etf_loc(:, ik), chw, cfac, dims)
ENDDO
CALL poolgather ( nbndsub, nkqtotf/2, nkf, etf_loc, etf_all )
!
! In case of k-point symmetry
IF (mp_mesh_k) THEN
BZtoIBZ(:) = 0
s_BZtoIBZ(:,:,:) = 0
!
IF ( mpime == ionode_id ) THEN
!
CALL set_sym_bl( )
!
! What we get from this call is BZtoIBZ
CALL kpoint_grid_epw ( nrot, time_reversal, .false., s, t_rev, bg, nkf1*nkf2*nkf3, &
nkf1,nkf2,nkf3, nkqtotf_tmp, xkf_tmp, wkf_tmp,BZtoIBZ,s_BZtoIBZ)
!
IF (iterative_bte) THEN
BZtoIBZ_tmp(:) = 0
DO ikbz=1, nkf1*nkf2*nkf3
BZtoIBZ_tmp(ikbz) = map_rebal( BZtoIBZ( ikbz ) )
ENDDO
BZtoIBZ(:) = BZtoIBZ_tmp(:)
ENDIF
!
ENDIF ! mpime
CALL mp_bcast( BZtoIBZ, ionode_id, inter_pool_comm )
!
ENDIF ! mp_mesh_k
!
DO iq=1, nqf
xxq = xqf (:, iq)
!
found(:) = 0
DO ik = 1, nkf
ikk = 2 * ik - 1
xkk = xkf(:, ikk)
xkq = xkk + xxq
!
CALL kpmq_map( xkk, (/0d0,0d0,0d0/), 1, ind1 )
CALL kpmq_map( xkk, xxq, 1, ind2 )
!
! Use k-point symmetry
IF (mp_mesh_k) THEN
IF ( (( minval ( abs(etf_all(:, BZtoIBZ(ind1)) - ef) ) < fsthick ) .and. &
( minval ( abs(etf_all(:, BZtoIBZ(ind2)) - ef) ) < fsthick )) ) THEN
found(my_pool_id+1) = 1
EXIT ! exit the loop
ENDIF
ELSE
IF ( (( minval ( abs(etf_all(:, ind1) - ef) ) < fsthick ) .and. &
( minval ( abs(etf_all(:, ind2) - ef) ) < fsthick )) ) THEN
found(my_pool_id+1) = 1
EXIT ! exit the loop
ENDIF
ENDIF
!
ENDDO ! k-loop
! If found on any k-point from the pools
CALL mp_sum(found, world_comm)
!
IF (SUM(found) > 0) THEN
totq = totq + 1
selecq(totq) = iq
!
IF (MOD(totq,500) == 0) THEN
WRITE(stdout,'(5x,a,i8,i8)')'Number selected, total',totq,iq
ENDIF
ENDIF
ENDDO ! iq
ELSE ! homogeneous
DO iq=1, nqf
xxq = xqf (:, iq)
found(:) = 0
DO ik = 1, nkf
ikk = 2 * ik - 1
xkk = xkf(:, ikk)
xkq = xkk + xxq
!
CALL dgemv('t', 3, nrr_k, twopi, irvec_r, 3, xkk, 1, 0.0_DP, rdotk, 1 )
CALL dgemv('t', 3, nrr_k, twopi, irvec_r, 3, xkq, 1, 0.0_DP, rdotk2, 1 )
IF (use_ws) THEN
DO iw=1, dims
DO iw2=1, dims
DO ir = 1, nrr_k
IF (ndegen_k(ir,iw2,iw) > 0) THEN
cfac(ir,iw2,iw) = exp( ci*rdotk(ir) ) / ndegen_k(ir,iw2,iw)
cfacq(ir,iw2,iw) = exp( ci*rdotk2(ir) ) / ndegen_k(ir,iw2,iw)
ENDIF
ENDDO
ENDDO
ENDDO
ELSE
cfac(:,1,1) = exp( ci*rdotk(:) ) / ndegen_k(:,1,1)
cfacq(:,1,1) = exp( ci*rdotk2(:) ) / ndegen_k(:,1,1)
ENDIF
CALL hamwan2bloch ( nbndsub, nrr_k, cufkk, etf_loc(:, ik), chw, cfac, dims)
CALL hamwan2bloch ( nbndsub, nrr_k, cufkq, etf_locq(:, ik), chw, cfacq, dims)
IF ( (( minval ( abs(etf_loc(:, ik) - ef) ) < fsthick ) .and. &
( minval ( abs(etf_locq(:, ik) - ef) ) < fsthick )) ) THEN
found(my_pool_id+1) = 1
EXIT ! exit the loop
ENDIF
ENDIF
!
ENDDO ! k-loop
! If found on any k-point from the pools
CALL mp_sum(found, world_comm)
!
IF (SUM(found) > 0) THEN
totq = totq + 1
selecq(totq) = iq
ENDDO ! ik
! If found on any k-point from the pools
CALL mp_sum(found, world_comm)
!
IF (MOD(totq,500) == 0) THEN
WRITE(stdout,'(5x,a,i8,i8)')'Number selected, total',totq,iq
IF (SUM(found) > 0) THEN
totq = totq + 1
selecq(totq) = iq
!
IF (MOD(totq,500) == 0) THEN
WRITE(stdout,'(5x,a,i8,i8)')'Number selected, total',totq,iq
ENDIF
ENDIF
ENDIF
!
ENDDO
ENDDO ! iq
ENDIF ! homogeneous
!
IF (mpime == ionode_id) THEN
OPEN(unit=iunselecq, file='selecq.fmt', action='write')
WRITE (iunselecq,*) totq ! Selected number of q-points

View File

@ -1 +1,2 @@
diam
n

View File

@ -60,7 +60,7 @@
phonselfen = .false.
a2f = .false.
fsthick = 1.2 ! eV
fsthick = 2.0 ! eV
eptemp = 1 ! K
degaussw = 0.05 ! eV

View File

@ -60,7 +60,7 @@
phonselfen = .false.
a2f = .false.
fsthick = 1.2 ! eV
fsthick = 2.0 ! eV
eptemp = 1 ! K
degaussw = 0.05 ! eV

View File

@ -1 +1,2 @@
si
n

View File

@ -1 +1,2 @@
si
n

View File

@ -1 +1,2 @@
si
n

View File

@ -1 +1,2 @@
si
n

View File

@ -1 +1,2 @@
sic
n

View File

@ -1 +1,2 @@
MgB2
n

View File

@ -1 +1,2 @@
sic
n