Commit 9123f441 authored by giannozz's avatar giannozz

Merge branch 'develop' into 'develop'

Select q-point extension to random grids.

See merge request !211
parents 50ed8228 145fc74d
#!/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*' )
......@@ -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.
! -----------------------------------------------------------------------
!
! 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 ( (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)
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)
......
......@@ -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)
......
......@@ -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
IF (homogeneous) THEN
! 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
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
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 )
!
IF ( mpime == ionode_id ) THEN
! In case of k-point symmetry
IF (mp_mesh_k) THEN
BZtoIBZ(:) = 0
s_BZtoIBZ(:,:,:) = 0
!
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 ( 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 )
!
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 ! mp_mesh_k
!
DO iq=1, nqf
xxq = xqf (:, iq)
!
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 )
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)
!
! 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
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
ELSE
IF ( (( minval ( abs(etf_all(:, ind1) - ef) ) < fsthick ) .and. &
( minval ( abs(etf_all(:, ind2) - ef) ) < fsthick )) ) THEN
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