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

Restart with arbitray number of cores after epmatwp has been written to file 

Some Ford doc added. 
Some unused variables cleaned. 



git-svn-id: http://qeforge.qe-forge.org/svn/q-e/trunk/espresso@12887 c92efa57-630b-4861-b058-cf58834340f0
This commit is contained in:
sponce 2016-08-30 17:23:01 +00:00
parent 95cf2364ef
commit cbd08e75df
12 changed files with 512 additions and 403 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

11
EPW/src/allocate_epwq.f90
View File

@ -26,14 +26,13 @@
USE spin_orb, ONLY : lspinorb
USE control_lr, ONLY : lgamma, nbnd_occ
USE phcom, ONLY : evq, dvpsi, dpsi, vlocq,&
dmuxc, w2, dyn00, t, tmq,&
npertx
dmuxc, w2, npertx
USE phus, ONLY : int1, int1_nc, int2, int2_so, &
int4, int4_nc, int5, int5_so, becsum_nc, &
alphasum, alphasum_nc, alphap
USE lr_symm_base, ONLY : rtau
USE qpoint, ONLY : eigqts
USE lrus, ONLY : becp1, int3, int3_nc, dpqq, dpqq_so
USE lrus, ONLY : becp1, int3, int3_nc
USE elph2, ONLY : elph, el_ph_mat
USE becmod, ONLY : becp, allocate_bec_type
USE uspp_param, ONLY : nhm
@ -81,9 +80,7 @@
ALLOCATE (rtau ( 3, 48, nat))
ALLOCATE (u ( 3 * nat, 3 * nat))
ALLOCATE (w2 ( 3 * nat))
ALLOCATE (dyn00 ( 3 * nat, 3 * nat))
ALLOCATE (t (npertx, npertx, 48,3 * nat))
ALLOCATE (tmq (npertx, npertx, 3 * nat))
! ALLOCATE (t (npertx, npertx, 48,3 * nat))
allocate (name_rap_mode( 3 * nat))
allocate (num_rap_mode( 3 * nat ))
ALLOCATE (npert ( 3 * nat))
@ -93,7 +90,6 @@
ALLOCATE (int3 ( nhm, nhm, npertx, nat, nspin))
ALLOCATE (int4 ( nhm * (nhm + 1)/2, 3, 3, nat, nspin))
ALLOCATE (int5 ( nhm * (nhm + 1)/2, 3, 3, nat , nat))
ALLOCATE (dpqq( nhm, nhm, 3, ntyp))
IF (noncolin) THEN
ALLOCATE(int1_nc( nhm, nhm, 3, nat, nspin))
ALLOCATE(int3_nc( nhm, nhm, npertx, nat, nspin))
@ -103,7 +99,6 @@
IF (lspinorb) THEN
ALLOCATE(int2_so( nhm, nhm, 3, nat, nat, nspin))
ALLOCATE(int5_so( nhm, nhm, 3, 3, nat, nat, nspin))
ALLOCATE(dpqq_so( nhm, nhm, nspin, 3, ntyp))
END IF
END IF
ALLOCATE (alphasum ( nhm * (nhm + 1)/2, 3, nat , nspin))

41
EPW/src/deallocate_epw.f90
View File

