scnc
/
quantum-espresso
mirror of https://gitlab.com/QEF/q-e.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Improved robustness of EPW (S. Ponce and C. Verdi) 

Explaintion:
The Wigner-Seitz construction in EPW was done by constructing a cell
centred unit cell. This is fine for electronic properties (this is what is done in wannier90).
However for phonon or electron-phonon properties, one can have issues when the cell
is tilded for example.
The proper way is to construct a set of WS vectors centred on pairs of atoms (phonons)
or atoms (el-ph).
In the matdyn code, a FT grid is constructed with weigths centred on pairs of atoms
and zeros everywhere else.
EPW now reproduced exactly the results of matdyn for the interpolated phonons at a
lower computation cost. Indeed we minimize the number of zeros by keeping the union
of values between all the cells.
In both cases this is very fast anyway but is important for el-ph properties.

A Wigner module was created to deal with this.

In addition, the subroutine cdiagh2 from PHonon/PH/rigid.f90 has been made public to
avoid code duplication in EPW
This commit is contained in:
Samuel Ponce 2018-07-01 12:12:38 +01:00
parent d0466b4b77
commit 73371a59c5
11 changed files with 1133 additions and 598 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
EPW/src/Makefile
View File

@ -29,6 +29,7 @@ superconductivity.o \
transportcom.o \
transport.o \
transport_iter.o \
wigner.o \
a2f.o \
allocate_epwq.o \
bcast_epw_input.o \
@ -93,8 +94,6 @@ wannierEPW.o \
wannierize.o \
pw2wan90epw.o \
wan2bloch.o \
wigner_seitz2.o \
wigner_seitz.o \
io_scattering.o \
system_mem_usage.o

116
EPW/src/bloch2wan.f90
View File

