From 145fc74df10e0351095354907fe302be9f2d6f54 Mon Sep 17 00:00:00 2001 From: Samuel Ponce Date: Sat, 17 Nov 2018 16:45:36 +0000 Subject: [PATCH] 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. --- EPW/bin/pp.py | 199 +++++++++++++++++++++++++- EPW/src/ephwann_shuffle.f90 | 73 +++++----- EPW/src/ephwann_shuffle_mem.f90 | 71 +++++----- EPW/src/transport.f90 | 222 +++++++++++++++++++----------- test-suite/epw_base/pp.in | 1 + test-suite/epw_mob/epw2.in | 2 +- test-suite/epw_mob/epw3.in | 2 +- test-suite/epw_mob/pp.in | 1 + test-suite/epw_mob_ibte/pp.in | 1 + test-suite/epw_mob_ibte_sym/pp.in | 1 + test-suite/epw_pl/pp.in | 1 + test-suite/epw_polar/pp.in | 1 + test-suite/epw_super/pp.in | 1 + test-suite/epw_trev/pp.in | 1 + 14 files changed, 420 insertions(+), 157 deletions(-) diff --git a/EPW/bin/pp.py b/EPW/bin/pp.py index 525f0b822..0864e4f4c 100644 --- a/EPW/bin/pp.py +++ b/EPW/bin/pp.py @@ -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*' ) + diff --git a/EPW/src/ephwann_shuffle.f90 b/EPW/src/ephwann_shuffle.f90 index 3389e739f..a8599109c 100644 --- a/EPW/src/ephwann_shuffle.f90 +++ b/EPW/src/ephwann_shuffle.f90 @@ -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) diff --git a/EPW/src/ephwann_shuffle_mem.f90 b/EPW/src/ephwann_shuffle_mem.f90 index 0ea56084b..963327edd 100644 --- a/EPW/src/ephwann_shuffle_mem.f90 +++ b/EPW/src/ephwann_shuffle_mem.f90 @@ -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) diff --git a/EPW/src/transport.f90 b/EPW/src/transport.f90 index 78c01e541..3ddf419dc 100644 --- a/EPW/src/transport.f90 +++ b/EPW/src/transport.f90 @@ -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 diff --git a/test-suite/epw_base/pp.in b/test-suite/epw_base/pp.in index 77f060859..c6caefcc3 100644 --- a/test-suite/epw_base/pp.in +++ b/test-suite/epw_base/pp.in @@ -1 +1,2 @@ diam +n diff --git a/test-suite/epw_mob/epw2.in b/test-suite/epw_mob/epw2.in index 2e1af8fd7..9e6226682 100644 --- a/test-suite/epw_mob/epw2.in +++ b/test-suite/epw_mob/epw2.in @@ -60,7 +60,7 @@ phonselfen = .false. a2f = .false. - fsthick = 1.2 ! eV + fsthick = 2.0 ! eV eptemp = 1 ! K degaussw = 0.05 ! eV diff --git a/test-suite/epw_mob/epw3.in b/test-suite/epw_mob/epw3.in index 1a1552d7c..b91713011 100644 --- a/test-suite/epw_mob/epw3.in +++ b/test-suite/epw_mob/epw3.in @@ -60,7 +60,7 @@ phonselfen = .false. a2f = .false. - fsthick = 1.2 ! eV + fsthick = 2.0 ! eV eptemp = 1 ! K degaussw = 0.05 ! eV diff --git a/test-suite/epw_mob/pp.in b/test-suite/epw_mob/pp.in index d58123204..384ad09f3 100644 --- a/test-suite/epw_mob/pp.in +++ b/test-suite/epw_mob/pp.in @@ -1 +1,2 @@ si +n diff --git a/test-suite/epw_mob_ibte/pp.in b/test-suite/epw_mob_ibte/pp.in index d58123204..384ad09f3 100644 --- a/test-suite/epw_mob_ibte/pp.in +++ b/test-suite/epw_mob_ibte/pp.in @@ -1 +1,2 @@ si +n diff --git a/test-suite/epw_mob_ibte_sym/pp.in b/test-suite/epw_mob_ibte_sym/pp.in index d58123204..384ad09f3 100644 --- a/test-suite/epw_mob_ibte_sym/pp.in +++ b/test-suite/epw_mob_ibte_sym/pp.in @@ -1 +1,2 @@ si +n diff --git a/test-suite/epw_pl/pp.in b/test-suite/epw_pl/pp.in index d58123204..384ad09f3 100644 --- a/test-suite/epw_pl/pp.in +++ b/test-suite/epw_pl/pp.in @@ -1 +1,2 @@ si +n diff --git a/test-suite/epw_polar/pp.in b/test-suite/epw_polar/pp.in index 12215da83..5f626e421 100644 --- a/test-suite/epw_polar/pp.in +++ b/test-suite/epw_polar/pp.in @@ -1 +1,2 @@ sic +n diff --git a/test-suite/epw_super/pp.in b/test-suite/epw_super/pp.in index cd6dc11ad..ec90b4406 100644 --- a/test-suite/epw_super/pp.in +++ b/test-suite/epw_super/pp.in @@ -1 +1,2 @@ MgB2 +n diff --git a/test-suite/epw_trev/pp.in b/test-suite/epw_trev/pp.in index 12215da83..5f626e421 100644 --- a/test-suite/epw_trev/pp.in +++ b/test-suite/epw_trev/pp.in @@ -1 +1,2 @@ sic +n