@ -23,10 +23,9 @@
!----------------------------------------------------------------------
USE phcom, ONLY : alphap, alphasum, alphasum_nc, &
becsum_nc, dmuxc, dpsi,&
drc, dpsi, dyn, dyn00, evq, w2, dvpsi,&
drc, dpsi, dyn, evq, w2, dvpsi,&
int5, vlocq, int2_so, int5_so
USE gc_lr, ONLY : grho, dvxc_rr, dvxc_s, dvxc_ss, dvxc_sr
USE lrus, ONLY : becp1, int3, int3_nc, dpqq, dpqq_so
USE lrus, ONLY : becp1, int3, int3_nc
USE phus, ONLY : int1, int1_nc, int2, int4, int4_nc
USE lr_symm_base, ONLY : rtau
USE noncollin_module, ONLY : m_loc
@ -35,12 +34,33 @@
USE elph2, ONLY : el_ph_mat, epf17, epsi, etf,&
etq, et_all, wf, wkf, wqf, wslen,&
xkq, xk_all, zstar, xkf, xqf, epmatwp
USE modes, ONLY : tmq, t, npert, u, name_rap_mode, num_rap_mode
USE epwcom, ONLY : epbread, epwread
USE modes, ONLY : t, npert, u, name_rap_mode, num_rap_mode
USE qpoint, ONLY : eigqts, igkq
USE klist, ONLY : nks
!
IMPLICIT NONE
INTEGER :: ik, ipol
!
IF ( epwread .and. .not. epbread ) THEN
! EPW variables ONLY
!
IF(ALLOCATED(el_ph_mat)) DEALLOCATE (el_ph_mat)
IF(ALLOCATED(epmatwp)) DEALLOCATE (epmatwp)
IF(ALLOCATED(epf17)) DEALLOCATE (epf17)
IF(ALLOCATED(etq)) DEALLOCATE (etq)
IF(ALLOCATED(etf)) DEALLOCATE (etf)
IF(ALLOCATED(wf)) DEALLOCATE (wf)
IF(ALLOCATED(xkq)) DEALLOCATE (xkq)
IF(ALLOCATED(xkf)) DEALLOCATE (xkf)
IF(ALLOCATED(wkf)) DEALLOCATE (wkf)
IF(ALLOCATED(xqf)) DEALLOCATE (xqf)
IF(ALLOCATED(wqf)) DEALLOCATE (wqf)
IF(ALLOCATED(xk_all)) DEALLOCATE (xk_all)
IF(ALLOCATED(et_all)) DEALLOCATE (et_all)
IF(ALLOCATED(wslen)) DEALLOCATE (wslen)
!
ELSE
!
IF (lgamma) THEN
IF(ASSOCIATED(evq)) NULLIFY(evq)
@ -62,10 +82,8 @@
if(allocated(name_rap_mode)) deallocate (name_rap_mode)
if(allocated(num_rap_mode)) deallocate (num_rap_mode)
IF(ALLOCATED(dyn)) DEALLOCATE (dyn)
IF(ALLOCATED(dyn00)) DEALLOCATE (dyn00)
IF(ALLOCATED(w2)) DEALLOCATE (w2)
IF(ASSOCIATED(t)) DEALLOCATE (t)
IF(ASSOCIATED(tmq)) DEALLOCATE (tmq)
!IF(ASSOCIATED(t)) DEALLOCATE (t)
IF(ALLOCATED(epsi)) DEALLOCATE (epsi)
IF(ALLOCATED(zstar)) DEALLOCATE (zstar)
!
@ -76,7 +94,6 @@
IF(ALLOCATED(int3)) DEALLOCATE (int3)
IF(ALLOCATED(int4)) DEALLOCATE (int4)
IF(ALLOCATED(int5)) DEALLOCATE (int5)
IF(ALLOCATED(dpqq)) DEALLOCATE (dpqq)
IF(ALLOCATED(int1_nc)) DEALLOCATE(int1_nc)
IF(ALLOCATED(int3_nc)) DEALLOCATE(int3_nc)
IF(ALLOCATED(int4_nc)) DEALLOCATE(int4_nc)
@ -84,7 +101,6 @@
IF(ALLOCATED(alphasum_nc)) DEALLOCATE(alphasum_nc)
IF(ALLOCATED(int2_so)) DEALLOCATE(int2_so)
IF(ALLOCATED(int5_so)) DEALLOCATE(int5_so)
IF(ALLOCATED(dpqq_so)) DEALLOCATE(dpqq_so)
IF(ALLOCATED(alphasum)) DEALLOCATE (alphasum)
!
if(allocated(alphap)) then
@ -108,12 +124,6 @@
!
IF(ALLOCATED(drc)) DEALLOCATE(drc)
!
IF(ALLOCATED(dvxc_rr)) DEALLOCATE (dvxc_rr)
IF(ALLOCATED(dvxc_sr)) DEALLOCATE (dvxc_sr)
IF(ALLOCATED(dvxc_ss)) DEALLOCATE (dvxc_ss)
IF(ALLOCATED(dvxc_s)) DEALLOCATE (dvxc_s)
IF(ALLOCATED(grho)) DEALLOCATE (grho)
!
! EPW variables
!
IF(ALLOCATED(el_ph_mat)) DEALLOCATE (el_ph_mat)
@ -130,5 +140,6 @@
IF(ALLOCATED(xk_all)) DEALLOCATE (xk_all)
IF(ALLOCATED(et_all)) DEALLOCATE (et_all)
IF(ALLOCATED(wslen)) DEALLOCATE (wslen)
ENDIF ! epwread .and. .not. epbread
!
END SUBROUTINE deallocate_epw