@ -1,5 +1,5 @@
!
! Copyright (C) 2010-2016 Samuel Ponce', Roxana Margine, Carla Verdi, Feliciano Giustino
!k Copyright (C) 2010-2016 Samuel Ponce', Roxana Margine, Carla Verdi, Feliciano Giustino
!
! This file is distributed under the terms of the GNU General Public
! License. See the file `LICENSE' in the root directory of the
@ -1187,8 +1187,8 @@
END SUBROUTINE ephbloch2wane
!
!--------------------------------------------------------------------------
SUBROUTINE ephbloch2wanp ( nbnd, nmodes, xk, nq, irvec, &
nrk, nrr, epmatwe )
SUBROUTINE ephbloch2wanp ( nbnd, nmodes, xk, nq, irvec_k, irvec_g, &
nrr_k, nrr_g, epmatwe )
!--------------------------------------------------------------------------
!!
!! From the EP Matrix in Electron Bloch representation (coarse mesh),
@ -1204,28 +1204,30 @@
USE io_global, ONLY : ionode_id
USE mp, ONLY : mp_barrier
USE mp_world, ONLY : mpime
!
implicit none
!
! Input variables - note irvec is dimensioned with nrr_k
! (which is assumed to be larger than nrr_q)
! Input variables
!
INTEGER, INTENT(in) :: nbnd
!! Number of electronic bands
INTEGER, INTENT(in) :: nrk
!! number of electronic WS points
INTEGER, INTENT(in) :: nmodes
!! number of branches
INTEGER, INTENT(in) :: nq
!! number of qpoints
INTEGER, INTENT(in) :: nrr
!! number of phononic WS points
INTEGER, INTENT(in) :: irvec(3, nrk)
!! Coordinates of real space vector (irvec is dimensioned with nrk)
INTEGER, INTENT(in) :: nrr_k
!! number of electronic WS points
INTEGER, INTENT(in) :: nrr_g
!! number of el-ph WS points
INTEGER, INTENT(in) :: irvec_k(3, nrr_k)
!! Coordinates of real space vector for electrons
INTEGER, INTENT(in) :: irvec_g(3, nrr_g)
!! Coordinates of real space vector for electron-phonon
!
REAL(kind=DP), INTENT(in) :: xk(3, nq)
!! Kpoint coordinates (cartesian in units of 2piba)
!
COMPLEX(kind=DP), INTENT(in) :: epmatwe(nbnd, nbnd, nrk, nmodes, nq)
COMPLEX(kind=DP), INTENT(in) :: epmatwe(nbnd, nbnd, nrr_k, nmodes, nq)
!! EP matrix in electron-wannier representation and phonon bloch representation
!! (Cartesian coordinates)
!
@ -1233,8 +1235,8 @@
!
! Work variables
!
INTEGER :: ik
!! Counter on k-point
INTEGER :: iq
!! Counter on q-point
INTEGER :: ir
!! Counter on WS points
INTEGER :: ire
@ -1263,14 +1265,14 @@
!
epmatwp = czero
!
DO ir = 1, nrr
DO ir = 1, nrr_g
!
DO ik = 1, nq
!
rdotk = twopi * dot_product( xk( :, ik), dble(irvec( :, ir) ))
cfac = exp( -ci*rdotk ) / dble(nq)
epmatwp(:,:,:,:,ir) = epmatwp(:,:,:,:,ir) + cfac * epmatwe(:,:,:,:,ik)
!
DO iq = 1, nq
!
rdotk = twopi * dot_product( xk( :, iq), dble(irvec_g( :, ir) ))
cfac = exp( -ci*rdotk ) / dble(nq)
epmatwp(:,:,:,:,ir) = epmatwp(:,:,:,:,ir) + cfac * epmatwe(:,:,:,:,iq)
!
ENDDO
!
! check spatial decay of e-p matrix elements in wannier basis - electrons
@ -1281,14 +1283,14 @@
IF (mpime.eq.ionode_id) THEN
IF (ir.eq.1) open(unit=iuwanep,file='decay.epmat_wanep',status='unknown')
IF (ir.eq.1) WRITE(iuwanep, '(a)') '# R_e, R_p, max_{m,n,nu} |g(m,n,nu;R_e,R_p)| '
DO ire = 1, nrk
DO ire = 1, nrr_k
!
rvec1 = dble(irvec(1,ire))*at(:,1) + &
dble(irvec(2,ire))*at(:,2) + &
dble(irvec(3,ire))*at(:,3)
rvec2 = dble(irvec(1,ir))*at(:,1) + &
dble(irvec(2,ir))*at(:,2) + &
dble(irvec(3,ir))*at(:,3)
rvec1 = dble(irvec_k(1,ire))*at(:,1) + &
dble(irvec_k(2,ire))*at(:,2) + &
dble(irvec_k(3,ire))*at(:,3)
rvec2 = dble(irvec_g(1,ir))*at(:,1) + &
dble(irvec_g(2,ir))*at(:,2) + &
dble(irvec_g(3,ir))*at(:,3)
len1 = sqrt(rvec1(1)**2.d0+rvec1(2)**2.d0+rvec1(3)**2.d0)
len2 = sqrt(rvec2(1)**2.d0+rvec2(2)**2.d0+rvec2(3)**2.d0)
tmp = maxval ( abs( epmatwp (:, :, ire, :, ir) ) )
@ -1299,7 +1301,7 @@
WRITE(iuwanep, '(5f15.10)') len1 * celldm(1) * bohr2ang, &
len2 * celldm(1) * bohr2ang, tmp
ENDDO
IF (ir.eq.nrr) CLOSE(iuwanep)
IF (ir.eq.nrr_g) CLOSE(iuwanep)
ENDIF
!
ENDDO
@ -1312,8 +1314,8 @@
!
! -----------------------------------------------------------
!--------------------------------------------------------------------------
SUBROUTINE ephbloch2wanp_mem ( nbnd, nmodes, xk, nq, irvec, &
nrk, nrr, epmatwe )
SUBROUTINE ephbloch2wanp_mem ( nbnd, nmodes, xk, nq, irvec_k, irvec_g, &
nrr_k, nrr_g, epmatwe )
!--------------------------------------------------------------------------
!
! From the EP Matrix in Electron Bloch representation (coarse mesh),
@ -1335,21 +1337,23 @@
!
INTEGER, INTENT(in) :: nbnd
!! Number of electronic bands
INTEGER, INTENT(in) :: nrk
INTEGER, INTENT(in) :: nrr_k
!! number of electronic WS points
INTEGER, INTENT(in) :: nrr_g
!! number of el-h WS points
INTEGER, INTENT(in) :: nmodes
!! number of branches
INTEGER, INTENT(in) :: nq
!! number of qpoints
INTEGER, INTENT(in) :: nrr
!! number of phononic WS points
INTEGER, INTENT(in) :: irvec (3, nrk)
!! Coordinates of real space vector (irvec is dimensioned with nrk)
INTEGER, INTENT(in) :: irvec_k (3, nrr_k)
!! Coordinates of real space vector
INTEGER, INTENT(in) :: irvec_g (3, nrr_g)
!! Coordinates of real space vector
!
REAL(kind=DP), INTENT(in) :: xk(3, nq)
!! Kpoint coordinates (cartesian in units of 2piba)
!
COMPLEX(kind=DP), INTENT(in) :: epmatwe (nbnd, nbnd, nrk, nmodes)
COMPLEX(kind=DP), INTENT(in) :: epmatwe (nbnd, nbnd, nrr_k, nmodes)
!! EP matrix in electron-wannier representation and phonon bloch representation
!! (Cartesian coordinates)
!
@ -1357,8 +1361,8 @@
!
! work variables
!
INTEGER :: ik
!! Counter on k-point
INTEGER :: iq
!! Counter on q-point
INTEGER :: ir
!! Counter on WS points
INTEGER :: ire
@ -1377,7 +1381,7 @@
COMPLEX(KIND=DP), ALLOCATABLE :: epmatwp_mem(:,:,:,:)
!! e-p matrix in Wannier basis
!
ALLOCATE (epmatwp_mem( nbnd, nbnd, nrk, nmodes))
ALLOCATE (epmatwp_mem( nbnd, nbnd, nrr_k, nmodes))
!
!----------------------------------------------------------
! Fourier transform to go into Wannier basis
@ -1390,39 +1394,39 @@
!
CALL cryst_to_cart (nq, xk, at, -1)
!
DO ir = 1, nrr
DO ir = 1, nrr_g
!
epmatwp_mem = czero
!
DO ik = 1, nq
DO iq = 1, nq
!
! direct read of epmatwe for this iq
CALL rwepmatw ( epmatwe, nbnd, nrk, nmodes, ik, iunepmatwe, -1)
CALL rwepmatw ( epmatwe, nbnd, nrr_k, nmodes, iq, iunepmatwe, -1)
!
rdotk = twopi * dot_product( xk( :, ik), dble(irvec( :, ir) ))
rdotk = twopi * dot_product( xk( :, iq), dble(irvec_g( :, ir) ))
cfac = exp( -ci*rdotk ) / dble(nq)
epmatwp_mem = epmatwp_mem + cfac * epmatwe
!
ENDDO
!
! direct write of epmatwp_mem for this ir
CALL rwepmatw (epmatwp_mem, nbnd, nrk, nmodes, ir, iunepmatwp, +1)
CALL rwepmatw (epmatwp_mem, nbnd, nrr_k, nmodes, ir, iunepmatwp, +1)
! check spatial decay of e-p matrix elements in wannier basis - electrons
! + phonons
!
! we plot: R_e, R_p, max_{m,n,nu} |g(m,n,nu;R_e,R_p)|
!
IF (mpime.eq.ionode_id) THEN
IF (ir.eq.1) open(unit=iuwanep,file='decay.epmat_wanep',status='unknown')
IF (ir.eq.1) WRITE(iuwanep, '(a)') '# R_e, R_p, max_{m,n,nu} |g(m,n,nu;R_e,R_p)| '
DO ire = 1, nrk
IF (mpime == ionode_id) THEN
IF (ir == 1) OPEN(unit=iuwanep, file='decay.epmat_wanep', status='unknown')
IF (ir == 1) WRITE(iuwanep, '(a)') '# R_e, R_p, max_{m,n,nu} |g(m,n,nu;R_e,R_p)| '
DO ire = 1, nrr_k
!
rvec1 = dble(irvec(1,ire))*at(:,1) + &
dble(irvec(2,ire))*at(:,2) + &
dble(irvec(3,ire))*at(:,3)
rvec2 = dble(irvec(1,ir))*at(:,1) + &
dble(irvec(2,ir))*at(:,2) + &
dble(irvec(3,ir))*at(:,3)
rvec1 = dble(irvec_k(1,ire))*at(:,1) + &
dble(irvec_k(2,ire))*at(:,2) + &
dble(irvec_k(3,ire))*at(:,3)
rvec2 = dble(irvec_g(1,ir))*at(:,1) + &
dble(irvec_g(2,ir))*at(:,2) + &
dble(irvec_g(3,ir))*at(:,3)
len1 = sqrt(rvec1(1)**2.d0+rvec1(2)**2.d0+rvec1(3)**2.d0)
len2 = sqrt(rvec2(1)**2.d0+rvec2(2)**2.d0+rvec2(3)**2.d0)
tmp = maxval ( abs( epmatwp_mem(:, :, ire, :) ) )
@ -1433,7 +1437,7 @@
WRITE(iuwanep, '(5f15.10)') len1 * celldm(1) * bohr2ang, &
len2 * celldm(1) * bohr2ang, tmp
ENDDO
IF (ir.eq.nrr) CLOSE(iuwanep)
IF (ir == nrr_g) CLOSE(iuwanep)
ENDIF
!
ENDDO

4
EPW/src/deallocate_epw.f90
View File

@ -32,7 +32,7 @@
USE control_lr, ONLY : nbnd_occ
USE becmod, ONLY : becp, deallocate_bec_type
USE elph2, ONLY : el_ph_mat, epf17, epsi, etf,&
etq, et_all, wf, wkf, wqf, wslen,&
etq, et_all, wf, wkf, wqf, &
xkq, xk_all, zstar, xkf, xqf, epmatwp, eps_rpa
USE epwcom, ONLY : epbread, epwread
USE modes, ONLY : npert, u, name_rap_mode, num_rap_mode
@ -58,7 +58,6 @@
IF(ALLOCATED(wqf)) DEALLOCATE (wqf)
IF(ALLOCATED(xk_all)) DEALLOCATE (xk_all)
IF(ALLOCATED(et_all)) DEALLOCATE (et_all)
IF(ALLOCATED(wslen)) DEALLOCATE (wslen)
IF(ALLOCATED(eps_rpa)) DEALLOCATE (eps_rpa)
IF(ALLOCATED(eps_rpa)) DEALLOCATE (eps_rpa)
!
@ -134,7 +133,6 @@
IF(ALLOCATED(wqf)) DEALLOCATE (wqf)
IF(ALLOCATED(xk_all)) DEALLOCATE (xk_all)
IF(ALLOCATED(et_all)) DEALLOCATE (et_all)
IF(ALLOCATED(wslen)) DEALLOCATE (wslen)
IF(ALLOCATED(eps_rpa)) DEALLOCATE (eps_rpa)
ENDIF ! epwread .and. .not. epbread
!

7
EPW/src/elph2.f90
View File

@ -51,7 +51,6 @@
etf_k(:,:), &! Saved interpolated KS eigenenergies for later used in q-parallelization (nbnd, nkqf)
etf_ks(:,:), &! interpolated eigenvalues (nbnd, nkqf) KS eigenvalues in the case of eig_read
wf(:,:), &! interpolated eigenfrequencies
wslen(:), &! length of the wigner seitz points in units of alat
gamma_all(:,:,:), &!
gamma_nest(:,:), &! Nesting function in the case of q-parallelization
gamma_v_all(:,:,:), &!
@ -92,17 +91,11 @@
nkqtotf, &! total number of k+q points (fine grid)
nqtotf, &! total number of q points (fine grid)
nrr, &! number of wigner-seitz points (elec interp only)
nrr_k, &! number of wigner-seitz points for electrons
nrr_q, &! number of wigner-seitz points for phonons
ibndmin, &! band bounds for slimming down electron-phonon matrix
ibndmax, &!
lower_band, &! Lower band index for image (band) parallelization
upper_band ! Upper band index for image (band) parallelization
INTEGER, ALLOCATABLE :: &
irvec(:,:), &! crys coordinates of wigner-seitz vectors (both elec and phon)
ndegen(:), &! corresponding degeneragy, electrons (old version)
ndegen_k(:), &! corresponding degeneragy, electrons
ndegen_q(:), &! corresponding degeneragy, phonons
igk(:), &! Index for k+G vector
igkq(:), &! Index for k+q+G vector
igk_k_all(:,:), &! Global index (in case of parallel)

214
EPW/src/ephwann_shuffle.f90
View File

@ -48,8 +48,8 @@
USE io_global, ONLY : stdout, ionode
USE io_epw, ONLY : lambda_phself, linewidth_phself, iunepmatwe, &
iunepmatwp, crystal, iunepmatwp2
USE elph2, ONLY : nrr_k, nrr_q, cu, cuq, lwin, lwinq, irvec, ndegen_k,&
ndegen_q, wslen, chw, chw_ks, cvmew, cdmew, rdw, &
USE elph2, ONLY : cu, cuq, lwin, lwinq, &
chw, chw_ks, cvmew, cdmew, rdw, &
epmatwp, epmatq, wf, etf, etf_k, etf_ks, xqf, xkf, &
wkf, dynq, nqtotf, nkqf, epf17, nkf, nqf, et_ks, &
ibndmin, ibndmax, lambda_all, dmec, dmef, vmef, &
@ -66,6 +66,7 @@
USE bloch2wan, ONLY : hambloch2wan, dmebloch2wan, dynbloch2wan, &
vmebloch2wan, ephbloch2wane, ephbloch2wanp, &
ephbloch2wanp_mem
USE wigner, ONLY : wigner_seitz_wrap
USE superconductivity, ONLY : write_ephmat, count_kpoints, kmesh_fine, &
kqmap_fine
USE transport, ONLY : transport_coeffs, scattering_rate_q
@ -104,14 +105,6 @@
CHARACTER (len=30) :: myfmt
!! Variable used for formatting output
!
INTEGER, ALLOCATABLE :: ndegen_kk(:)
!! Number of degeneresence with nrr_k bounds.
INTEGER, ALLOCATABLE :: ndegen_qq(:)
!! Number of degeneresence with nrr_q bounds.
INTEGER, ALLOCATABLE :: irvec_kk(:,:)
!! Irevec on the nrr_k bounds
INTEGER, ALLOCATABLE :: irvec_qq(:,:)
!! Irevec on the nrr_k bounds
INTEGER :: ios
!! integer variable for I/O control
INTEGER :: iq
@ -158,6 +151,27 @@
!! Return virtual and resisdent memory from system
INTEGER :: ierr
!! Error status
INTEGER :: nrr_k
!! Number of WS points for electrons
INTEGER :: nrr_q
!! Number of WS points for phonons
INTEGER :: nrr_g
!! Number of WS points for electron-phonons
INTEGER, ALLOCATABLE :: irvec_k(:,:)
!! integer components of the ir-th Wigner-Seitz grid point in the basis
!! of the lattice vectors for electrons
INTEGER, ALLOCATABLE :: irvec_q(:,:)
!! integer components of the ir-th Wigner-Seitz grid point for phonons
INTEGER, ALLOCATABLE :: irvec_g(:,:)
!! integer components of the ir-th Wigner-Seitz grid point for electron-phonon
INTEGER, ALLOCATABLE :: ndegen_k (:)
!! Wigner-Seitz number of degenerescence (weights) for the electrons grid
INTEGER, ALLOCATABLE :: ndegen_q (:,:,:)
!! Wigner-Seitz weights for the phonon grid that depend on
!! atomic positions $R + \tau(nb) - \tau(na)$
INTEGER, ALLOCATABLE :: ndegen_g (:,:)
!! Wigner-Seitz weights for the electron-phonon grid that depend on
!! atomic positions $R - \tau(na)$
INTEGER, PARAMETER :: nrwsx=200
!! Maximum number of real-space Wigner-Seitz
!
@ -187,6 +201,12 @@
!! Wigner-Size supercell vectors, store in real instead of integer
REAL(kind=DP), ALLOCATABLE :: rdotk(:)
!! $r\cdot k$
REAL(kind=DP), ALLOCATABLE :: wslen_k(:)
!! real-space length for electrons, in units of alat
REAL(kind=DP), ALLOCATABLE :: wslen_q(:)
!! real-space length for phonons, in units of alat
REAL(kind=DP), ALLOCATABLE :: wslen_g(:)
!! real-space length for electron-phonons, in units of alat
!
COMPLEX(kind=DP) :: tablex (4*nk1+1,nkf1)
!! Look-up table for the exponential (speed optimization) in the case of
@ -259,16 +279,10 @@
cuq ( nbnd, nbndsub, nks), &
lwin ( nbnd, nks ), &
lwinq ( nbnd, nks ), &
exband ( nbnd ), &
irvec (3, 20*nk1*nk2*nk3), &
ndegen_k (20*nk1*nk2*nk3), &
ndegen_q (20*nq1*nq2*nq3), &
wslen(20*nk1*nk2*nk3) )
exband ( nbnd ) )
!
CALL start_clock ( 'ephwann' )
!
! DBSP
! HERE loadkmesh
IF ( epwread ) THEN
!
! Might have been pre-allocate depending of the restart configuration
@ -333,27 +347,17 @@
!
ALLOCATE( w2( 3*nat) )
!
! determine Wigner-Seitz points
! Determine Wigner-Seitz points
!
CALL wigner_seitz2 &
( nk1, nk2, nk3, nq1, nq2, nq3, nrr_k, nrr_q, irvec, wslen, ndegen_k, ndegen_q )
! Inside we allocate irvec_k, irvec_q, irvec_g, ndegen_k, ndegen_q, ndegen_g,
! wslen_k, wslen_q, wslen_g
CALL wigner_seitz_wrap ( nk1, nk2, nk3, nq1, nq2, nq3, irvec_k, irvec_q, irvec_g, &
ndegen_k, ndegen_q, ndegen_g, wslen_k, wslen_q, wslen_g )
!
ALLOCATE ( ndegen_kk(nrr_k) )
ALLOCATE ( ndegen_qq(nrr_q) )
ALLOCATE ( irvec_kk(3,nrr_k) )
ALLOCATE ( irvec_qq(3,nrr_q) )
DO ir = 1, nrr_k
ndegen_kk(ir) = ndegen_k(ir)
irvec_kk(:,ir) = irvec(:,ir)
ENDDO
DO ir = 1, nrr_q
ndegen_qq(ir) = ndegen_q(ir)
irvec_qq(:,ir) = irvec(:,ir)
ENDDO
DEALLOCATE(ndegen_k)
DEALLOCATE(ndegen_q)
DEALLOCATE(irvec)
! Determine the size of the respective WS sets based on the length of the matrices
nrr_k = SIZE(irvec_k(1,:))
nrr_q = SIZE(irvec_q(1,:))
nrr_g = SIZE(irvec_g(1,:))
!
#ifndef __MPI
! Open like this only in sequential. Otherwize open with MPI-open
@ -372,7 +376,7 @@
! read all quantities in Wannier representation from file
! in parallel case all pools read the same file
!
CALL epw_read
CALL epw_read(nrr_k, nrr_q, nrr_g)
!
ELSE !if not epwread (i.e. need to calculate fmt file)
!
@ -406,7 +410,7 @@
! SP : Let the user chose. If false use files on disk
IF (etf_mem == 0) THEN
ALLOCATE(epmatwe ( nbndsub, nbndsub, nrr_k, nmodes, nqc))
ALLOCATE (epmatwp ( nbndsub, nbndsub, nrr_k, nmodes, nrr_q))
ALLOCATE (epmatwp ( nbndsub, nbndsub, nrr_k, nmodes, nrr_g))
ELSE
ALLOCATE(epmatwe_mem ( nbndsub, nbndsub, nrr_k, nmodes))
epmatwe_mem(:,:,:,:) = czero
@ -415,31 +419,31 @@
! Hamiltonian
!
CALL hambloch2wan &
( nbnd, nbndsub, nks, nkstot, et, xk, cu, lwin, exband, nrr_k, irvec_kk, wslen, chw )
( nbnd, nbndsub, nks, nkstot, et, xk, cu, lwin, exband, nrr_k, irvec_k, wslen_k, chw )
!
! Kohn-Sham eigenvalues
!
IF (eig_read) THEN
WRITE (stdout,'(5x,a)') "Interpolating MB and KS eigenvalues"
CALL hambloch2wan &
( nbnd, nbndsub, nks, nkstot, et_ks, xk, cu, lwin, exband, nrr_k, irvec_kk, wslen, chw_ks )
( nbnd, nbndsub, nks, nkstot, et_ks, xk, cu, lwin, exband, nrr_k, irvec_k, wslen_k, chw_ks )
ENDIF
!
IF (vme) THEN
! Transform of position matrix elements
! PRB 74 195118 (2006)
CALL vmebloch2wan &
( nbnd, nbndsub, nks, nkstot, xk, cu, nrr_k, irvec_kk, wslen, lwin, exband )
( nbnd, nbndsub, nks, nkstot, xk, cu, nrr_k, irvec_k, wslen_k, lwin, exband )
ELSE
! Dipole
CALL dmebloch2wan &
( nbnd, nbndsub, nks, nkstot, dmec, xk, cu, nrr_k, irvec_kk, wslen, lwin, exband )
( nbnd, nbndsub, nks, nkstot, dmec, xk, cu, nrr_k, irvec_k, wslen_k, lwin, exband )
ENDIF
!
! Dynamical Matrix
!
IF (.not. lifc) CALL dynbloch2wan &
( nmodes, nqc, xqc, dynq, nrr_q, irvec_qq, wslen )
( nmodes, nqc, xqc, dynq, nrr_q, irvec_q, wslen_q )
!
!
! Electron-Phonon vertex (Bloch el and Bloch ph -> Wannier el and Bloch ph)
@ -457,11 +461,11 @@
IF (etf_mem == 0) THEN
CALL ephbloch2wane &
( nbnd, nbndsub, nks, nkstot, xk, cu, cuq, &
epmatq(:,:,:,imode,iq), nrr_k, irvec_kk, wslen, epmatwe(:,:,:,imode,iq) )
epmatq(:,:,:,imode,iq), nrr_k, irvec_k, wslen_k, epmatwe(:,:,:,imode,iq) )
ELSE
CALL ephbloch2wane &
( nbnd, nbndsub, nks, nkstot, xk, cu, cuq, &
epmatq(:,:,:,imode,iq), nrr_k, irvec_kk, wslen, epmatwe_mem(:,:,:,imode) )
epmatq(:,:,:,imode,iq), nrr_k, irvec_k, wslen_k, epmatwe_mem(:,:,:,imode) )
!
ENDIF
!
@ -481,18 +485,18 @@
IF (ionode) THEN
IF (etf_mem == 0) THEN
CALL ephbloch2wanp &
( nbndsub, nmodes, xqc, nqc, irvec_kk, nrr_k, nrr_q, epmatwe )
( nbndsub, nmodes, xqc, nqc, irvec_k, irvec_g, nrr_k, nrr_g, epmatwe )
ELSE
CALL ephbloch2wanp_mem &
( nbndsub, nmodes, xqc, nqc, irvec_kk, nrr_k, nrr_q, epmatwe_mem )
( nbndsub, nmodes, xqc, nqc, irvec_k, irvec_g, nrr_k, nrr_g, epmatwe_mem )
ENDIF
ENDIF
!
CALL mp_barrier(inter_pool_comm)
!
IF ( epwwrite ) THEN
CALL epw_write
CALL epw_read
CALL epw_write(nrr_k, nrr_q, nrr_g)
CALL epw_read(nrr_k, nrr_q, nrr_g)
ENDIF
!
ENDIF
@ -561,7 +565,7 @@
ALLOCATE(rdotk(nrr_k))
! This is simply because dgemv take only real number (not integer)
ALLOCATE(irvec_r(3,nrr_k))
irvec_r = REAL(irvec_kk,KIND=dp)
irvec_r = REAL(irvec_k,KIND=dp)
!
! SP: Create a look-up table for the exponential of the factor.
! This can only work with homogeneous fine grids.
@ -645,7 +649,7 @@
! SP: Compute the cfac only once here since the same are use in both hamwan2bloch and dmewan2bloch
! + optimize the 2\pi r\cdot k with Blas
CALL dgemv('t', 3, nrr_k, twopi, irvec_r, 3, xxk, 1, 0.0_DP, rdotk, 1 )
cfac(:) = exp( ci*rdotk ) / ndegen_kk(:)
cfac(:) = exp( ci*rdotk ) / ndegen_k(:)
!
CALL hamwan2bloch &
( nbndsub, nrr_k, cufkk, etf(:, ik), chw, cfac)
@ -934,7 +938,7 @@
!
IF (.not. lifc) THEN
CALL dynwan2bloch &
( nmodes, nrr_q, irvec_qq, ndegen_qq, xxq, uf, w2 )
( nmodes, nrr_q, irvec_q, ndegen_q, xxq, uf, w2 )
ELSE
CALL dynifc2blochf ( nmodes, rws, nrws, xxq, uf, w2 )
ENDIF
@ -966,7 +970,7 @@
!CALL start_clock ( 'cl2' )
IF (.NOT. longrange) THEN
CALL ephwan2blochp &
( nmodes, xxq, irvec_qq, ndegen_qq, nrr_q, uf, epmatwef, nbndsub, nrr_k )
( nmodes, xxq, irvec_g, ndegen_g, nrr_g, uf, epmatwef, nbndsub, nrr_k )
ENDIF
!CALL stop_clock ( 'cl2' )
!
@ -1005,10 +1009,10 @@
!
! SP: Look-up table is more effecient than calling the exp function.
DO ir = 1, nrr_k
cfac(ir) = ( tablex(irvec_kk(1,ir)+2*nk1+1,xkk1) *&
tabley(irvec_kk(2,ir)+2*nk2+1,xkk2) * tablez(irvec_kk(3,ir)+2*nk3+1,xkk3) ) / ndegen_kk(ir)
cfacq(ir) = ( tableqx(irvec_kk(1,ir)+2*nk1+1,xkq1) *&
tableqy(irvec_kk(2,ir)+2*nk2+1,xkq2) * tableqz(irvec_kk(3,ir)+2*nk3+1,xkq3) ) / ndegen_kk(ir)
cfac(ir) = ( tablex(irvec_k(1,ir)+2*nk1+1,xkk1) *&
tabley(irvec_k(2,ir)+2*nk2+1,xkk2) * tablez(irvec_k(3,ir)+2*nk3+1,xkk3) ) / ndegen_k(ir)
cfacq(ir) = ( tableqx(irvec_k(1,ir)+2*nk1+1,xkq1) *&
tableqy(irvec_k(2,ir)+2*nk2+1,xkq2) * tableqz(irvec_k(3,ir)+2*nk3+1,xkq3) ) / ndegen_k(ir)
ENDDO
!DBSP
!IF ( (iq == 1) .and. (ik ==12)) THEN
@ -1019,9 +1023,9 @@
!ENDIF
ELSE
CALL dgemv('t', 3, nrr_k, twopi, irvec_r, 3, xkk, 1, 0.0_DP, rdotk, 1 )
cfac(:) = exp( ci*rdotk(:) ) / ndegen_kk(:)
cfac(:) = exp( ci*rdotk(:) ) / ndegen_k(:)
CALL dgemv('t', 3, nrr_k, twopi, irvec_r, 3, xkq, 1, 0.0_DP, rdotk, 1 )
cfacq(:) = exp( ci*rdotk(:) ) / ndegen_kk(:)
cfacq(:) = exp( ci*rdotk(:) ) / ndegen_k(:)
ENDIF
!
! ------------------------------------------------------
@ -1049,14 +1053,14 @@
!
IF (eig_read) THEN
CALL vmewan2bloch &
( nbndsub, nrr_k, irvec_kk, cufkk, vmef(:,:,:, ikk), etf(:,ikk), etf_ks(:,ikk), chw_ks, cfac )
( nbndsub, nrr_k, irvec_k, cufkk, vmef(:,:,:, ikk), etf(:,ikk), etf_ks(:,ikk), chw_ks, cfac )
CALL vmewan2bloch &
( nbndsub, nrr_k, irvec_kk, cufkq, vmef(:,:,:, ikq), etf(:,ikq), etf_ks(:,ikq), chw_ks, cfacq )
( nbndsub, nrr_k, irvec_k, cufkq, vmef(:,:,:, ikq), etf(:,ikq), etf_ks(:,ikq), chw_ks, cfacq )
ELSE
CALL vmewan2bloch &
( nbndsub, nrr_k, irvec_kk, cufkk, vmef(:,:,:, ikk), etf(:,ikk), etf_ks(:,ikk), chw, cfac )
( nbndsub, nrr_k, irvec_k, cufkk, vmef(:,:,:, ikk), etf(:,ikk), etf_ks(:,ikk), chw, cfac )
CALL vmewan2bloch &
( nbndsub, nrr_k, irvec_kk, cufkq, vmef(:,:,:, ikq), etf(:,ikq), etf_ks(:,ikq), chw, cfacq )
( nbndsub, nrr_k, irvec_k, cufkq, vmef(:,:,:, ikq), etf(:,ikq), etf_ks(:,ikq), chw, cfacq )
ENDIF
ELSE
!
@ -1414,8 +1418,15 @@
DEALLOCATE(cfacq)
DEALLOCATE(rdotk)
DEALLOCATE(irvec_r)
DEALLOCATE(irvec_kk)
DEALLOCATE(irvec_qq)
DEALLOCATE(irvec_k)
DEALLOCATE(irvec_q)
DEALLOCATE(irvec_g)
DEALLOCATE(ndegen_k)
DEALLOCATE(ndegen_q)
DEALLOCATE(ndegen_g)
DEALLOCATE(wslen_k)
DEALLOCATE(wslen_q)
DEALLOCATE(wslen_g)
IF ( ALLOCATED(inv_tau_all) ) DEALLOCATE( inv_tau_all )
IF ( ALLOCATED(inv_tau_allcb) ) DEALLOCATE( inv_tau_allcb )
IF ( ALLOCATED(zi_allvb) ) DEALLOCATE( zi_allvb )
@ -1426,12 +1437,12 @@
END SUBROUTINE ephwann_shuffle
!
!-------------------------------------------
SUBROUTINE epw_write
SUBROUTINE epw_write(nrr_k, nrr_q, nrr_g)
!-------------------------------------------
!
USE epwcom, ONLY : nbndsub, vme, eig_read, etf_mem
USE pwcom, ONLY : ef, nelec, isk
USE elph2, ONLY : nrr_k, nrr_q, chw, rdw, cdmew, cvmew, chw_ks, &
USE elph2, ONLY : chw, rdw, cdmew, cvmew, chw_ks, &
zstar, epsi, epmatwp
USE ions_base, ONLY : amass, ityp, nat, tau
USE cell_base, ONLY : at, bg, omega, alat
@ -1446,9 +1457,27 @@
USE io_global, ONLY : ionode_id
!
implicit none
LOGICAL :: exst
INTEGER :: ibnd, jbnd, jmode, imode, irk, irq, ipol, lrepmatw
!
INTEGER, INTENT(in) :: nrr_k
!! Number of WS vectors for the electrons
INTEGER, INTENT(in) :: nrr_q
!! Number of WS vectors for the phonons
INTEGER, INTENT(in) :: nrr_g
!! Number of WS vectors for the electron-phonons
INTEGER :: ibnd, jbnd
!! Band index
INTEGER :: jmode, imode
!! Mode index
INTEGER :: irk, irq, irg
!! WS vector looping index on electron, phonons and el-ph
INTEGER :: ipol
!! Cartesian direction (polarison direction)
INTEGER :: lrepmatw
!! Record length
CHARACTER (len=256) :: filint
!! Name of the file
LOGICAL :: exst
!! The file exists
!
WRITE(6,'(/5x,"Writing Hamiltonian, Dynamical matrix and EP vertex in Wann rep to file"/)')
!
@ -1475,7 +1504,7 @@
WRITE (crystal,*) isk
WRITE (crystal,*) noncolin
WRITE (epwdata,*) ef
WRITE (epwdata,*) nbndsub, nrr_k, nmodes, nrr_q
WRITE (epwdata,*) nbndsub, nrr_k, nmodes, nrr_q, nrr_g
WRITE (epwdata,*) zstar, epsi
!
DO ibnd = 1, nbndsub
@ -1511,8 +1540,8 @@
lrepmatw = 2 * nbndsub * nbndsub * nrr_k * nmodes
filint = trim(prefix)//'.epmatwp'
CALL diropn (iunepmatwp, 'epmatwp', lrepmatw, exst)
DO irq = 1, nrr_q
CALL davcio ( epmatwp(:,:,:,:,irq), lrepmatw, iunepmatwp, irq, +1 )
DO irg = 1, nrr_g
CALL davcio ( epmatwp(:,:,:,:,irg), lrepmatw, iunepmatwp, irg, +1 )
ENDDO
!
CLOSE(iunepmatwp)
@ -1534,11 +1563,11 @@
END SUBROUTINE epw_write
!---------------------------------
!---------------------------------
SUBROUTINE epw_read()
SUBROUTINE epw_read(nrr_k, nrr_q, nrr_g)
!---------------------------------
USE epwcom, ONLY : nbndsub, vme, eig_read, etf_mem, lifc
USE pwcom, ONLY : ef
USE elph2, ONLY : nrr_k, nrr_q, chw, rdw, epmatwp, &
USE elph2, ONLY : chw, rdw, epmatwp, &
cdmew, cvmew, chw_ks, zstar, epsi
USE ions_base, ONLY : nat
USE phcom, ONLY : nmodes
@ -1556,11 +1585,29 @@
!
implicit none
!
LOGICAL :: exst
INTEGER, INTENT (out) :: nrr_k
!! Number of WS vectors for the electrons
INTEGER, INTENT (out) :: nrr_q
!! Number of WS vectors for the phonons
INTEGER, INTENT (out) :: nrr_g
!! Number of WS vectors for the electron-phonons
INTEGER :: ibnd, jbnd
!! Band index
INTEGER :: jmode, imode
!! Mode index
INTEGER :: irk, irq, irg
!! WS vector looping index on electron, phonons and el-ph
INTEGER :: ipol
!! Cartesian direction (polarison direction)
INTEGER :: lrepmatw
!! Record length
INTEGER :: ios
!! Status of files
CHARACTER (len=256) :: filint
INTEGER :: ibnd, jbnd, jmode, imode, irk, irq, &
ipol, ios, lrepmatw
!
!! Name of the file
LOGICAL :: exst
!! The file exists
!
WRITE(stdout,'(/5x,"Reading Hamiltonian, Dynamical matrix and EP vertex in Wann rep from file"/)')
call flush(6)
!
@ -1582,7 +1629,7 @@
IF (ios /= 0) call errore ('ephwann_shuffle', 'error opening dmedata.fmt',iundmedata)
ENDIF
READ (epwdata,*) ef
READ (epwdata,*) nbndsub, nrr_k, nmodes, nrr_q
READ (epwdata,*) nbndsub, nrr_k, nmodes, nrr_q, nrr_g
READ (epwdata,*) zstar, epsi
!
ENDIF
@ -1600,6 +1647,9 @@
!
CALL mp_bcast (nrr_q, ionode_id, inter_pool_comm)
CALL mp_bcast (nrr_q, root_pool, intra_pool_comm)
!
CALL mp_bcast (nrr_g, ionode_id, inter_pool_comm)
CALL mp_bcast (nrr_g, root_pool, intra_pool_comm)
!
CALL mp_bcast (zstar, ionode_id, inter_pool_comm)
CALL mp_bcast (zstar, root_pool, intra_pool_comm)
@ -1666,7 +1716,7 @@
IF (lifc) CALL read_ifc
!
IF (etf_mem == 0) then
IF (.not. ALLOCATED(epmatwp)) ALLOCATE ( epmatwp ( nbndsub, nbndsub, nrr_k, nmodes, nrr_q) )
IF (.not. ALLOCATED(epmatwp)) ALLOCATE ( epmatwp ( nbndsub, nbndsub, nrr_k, nmodes, nrr_g) )
epmatwp = czero
IF (mpime.eq.ionode_id) THEN
! SP: The call to epmatwp is now inside the loop
@ -1677,8 +1727,8 @@
lrepmatw = 2 * nbndsub * nbndsub * nrr_k * nmodes
filint = trim(prefix)//'.epmatwp'
CALL diropn (iunepmatwp, 'epmatwp', lrepmatw, exst)
DO irq = 1, nrr_q
CALL davcio ( epmatwp(:,:,:,:,irq), lrepmatw, iunepmatwp, irq, -1 )
DO irg = 1, nrr_g
CALL davcio ( epmatwp(:,:,:,:,irg), lrepmatw, iunepmatwp, irg, -1 )
ENDDO
!
CLOSE(iunepmatwp)

146
EPW/src/ephwann_shuffle_mem.f90
View File

@ -50,8 +50,8 @@
USE io_global, ONLY : stdout, ionode
USE io_epw, ONLY : lambda_phself, linewidth_phself, iunepmatwe, &
iunepmatwp, crystal, iunepmatwp2
USE elph2, ONLY : nrr_k, nrr_q, cu, cuq, lwin, lwinq, irvec, ndegen_k,&
ndegen_q, wslen, chw, chw_ks, cvmew, cdmew, rdw, &
USE elph2, ONLY : cu, cuq, lwin, lwinq, &
chw, chw_ks, cvmew, cdmew, rdw, &
epmatq, wf, etf, etf_k, etf_ks, xqf, xkf, &
wkf, dynq, nqtotf, nkqf, epf17, nkf, nqf, et_ks, &
ibndmin, ibndmax, lambda_all, dmec, dmef, vmef, &
@ -60,6 +60,7 @@
nkqtotf, sigmar_all, zi_allvb, inv_tau_all, Fi_all, &
F_current, F_SERTA, inv_tau_allcb, zi_allcb, exband,&
F_currentcb, F_SERTAcb, Fi_allcb
USE wigner, ONLY : wigner_seitz_wrap
USE transportcom, ONLY : transp_temp, mobilityh_save, mobilityel_save, lower_bnd, &
upper_bnd, ixkqf_tr, s_BZtoIBZ_full
USE wan2bloch, ONLY : dmewan2bloch, hamwan2bloch, dynwan2bloch, &
@ -105,14 +106,6 @@
CHARACTER (len=30) :: myfmt
!! Variable used for formatting output
!
INTEGER, ALLOCATABLE :: ndegen_kk(:)
!! Number of degeneresence with nrr_k bounds.
INTEGER, ALLOCATABLE :: ndegen_qq(:)
!! Number of degeneresence with nrr_q bounds.
INTEGER, ALLOCATABLE :: irvec_kk(:,:)
!! Irevec on the nrr_k bounds
INTEGER, ALLOCATABLE :: irvec_qq(:,:)
!! Irevec on the nrr_k bounds
INTEGER :: ios
!! integer variable for I/O control
INTEGER :: iq
@ -159,6 +152,27 @@
!! Return virtual and resisdent memory from system
INTEGER :: ierr
!! Error status
INTEGER :: nrr_k
!! Number of WS points for electrons
INTEGER :: nrr_q
!! Number of WS points for phonons
INTEGER :: nrr_g
!! Number of WS points for electron-phonons
INTEGER, ALLOCATABLE :: irvec_k(:,:)
!! integer components of the ir-th Wigner-Seitz grid point in the basis
!! of the lattice vectors for electrons
INTEGER, ALLOCATABLE :: irvec_q(:,:)
!! integer components of the ir-th Wigner-Seitz grid point for phonons
INTEGER, ALLOCATABLE :: irvec_g(:,:)
!! integer components of the ir-th Wigner-Seitz grid point for electron-phonon
INTEGER, ALLOCATABLE :: ndegen_k (:)
!! Wigner-Seitz number of degenerescence (weights) for the electrons grid
INTEGER, ALLOCATABLE :: ndegen_q (:,:,:)
!! Wigner-Seitz weights for the phonon grid that depend on
!! atomic positions $R + \tau(nb) - \tau(na)$
INTEGER, ALLOCATABLE :: ndegen_g (:,:)
!! Wigner-Seitz weights for the electron-phonon grid that depend on
!! atomic positions $R - \tau(na)$
INTEGER, PARAMETER :: nrwsx=200
!! Maximum number of real-space Wigner-Seitz
!
@ -188,6 +202,12 @@
!! Wigner-Size supercell vectors, store in real instead of integer
REAL(kind=DP), ALLOCATABLE :: rdotk(:)
!! $r\cdot k$
REAL(kind=DP), ALLOCATABLE :: wslen_k(:)
!! real-space length for electrons, in units of alat
REAL(kind=DP), ALLOCATABLE :: wslen_q(:)
!! real-space length for phonons, in units of alat
REAL(kind=DP), ALLOCATABLE :: wslen_g(:)
!! real-space length for electron-phonons, in units of alat
!
COMPLEX(kind=DP) :: tablex (4*nk1+1,nkf1)
!! Look-up table for the exponential (speed optimization) in the case of
@ -262,16 +282,10 @@
cuq ( nbnd, nbndsub, nks), &
lwin ( nbnd, nks ), &
lwinq ( nbnd, nks ), &
exband ( nbnd ), &
irvec (3, 20*nk1*nk2*nk3), &
ndegen_k (20*nk1*nk2*nk3), &
ndegen_q (20*nq1*nq2*nq3), &
wslen(20*nk1*nk2*nk3) )
exband ( nbnd ) )
!
CALL start_clock ( 'ephwann' )
!
! DBSP
! HERE loadkmesh
IF ( epwread ) THEN
!
! Might have been pre-allocate depending of the restart configuration
@ -336,28 +350,17 @@
!
ALLOCATE( w2( 3*nat) )
!
! determine Wigner-Seitz points
! Determine Wigner-Seitz points
!
CALL wigner_seitz2 &
( nk1, nk2, nk3, nq1, nq2, nq3, nrr_k, nrr_q, irvec, wslen, ndegen_k, ndegen_q )
! Inside we allocate irvec_k, irvec_q, irvec_g, ndegen_k, ndegen_q, ndegen_g,
! wslen_k, wslen_q, wslen_g
CALL wigner_seitz_wrap ( nk1, nk2, nk3, nq1, nq2, nq3, irvec_k, irvec_q, irvec_g, &
ndegen_k, ndegen_q, ndegen_g, wslen_k, wslen_q, wslen_g )
!
ALLOCATE ( ndegen_kk(nrr_k) )
ALLOCATE ( ndegen_qq(nrr_q) )
ALLOCATE ( irvec_kk(3,nrr_k) )
ALLOCATE ( irvec_qq(3,nrr_q) )
DO ir = 1, nrr_k
ndegen_kk(ir) = ndegen_k(ir)
irvec_kk(:,ir) = irvec(:,ir)
ENDDO
DO ir = 1, nrr_q
ndegen_qq(ir) = ndegen_q(ir)
irvec_qq(:,ir) = irvec(:,ir)
ENDDO
DEALLOCATE(ndegen_k)
DEALLOCATE(ndegen_q)
DEALLOCATE(irvec)
! Determine the size of the respective WS sets based on the length of the matrices
nrr_k = SIZE(irvec_k(1,:))
nrr_q = SIZE(irvec_q(1,:))
nrr_g = SIZE(irvec_g(1,:))
!
#ifndef __MPI
! Open like this only in sequential. Otherwize open with MPI-open
@ -413,31 +416,32 @@
! Hamiltonian
!
CALL hambloch2wan &
( nbnd, nbndsub, nks, nkstot, et, xk, cu, lwin, exband, nrr_k, irvec_kk, wslen, chw )
( nbnd, nbndsub, nks, nkstot, et, xk, cu, lwin, exband, nrr_k, irvec_k, wslen_k, chw )
!
! Kohn-Sham eigenvalues
!
IF (eig_read) THEN
WRITE (stdout,'(5x,a)') "Interpolating MB and KS eigenvalues"
CALL hambloch2wan &
( nbnd, nbndsub, nks, nkstot, et_ks, xk, cu, lwin, exband, nrr_k, irvec_kk, wslen, chw_ks )
( nbnd, nbndsub, nks, nkstot, et_ks, xk, cu, lwin, exband, nrr_k, irvec_k, wslen_k, chw_ks )
ENDIF
!
IF (vme) THEN
! Transform of position matrix elements
! PRB 74 195118 (2006)
CALL vmebloch2wan &
( nbnd, nbndsub, nks, nkstot, xk, cu, nrr_k, irvec_kk, wslen, lwin, exband )
( nbnd, nbndsub, nks, nkstot, xk, cu, nrr_k, irvec_k, wslen_k, lwin, exband )
ELSE
! Dipole
CALL dmebloch2wan &
( nbnd, nbndsub, nks, nkstot, dmec, xk, cu, nrr_k, irvec_kk, wslen, lwin, exband )
( nbnd, nbndsub, nks, nkstot, dmec, xk, cu, nrr_k, irvec_k, wslen_k, lwin, exband )
ENDIF
!
! Dynamical Matrix
!
IF (.not. lifc) CALL dynbloch2wan &
( nmodes, nqc, xqc, dynq, nrr_q, irvec_qq, wslen )
( nmodes, nqc, xqc, dynq, nrr_q, irvec_q, wslen_q )
!
!
! Electron-Phonon vertex (Bloch el and Bloch ph -> Wannier el and Bloch ph)
!
@ -453,7 +457,7 @@
!
CALL ephbloch2wane &
( nbnd, nbndsub, nks, nkstot, xk, cu, cuq, &
epmatq (:,:,:,imode,iq), nrr_k, irvec_kk, wslen, epmatwe_mem(:,:,:,imode) )
epmatq (:,:,:,imode,iq), nrr_k, irvec_k, wslen_k, epmatwe_mem(:,:,:,imode) )
!
ENDDO
! Only the master node writes
@ -470,7 +474,7 @@
! Only master perform this task. Need to be parallelize in the future (SP)
IF (ionode) THEN
CALL ephbloch2wanp_mem &
( nbndsub, nmodes, xqc, nqc, irvec_kk, nrr_k, nrr_q, epmatwe_mem )
( nbndsub, nmodes, xqc, nqc, irvec_k, irvec_g, nrr_k, nrr_g, epmatwe_mem )
ENDIF
!
CALL mp_barrier(inter_pool_comm)
@ -539,7 +543,7 @@
ALLOCATE(rdotk(nrr_k))
! This is simply because dgemv take only real number (not integer)
ALLOCATE(irvec_r(3,nrr_k))
irvec_r = REAL(irvec_kk,KIND=dp)
irvec_r = REAL(irvec_k,KIND=dp)
!
! SP: Create a look-up table for the exponential of the factor.
! This can only work with homogeneous fine grids.
@ -622,7 +626,7 @@
! SP: Compute the cfac only once here since the same are use in both hamwan2bloch and dmewan2bloch
! + optimize the 2\pi r\cdot k with Blas
CALL dgemv('t', 3, nrr_k, twopi, irvec_r, 3, xxk, 1, 0.0_DP, rdotk, 1 )
cfac(:) = exp( ci*rdotk ) / ndegen_kk(:)
cfac(:) = exp( ci*rdotk ) / ndegen_k(:)
!
CALL hamwan2bloch &
( nbndsub, nrr_k, cufkk, etf (:, ik), chw, cfac)
@ -824,18 +828,22 @@
ENDIF
IF ( iterative_bte ) THEN
!
CALL F_read(iter, iq_restart, nqf,nkqtotf/2, error_h, error_el)
IF (int_mob .AND. carrier) THEN
CALL F_read(iter, iq_restart, nqf, nkqtotf/2, error_h, error_el, .TRUE.)
ELSE
CALL F_read(iter, iq_restart, nqf, nkqtotf/2, error_h, error_el, .FALSE.)
ENDIF
!
IF (int_mob .OR. (ncarrier < 1E5)) THEN
IF ( error_h < eps2 ) WRITE(stdout,'(5x,a)') repeat('=',67)
IF ( error_h < eps2 ) &
WRITE( stdout,'(5x,"IBTE is converged with value for hole mobility of",1E18.6," "/)') mobilityh_save
WRITE( stdout,'(5x,"IBTE is converged with value for hole mobility of",1E18.6," "/)') MAXVAL(mobilityh_save(:))
IF ( error_h < eps2 ) WRITE(stdout,'(5x,a)') repeat('=',67)
ENDIF
IF (int_mob .OR. (ncarrier > 1E5)) THEN
IF ( error_el < eps2 ) WRITE(stdout,'(5x,a)') repeat('=',67)
IF ( error_el < eps2 ) &
WRITE( stdout,'(5x,"IBTE is converged with value for electron mobility of",1E18.6," "/)') mobilityel_save
WRITE( stdout,'(5x,"IBTE is converged with value for electron mobility of",1E18.6," "/)') MAXVAL(mobilityel_save(:))
IF ( error_el < eps2 ) WRITE(stdout,'(5x,a)') repeat('=',67)
ENDIF
!
@ -862,8 +870,8 @@
ENDIF
WRITE(stdout,'(/5x,"Iteration number:", i10," "/)') iter
!
IF (nstemp > 1) CALL errore('ephwann_shuffle', &
'Iterative BTE can only be done at 1 temperature, nstemp = 1.',1)
!IF (nstemp > 1) CALL errore('ephwann_shuffle', &
! 'Iterative BTE can only be done at 1 temperature, nstemp = 1.',1)
!
IF (iter > maxiter) CALL errore('ephwann_shuffle', &
'The iteration reached the maximum but did not converge. ',1)
@ -875,7 +883,7 @@
!
!iter = iter +1
!
DO iq=iq_restart, nqf
DO iq = iq_restart, nqf
!
CALL start_clock ( 'ep-interp' )
!
@ -896,7 +904,7 @@
!
IF (.not. lifc) THEN
CALL dynwan2bloch &
( nmodes, nrr_q, irvec_qq, ndegen_qq, xxq, uf, w2 )
( nmodes, nrr_q, irvec_q, ndegen_q, xxq, uf, w2 )
ELSE
CALL dynifc2blochf ( nmodes, rws, nrws, xxq, uf, w2 )
ENDIF
@ -935,7 +943,7 @@
!CALL start_clock ( 'cl2' )
IF (.NOT. longrange) THEN
CALL ephwan2blochp_mem &
(imode, nmodes, xxq, irvec_qq, ndegen_qq, nrr_q, epmatwef, nbndsub, nrr_k )
(imode, nmodes, xxq, irvec_g, ndegen_g, nrr_g, epmatwef, nbndsub, nrr_k )
ENDIF
!CALL stop_clock ( 'cl2' )
!
@ -974,10 +982,10 @@
!
! SP: Look-up table is more effecient than calling the exp function.
DO ir = 1, nrr_k
cfac(ir) = ( tablex(irvec_kk(1,ir)+2*nk1+1,xkk1) *&
tabley(irvec_kk(2,ir)+2*nk2+1,xkk2) * tablez(irvec_kk(3,ir)+2*nk3+1,xkk3) ) / ndegen_kk(ir)
cfacq(ir) = ( tableqx(irvec_kk(1,ir)+2*nk1+1,xkq1) *&
tableqy(irvec_kk(2,ir)+2*nk2+1,xkq2) * tableqz(irvec_kk(3,ir)+2*nk3+1,xkq3) ) / ndegen_kk(ir)
cfac(ir) = ( tablex(irvec_k(1,ir)+2*nk1+1,xkk1) *&
tabley(irvec_k(2,ir)+2*nk2+1,xkk2) * tablez(irvec_k(3,ir)+2*nk3+1,xkk3) ) / ndegen_k(ir)
cfacq(ir) = ( tableqx(irvec_k(1,ir)+2*nk1+1,xkq1) *&
tableqy(irvec_k(2,ir)+2*nk2+1,xkq2) * tableqz(irvec_k(3,ir)+2*nk3+1,xkq3) ) / ndegen_k(ir)
ENDDO
!DBSP
!IF ( (iq == 1) .and. (ik ==12)) THEN
@ -988,9 +996,9 @@
!ENDIF
ELSE
CALL dgemv('t', 3, nrr_k, twopi, irvec_r, 3, xkk, 1, 0.0_DP, rdotk, 1 )
cfac(:) = exp( ci*rdotk(:) ) / ndegen_kk(:)
cfac(:) = exp( ci*rdotk(:) ) / ndegen_k(:)
CALL dgemv('t', 3, nrr_k, twopi, irvec_r, 3, xkq, 1, 0.0_DP, rdotk, 1 )
cfacq(:) = exp( ci*rdotk(:) ) / ndegen_kk(:)
cfacq(:) = exp( ci*rdotk(:) ) / ndegen_k(:)
ENDIF
!
! ------------------------------------------------------
@ -1018,14 +1026,14 @@
!
IF (eig_read) THEN
CALL vmewan2bloch &
( nbndsub, nrr_k, irvec_kk, cufkk, vmef(:,:,:, ikk), etf(:,ikk), etf_ks(:,ikk), chw_ks, cfac )
( nbndsub, nrr_k, irvec_k, cufkk, vmef(:,:,:, ikk), etf(:,ikk), etf_ks(:,ikk), chw_ks, cfac )
CALL vmewan2bloch &
( nbndsub, nrr_k, irvec_kk, cufkq, vmef(:,:,:, ikq), etf(:,ikq), etf_ks(:,ikq), chw_ks, cfacq )
( nbndsub, nrr_k, irvec_k, cufkq, vmef(:,:,:, ikq), etf(:,ikq), etf_ks(:,ikq), chw_ks, cfacq )
ELSE
CALL vmewan2bloch &
( nbndsub, nrr_k, irvec_kk, cufkk, vmef(:,:,:, ikk), etf(:,ikk), etf_ks(:,ikk), chw, cfac )
( nbndsub, nrr_k, irvec_k, cufkk, vmef(:,:,:, ikk), etf(:,ikk), etf_ks(:,ikk), chw, cfac )
CALL vmewan2bloch &
( nbndsub, nrr_k, irvec_kk, cufkq, vmef(:,:,:, ikq), etf(:,ikq), etf_ks(:,ikq), chw, cfacq )
( nbndsub, nrr_k, irvec_k, cufkq, vmef(:,:,:, ikq), etf(:,ikq), etf_ks(:,ikq), chw, cfacq )
ENDIF
ELSE
!
@ -1386,8 +1394,14 @@
DEALLOCATE(cfacq)
DEALLOCATE(rdotk)
DEALLOCATE(irvec_r)
DEALLOCATE(irvec_kk)
DEALLOCATE(irvec_qq)
DEALLOCATE(irvec_k)
DEALLOCATE(irvec_q)
DEALLOCATE(ndegen_k)
DEALLOCATE(ndegen_q)
DEALLOCATE(ndegen_g)
DEALLOCATE(wslen_k)
DEALLOCATE(wslen_q)
DEALLOCATE(wslen_g)
!
IF (.not. iterative_bte) CALL transport_coeffs (ef0,efcb)
!

136
EPW/src/wan2bloch.f90
View File

@ -160,7 +160,7 @@
END SUBROUTINE hamwan2bloch
!
!--------------------------------------------------------------------------
SUBROUTINE dynwan2bloch ( nmodes, nrr, irvec, ndegen, xxq, cuf, eig)
SUBROUTINE dynwan2bloch ( nmodes, nrr_q, irvec_q, ndegen_q, xxq, cuf, eig)
!--------------------------------------------------------------------------
!!
!!
@ -184,16 +184,17 @@
USE elph2, ONLY : rdw, epsi, zstar
USE epwcom, ONLY : lpolar, lphase
USE constants_epw, ONLY : twopi, ci, czero, zero, one, eps12
USE rigid, ONLY : cdiagh2
!
implicit none
!
INTEGER, INTENT (in) :: nmodes
!! number of modes (possibly of the optimal subspace)
INTEGER, INTENT (in) :: nrr
INTEGER, INTENT (in) :: nrr_q
!! number of WS points
INTEGER, INTENT (in) :: irvec(3, nrr)
INTEGER, INTENT (in) :: irvec_q(3, nrr_q)
!! coordinates of phononic WS points
INTEGER, INTENT (in) :: ndegen(nrr)
INTEGER, INTENT (in) :: ndegen_q(nrr_q, nat, nat)
!! degeneracy of WS points
!
REAL(kind=DP), INTENT (in) :: xxq(3)
@ -251,11 +252,21 @@
xq = xxq
chf(:,:) = czero
!
DO ir = 1, nrr
DO ir = 1, nrr_q
!
rdotk = twopi * dot_product( xq, dble(irvec( :, ir) ))
cfac = exp( ci*rdotk ) / dble( ndegen(ir) )
chf = chf + cfac * rdw(:,:, ir )
rdotk = twopi * dot_product( xq, dble(irvec_q( :, ir) ))
!
DO na = 1, nat
DO nb = 1, nat
IF ( ndegen_q(ir, na, nb) > 0 ) THEN
cfac = exp( ci*rdotk ) / dble( ndegen_q(ir,na,nb) )
! To map atom coordinate to mode basis.
chf(3*(na-1)+1:3*na, 3*(nb-1)+1:3*nb) = &
chf(3*(na-1)+1:3*na, 3*(nb-1)+1:3*nb) &
+ cfac * rdw(3*(na-1)+1:3*na, 3*(nb-1)+1:3*nb, ir )
ENDIF
ENDDO
ENDDO
!
ENDDO
!
@ -298,9 +309,10 @@
ENDDO
ENDDO
!
CALL zhpevx ('V', 'A', 'U', nmodes, champ, zero, zero, &
0, 0, -one, neig, w, cz, nmodes, cwork, &
rwork, iwork, ifail, info)
!CALL zhpevx ('V', 'A', 'U', nmodes, champ, zero, zero, &
! 0, 0, -one, neig, w, cz, nmodes, cwork, &
! rwork, iwork, ifail, info)
CALL cdiagh2(nmodes,chf,nmodes,w,cz)
!
! clean noise
DO jmode=1,nmodes
@ -1012,31 +1024,33 @@
END SUBROUTINE vmewan2bloch
!
!---------------------------------------------------------------------------
SUBROUTINE ephwan2blochp ( nmodes, xxq, irvec, ndegen, nrr_q, cuf, epmatf, nbnd, nrr_k )
SUBROUTINE ephwan2blochp ( nmodes, xxq, irvec_g, ndegen_g, nrr_g, cuf, epmatf, nbnd, nrr_k )
!---------------------------------------------------------------------------
!!
!! even though this is for phonons, I use the same notations
!! adopted for the electronic case (nmodes->nmodes etc)
!!
USE kinds, only : DP
USE epwcom, only : etf_mem
USE elph2, only : epmatwp
USE constants_epw, ONLY : twopi, ci, czero, cone
USE io_epw, ONLY : iunepmatwp, iunepmatwp2
USE mp_global, ONLY : mp_sum
USE mp_world, ONLY : world_comm
USE parallel_include
USE kinds, ONLY : DP
USE epwcom, ONLY : etf_mem
USE elph2, ONLY : epmatwp
USE constants_epw, ONLY : twopi, ci, czero, cone
USE io_epw, ONLY : iunepmatwp, iunepmatwp2
USE mp_global, ONLY : mp_sum
USE mp_world, ONLY : world_comm
USE parallel_include, ONLY : MPI_OFFSET_KIND, MPI_SEEK_SET, &
MPI_DOUBLE_PRECISION, MPI_STATUS_IGNORE
USE ions_base, ONLY : nat
implicit none
!
! input variables
!
INTEGER, INTENT (in) :: nmodes
!! Total number of modes
INTEGER, INTENT (in) :: nrr_q
INTEGER, INTENT (in) :: nrr_g
!! Number of phononic WS points
INTEGER, INTENT (in) :: irvec ( 3, nrr_q)
INTEGER, INTENT (in) :: irvec_g ( 3, nrr_g)
!! Coordinates of WS points
INTEGER, INTENT (in) :: ndegen (nrr_q)
INTEGER, INTENT (in) :: ndegen_g (nrr_g, nat)
!! Number of degeneracy of WS points
INTEGER, INTENT (in) :: nbnd
!! Number of bands
@ -1059,6 +1073,8 @@
!! Ending ir for this pool
INTEGER :: ierr
!! Return if there is an error
INTEGER :: na
!! Atom index
#if defined(__MPI)
INTEGER (kind=MPI_OFFSET_KIND) :: lrepmatw
!! Offset to tell where to start reading the file
@ -1076,7 +1092,7 @@
!
COMPLEX(kind=DP) :: eptmp( nbnd, nbnd, nrr_k, nmodes)
!! Temporary matrix to store el-ph
COMPLEX(kind=DP) :: cfac(nrr_q)
COMPLEX(kind=DP) :: cfac(nat, nrr_g)
!! Factor for the FT
COMPLEX(kind=DP), ALLOCATABLE :: epmatw( :,:,:,:)
!! El-ph matrix elements
@ -1093,24 +1109,29 @@
! g~(R_e,q') is epmatf(nmodes, nmodes, ik )
! every pool works with its own subset of k points on the fine grid
!
CALL para_bounds(ir_start, ir_stop, nrr_q)
CALL para_bounds(ir_start, ir_stop, nrr_g)
!
eptmp = czero
cfac(:) = czero
eptmp(:,:,:,:) = czero
cfac(:,:) = czero
!
DO ir = ir_start, ir_stop
!
! note xxq is assumed to be already in cryst coord
!
rdotk = twopi * dot_product ( xxq, dble(irvec(:, ir)) )
cfac(ir) = exp( ci*rdotk ) / dble( ndegen(ir) )
!
! note xxq is assumed to be already in cryst coord
!
rdotk = twopi * dot_product ( xxq, dble(irvec_g(:, ir)) )
DO na = 1, nat
IF (ndegen_g(ir,na) > 0) &
cfac(na,ir) = exp( ci*rdotk ) / dble( ndegen_g(ir,na) )
ENDDO
ENDDO
!
IF (etf_mem == 0) then
!
! SP: This is faster by 20 %
Call zgemv( 'n', nbnd * nbnd * nrr_k * nmodes, ir_stop - ir_start + 1, cone, &
epmatwp(1,1,1,1,ir_start), nbnd * nbnd * nrr_k * nmodes, cfac(ir_start), 1, czero, eptmp, 1 )
DO na = 1, nat
Call zgemv( 'n', nbnd * nbnd * nrr_k * 3, ir_stop - ir_start + 1, cone, &
epmatwp(:,:,:,3*(na-1)+1:3*na,ir_start:ir_stop), nbnd * nbnd * nrr_k * 3, &
cfac(na,ir_start:ir_stop), 1, czero, eptmp(:,:,:,3*(na-1)+1:3*na), 1 )
ENDDO
!
ELSE
!
@ -1159,7 +1180,10 @@
CALL rwepmatw ( epmatw, nbnd, nrr_k, nmodes, ir, iunepmatwp, -1)
#endif
!
CALL ZAXPY(nbnd * nbnd * nrr_k * nmodes, cfac(ir), epmatw, 1, eptmp, 1)
DO na = 1, nat
CALL ZAXPY(nbnd * nbnd * nrr_k * 3, cfac(na,ir), epmatw(:,:,:,3*(na-1)+1:3*na), 1, &
eptmp(:,:,:,3*(na-1)+1:3*na), 1)
ENDDO
!
ENDDO
DEALLOCATE(epmatw)
@ -1359,18 +1383,22 @@
END SUBROUTINE ephwan2bloch_mem
!
!---------------------------------------------------------------------------
SUBROUTINE ephwan2blochp_mem (imode, nmodes, xxq, irvec, ndegen, nrr_q, epmatf, nbnd, nrr_k )
SUBROUTINE ephwan2blochp_mem (imode, nmodes, xxq, irvec_g, ndegen_g, nrr_g, epmatf, nbnd, nrr_k )
!---------------------------------------------------------------------------
!!
!! Even though this is for phonons, I use the same notations
!! adopted for the electronic case (nmodes->nmodes etc)
!!
USE kinds, only : DP
USE constants_epw, ONLY : twopi, ci, czero
USE io_files, ONLY : prefix, tmp_dir
USE mp_global, ONLY : mp_sum
USE mp_world, ONLY : world_comm
USE parallel_include
USE kinds, ONLY : DP
USE constants_epw, ONLY : twopi, ci, czero
USE io_files, ONLY : prefix, tmp_dir
USE mp_global, ONLY : mp_sum
USE mp_world, ONLY : world_comm
USE parallel_include, ONLY : MPI_OFFSET_KIND, MPI_SEEK_SET, &
MPI_DOUBLE_PRECISION, MPI_STATUS_IGNORE, &
MPI_MODE_RDONLY,MPI_INFO_NULL
USE ions_base, ONLY : nat
implicit none
!
! input variables
@ -1379,11 +1407,11 @@
!! Current mode
INTEGER, INTENT (in) :: nmodes
!! Total number of modes
INTEGER, INTENT (in) :: nrr_q
INTEGER, INTENT (in) :: nrr_g
!! Number of phononic WS points
INTEGER, INTENT (in) :: irvec( 3, nrr_q)
INTEGER, INTENT (in) :: irvec_g( 3, nrr_g)
!! Coordinates of WS points
INTEGER, INTENT (in) :: ndegen(nrr_q)
INTEGER, INTENT (in) :: ndegen_g(nrr_g, nat)
!! Number of degeneracy of WS points
INTEGER, INTENT (in) :: nbnd
!! Number of bands
@ -1409,6 +1437,8 @@
!! Return the file unit
INTEGER :: ierr
!! Return if there is an error
INTEGER :: na
!! Index on atom
#if defined(__MPI)
INTEGER (kind=MPI_OFFSET_KIND) :: lrepmatw
!! Offset to tell where to start reading the file
@ -1419,7 +1449,7 @@
REAL(kind=DP) :: rdotk
!! Exponential for the FT
!
COMPLEX(kind=DP) :: cfac(nrr_q)
COMPLEX(kind=DP) :: cfac(nrr_g)
!! Factor for the FT
COMPLEX(kind=DP), ALLOCATABLE :: epmatw( :,:,:)
!! El-ph matrix elements
@ -1435,7 +1465,7 @@
! g~(R_e,q') is epmatf(nmodes, nmodes, ik )
! every pool works with its own subset of k points on the fine grid
!
CALL para_bounds(ir_start, ir_stop, nrr_q)
CALL para_bounds(ir_start, ir_stop, nrr_g)
!
#if defined(__MPI)
filint = trim(tmp_dir)//trim(prefix)//'.epmatwp1'
@ -1446,10 +1476,12 @@
cfac(:) = czero
!
DO ir = ir_start, ir_stop
!
! note xxq is assumed to be already in cryst coord
rdotk = twopi * dot_product ( xxq, dble(irvec(:, ir)) )
cfac(ir) = exp( ci*rdotk ) / dble( ndegen(ir) )
!
! note xxq is assumed to be already in cryst coord
rdotk = twopi * dot_product ( xxq, dble(irvec_g(:, ir)) )
na = (imode - 1) / 3 + 1
IF (ndegen_g(ir, na) > 0) &
cfac(ir) = exp( ci*rdotk ) / dble( ndegen_g(ir, na) )
ENDDO
!
ALLOCATE(epmatw( nbnd, nbnd, nrr_k))

774
EPW/src/wigner.f90 Normal file
View File

@ -0,0 +1,774 @@
!
! Copyright (C) 2010-2019 Samuel Ponce', Roxana Margine, Carla Verdi, Feliciano Giustino
! Copyright (C) 2007-2009 Jesse Noffsinger, Brad Malone, Feliciano Giustino
!
! This file is distributed under the terms of the GNU General Public
! License. See the file `LICENSE' in the root directory of the
! present distribution, or http://www.gnu.org/copyleft.gpl.txt .
!
!----------------------------------------------------------------------
MODULE wigner
!----------------------------------------------------------------------
!!
!! This module contains all the subroutine linked creation of Wigner-Seitz cell
!!
IMPLICIT NONE
!
CONTAINS
!-----------------------------------------------------------------
SUBROUTINE wigner_seitz_wrap (nk1, nk2, nk3, nq1, nq2, nq3, &
irvec_k, irvec_q, irvec_g, &
ndegen_k, ndegen_q, ndegen_g, &
wslen_k, wslen_q, wslen_g )
!-----------------------------------------------------------------
!!
!! June 2018 - SP - CV
!!
!! This routine wrap the call to three Wigner-Seitz routines:
!! wigner_seitzk : Contruct a grid of points that fall inside of (and eventually
!! on the surface of) the Wigner-Seitz supercell centered on the
!! origin of the Bravais lattice. Use for electronic properties.
!! wigner_seitzq : Creates a set of WS vectors for each pair of atoms tau(nb)-tau(na)
!! On exiting, ndegen_q contains the degeneracies of each pairs
!! of atoms while irvec_q contains the minimal communal sets of WS vectors.
!! Used for phonon properties
!! wigner_seitzg : Creates a set of WS vector for each atoms tau(na).
!! On exiting, ndegen_g contains the degeneracies of each atoms
!! while irvec_g contains the minimal communal sets of WS vectors.
!! Used for electron-phonon properties.
!!
!! Note 1: ndegen_k is always > 0 while ndegen_q and ndegen_g might contains 0 weigths.
!! Note 2: No sorting of vectors is needed anymore
!! Note 3: The dimension 20*nk1*nk2*nk3 should be safe enough.
!! Note 4: The Wigner-Seitz construction in EPW was done by constructing a cell
!! centred unit cell. This is fine for electronic properties (this is what is done in wannier90).
!! However for phonon or electron-phonon properties, one can have issues when the cell
!! is tilted for example.
!! The proper way is to construct a set of WS vectors centred on pairs of atoms (phonons)
!! or atoms (el-ph).
!! In the matdyn code, a FT grid is constructed with weights centred on pairs of atoms
!! and zeros everywhere else.
!! EPW now reproduced exactly the results of matdyn for the interpolated phonons at a
!! lower computation cost. Indeed we minimize the number of zeros by keeping the union
!! of values between all the cells.
!! In both cases this is very fast anyway but is important for el-ph properties.
!!
!-----------------------------------------------------------------
USE kinds, ONLY : DP
USE ions_base, ONLY : nat
!
implicit none
!
INTEGER, INTENT (in) :: nk1
!! size of the uniform k mesh
INTEGER, INTENT (in) :: nk2
!! size of the uniform k mesh
INTEGER, INTENT (in) :: nk3
!! size of the uniform k mesh
INTEGER, INTENT (in) :: nq1
!! size of the uniform q mesh
INTEGER, INTENT (in) :: nq2
!! size of the uniform q mesh
INTEGER, INTENT (in) :: nq3
!! size of the uniform q mesh
INTEGER, ALLOCATABLE, INTENT (out) :: irvec_k(:,:)
!! integer components of the ir-th Wigner-Seitz grid point in the basis
!! of the lattice vectors for electrons
INTEGER, ALLOCATABLE, INTENT (out) :: irvec_q(:,:)
!! integer components of the ir-th Wigner-Seitz grid point for phonons
INTEGER, ALLOCATABLE, INTENT (out) :: irvec_g(:,:)
!! integer components of the ir-th Wigner-Seitz grid point for electron-phonon
INTEGER, ALLOCATABLE, INTENT (out) :: ndegen_k (:)
!! Wigner-Seitz number of degenerescence (weights) for the electrons grid
INTEGER, ALLOCATABLE, INTENT (out) :: ndegen_q (:,:,:)
!! Wigner-Seitz weights for the phonon grid that depend on
!! atomic positions $R + \tau(nb) - \tau(na)$
INTEGER, ALLOCATABLE, INTENT (out) :: ndegen_g (:,:)
!! Wigner-Seitz weights for the electron-phonon grid that depend on
!! atomic positions $R - \tau(na)$
REAL(kind=DP), ALLOCATABLE, INTENT (out) :: wslen_k(:)
!! real-space length for electrons, in units of alat
REAL(kind=DP), ALLOCATABLE, INTENT (out) :: wslen_q(:)
!! real-space length for phonons, in units of alat
REAL(kind=DP), ALLOCATABLE, INTENT (out) :: wslen_g(:)
!! real-space length for electron-phonons, in units of alat
!
! Work Variables
INTEGER :: ir
!! Index for WS vectors
INTEGER :: nrr_k
!! maximum number of WS vectors for the electrons
INTEGER :: nrr_q
!! maximum number of WS vectors for the phonons
INTEGER :: nrr_g
!! maximum number of WS vectors for the electron-phonon
INTEGER :: irvec_kk (3,20*nk1*nk2*nk3)
!! local integer components of the ir-th Wigner-Seitz grid point
!! in the basis of the lattice vectors for electrons
INTEGER :: irvec_qq (3,20*nq1*nq2*nq3)
!! local integer components of the ir-th Wigner-Seitz grid point for phonons
INTEGER :: irvec_gg (3,20*nk1*nk2*nk3)
!! local integer components of the ir-th Wigner-Seitz grid point for electron-phonons
!! We use nk1 instead of nq1 because the k-grid is always larger or equal to q-grid.
INTEGER :: ndegen_kk (20*nk1*nk2*nk3)
!! local Wigner-Seitz number of degenerescence (weights) for the electrons grid
INTEGER :: ndegen_qq (20*nq1*nq2*nq3, nat, nat)
!! local Wigner-Seitz number of degenerescence (weights) for the phonons grid
INTEGER :: ndegen_gg (20*nk1*nk2*nk3, nat)
!! local Wigner-Seitz number of degenerescence (weights) for the electron-phonons grid
REAL(kind=DP) :: wslen_kk (20*nk1*nk2*nk3)
!! local real-space length for electrons, in units of alat
REAL(kind=DP) :: wslen_qq (20*nq1*nq2*nq3)
!! local real-space length for phonons, in units of alat
REAL(kind=DP) :: wslen_gg (20*nk1*nk2*nk3)
!! local real-space length for electron-phonon, in units of alat
!
! Check the bounds
IF ( nq1 > nk1 .OR. nq2 > nk2 .OR. nq3 > nk3 ) call errore &
('wigner_seitz_wrap',' the phonon grid should be smaller than electron grid',1)
!
! Now generated the un-sorted points for the electrons, phonons and electron-phonon
!
CALL wigner_seitzk ( nk1, nk2, nk3, irvec_kk, ndegen_kk, wslen_kk, nrr_k)
CALL wigner_seitzq ( nq1, nq2, nq3, irvec_qq, ndegen_qq, wslen_qq, nrr_q)
CALL wigner_seitzg ( nq1, nq2, nq3, irvec_gg, ndegen_gg, wslen_gg, nrr_g)
!
ALLOCATE ( irvec_k(3,nrr_k) )
ALLOCATE ( irvec_q(3,nrr_q) )
ALLOCATE ( irvec_g(3,nrr_g) )
ALLOCATE ( ndegen_k(nrr_k) )
ALLOCATE ( ndegen_q(nrr_q, nat, nat) )
ALLOCATE ( ndegen_g(nrr_g, nat) )
ALLOCATE ( wslen_k(nrr_k) )
ALLOCATE ( wslen_q(nrr_q) )
ALLOCATE ( wslen_g(nrr_g) )
!
! Create vectors with correct size.
DO ir = 1, nrr_k
ndegen_k(ir) = ndegen_kk(ir)
irvec_k(:,ir) = irvec_kk(:,ir)
wslen_k(ir) = wslen_kk(ir)
ENDDO
DO ir = 1, nrr_q
ndegen_q(ir,:,:) = ndegen_qq(ir,:,:)
irvec_q(:,ir) = irvec_qq(:,ir)
wslen_q(ir) = wslen_qq(ir)
ENDDO
DO ir = 1, nrr_g
ndegen_g(ir,:) = ndegen_gg(ir,:)
irvec_g(:,ir) = irvec_gg(:,ir)
wslen_g(ir) = wslen_gg(ir)
ENDDO
!
END SUBROUTINE wigner_seitz_wrap
!-----------------------------------------------------------------------------
!
!-----------------------------------------------------------------------------
SUBROUTINE wigner_seitzk (nk1, nk2, nk3, irvec_kk, ndegen_kk, wslen_kk, nrr_k)
!-----------------------------------------------------------------------------
!!
!! Calculates a grid of points that fall inside of (and eventually
!! on the surface of) the Wigner-Seitz supercell centered on the
!! origin of the Bravais lattice with primitive translations
!! nk1*a_1+nk2*a_2+nk3*a_3
!!
!-----------------------------------------------------------------
USE kinds, ONLY : DP
USE cell_base, ONLY : at
USE constants_epw, ONLY : eps8
!
implicit none
!
INTEGER, INTENT (in) :: nk1
!! size of the uniform k mesh
INTEGER, INTENT (in) :: nk2
!! size of the uniform k mesh
INTEGER, INTENT (in) :: nk3
!! size of the uniform k mesh
INTEGER, INTENT (out) :: irvec_kk(3,20*nk1*nk2*nk3)
!! integer components of the ir-th Wigner-Seitz grid point in the basis of the lattice vectors
INTEGER, INTENT (out) :: ndegen_kk(20*nk1*nk2*nk3)
!! Number of degeneracies
INTEGER, INTENT (out) :: nrr_k
!! number of Wigner-Seitz grid points
!
REAL(kind=DP), INTENT (out) :: wslen_kk(20*nk1*nk2*nk3)
!! real-space length, in units of alat
!
! Work variables
INTEGER :: n1, n2, n3
!! Index for the larger supercell
INTEGER :: i1, i2, i3
!! Index to compute |r-R| distance
INTEGER :: i
!! Iterative index
INTEGER :: ndiff(3)
!! Distances
INTEGER :: ipol, jpol
!! Cartesian direction
INTEGER, ALLOCATABLE :: ind(:)
!! Index of sorting
INTEGER :: nind
!! The metric tensor
!
REAL(kind=DP) :: adot(3,3)
!! Dot product between lattice vector
REAL(kind=DP) :: dist(125)
!! Contains the distance squared |r-R|^2
REAL(kind=DP) :: mindist
!! Minimum distance
REAL(kind=DP) :: tot
!! Sum of all the weigths
!
LOGICAL :: found
!! True if the vector has been found
!
nind = 20*nk1*nk2*nk3
IF (nind < 125) then
nind = 125
ENDIF
ALLOCATE (ind(nind))
!
DO ipol = 1, 3
DO jpol = 1, 3
adot (ipol, jpol) = dot_product ( at(:,ipol), at(:,jpol) )
ENDDO
ENDDO
!
nrr_k = 0
DO n1 = -2*nk1, 2*nk1
DO n2 = -2*nk2, 2*nk2
DO n3 = -2*nk3, 2*nk3
! Loop over the 5^3 = 125 points R. R=0 corresponds to i1=i2=i3=2, or icnt=63
i = 0
dist(:) = 0.d0
DO i1 = -2, 2
DO i2 = -2, 2
DO i3 = -2, 2
i = i + 1
! Calculate distance squared |r-R|^2
ndiff(1) = n1 - i1*nk1
ndiff(2) = n2 - i2*nk2
ndiff(3) = n3 - i3*nk3
DO ipol = 1, 3
DO jpol = 1, 3
dist(i) = dist(i) + dble(ndiff(ipol))*adot(ipol,jpol)*dble(ndiff(jpol))
ENDDO
ENDDO
!
ENDDO ! i3
ENDDO ! i2
ENDDO ! i1
!
! Sort the 125 vectors R by increasing value of |r-R|^2
ind(1) = 0 ! required for hpsort_eps (see the subroutine)
CALL hpsort_eps_epw( 125, dist, ind, eps8)
!
! Find all the vectors R with the (same) smallest |r-R|^2;
! if R=0 is one of them, then the current point r belongs to
! Wignez-Seitz cell => set found to true
!
found = .false.
i = 1
mindist = dist(1)
DO WHILE ( abs(dist(i)-mindist) < eps8 .and. i < 125 )
IF (ind(i) == 63) found = .true.
i = i + 1
ENDDO
!
IF (found) THEN
nrr_k = nrr_k + 1
ndegen_kk (nrr_k) = i - 1
irvec_kk (:, nrr_k) = (/ n1, n2, n3 /)
ENDIF
ENDDO ! n3
ENDDO ! n2
ENDDO ! n1
!
wslen_kk(:) = 0.d0
DO i = 1, nrr_k
DO ipol = 1, 3
DO jpol = 1, 3
wslen_kk(i) = wslen_kk(i) + dble(irvec_kk(ipol,i))*adot(ipol,jpol)*dble(irvec_kk(jpol,i))
ENDDO
ENDDO
wslen_kk(i) = sqrt(wslen_kk(i))
ENDDO
!
! Check the "sum rule"
!
tot = 0.d0
DO i = 1, nrr_k
tot = tot + 1.d0 / dble (ndegen_kk(i))
ENDDO
!
IF(abs(tot-dble(nk1*nk2*nk3)) > eps8) call errore &
('wigner_seitzk',' weights do not add up to nk1*nk2*nk3',1)
!
IF(nrr_k > 20*nk1*nk2*nk3) call errore &
('wigner_seitzk',' too many WS points, try to increase the hard bound 20*nk1*nk2*nk3',1)
!
DEALLOCATE(ind)
!
!-----------------------------------------------------------------------------
END SUBROUTINE wigner_seitzk
!-----------------------------------------------------------------------------
!
!-----------------------------------------------------------------
SUBROUTINE wigner_seitzq (nq1, nq2, nq3, irvec_qq, ndegen_qq, wslen_qq, nrr_q)
!-----------------------------------------------------------------
!!
!! Calculates a grid of points that fall inside of (and eventually
!! on the surface of) the Wigner-Seitz supercell centered on
!! each pair of atoms.
!! Follows Eq. 66 of PRB 55, 10355 (1997).
!! We are part of the WS if $R_b + \tau_{\kappa'} - \tau_\kappa$ is inside the
!! supercell.
!!
!-----------------------------------------------------------------
USE kinds, ONLY : DP
USE cell_base, ONLY : at, bg
USE ions_base, ONLY : nat, tau
USE constants_epw, ONLY : eps8, eps6
!
implicit none
!
INTEGER, INTENT (in) :: nq1
!! size of the uniform k mesh
INTEGER, INTENT (in) :: nq2
!! size of the uniform k mesh
INTEGER, INTENT (in) :: nq3
!! size of the uniform k mesh
INTEGER, INTENT (out) :: irvec_qq(3,20*nq1*nq2*nq3)
!! integer components of the ir-th Wigner-Seitz grid point in the basis of the lattice vectors
INTEGER, INTENT (out) :: ndegen_qq(20*nq1*nq2*nq3,nat,nat)
!! Number of degeneracies
INTEGER, INTENT (out) :: nrr_q
!! number of Wigner-Seitz grid points
REAL(kind=DP), INTENT (out) :: wslen_qq(20*nq1*nq2*nq3)
!! real-space length, in units of alat
!
! work variables
INTEGER :: n1, n2, n3
!! Index for the larger supercell
INTEGER :: i1, i2, i3
!! Index to compute |r-R| distance
INTEGER :: na, nb
!! Atom index
INTEGER :: i
!! Iterative index
INTEGER :: ir
!! Iterative index on the pair of atoms
INTEGER :: irtot
!! Iterative index on the combined pair of atoms
INTEGER :: ipol, jpol
!! Cartesian direction
INTEGER, ALLOCATABLE :: ind(:)
!! Index of sorting
INTEGER :: nind
!! The metric tensor
INTEGER :: nrr_tmp(nat,nat)
!! Temporary array that contains the max number of WS vectors
!! for a pair of atoms.
INTEGER :: irvec_tmp(3,20*nq1*nq2*nq3,nat,nat)
!! Temporary WS vectors for each atoms pair
INTEGER :: ndegen_tmp(20*nq1*nq2*nq3,nat,nat)
!! Temporary WS vectors weight for each atoms pair
!
REAL(kind=DP) :: adot(3,3)
!! Dot product between lattice vector
REAL(kind=DP) :: dist(125)
!! Contains the distance squared |r-R|^2
REAL(kind=DP) :: mindist
!! Minimum distance
REAL(kind=DP) :: tot, tot2
!! Sum of all the weigths
REAL(kind=DP) :: ndiff(3)
!! Distances. Must be now real because of fractional atoms
!
LOGICAL :: found
!! True if the vector has been found
!
nind = 20*nq1*nq2*nq3
IF (nind .lt. 125) THEN
nind = 125
ENDIF
ALLOCATE (ind(nind))
!
DO ipol = 1, 3
DO jpol = 1, 3
adot (ipol, jpol) = dot_product ( at(:,ipol), at(:,jpol) )
ENDDO
ENDDO
!
CALL cryst_to_cart(nat,tau(:,:),bg,-1)
!
! Loop over grid points r on a unit cell that is 8 times larger than a
! primitive supercell. In the end nrr contains the total number of grids
! points that have been found in the Wigner-Seitz cell
!
nrr_tmp(:,:) = 0
DO na=1, nat
DO nb=1, nat
DO n1 = -2*nq1, 2*nq1
DO n2 = -2*nq2, 2*nq2
DO n3 = -2*nq3, 2*nq3
!
! Loop over the 5^3 = 125 points R. R=0 corresponds to i1=i2=i3=2, or icnt=63
!
i = 0
dist(:) = 0.d0
DO i1 = -2, 2
DO i2 = -2, 2
DO i3 = -2, 2
i = i + 1
!
! Calculate distance squared |r-R|^2
!
ndiff(1) = n1 - i1*nq1 + tau(1,nb)-tau(1,na)
ndiff(2) = n2 - i2*nq2 + tau(2,nb)-tau(2,na)
ndiff(3) = n3 - i3*nq3 + tau(3,nb)-tau(3,na)
DO ipol = 1, 3
DO jpol = 1, 3
dist(i) = dist(i) + dble(ndiff(ipol))*adot(ipol,jpol)*dble(ndiff(jpol))
ENDDO
ENDDO
!
ENDDO
ENDDO
ENDDO
!
! Sort the 125 vectors R by increasing value of |r-R|^2
ind(1) = 0 ! required for hpsort_eps (see the subroutine)
CALL hpsort_eps_epw( 125, dist, ind, eps8)
!
! Find all the vectors R with the (same) smallest |r-R|^2;
! if R=0 is one of them, then the current point r belongs to
! Wignez-Seitz cell => set found to true
!
found = .false.
i = 1
mindist = dist(1)
DO WHILE ( abs(dist(i)-mindist) < eps8 .and. i < 125 )
IF (ind(i) == 63) found = .true.
i = i + 1
ENDDO
!
IF (found) THEN
nrr_tmp(na, nb) = nrr_tmp(na, nb) + 1
ndegen_tmp (nrr_tmp(na,nb), na, nb) = i - 1
irvec_tmp (:, nrr_tmp(na,nb), na, nb) = (/ n1, n2, n3 /)
ENDIF
!
ENDDO ! n3
ENDDO ! n2
ENDDO ! n3
ENDDO ! nb
ENDDO ! na
!
! Now creates a global set of WS vectors from all the atoms pair.
! Also remove the duplicated ones.
nrr_q = nrr_tmp(1, 1)
irvec_qq(:,:) = irvec_tmp(:, :, 1,1)
DO na = 1, nat
DO nb = 1, nat
DO ir=1, nrr_tmp(na, nb)
found = .false.
DO irtot = 1, nrr_q
IF ( ALL(irvec_tmp(:, ir, na, nb) == irvec_qq(:, irtot)) ) THEN
found = .true.
ENDIF
ENDDO !nrr
IF(.not. found) THEN
nrr_q = nrr_q + 1
irvec_qq(:, nrr_q) = irvec_tmp(:, ir, na, nb)
ENDIF
ENDDO ! ir
ENDDO ! nb
ENDDO ! na
!
! Creates a pair of atoms depended degeneracy array but with a number of WS
! vectors per pair that is equal to the global set. Populate with zero weights
! the one that are not part of that pair set.
ndegen_qq(:,:,:) = 0
DO na = 1, nat
DO nb = 1, nat
DO ir=1, nrr_tmp(na,nb)
DO irtot = 1, nrr_q
IF ( ALL(irvec_qq(:, irtot) == irvec_tmp(:, ir, na, nb)) ) THEN
ndegen_qq(irtot, na, nb) = ndegen_tmp(ir, na, nb)
ENDIF
ENDDO
ENDDO
ENDDO
ENDDO
!
DO na = 1, nat
DO nb = 1, nat
tot = 0.d0
tot2 = 0.d0
DO i = 1, nrr_q
IF (ndegen_qq(i, na, nb) > 0 ) THEN
tot2 = tot2 + 1.d0 / dble(ndegen_qq(i, na, nb))
ENDIF
ENDDO
DO i = 1, nrr_tmp(na, nb)
tot = tot + 1.d0 / dble(ndegen_tmp(i, na, nb))
ENDDO
!
!print*,'na, nb, tot tot2 ',na,nb,tot,tot2,nq1,nq2,nq3,dble(nq1*nq2*nq3)
IF(abs(tot-dble(nq1*nq2*nq3)) > eps8) call errore &
('wigner_seitzq',' weights do not add up to nq1*nq2*nq3',1)
IF(abs(tot-tot2) > eps8) call errore &
('wigner_seitzq',' weigths of pair of atoms is not equal to global weights',1)
ENDDO
ENDDO
!
IF(nrr_q > 20*nq1*nq2*nq3) call errore &
('wigner_seitzq','too many WS points, try to increase the bound 20*nq1*nq2*nq3',1)
!
! Now we compute the WS distances
wslen_qq(:) = 0.d0
DO i = 1, nrr_q
DO ipol = 1, 3
DO jpol = 1, 3
wslen_qq(i) = wslen_qq(i) + dble(irvec_qq(ipol,i))*adot(ipol,jpol)*dble(irvec_qq(jpol,i))
ENDDO
ENDDO
wslen_qq(i) = sqrt(wslen_qq(i))
ENDDO
!
CALL cryst_to_cart(nat,tau(:,:),at,1)
!
DEALLOCATE(ind)
!
! -----------------------------------------------------------------------------------------
END SUBROUTINE wigner_seitzq
! -----------------------------------------------------------------------------------------
!
!-----------------------------------------------------------------
SUBROUTINE wigner_seitzg (nq1, nq2, nq3, irvec_gg, ndegen_gg, wslen_gg, nrr_g)
!-----------------------------------------------------------------
!!
!! Calculates a grid of points that fall inside of (and eventually
!! on the surface of) the Wigner-Seitz supercell centered on
!! each atoms.
!! Follows Eq. 66 of PRB 55, 10355 (1997).
!! We are part of the WS if $R_b - \tau_{\kappa}$ is inside the
!! supercell.
!!
!-----------------------------------------------------------------
USE kinds, ONLY : DP
USE cell_base, ONLY : at, bg
USE ions_base, ONLY : nat, tau
USE constants_epw, ONLY : eps8
!
implicit none
!
INTEGER, INTENT (in) :: nq1
!! size of the uniform k mesh
INTEGER, INTENT (in) :: nq2
!! size of the uniform k mesh
INTEGER, INTENT (in) :: nq3
!! size of the uniform k mesh
INTEGER, INTENT (out) :: irvec_gg(3,20*nq1*nq2*nq3)
!! integer components of the ir-th Wigner-Seitz grid point in the basis of the lattice vectors
INTEGER, INTENT (out) :: ndegen_gg(20*nq1*nq2*nq3,nat)
!! Number of degeneracies
INTEGER, INTENT (out) :: nrr_g
!! number of Wigner-Seitz grid points
REAL(kind=DP), INTENT (out) :: wslen_gg(20*nq1*nq2*nq3)
!! real-space length, in units of alat
!
! work variables
INTEGER :: n1, n2, n3
!! Index for the larger supercell
INTEGER :: i1, i2, i3
!! Index to compute |r-R| distance
INTEGER :: na
!! Atom index
INTEGER :: i
!! Iterative index
INTEGER :: ir
!! Iterative index on the pair of atoms
INTEGER :: irtot
!! Iterative index on the combined pair of atoms
INTEGER :: ipol, jpol
!! Cartesian direction
INTEGER, ALLOCATABLE :: ind(:)
!! Index of sorting
INTEGER :: nind
!! The metric tensor
INTEGER :: nrr_tmp(nat)
!! Temporary array that contains the max number of WS vectors
!! for a pair of atoms.
INTEGER :: irvec_tmp(3,20*nq1*nq2*nq3,nat)
!! Temporary WS vectors for each atoms
INTEGER :: ndegen_tmp(20*nq1*nq2*nq3,nat)
!! Temporary WS vectors weigths for each atoms
!
REAL(kind=DP) :: adot(3,3)
!! Dot product between lattice vector
REAL(kind=DP) :: dist(125)
!! Contains the distance squared |r-R|^2
REAL(kind=DP) :: mindist
!! Minimum distance
REAL(kind=DP) :: tot, tot2
!! Sum of all the weigths
REAL(kind=DP) :: ndiff(3)
!! Distances. Must be now real because of fractional atoms
!
LOGICAL :: found
!! True if the vector has been found
!
nind = 20*nq1*nq2*nq3
IF (nind .lt. 125) THEN
nind = 125
ENDIF
ALLOCATE (ind(nind))
!
DO ipol = 1, 3
DO jpol = 1, 3
adot (ipol, jpol) = dot_product ( at(:,ipol), at(:,jpol) )
ENDDO
ENDDO
!
CALL cryst_to_cart(nat,tau(:,:),bg,-1)
!
! Loop over grid points r on a unit cell that is 8 times larger than a
! primitive supercell. In the end nrr contains the total number of grids
! points that have been found in the Wigner-Seitz cell
!
nrr_tmp(:) = 0
DO na=1, nat
DO n1 = -2*nq1, 2*nq1
DO n2 = -2*nq2, 2*nq2
DO n3 = -2*nq3, 2*nq3
!
! Loop over the 5^3 = 125 points R. R=0 corresponds to i1=i2=i3=2, or icnt=63
!
i = 0
dist(:) = 0.d0
DO i1 = -2, 2
DO i2 = -2, 2
DO i3 = -2, 2
i = i + 1
!
! Calculate distance squared |r-R|^2
!
ndiff(1) = n1 - i1*nq1 -tau(1,na)
ndiff(2) = n2 - i2*nq2 -tau(2,na)
ndiff(3) = n3 - i3*nq3 -tau(3,na)
DO ipol = 1, 3
DO jpol = 1, 3
dist(i) = dist(i) + dble(ndiff(ipol))*adot(ipol,jpol)*dble(ndiff(jpol))
ENDDO
ENDDO
!
ENDDO
ENDDO
ENDDO
!
! Sort the 125 vectors R by increasing value of |r-R|^2
ind(1) = 0 ! required for hpsort_eps (see the subroutine)
CALL hpsort_eps_epw( 125, dist, ind, eps8)
!
! Find all the vectors R with the (same) smallest |r-R|^2;
! if R=0 is one of them, then the current point r belongs to
! Wignez-Seitz cell => set found to true
!
found = .false.
i = 1
mindist = dist(1)
DO WHILE ( abs(dist(i)-mindist) < eps8 .and. i < 125 )
IF (ind(i) == 63) found = .true.
i = i + 1
ENDDO
!
IF (found) THEN
nrr_tmp(na) = nrr_tmp(na) + 1
ndegen_tmp (nrr_tmp(na), na) = i - 1
irvec_tmp (:, nrr_tmp(na), na) = (/ n1, n2, n3 /)
ENDIF
!
ENDDO ! n3
ENDDO ! n2
ENDDO ! n3
ENDDO ! na
!
! Now creates a global set of WS vectors from all the atoms pair.
! Also remove the duplicated ones.
nrr_g = nrr_tmp(1)
irvec_gg(:,:) = irvec_tmp(:, :, 1)
DO na = 1, nat
DO ir=1, nrr_tmp(na)
found = .false.
DO irtot = 1, nrr_g
IF ( ALL(irvec_tmp(:, ir, na) == irvec_gg(:, irtot)) ) THEN
found = .true.
ENDIF
ENDDO !nrr
IF(.not. found) THEN
nrr_g = nrr_g + 1
irvec_gg(:, nrr_g) = irvec_tmp(:, ir, na)
ENDIF
ENDDO ! ir
ENDDO ! na
!
! Creates a pair of atoms depended degeneracy array but with a number of WS
! vectors per pair that is equal to the global set. Populate with zero weights
! the one that are not part of that pair set.
ndegen_gg(:,:) = 0
DO na = 1, nat
DO ir=1, nrr_tmp(na)
DO irtot = 1, nrr_g
IF ( ALL(irvec_gg(:, irtot) == irvec_tmp(:, ir, na)) ) THEN
ndegen_gg(irtot, na) = ndegen_tmp(ir, na)
ENDIF
ENDDO
ENDDO
ENDDO
!
DO na = 1, nat
tot = 0.d0
tot2 = 0.d0
DO i = 1, nrr_g
IF (ndegen_gg(i, na) > 0 ) THEN
tot2 = tot2 + 1.d0 / dble(ndegen_gg(i, na))
ENDIF
ENDDO
DO i = 1, nrr_tmp(na)
tot = tot + 1.d0 / dble(ndegen_tmp(i, na))
ENDDO
!
IF(abs(tot-dble(nq1*nq2*nq3)) > eps8) call errore &
('wigner_seitzg',' weights do not add up to nq1*nq2*nq3',1)
IF(abs(tot-tot2) > eps8) call errore &
('wigner_seitzg',' weigths of pair of atoms is not equal to global weights',1)
ENDDO
!
IF(nrr_g > 20*nq1*nq2*nq3) call errore &
('wigner_seitzq','too many WS points, try to increase the bound 20*nq1*nq2*nq3',1)
!
! Now we compute the WS distances
wslen_gg(:) = 0.d0
DO i = 1, nrr_g
DO ipol = 1, 3
DO jpol = 1, 3
wslen_gg(i) = wslen_gg(i) + dble(irvec_gg(ipol,i))*adot(ipol,jpol)*dble(irvec_gg(jpol,i))
ENDDO
ENDDO
wslen_gg(i) = sqrt(wslen_gg(i))
ENDDO
!
CALL cryst_to_cart(nat,tau(:,:),at,1)
!
DEALLOCATE(ind)
!
! -----------------------------------------------------------------------------------------
END SUBROUTINE wigner_seitzg
! -----------------------------------------------------------------------------------------
!
END MODULE wigner

212
EPW/src/wigner_seitz.f90
View File

@ -1,212 +0,0 @@
!
! Copyright (C) 2010-2016 Samuel Ponce', Roxana Margine, Carla Verdi, Feliciano Giustino
! Copyright (C) 2007-2009 Jesse Noffsinger, Brad Malone, Feliciano Giustino
!
! This file is distributed under the terms of the GNU General Public
! License. See the file `LICENSE' in the root directory of the
! present distribution, or http://www.gnu.org/copyleft.gpl.txt .
!
! Adapted from the original f77 Wannier code of Marzari, Vanderbilt,
! and Souza
!
!-----------------------------------------------------------------
subroutine wigner_seitz (nk1, nk2, nk3, irvec, nrr, ndegen, wslen)
!-----------------------------------------------------------------
!!
!! Calculates a grid of points that fall inside of (and eventually
!! on the surface of) the Wigner-Seitz supercell centered on the
!! origin of the Bravais lattice with primitive translations
!! nk1*a_1+nk2*a_2+nk3*a_3
!!
!!
!! w.r.t. the original version the wigner-seitz vectors are sorted
!! by increasing lenght in output. In this way the electron and
!! phonon wigner-seitz vectors should always be the same (even though
!! the number of them may differ)
!!
!! BUG FIX: in the case of the tetragonal cell with c>a (LSCO)
!! the WS points are correctly determined, but there are a few points
!! with the wrong degeneracies. To avoid this I search for points
!! in the -2:2 replicas (5^2 replicas). I had the same problem in
!! createkmap for the g0vec shift, and also there I have fixed it
!! by extending the replicas to -2:2 instead of -1:1. FG May 07
!!
!! 02/2017, SP: I modified the routine to be closer to the current
!! Wannier90 routine hamiltonian.F90
!-----------------------------------------------------------------
USE kinds, ONLY : DP
USE cell_base, ONLY : at
!
implicit none
!
INTEGER, INTENT (in) :: nk1
!! size of the uniform k mesh
INTEGER, INTENT (in) :: nk2
!! size of the uniform k mesh
INTEGER, INTENT (in) :: nk3
!! size of the uniform k mesh
INTEGER, INTENT (out) :: irvec(3,20*nk1*nk2*nk3)
!! integer components of the ir-th Wigner-Seitz grid point in the basis of the lattice vectors
INTEGER, INTENT (out) :: ndegen(20*nk1*nk2*nk3)
!! Number of degeneracies
INTEGER, INTENT (out) :: nrr
!! number of Wigner-Seitz grid points
REAL(kind=DP), INTENT (out) :: wslen(20*nk1*nk2*nk3)
!! real-space length, in units of alat
!
! work variables
integer :: irvec_ (3,20*nk1*nk2*nk3), ndegen_ (20*nk1*nk2*nk3)
real(kind=DP), parameter :: eps = 1.d-8
integer :: n1, n2, n3, i1, i2, i3, i, ipol, jpol, ndiff(3)!, ind(2*nk1*nk2*nk3)
integer, allocatable :: ind(:)
real(kind=DP) :: tot, mindist, adot(3,3), dist(125)
logical :: found
!
! the metric tensor
!
INTEGER :: nind
!
nind = 20*nk1*nk2*nk3
IF (nind .lt. 125) then
nind = 125
ENDIF
IF (allocated(ind)) deallocate (ind)
allocate (ind(nind))
!
DO ipol = 1, 3
DO jpol = 1, 3
adot (ipol, jpol) = dot_product ( at(:,ipol), at(:,jpol) )
ENDDO
ENDDO
!
! Loop over grid points r on a unit cell that is 8 times larger than a
! primitive supercell. In the end nrr contains the total number of grids
! points that have been found in the Wigner-Seitz cell
!
! SP: I modified the routine to be closer to Wannier90.
! If problem happens, try reverting to old (more time consuming but could
! be more accurate in weird shaped primitive cells)
nrr = 0
!DO n1 = -nk1, nk1
! DO n2 = -nk2, nk2
! DO n3 = -nk3, nk3
DO n1 = 0, 4*nk1
DO n2 = 0, 4*nk2
DO n3 = 0, 4*nk3
!
! Loop over the 5^3 = 125 points R. R=0 corresponds to i1=i2=i3=2, or icnt=63
!
i = 0
!DO i1 = -2, 2
! DO i2 = -2, 2
! DO i3 = -2, 2
DO i1 = 0, 4
DO i2 = 0, 4
DO i3 = 0, 4
i = i + 1
!
! Calculate distance squared |r-R|^2
!
ndiff(1) = n1 - i1*nk1
ndiff(2) = n2 - i2*nk2
ndiff(3) = n3 - i3*nk3
dist(i) = 0.d0
DO ipol = 1, 3
DO jpol = 1, 3
!dist(i) = dist(i) + real(ndiff(ipol),kind=DP)*adot(ipol,jpol)*real(ndiff(jpol),kind=DP)
dist(i) = dist(i) + dble(ndiff(ipol))*adot(ipol,jpol)*dble(ndiff(jpol))
ENDDO
ENDDO
!
ENDDO
ENDDO
ENDDO
!
! Sort the 125 vectors R by increasing value of |r-R|^2
ind(1) = 0 ! required for hpsort_eps (see the subroutine)
CALL hpsort_eps_epw( 125, dist, ind, eps)
!
! Find all the vectors R with the (same) smallest |r-R|^2;
! if R=0 is one of them, then the current point r belongs to
! Wignez-Seitz cell => set found to true
!
found = .false.
i = 1
mindist = dist(1)
DO while ( abs(dist(i)-mindist).le.eps .and. i.lt.125 )
IF (ind(i).eq.63) found = .true.
i = i + 1
ENDDO
!
IF (found) then
nrr = nrr + 1
ndegen (nrr) = i - 1
irvec (1, nrr) = n1 - 2*nk1
irvec (2, nrr) = n2 - 2*nk2
irvec (3, nrr) = n3 - 2*nk3
ENDIF
!
!mindist=minval(dist)
!if (abs(dist(63) - mindist ) .lt. eps ) then
! nrr = nrr + 1
! ndegen(nrr)=0
! do i=1,125
! if (abs (dist (i) - mindist) .lt. eps ) ndegen(nrr)=ndegen(nrr)+1
! end do
! irvec(1, nrr) = n1
! irvec(2, nrr) = n2
! irvec(3, nrr) = n3
! !
!end if
!
ENDDO
ENDDO
ENDDO
!
! Check the "sum rule"
!
tot = 0.d0
DO i = 1, nrr
tot = tot + 1.d0 / dble (ndegen(i))
ENDDO
!
IF(abs(tot-dble(nk1*nk2*nk3)).gt.eps) call errore &
('wigner_seitz','weights do not add up to nk1*nk2*nk3',1)
!
!@ JN it happens in 2d and 1d systems with small course grids. I've changed to 20**3
!@ could calculate the max number of elements at the beginning
! Hopefully this will never happen, i.e., I think 2*nk1*nk2*nk3 is
! an upper bound to the number of lattice points in (or on
! the surface of) the Wigner-Seitz supercell
!
IF(nrr.gt.20*nk1*nk2*nk3) call errore &
('wigner_seitz','too many wigseit points, try to increase the bound 20*nk1*nk2*nk3',1)
!
! Now sort the wigner-seitz vectors by increasing magnitude
!
DO i = 1, nrr
wslen(i) = 0.d0
DO ipol = 1, 3
DO jpol = 1, 3
wslen(i) = wslen(i) + dble(irvec(ipol,i))*adot(ipol,jpol)*dble(irvec(jpol,i))
ENDDO
ENDDO
wslen(i) = sqrt(wslen(i))
ENDDO
!
ind(1) = 0 ! required for hpsort_eps (see the subroutine)
CALL hpsort_eps( nrr, wslen, ind, eps)
!
! now wslen is already sorted, but we still have to sort
! irvec and ndegen
!
DO i = 1, nrr
ndegen_ (i) = ndegen(ind(i))
irvec_ (:,i) = irvec(:,ind(i))
ENDDO
ndegen = ndegen_
irvec = irvec_
!
end subroutine wigner_seitz

117
EPW/src/wigner_seitz2.f90
View File

@ -1,117 +0,0 @@
!
! Copyright (C) 2010-2016 Samuel Ponce', Roxana Margine, Carla Verdi, Feliciano Giustino
! Copyright (C) 2007-2009 Jesse Noffsinger, Brad Malone, Feliciano Giustino
!
! This file is distributed under the terms of the GNU General Public
! License. See the file `LICENSE' in the root directory of the
! present distribution, or http://www.gnu.org/copyleft.gpl.txt .
!
!-----------------------------------------------------------------
subroutine wigner_seitz2 (nk1, nk2, nk3, nq1, nq2, nq3,&
nrr_k, nrr_q, irvec, wslen, ndegen_k, ndegen_q)
!-----------------------------------------------------------------
!!
!! We have nk1*nk2*nk3 electron points and nq1*nq2*nq3 phonon points
!! on the same grid. Assuming nq_i <= nk_i, i=1..3 we sort the corresponding
!! wigner-seitz points in such a way that a subset 1...nrr_q < nrr_k gives
!! the WS points for the phonons, while the full set 1..nrr_k gives the
!! WS points for the electrons
!!
!! the unsorted electron and phonon grids are obtained by calling
!! wigner_seitz.f90
!!
!! On exit, we have the same irvec for electrons and phonons, but
!! different ndegen
!!
!-----------------------------------------------------------------
USE kinds, ONLY : DP
!
implicit none
!
INTEGER, INTENT (in) :: nk1
!! size of the uniform k mesh
INTEGER, INTENT (in) :: nk2
!! size of the uniform k mesh
INTEGER, INTENT (in) :: nk3
!! size of the uniform k mesh
INTEGER, INTENT (in) :: nq1
!! size of the uniform q mesh
INTEGER, INTENT (in) :: nq2
!! size of the uniform q mesh
INTEGER, INTENT (in) :: nq3
!! size of the uniform q mesh
INTEGER, INTENT (out) :: irvec(3,20*nk1*nk2*nk3)
!! integer components of the ir-th Wigner-Seitz grid point in the basis of the lattice vectors
INTEGER, INTENT (out) :: nrr_k
!! number of Wigner-Seitz grid points for electrons
INTEGER, INTENT (out) :: nrr_q
!! number of Wigner-Seitz grid points for electrons
!
REAL(kind=DP), INTENT (out) :: wslen(2*nk1*nk2*nk3)
!! real-space length, in units of alat
!
! work variables
INTEGER :: irvec_k (3,20*nk1*nk2*nk3), ndegen_k (20*nk1*nk2*nk3), &
irvec_q (3,20*nk1*nk2*nk3), ndegen_q (20*nk1*nk2*nk3), &
ind2(20*nk1*nk2*nk3), ire, ir, nind
INTEGER, ALLOCATABLE :: ind(:)
REAL(kind=DP) :: wslen_k (20*nk1*nk2*nk3), wslen_q (20*nk1*nk2*nk3)
!
! The allocation of the sorting arrays is not very clean. However,
! for the moment it works.
!
nind = 20*nk1*nk2*nk3
IF (nind .lt. 125) then
nind = 125
ENDIF
ALLOCATE (ind(nind))
!
! initialization for ihpsort (not to be removed!)
!
ind = 0
ind2 = 0
!
! check the bounds
!
IF ( nq1.gt.nk1 .or. nq2.gt.nk2 .or. nq3.gt.nk3 ) call errore &
('wigner_seitz2','phonon grid should be smaller than electron grid',1)
!
! now generated the un-sorted points for both electrons and phonons
!
CALL wigner_seitz ( nk1, nk2, nk3, irvec_k, nrr_k, ndegen_k, wslen_k)
CALL wigner_seitz ( nq1, nq2, nq3, irvec_q, nrr_q, ndegen_q, wslen_q)
!
! loop on phonon points and find the match in the corresponding electronic list
!
DO ir = 1, nrr_k
!
! fake index which is useful in ihpsort
! (I need to split the two subsets 1...nrr_q and nrr_q+1...nrr_k:
! here below ind() will be between 1 and nrr_q. Therefore, if ind()
! is larger than nrr_q it must belong to the second subset)
!
ind (ir) = nrr_q + ir
ENDDO
DO ire = 1, nrr_k
ir = 1
DO while ( ( irvec_k(1,ire).ne.irvec_q(1,ir) .or. &
irvec_k(2,ire).ne.irvec_q(2,ir) .or. &
irvec_k(3,ire).ne.irvec_q(3,ir) ) .and. &
ir.le.nrr_q )
ir = ir + 1
ENDDO
IF ( ir.le.nrr_q ) ind (ire) = ir - 1
ENDDO
!
! Sort the electronic points accordingly
!
CALL ihpsort ( nrr_k, ind, ind2 )
!
DO ir = 1, nrr_k
irvec (:, ir) = irvec_k (:, ind2(ir) )
ndegen_k (ir) = ndegen_k ( ind2(ir) )
wslen (ir) = wslen_k ( ind2(ir) )
ENDDO
!
end subroutine wigner_seitz2

2
PHonon/PH/rigid.f90
View File

@ -6,7 +6,7 @@
! or http://www.gnu.org/copyleft/gpl.txt .
!
MODULE rigid
PUBLIC :: rgd_blk, dyndiag, nonanal, nonanal_ifc
PUBLIC :: rgd_blk, dyndiag, nonanal, nonanal_ifc, cdiagh2
PRIVATE
CONTAINS
!