41
EPW/src/dynwan2bloch.f90
View File

@ -7,7 +7,7 @@
! present distribution, or http://www.gnu.org/copyleft.gpl.txt .
!
!--------------------------------------------------------------------------
SUBROUTINE dynwan2bloch ( nbnd, nrr, irvec, ndegen, xq, cuf, eig)
SUBROUTINE dynwan2bloch ( nbnd, nrr, irvec, ndegen, xxq, cuf, eig)
!--------------------------------------------------------------------------
!!
!!
@ -40,37 +40,33 @@
USE elph2, ONLY : rdw, epsi, zstar
USE epwcom, ONLY : lpolar
USE constants_epw, ONLY : twopi, ci, czero
!
implicit none
!
! input variables
INTEGER, INTENT (in) :: nbnd
!! number of bands (possibly of the optimal subspace)
INTEGER, INTENT (in) :: nrr
!! kpoint number for the interpolation
INTEGER, INTENT (in) :: irvec (3, nrr)
!! record length and unit for direct write of rotation matrix
INTEGER, INTENT (in) :: ndegen (nrr)
!! number of WS points, crystal coordinates, degeneracy
!
integer :: nbnd, nrr, irvec (3, nrr), ndegen (nrr)
! number of bands (possibly of the optimal subspace)
! kpoint number for the interpolation
! record length and unit for direct write of rotation matrix
! number of WS points, crystal coordinates, degeneracy
!
! Hamiltonian in wannier basis
!
real(kind=DP) :: xq (3)
! kpoint coordinates for the interpolation
!
! output variables
!
real(kind=DP) :: eig (nbnd)
! interpolated hamiltonian eigenvalues for this kpoint
complex(kind=DP) :: cuf(nbnd, nbnd)
! Rotation matrix, fine mesh
REAL(kind=DP), INTENT (in) :: xxq (3)
!! kpoint coordinates for the interpolation
REAL(kind=DP), INTENT (out) :: eig (nbnd)
!! interpolated hamiltonian eigenvalues for this kpoint
COMPLEX(kind=DP), INTENT (out) :: cuf(nbnd, nbnd)
!! Rotation matrix, fine mesh
!
! work variables
! variables for lapack ZHPEVX
!
integer :: neig, info, ifail( nbnd ), iwork( 5*nbnd )
real(kind=DP) :: w( nbnd ), rwork( 7*nbnd )
complex(kind=DP) :: champ( nbnd*(nbnd+1)/2 ), &
cwork( 2*nbnd ), cz( nbnd, nbnd)
!
! work variables
!
real(kind=DP) :: xq (3)
complex(kind=DP) :: chf(nbnd, nbnd)
! Hamiltonian in Bloch basis, fine mesh
integer :: ibnd, jbnd, ir, na, nb
@ -88,6 +84,7 @@
! H~_k is chf ( nbnd, nbnd, 2*ik-1 )
! H~_k+q is chf ( nbnd, nbnd, 2*ik )
!
xq = xxq
chf (:,:) = czero
!
DO ir = 1, nrr

18
EPW/src/elphon_shuffle_wrap.f90
View File

@ -215,8 +215,14 @@
IF (allocated(qrad)) deallocate(qrad)
allocate (qrad (maxvalue, nbetam*(nbetam+1)/2,lmaxq, nsp))
ENDIF
!
! DO not perform the check if restart
IF ( epwread .and. .not. epbread ) then
continue
ELSE
IF (nkstot .ne. nk1*nk2*nk3 ) &
CALL errore('elphon_shuffle_wrap','nscf run inconsistent with epw input',1)
ENDIF
!
! READ in external electronic eigenvalues. e.g. GW
!
@ -251,8 +257,14 @@
et_mb(:,:) = et(:,1:nks)
ENDIF
!
! Do not recompute dipole matrix elements
IF ( epwread .and. .not. epbread ) then
continue
ELSE
! compute coarse grid dipole matrix elements. Very fast
CALL compute_pmn_para
ENDIF
!CALL compute_pmn_para
!
! gather electronic eigenvalues for subsequent shuffle
!
@ -277,6 +289,11 @@
ENDIF
!
CALL mp_barrier(inter_pool_comm)
!
! Do not do symmetry stuff
IF ( epwread .and. .not. epbread ) then
CONTINUE
ELSE
!
! allocate dynamical matrix and ep matrix for all q's
!
@ -310,6 +327,7 @@
CALL sgam_ph_new (at, bg, nsym, s, irt, tau, rtau, nat)
CALL set_giq (xq,s,nsymq,nsym,irotmq,minus_q,gi,gimq)
ENDDO
ENDIF ! epwread .and. .not. epbread
!
! CV: if we read the .fmt files we don't need to read the .epb anymore
!

78
EPW/src/ephwann_shuffle.f90
View File

@ -26,9 +26,10 @@
!
USE kinds, ONLY : DP
USE pwcom, ONLY : nbnd, nks, nkstot, isk, &
et, xk, at, bg, ef, nelec
et, xk, ef, nelec
USE cell_base, ONLY : at, bg
USE start_k, ONLY : nk1, nk2, nk3
USE ions_base, ONLY : amass, ityp
USE ions_base, ONLY : nat, amass, ityp
USE phcom, ONLY : nq1, nq2, nq3, nmodes, w2
USE epwcom, ONLY : nbndsub, lrepmatf, fsthick, epwread, &
epwwrite, ngaussw, degaussw, lpolar, &
@ -42,7 +43,7 @@
USE io_files, ONLY : prefix, diropn
USE io_global, ONLY : stdout, ionode
USE io_epw, ONLY : lambda_phself, linewidth_phself, iunepmatf, &
iunepmatwe, iunepmatwp
iunepmatwe, iunepmatwp, crystal
USE elph2, ONLY : nrr_k, nrr_q, cu, cuq, lwin, lwinq, irvec, ndegen_k, ndegen_q, &
wslen, chw, chw_ks, cvmew, cdmew, rdw, epmatwp, epmatq, &
wf, etf, etf_k, etf_ks, xqf, xkf, wkf, &
@ -54,7 +55,7 @@
#endif
USE mp, ONLY : mp_barrier, mp_bcast, mp_sum
USE io_global, ONLY : ionode_id
USE mp_global, ONLY : inter_pool_comm
USE mp_global, ONLY : inter_pool_comm, intra_pool_comm, root_pool
USE mp_world, ONLY : mpime
!
implicit none
@ -157,6 +158,45 @@
!
! DBSP
! HERE loadkmesh
IF ( epwread ) THEN
!
! We need some crystal info
IF (mpime.eq.ionode_id) THEN
!
OPEN(unit=crystal,file='crystal.fmt',status='old',iostat=ios)
READ (crystal,*) nat
READ (crystal,*) nmodes
READ (crystal,*) nelec
READ (crystal,*) at
READ (crystal,*) bg
READ (crystal,*) amass
ALLOCATE( ityp( nat ) )
READ (crystal,*) ityp
!
ENDIF
CALL mp_bcast (nat, ionode_id, inter_pool_comm)
CALL mp_bcast (nat, root_pool, intra_pool_comm)
IF (mpime /= ionode_id) ALLOCATE( ityp( nat ) )
CALL mp_bcast (nmodes, ionode_id, inter_pool_comm)
CALL mp_bcast (nmodes, root_pool, intra_pool_comm)
CALL mp_bcast (nelec, ionode_id, inter_pool_comm)
CALL mp_bcast (nelec, root_pool, intra_pool_comm)
CALL mp_bcast (at, ionode_id, inter_pool_comm)
CALL mp_bcast (at, root_pool, intra_pool_comm)
CALL mp_bcast (bg, ionode_id, inter_pool_comm)
CALL mp_bcast (bg, root_pool, intra_pool_comm)
CALL mp_bcast (amass, ionode_id, inter_pool_comm)
CALL mp_bcast (amass, root_pool, intra_pool_comm)
CALL mp_bcast (ityp, ionode_id, inter_pool_comm)
CALL mp_bcast (ityp, root_pool, intra_pool_comm)
IF (mpime.eq.ionode_id) THEN
CLOSE(crystal)
ENDIF
CALL mp_barrier(inter_pool_comm)
!
ELSE
continue
ENDIF
!
! determine Wigner-Seitz points
!
@ -182,6 +222,7 @@
CALL epw_read
!
ELSE !if not epwread (i.e. need to calculate fmt file)
!
!
xxq = 0.d0
CALL loadumat &
@ -970,11 +1011,14 @@ SUBROUTINE epw_write
!
USE kinds, ONLY : DP
USE epwcom, ONLY : nbndsub, vme, eig_read, etf_mem
USE pwcom, ONLY : ef
USE pwcom, ONLY : ef, nelec
USE elph2, ONLY : nrr_k, nrr_q, chw, rdw, cdmew, cvmew, chw_ks, &
zstar, epsi, epmatwp
USE ions_base, ONLY : amass, ityp, nat
USE cell_base, ONLY : at, bg
USE phcom, ONLY : nmodes
USE io_epw, ONLY : epwdata, iundmedata, iunvmedata, iunksdata, iunepmatwp
USE io_epw, ONLY : epwdata, iundmedata, iunvmedata, iunksdata, iunepmatwp, &
crystal
USE io_files, ONLY : prefix, diropn
USE mp, ONLY : mp_barrier
USE mp_global, ONLY : inter_pool_comm
@ -992,9 +1036,17 @@ SUBROUTINE epw_write
IF (mpime.eq.ionode_id) THEN
!
OPEN(unit=epwdata,file='epwdata.fmt')
OPEN(unit=crystal,file='crystal.fmt')
OPEN(unit=iundmedata,file='dmedata.fmt')
IF (vme) OPEN(unit=iunvmedata,file='vmedata.fmt')
IF (eig_read) OPEN(unit=iunksdata,file='ksdata.fmt')
WRITE (crystal,*) nat
WRITE (crystal,*) nmodes
WRITE (crystal,*) nelec
WRITE (crystal,*) at
WRITE (crystal,*) bg
WRITE (crystal,*) amass
WRITE (crystal,*) ityp
WRITE (epwdata,*) ef
WRITE (epwdata,*) nbndsub, nrr_k, nmodes, nrr_q
WRITE (epwdata,*) zstar, epsi
@ -1036,7 +1088,7 @@ SUBROUTINE epw_write
ENDDO
ENDDO
!DBSP
!filint = trim(prefix)//'.epmatwp'
filint = trim(prefix)//'.epmatwp'
CALL diropn (iunepmatwp, 'epmatwp', lrepmatw, exst)
CALL davcio ( aux, lrepmatw, iunepmatwp, 1, +1 )
CLOSE(iunepmatwp)
@ -1044,6 +1096,7 @@ SUBROUTINE epw_write
ENDIF
!
CLOSE(epwdata)
CLOSE(crystal)
CLOSE(iundmedata)
IF (vme) CLOSE(iunvmedata)
IF (eig_read) CLOSE(iunksdata)
@ -1061,6 +1114,7 @@ SUBROUTINE epw_read()
USE pwcom, ONLY : ef
USE elph2, ONLY : nrr_k, nrr_q, chw, rdw, epmatwp, &
cdmew, cvmew, chw_ks, zstar, epsi
USE ions_base, ONLY : nat
USE phcom, ONLY : nmodes
USE io_global, ONLY : stdout
USE io_files, ONLY : prefix, diropn
@ -1084,6 +1138,11 @@ SUBROUTINE epw_read()
!
WRITE(stdout,'(/5x,"Reading Hamiltonian, Dynamical matrix and EP vertex in Wann rep from file"/)')
call flush(6)
!
! This is important in restart mode as zstar etc has not been allocated
IF (.NOT. ALLOCATED (zstar) ) ALLOCATE( zstar(3,3,nat) )
IF (.NOT. ALLOCATED (epsi) ) ALLOCATE( epsi(3,3) )
IF (mpime.eq.ionode_id) THEN
!
OPEN(unit=epwdata,file='epwdata.fmt',status='old',iostat=ios)
@ -1170,8 +1229,13 @@ SUBROUTINE epw_read()
!
lrepmatw = 2 * nbndsub * nbndsub * nrr_k * nmodes * nrr_q
filint = trim(prefix)//'.epmatwp'
!CALL diropn (iunepmatwp, filint, lrepmatw, exst)
CALL diropn (iunepmatwp, 'epmatwp', lrepmatw, exst)
CALL davcio ( aux, lrepmatw, iunepmatwp, 1, -1 )
!READ( UNIT = iunepmatwp, REC = 1, IOSTAT = ios ) aux
i = 0
DO irq = 1, nrr_q
DO imode = 1, nmodes

22
EPW/src/epw.f90
View File

@ -27,7 +27,7 @@
USE control_flags, ONLY : gamma_only
USE control_epw, ONLY : wannierize
USE global_version, ONLY : version_number
USE epwcom, ONLY : filukk, eliashberg, ep_coupling
USE epwcom, ONLY : filukk, eliashberg, ep_coupling, epwread, epbread
USE environment, ONLY : environment_start
USE elph2, ONLY : elph
! Flag to perform an electron-phonon calculation. If .true.
@ -90,7 +90,18 @@ write(stdout,'(a)') "
CALL epw_readin
!
CALL allocate_epwq
IF ( epwread .AND. .NOT. epbread ) THEN
write(stdout,'(a)') " "
write(stdout,'(a)') " ------------------------------------------------------------------------ "
write(stdout,'(a)') " RESTART - RESTART - RESTART - RESTART "
write(stdout,'(a)') " Restart is done without reading PWSCF save file. "
write(stdout,'(a)') " Be aware that some consistency checks are therefore not done. "
write(stdout,'(a)') " ------------------------------------------------------------------------ "
write(stdout,'(a)') " "
ELSE
CALL epw_setup
ENDIF
!
! Print run info to stdout
!
@ -98,11 +109,20 @@ write(stdout,'(a)') "
!
IF ( ep_coupling ) THEN
!
! In case of restart with arbitrary number of cores.
IF ( epwread .and. .not. epbread ) THEN
continue
ELSE
CALL openfilepw
ENDIF
!
CALL print_clock( 'EPW' )
!
IF ( epwread .and. .not. epbread ) THEN
continue
ELSE
CALL epw_init(.true.)
ENDIF
!
CALL print_clock( 'EPW' )
!

4
EPW/src/epw_readin.f90
View File

@ -593,7 +593,11 @@
modenum_aux = modenum
!
! SP: This initialized nspin and nspin_mag
IF ( epwread .and. .not. epbread ) then
continue
ELSE
CALL read_file
ENDIF
!
! nbnd comes out of readfile
IF (nbndsub.eq.0) nbndsub = nbnd

3
EPW/src/io_epw.f90
View File

@ -28,7 +28,7 @@
PUBLIC :: epwdata, iundmedata, iunvmedata, iunksdata, iudyn, iukgmap, iuepb,&
iunepmatf, iurecover, iufilfreq, iufilegnv, iufileph, iufilkqmap, &
iufilikmap, iueig, iunepmatwp, iunepmatwe, iunkf, iunqf, &
iufileig, iukmap
iufileig, iukmap, crystal
PUBLIC :: iuwinfil, iun_plot, iuukk, iuprojfil !, iummn
!
! Output of physically relevant quantities (60-100)
@ -93,6 +93,7 @@
INTEGER :: iunqf = 118 ! The unit with the fine q-point mesh in crystal coord.
INTEGER :: iufileig = 119 ! The unit with eigenenergies [band.eig]
INTEGER :: iukmap = 120 ! Unit for the k-point map generation
INTEGER :: crystal = 121 ! Unit for crystal data
!
! Output quantites related to Wannier (201-250)

3
EPW/src/rigid_epw.f90
View File

@ -227,8 +227,7 @@ SUBROUTINE rgd_blk_epw(nq1,nq2,nq3,q,uq,epmat,nmodes,epsil,zeu,bmat,signe)
!
! work variables
!
real(DP) :: &
qeq, &! <q+G| epsil | q+G>
real(DP) :: qeq, &! <q+G| epsil | q+G>
arg, zaq, &
g1, g2, g3, gmax, alph, geg
integer :: na, ipol, im, m1,m2,m3, nrx1,nrx2,nrx3

2
EPW/src/sort.f90
View File

@ -28,7 +28,7 @@
!
! adapted from Numerical Recipes pg. 329 (new edition)
!
use kinds
use kinds, ONLY : DP
implicit none
!-input/output variables
integer, intent(in) :: n

83
EPW/src/wigner_seitz.f90
View File

@ -12,48 +12,52 @@
!-----------------------------------------------------------------
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
!
!!
!! 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
!!
!-----------------------------------------------------------------
USE kinds, only : DP
USE kinds, ONLY : DP
USE cell_base, ONLY : at
!
implicit none
!
integer :: nk1, nk2, nk3, irvec(3,20*nk1*nk2*nk3), ndegen(20*nk1*nk2*nk3), nrr
!@ integer :: nk1, nk2, nk3, irvec(3,2*nk1*nk2*nk3), ndegen(2*nk1*nk2*nk3), nrr
real(kind=DP) :: wslen(20*nk1*nk2*nk3)
!@ real(kind=DP) :: wslen(2*nk1*nk2*nk3)
!
! size of the uniform k mesh (input)
! integer components of the ir-th Wigner-Seitz grid point
! in the basis of the lattice vectors (output)
! number of Wigner-Seitz grid points (output)
! real-space length, in units of alat
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)
!@ integer :: irvec_ (3,2*nk1*nk2*nk3), ndegen_ (2*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(27)
real(kind=DP) :: tot, mindist, adot(3,3), dist(125)
logical :: found
!
@ -79,20 +83,13 @@
! points that have been found in the Wigner-Seitz cell
!
nrr = 0
! do n1 = 0, 2*nk1
! do n2 = 0, 2*nk2
! do n3 = 0, 2*nk3
DO n1 = 0, 4*nk1
DO n2 = 0, 4*nk2
DO n3 = 0, 4*nk3
!
! ! Loop over the 27 points R. R=0 corresponds to i1=i2=i3=1, or icnt=14
! Loop over the 5^3 = 125 points R. R=0 corresponds to i1=i2=i3=2, or icnt=63
!
i = 0
! do i1 = 0, 2
! do i2 = 0, 2
! do i3 = 0, 2
DO i1 = 0, 4
DO i2 = 0, 4
DO i3 = 0, 4
@ -114,7 +111,6 @@
ENDDO
ENDDO
!
! ! Sort the 27 vectors R by increasing value of |r-R|^2
! Sort the 125 vectors R by increasing value of |r-R|^2
!
! NOTA BENE: hpsort really sorts the dist vector
@ -123,8 +119,8 @@
! dist(i), while ind(i) is kept.
!
ind(1) = 0 ! required for hpsort_eps (see the subroutine)
! call hpsort_eps( 27, dist, ind, eps)
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
@ -134,7 +130,6 @@
i = 1
mindist = dist(1)
DO while ( abs(dist(i)-mindist).le.eps .and. i.lt.125 )
! if (ind(i).eq.14) found = .true.
IF (ind(i).eq.63) found = .true.
i = i + 1
ENDDO
@ -144,9 +139,6 @@
IF (found) then
nrr = nrr + 1
ndegen (nrr) = i - 1
! irvec (1, nrr) = n1 - nk1
! irvec (2, nrr) = n2 - nk2
! irvec (3, nrr) = n3 - nk3
irvec (1, nrr) = n1 - 2*nk1
irvec (2, nrr) = n2 - 2*nk2
irvec (3, nrr) = n3 - 2*nk3
@ -155,6 +147,7 @@
ENDDO
ENDDO
ENDDO
!
! Check the "sum rule"
!
tot = 0.d0

73
EPW/src/wigner_seitz2.f90
View File

@ -10,45 +10,52 @@
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
!
!!
!! 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
USE kinds, ONLY : DP
!
implicit none
!
integer :: nk1, nk2, nk3, nq1, nq2, nq3, &
irvec(3,20*nk1*nk2*nk3), nrr_k, nrr_q
!@ irvec(3,2*nk1*nk2*nk3), nrr_k, nrr_q
real(kind=DP) :: wslen(2*nk1*nk2*nk3)
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
!
! size of the uniform k mesh (input)
! size of the uniform q mesh (input)
! integer components of the ir-th Wigner-Seitz grid point
! in the basis of the lattice vectors (output)
! number of Wigner-Seitz grid points for electrons and for phonons (output)
! real-space length, in units of alat
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), &
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 :: irvec_k (3,2*nk1*nk2*nk3), ndegen_k (2*nk1*nk2*nk3), &
!@ irvec_q (3,2*nk1*nk2*nk3), ndegen_q (2*nk1*nk2*nk3), &
!@ ind2(2*nk1*nk2*nk3), ire, ir, nind
! ind(2*nk1*nk2*nk3), ind2(2*nk1*nk2*nk3), ire, ir
integer, allocatable :: ind(:)
real(kind=DP) :: wslen_k (20*nk1*nk2*nk3), wslen_q (20*nk1*nk2*nk3)
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.
@ -57,7 +64,7 @@
IF (nind .lt. 125) then
nind = 125
ENDIF
allocate (ind(nind))
ALLOCATE (ind(nind))
!
! initialization for ihpsort (not to be removed!)
!
@ -96,7 +103,7 @@
IF ( ir.le.nrr_q ) ind (ire) = ir - 1
ENDDO
!
! sort the electronic points accordingly
! Sort the electronic points accordingly
!
CALL ihpsort ( nrr_k, ind, ind2 )
!