mirror of https://gitlab.com/QEF/q-e.git
pw2wannier90: Merge SCDM for collinear and noncollinear. Pool implemented for both.
This commit is contained in:
Jae-Mo Lihm
parent
5ccf2b2b3f
commit
c824589e20
|
|
@ -611,11 +611,7 @@ PROGRAM pw2wannier90
|
|||
WRITE(stdout,*) ' *** Compute A with SCDM-k'
|
||||
WRITE(stdout,*) ' --------------------------'
|
||||
WRITE(stdout,*)
|
||||
if (noncolin) then
|
||||
CALL compute_amn_with_scdm_spinor
|
||||
else
|
||||
CALL compute_amn_with_scdm
|
||||
end if
|
||||
CALL compute_amn_with_scdm
|
||||
ELSE
|
||||
WRITE(stdout,*) ' --------------------------'
|
||||
WRITE(stdout,*) ' *** Compute A projections'
|
||||
|
|
@ -3778,7 +3774,7 @@ SUBROUTINE compute_amn_with_scdm
|
|||
USE wvfct, ONLY : nbnd, et, npwx
|
||||
USE gvecw, ONLY : gcutw
|
||||
USE control_flags, ONLY : gamma_only
|
||||
USE wavefunctions, ONLY : evc, psic
|
||||
USE wavefunctions, ONLY : evc, psic, psic_nc
|
||||
USE io_files, ONLY : nwordwfc, iunwfc
|
||||
USE wannier
|
||||
USE klist, ONLY : nkstot, xk, ngk, igk_k, nks
|
||||
|
|
@ -3800,14 +3796,19 @@ SUBROUTINE compute_amn_with_scdm
|
|||
qr_tau(:), cwork(:), cwork2(:), Umat(:,:), VTmat(:,:), Amat(:,:) ! vv: complex arrays for the SVD factorization
|
||||
COMPLEX(DP), ALLOCATABLE :: phase_g(:,:) ! jml
|
||||
REAL(DP), ALLOCATABLE :: rwork(:), rwork2(:), singval(:), rpos(:,:) ! vv: Real array for the QR factorization and SVD
|
||||
INTEGER, ALLOCATABLE :: piv(:) ! vv: Pivot array in the QR factorization
|
||||
INTEGER, ALLOCATABLE :: piv(:)
|
||||
!! vv: Pivot array in the QR factorization
|
||||
INTEGER, ALLOCATABLE :: piv_pos(:)
|
||||
!! Position of the pivot points
|
||||
INTEGER, ALLOCATABLE :: piv_spin(:)
|
||||
!! Spin index of the pivot points. 1 for spin up, 2 for spin down.
|
||||
COMPLEX(DP) :: tmp_cwork(2)
|
||||
COMPLEX(DP) :: nowfc_tmp ! jml
|
||||
REAL(DP):: norm_psi, focc, arg, tpi_r_dot_g, xk_cry(3), rpos_cart(3)
|
||||
INTEGER :: ik, npw, ibnd, iw, nrtot, info, lcwork, locibnd, &
|
||||
ib, gamma_idx, minmn, minmn2, maxmn2, &
|
||||
ig, ipool_gamma, ik_gamma_loc, i, j, k, ik_g_w90, nxxs ! jml
|
||||
COMPLEX(DP), ALLOCATABLE :: psic_all(:)
|
||||
ig, ipool_gamma, ik_gamma_loc, i, j, k, ik_g_w90, nxxs, count_piv_spin_up ! jml
|
||||
COMPLEX(DP), ALLOCATABLE :: psic_all(:, :)
|
||||
!
|
||||
INTEGER, EXTERNAL :: global_kpoint_index
|
||||
!
|
||||
|
|
@ -3838,21 +3839,23 @@ SUBROUTINE compute_amn_with_scdm
|
|||
nrtot = dffts%nr1*dffts%nr2*dffts%nr3
|
||||
nxxs = dffts%nr1x * dffts%nr2x * dffts%nr3x
|
||||
info = 0
|
||||
minmn = MIN(num_bands,nrtot)
|
||||
minmn = MIN(num_bands, npol*nrtot)
|
||||
ALLOCATE(qr_tau(2*minmn))
|
||||
ALLOCATE(piv(nrtot))
|
||||
ALLOCATE(piv(npol*nrtot))
|
||||
ALLOCATE(rwork(2*npol*nrtot))
|
||||
piv(:) = 0
|
||||
ALLOCATE(rwork(2*nrtot))
|
||||
rwork(:) = 0.0_DP
|
||||
|
||||
ALLOCATE(nowfc(n_wannier,num_bands))
|
||||
ALLOCATE(psi_gamma(nrtot,num_bands))
|
||||
ALLOCATE(psi_gamma(npol*nrtot,num_bands))
|
||||
minmn2 = MIN(num_bands,n_wannier)
|
||||
maxmn2 = MAX(num_bands,n_wannier)
|
||||
ALLOCATE(rwork2(5*minmn2))
|
||||
|
||||
ALLOCATE(psic_all(nxxs))
|
||||
ALLOCATE(piv_pos(n_wannier))
|
||||
ALLOCATE(piv_spin(n_wannier))
|
||||
ALLOCATE(rpos(3, n_wannier))
|
||||
ALLOCATE(psic_all(nxxs, npol))
|
||||
ALLOCATE(phase(n_wannier))
|
||||
ALLOCATE(singval(n_wannier))
|
||||
ALLOCATE(Umat(num_bands,n_wannier))
|
||||
|
|
@ -3912,33 +3915,58 @@ SUBROUTINE compute_amn_with_scdm
|
|||
END IF
|
||||
!
|
||||
npw = ngk(ik)
|
||||
! vv: Compute unk's on a real grid (the fft grid)
|
||||
psic(:) = (0.0_DP, 0.0_DP)
|
||||
psic(dffts%nl (igk_k (1:npw,ik) ) ) = evc (1:npw,ibnd)
|
||||
CALL invfft ('Wave', psic, dffts)
|
||||
!
|
||||
psic_all(:) = (0.0_DP, 0.0_DP)
|
||||
IF (noncolin) THEN
|
||||
psic_nc(:, :) = (0.0_DP, 0.0_DP)
|
||||
psic_nc( dffts%nl( igk_k(1:npw, ik) ), 1) = evc(1:npw, ibnd)
|
||||
psic_nc( dffts%nl( igk_k(1:npw, ik) ), 2) = evc(1+npwx:npw+npwx, ibnd)
|
||||
CALL invfft('Wave', psic_nc(:,1), dffts)
|
||||
CALL invfft('Wave', psic_nc(:,2), dffts)
|
||||
!
|
||||
#if defined(__MPI)
|
||||
CALL gather_grid(dffts, psic, psic_all)
|
||||
CALL gather_grid(dffts, psic_nc(:, 1), psic_all(:, 1))
|
||||
CALL gather_grid(dffts, psic_nc(:, 2), psic_all(:, 2))
|
||||
#else
|
||||
psic_all(1:nrtot) = psic(1:nrtot)
|
||||
psic_all = psic_nc
|
||||
#endif
|
||||
! vv: Gamma only
|
||||
! vv: Build Psi_k = Unk * focc
|
||||
norm_psi = SQRT(SUM( ABS(psic_all(1:nrtot))**2 ))
|
||||
psi_gamma(1:nrtot, locibnd) = psic_all(1:nrtot) * (focc / norm_psi)
|
||||
! vv: Gamma only
|
||||
! vv: Build Psi_k = Unk * focc
|
||||
norm_psi = SQRT(SUM(ABS(psic_all(1:nrtot, 1:2))**2))
|
||||
psi_gamma(1:nrtot, locibnd) = psic_all(1:nrtot, 1) * focc / norm_psi
|
||||
psi_gamma(1+nrtot:2*nrtot, locibnd) = psic_all(1:nrtot, 2) * focc / norm_psi
|
||||
ELSE
|
||||
! spin-collinear case
|
||||
! vv: Compute unk's on a real grid (the fft grid)
|
||||
psic(:) = (0.0_DP, 0.0_DP)
|
||||
psic( dffts%nl( igk_k(1:npw,ik) ) ) = evc(1:npw, ibnd)
|
||||
CALL invfft ('Wave', psic, dffts)
|
||||
!
|
||||
psic_all(:, 1) = (0.0_DP, 0.0_DP)
|
||||
#if defined(__MPI)
|
||||
CALL gather_grid(dffts, psic, psic_all(:, 1))
|
||||
#else
|
||||
psic_all(1:nrtot, 1) = psic(1:nrtot)
|
||||
#endif
|
||||
! vv: Gamma only
|
||||
! vv: Build Psi_k = Unk * focc
|
||||
norm_psi = SQRT(SUM( ABS(psic_all(1:nrtot, 1))**2 ))
|
||||
psi_gamma(1:nrtot, locibnd) = psic_all(1:nrtot, 1) * (focc / norm_psi)
|
||||
ENDIF
|
||||
!
|
||||
ENDDO
|
||||
!
|
||||
! vv: Perform QR factorization with pivoting on Psi_Gamma
|
||||
! vv: Preliminary call to define optimal values for lwork and cwork size
|
||||
! Perform QR factorization only in a single processer
|
||||
IF(me_pool == root_pool) THEN
|
||||
CALL ZGEQP3(num_bands,nrtot,TRANSPOSE(CONJG(psi_gamma)),num_bands,piv,qr_tau,tmp_cwork,-1,rwork,info)
|
||||
CALL ZGEQP3(num_bands, npol*nrtot, TRANSPOSE(CONJG(psi_gamma)), num_bands, &
|
||||
piv, qr_tau, tmp_cwork, -1, rwork, info)
|
||||
IF (info/=0) CALL errore('compute_amn', 'Error in priliminary call for the QR factorization', 1)
|
||||
lcwork = AINT(REAL(tmp_cwork(1)))
|
||||
piv(:) = 0
|
||||
ALLOCATE(cwork(lcwork))
|
||||
CALL ZGEQP3(num_bands,nrtot,TRANSPOSE(CONJG(psi_gamma)),num_bands,piv,qr_tau,cwork,lcwork,rwork,info)
|
||||
CALL ZGEQP3(num_bands, npol*nrtot, TRANSPOSE(CONJG(psi_gamma)), num_bands, &
|
||||
piv, qr_tau, cwork, lcwork, rwork, info)
|
||||
DEALLOCATE(cwork)
|
||||
ENDIF
|
||||
CALL mp_bcast(info, root_pool, intra_pool_comm)
|
||||
|
|
@ -3951,9 +3979,30 @@ SUBROUTINE compute_amn_with_scdm
|
|||
!
|
||||
CALL stop_clock('scdm_QRCP')
|
||||
!
|
||||
! noncollinear case: calculate position and spin part of piv
|
||||
IF (noncolin) THEN
|
||||
count_piv_spin_up = 0
|
||||
DO iw = 1, n_wannier
|
||||
IF (piv(iw) <= nrtot) THEN
|
||||
piv_pos(iw) = piv(iw)
|
||||
piv_spin(iw) = 1
|
||||
count_piv_spin_up = count_piv_spin_up + 1
|
||||
ELSE
|
||||
piv_pos(iw) = piv(iw) - nrtot
|
||||
piv_spin(iw) = -1
|
||||
ENDIF
|
||||
ENDDO
|
||||
WRITE(stdout, '(a,I5)') " Number of pivot points with spin up : ", count_piv_spin_up
|
||||
WRITE(stdout, '(a,I5)') " Number of pivot points with spin down: ", n_wannier - count_piv_spin_up
|
||||
ELSE
|
||||
piv_pos(:) = piv(:)
|
||||
piv_spin(:) = 0
|
||||
WRITE(stdout, '(a,i8)') ' Number of pivot points: ', n_wannier
|
||||
ENDIF
|
||||
!
|
||||
! vv: Compute the points in 3d grid
|
||||
DO iw = 1, n_wannier
|
||||
i = piv(iw) - 1
|
||||
i = piv_pos(iw) - 1
|
||||
k = i / (dffts%nr1 * dffts%nr2)
|
||||
i = i - (dffts%nr1 * dffts%nr2) * k
|
||||
j = i / dffts%nr1
|
||||
|
|
@ -3964,12 +4013,20 @@ SUBROUTINE compute_amn_with_scdm
|
|||
rpos(:, iw) = rpos(:, iw) - ANINT(rpos(:, iw))
|
||||
ENDDO
|
||||
!
|
||||
WRITE(stdout, '(a,i8)') ' Number of pivot points: ', n_wannier
|
||||
WRITE(stdout, '(a)') ' Pivot point positions (alat units):'
|
||||
! Print pivot positions (and spin indices)
|
||||
IF (noncolin) THEN
|
||||
WRITE(stdout, '(a)') ' Pivot point positions (alat units) and spin indices:'
|
||||
ELSE
|
||||
WRITE(stdout, '(a)') ' Pivot point positions (alat units):'
|
||||
ENDIF
|
||||
DO iw = 1, n_wannier
|
||||
rpos_cart(:) = rpos(:, iw)
|
||||
CALL cryst_to_cart(1, rpos_cart, at, +1)
|
||||
WRITE(stdout, '(I8,3F12.6)') iw, rpos_cart
|
||||
IF (noncolin) THEN
|
||||
WRITE(stdout, '(I8,3F12.6,I3)') iw, rpos_cart, piv_spin(iw)
|
||||
ELSE
|
||||
WRITE(stdout, '(I8,3F12.6)') iw, rpos_cart
|
||||
ENDIF
|
||||
ENDDO
|
||||
WRITE(stdout, *)
|
||||
!
|
||||
|
|
@ -4030,12 +4087,24 @@ SUBROUTINE compute_amn_with_scdm
|
|||
END IF
|
||||
!
|
||||
norm_psi = SUM( ABS(evc(1:npw, ibnd))**2 )
|
||||
IF (noncolin) norm_psi = norm_psi + SUM( ABS(evc(1+npwx:npw+npwx, ibnd))**2 )
|
||||
CALL mp_sum(norm_psi, intra_pool_comm)
|
||||
norm_psi = SQRT(norm_psi)
|
||||
!
|
||||
! jml: nowfc = sum_G (psi(G) * exp(i*G*r)) * focc * phase(iw) / norm_psi
|
||||
DO iw = 1, n_wannier
|
||||
nowfc_tmp = SUM( evc(1:npw, ibnd) * phase_g(1:npw, iw) )
|
||||
IF (noncolin) THEN
|
||||
IF (piv_spin(iw) == 1) THEN
|
||||
! spin up
|
||||
nowfc_tmp = SUM( evc(1:npw, ibnd) * phase_g(1:npw, iw) )
|
||||
ELSE
|
||||
! spin down
|
||||
nowfc_tmp = SUM( evc(1+npwx:npw+npwx, ibnd) * phase_g(1:npw, iw) )
|
||||
ENDIF
|
||||
ELSE
|
||||
! spin collinear
|
||||
nowfc_tmp = SUM( evc(1:npw, ibnd) * phase_g(1:npw, iw) )
|
||||
ENDIF
|
||||
nowfc(iw,locibnd) = nowfc_tmp * phase(iw) * focc / norm_psi
|
||||
ENDDO
|
||||
!
|
||||
|
|
@ -4079,6 +4148,8 @@ SUBROUTINE compute_amn_with_scdm
|
|||
DEALLOCATE(psi_gamma)
|
||||
DEALLOCATE(nowfc)
|
||||
DEALLOCATE(piv)
|
||||
DEALLOCATE(piv_pos)
|
||||
DEALLOCATE(piv_spin)
|
||||
DEALLOCATE(qr_tau)
|
||||
DEALLOCATE(rwork)
|
||||
DEALLOCATE(rwork2)
|
||||
|
|
@ -4096,394 +4167,6 @@ SUBROUTINE compute_amn_with_scdm
|
|||
!
|
||||
END SUBROUTINE compute_amn_with_scdm
|
||||
|
||||
|
||||
SUBROUTINE compute_amn_with_scdm_spinor
|
||||
!
|
||||
! jml: scdm for noncollinear case
|
||||
!
|
||||
USE constants, ONLY : rytoev, pi
|
||||
USE io_global, ONLY : stdout, ionode, ionode_id
|
||||
USE wvfct, ONLY : nbnd, et, npwx
|
||||
USE gvecw, ONLY : gcutw
|
||||
USE control_flags, ONLY : gamma_only
|
||||
USE wavefunctions, ONLY : evc, psic_nc
|
||||
USE io_files, ONLY : nwordwfc, iunwfc
|
||||
USE wannier
|
||||
USE klist, ONLY : nkstot, xk, ngk, igk_k
|
||||
USE gvect, ONLY : g, ngm, mill
|
||||
USE fft_base, ONLY : dffts !vv: unk for the SCDM-k algorithm
|
||||
USE scatter_mod, ONLY : gather_grid
|
||||
USE fft_interfaces, ONLY : invfft !vv: inverse fft transform for computing the unk's on a grid
|
||||
USE noncollin_module,ONLY : noncolin, npol
|
||||
USE mp, ONLY : mp_bcast, mp_barrier, mp_sum
|
||||
USE mp_world, ONLY : world_comm
|
||||
USE mp_pools, ONLY : intra_pool_comm
|
||||
USE cell_base, ONLY : at
|
||||
USE ions_base, ONLY : ntyp => nsp, tau
|
||||
USE uspp, ONLY : okvan
|
||||
|
||||
IMPLICIT NONE
|
||||
|
||||
INTEGER, EXTERNAL :: find_free_unit
|
||||
COMPLEX(DP), ALLOCATABLE :: phase(:), nowfc1(:,:), nowfc(:,:), psi_gamma(:,:), &
|
||||
qr_tau(:), cwork(:), cwork2(:), Umat(:,:), VTmat(:,:), Amat(:,:) ! vv: complex arrays for the SVD factorization
|
||||
COMPLEX(DP), ALLOCATABLE :: phase_g(:,:) ! jml
|
||||
REAL(DP), ALLOCATABLE :: focc(:), rwork(:), rwork2(:), singval(:), rpos(:,:), cpos(:,:) ! vv: Real array for the QR factorization and SVD
|
||||
INTEGER, ALLOCATABLE :: piv(:) ! vv: Pivot array in the QR factorization
|
||||
INTEGER, ALLOCATABLE :: piv_pos(:), piv_spin(:) ! jml: position and spin index of piv
|
||||
COMPLEX(DP) :: tmp_cwork(2)
|
||||
COMPLEX(DP) :: nowfc_tmp ! jml
|
||||
REAL(DP):: ddot, sumk, norm_psi, f_gamma, tpi_r_dot_g
|
||||
INTEGER :: ik, npw, ibnd, iw, ikevc, nrtot, ipt, info, lcwork, locibnd, &
|
||||
jpt,kpt,lpt, ib, istart, gamma_idx, minmn, minmn2, maxmn2, numbands, nbtot, &
|
||||
ig, ig_local, count_piv_spin, ispin ! jml
|
||||
CHARACTER (len=9) :: cdate,ctime
|
||||
CHARACTER (len=60) :: header
|
||||
LOGICAL :: found_gamma
|
||||
|
||||
#if defined(__MPI)
|
||||
INTEGER :: nxxs
|
||||
COMPLEX(DP),ALLOCATABLE :: psic_all(:,:)
|
||||
nxxs = dffts%nr1x * dffts%nr2x * dffts%nr3x
|
||||
ALLOCATE(psic_all(nxxs, 2) )
|
||||
#endif
|
||||
|
||||
! vv: Write info about SCDM in output
|
||||
IF (TRIM(scdm_entanglement) == 'isolated') THEN
|
||||
WRITE(stdout,'(1x,a,a/)') 'Case : ',trim(scdm_entanglement)
|
||||
ELSEIF (TRIM(scdm_entanglement) == 'erfc' .OR. &
|
||||
TRIM(scdm_entanglement) == 'gaussian') THEN
|
||||
WRITE(stdout,'(1x,a,a)') 'Case : ',trim(scdm_entanglement)
|
||||
WRITE(stdout,'(1x,a,f10.3,a/,1x,a,f10.3,a/)') 'mu = ', scdm_mu, ' eV', 'sigma =', scdm_sigma, ' eV'
|
||||
ENDIF
|
||||
|
||||
CALL start_clock( 'compute_amn' )
|
||||
|
||||
! vv: Error for using SCDM with Ultrasoft pseudopotentials
|
||||
!IF (okvan) THEN
|
||||
! call errore('pw2wannier90','The SCDM method does not work with Ultrasoft pseudopotential yet.',1)
|
||||
!ENDIF
|
||||
|
||||
! vv: Error for using SCDM with gamma_only
|
||||
IF (gamma_only) THEN
|
||||
call errore('pw2wannier90','The SCDM method does not work with gamma_only calculations.',1)
|
||||
ENDIF
|
||||
! vv: Allocate all the variables for the SCDM method:
|
||||
! 1)For the QR decomposition
|
||||
! 2)For the unk's on the real grid
|
||||
! 3)For the SVD
|
||||
IF(TRIM(scdm_entanglement) == 'isolated') THEN
|
||||
numbands=n_wannier
|
||||
nbtot=n_wannier + nexband
|
||||
ELSE
|
||||
numbands=nbnd-nexband
|
||||
nbtot=nbnd
|
||||
ENDIF
|
||||
nrtot = dffts%nr1*dffts%nr2*dffts%nr3
|
||||
info = 0
|
||||
minmn = MIN(numbands,nrtot*2) ! jml: spinor
|
||||
ALLOCATE(qr_tau(2*minmn))
|
||||
ALLOCATE(piv(nrtot*2)) ! jml: spinor
|
||||
ALLOCATE(piv_pos(n_wannier)) ! jml: spinor
|
||||
ALLOCATE(piv_spin(n_wannier)) ! jml: spinor
|
||||
piv(:) = 0
|
||||
ALLOCATE(rwork(2*nrtot*2)) ! jml: spinor
|
||||
rwork(:) = 0.0_DP
|
||||
|
||||
ALLOCATE(kpt_latt(3,iknum))
|
||||
ALLOCATE(nowfc1(n_wannier,numbands))
|
||||
ALLOCATE(nowfc(n_wannier,numbands))
|
||||
ALLOCATE(psi_gamma(nrtot*2,numbands)) ! jml: spinor
|
||||
ALLOCATE(focc(numbands))
|
||||
minmn2 = MIN(numbands,n_wannier)
|
||||
maxmn2 = MAX(numbands,n_wannier)
|
||||
ALLOCATE(rwork2(5*minmn2))
|
||||
|
||||
ALLOCATE(rpos(nrtot,3)) ! jml: spinor
|
||||
ALLOCATE(cpos(n_wannier,3))
|
||||
ALLOCATE(phase(n_wannier))
|
||||
ALLOCATE(singval(n_wannier))
|
||||
ALLOCATE(Umat(numbands,n_wannier))
|
||||
ALLOCATE(VTmat(n_wannier,n_wannier))
|
||||
ALLOCATE(Amat(numbands,n_wannier))
|
||||
|
||||
IF (wan_mode=='library') ALLOCATE(a_mat(num_bands,n_wannier,iknum))
|
||||
|
||||
IF (wan_mode=='standalone') THEN
|
||||
iun_amn = find_free_unit()
|
||||
IF (ionode) OPEN (unit=iun_amn, file=trim(seedname)//".amn",form='formatted')
|
||||
ENDIF
|
||||
|
||||
WRITE(stdout,'(a,i8)') ' AMN: iknum = ',iknum
|
||||
!
|
||||
IF (wan_mode=='standalone') THEN
|
||||
CALL date_and_tim( cdate, ctime )
|
||||
header='Created on '//cdate//' at '//ctime//' with SCDM '
|
||||
IF (ionode) THEN
|
||||
WRITE (iun_amn,*) header
|
||||
WRITE (iun_amn,'(3i8,xxx,2f10.6)') numbands, iknum, n_wannier, scdm_mu, scdm_sigma
|
||||
ENDIF
|
||||
ENDIF
|
||||
|
||||
!vv: Find Gamma-point index in the list of k-vectors
|
||||
ik = 0
|
||||
gamma_idx = 1
|
||||
sumk = -1.0_DP
|
||||
found_gamma = .false.
|
||||
kpt_latt(:,1:iknum)=xk(:,1:iknum)
|
||||
CALL cryst_to_cart(iknum,kpt_latt,at,-1)
|
||||
DO WHILE(sumk/=0.0_DP .and. ik < iknum)
|
||||
ik = ik + 1
|
||||
sumk = ABS(kpt_latt(1,ik)**2 + kpt_latt(2,ik)**2 + kpt_latt(3,ik)**2)
|
||||
IF (sumk==0.0_DP) THEN
|
||||
found_gamma = .true.
|
||||
gamma_idx = ik
|
||||
ENDIF
|
||||
END DO
|
||||
IF (.not. found_gamma) call errore('compute_amn','No Gamma point found.',1)
|
||||
|
||||
f_gamma = 0.0_DP
|
||||
ik = gamma_idx
|
||||
locibnd = 0
|
||||
CALL davcio (evc, 2*nwordwfc, iunwfc, ik, -1 )
|
||||
DO ibnd=1,nbtot
|
||||
IF(excluded_band(ibnd)) CYCLE
|
||||
locibnd = locibnd + 1
|
||||
! check locibnd <= numbands
|
||||
IF (locibnd > numbands) call errore('compute_amn','Something wrong with the number of bands. Check exclude_bands.')
|
||||
IF(TRIM(scdm_entanglement) == 'isolated') THEN
|
||||
f_gamma = 1.0_DP
|
||||
ELSEIF (TRIM(scdm_entanglement) == 'erfc') THEN
|
||||
f_gamma = 0.5_DP*ERFC((et(ibnd,ik)*rytoev - scdm_mu)/scdm_sigma)
|
||||
ELSEIF (TRIM(scdm_entanglement) == 'gaussian') THEN
|
||||
f_gamma = EXP(-1.0_DP*((et(ibnd,ik)*rytoev - scdm_mu)**2)/(scdm_sigma**2))
|
||||
ELSE
|
||||
call errore('compute_amn','scdm_entanglement value not recognized.',1)
|
||||
END IF
|
||||
npw = ngk(ik)
|
||||
! vv: Compute unk's on a real grid (the fft grid)
|
||||
psic_nc(:,:) = (0.D0,0.D0)
|
||||
psic_nc(dffts%nl (igk_k (1:npw,ik) ), 1) = evc (1:npw,ibnd)
|
||||
psic_nc(dffts%nl (igk_k (1:npw,ik) ), 2) = evc (1+npwx:npw+npwx,ibnd)
|
||||
CALL invfft ('Wave', psic_nc(:,1), dffts)
|
||||
CALL invfft ('Wave', psic_nc(:,2), dffts)
|
||||
|
||||
#if defined(__MPI)
|
||||
CALL gather_grid(dffts, psic_nc(:,1), psic_all(:,1))
|
||||
CALL gather_grid(dffts, psic_nc(:,2), psic_all(:,2))
|
||||
norm_psi = sqrt( real(sum(psic_all(1:nrtot, 1)*conjg(psic_all(1:nrtot, 1))),kind=DP) &
|
||||
+real(sum(psic_all(1:nrtot, 2)*conjg(psic_all(1:nrtot, 2))),kind=DP) )
|
||||
! vv: Gamma only
|
||||
! vv: Build Psi_k = Unk * focc
|
||||
psi_gamma(1:nrtot, locibnd) = psic_all(1:nrtot, 1) * f_gamma / norm_psi
|
||||
psi_gamma(1+nrtot:2*nrtot,locibnd) = psic_all(1:nrtot, 2) * f_gamma / norm_psi
|
||||
#else
|
||||
norm_psi = sqrt( real(sum(psic_nc(1:nrtot, 1)*conjg(psic_nc(1:nrtot, 1))),kind=DP) &
|
||||
+real(sum(psic_nc(1:nrtot, 2)*conjg(psic_nc(1:nrtot, 2))),kind=DP) )
|
||||
psi_gamma(1:nrtot, locibnd) = psic_nc(1:nrtot, 1) * f_gamma / norm_psi
|
||||
psi_gamma(1+nrtot:2*nrtot,locibnd) = psic_nc(1:nrtot, 2) * f_gamma / norm_psi
|
||||
#endif
|
||||
ENDDO
|
||||
|
||||
! vv: Perform QR factorization with pivoting on Psi_Gamma
|
||||
! vv: Preliminary call to define optimal values for lwork and cwork size
|
||||
CALL ZGEQP3(numbands,nrtot*2,TRANSPOSE(CONJG(psi_gamma)),numbands,piv,qr_tau,tmp_cwork,-1,rwork,info)
|
||||
IF(info/=0) call errore('compute_amn','Error in computing the QR factorization',1)
|
||||
lcwork = AINT(REAL(tmp_cwork(1)))
|
||||
tmp_cwork(:) = (0.0_DP,0.0_DP)
|
||||
piv(:) = 0
|
||||
rwork(:) = 0.0_DP
|
||||
ALLOCATE(cwork(lcwork))
|
||||
cwork(:) = (0.0_DP,0.0_DP)
|
||||
#if defined(__MPI)
|
||||
IF(ionode) THEN
|
||||
CALL ZGEQP3(numbands,nrtot*2,TRANSPOSE(CONJG(psi_gamma)),numbands,piv,qr_tau,cwork,lcwork,rwork,info)
|
||||
IF(info/=0) call errore('compute_amn','Error in computing the QR factorization',1)
|
||||
ENDIF
|
||||
CALL mp_bcast(piv,ionode_id,world_comm)
|
||||
#else
|
||||
! vv: Perform QR factorization with pivoting on Psi_Gamma
|
||||
CALL ZGEQP3(numbands,nrtot*2,TRANSPOSE(CONJG(psi_gamma)),numbands,piv,qr_tau,cwork,lcwork,rwork,info)
|
||||
IF(info/=0) call errore('compute_amn','Error in computing the QR factorization',1)
|
||||
#endif
|
||||
DEALLOCATE(cwork)
|
||||
tmp_cwork(:) = (0.0_DP,0.0_DP)
|
||||
|
||||
! jml: calculate position and spin part of piv
|
||||
count_piv_spin = 0
|
||||
DO iw = 1, n_wannier
|
||||
IF (piv(iw) .le. nrtot) then
|
||||
piv_pos(iw) = piv(iw)
|
||||
piv_spin(iw) = 1
|
||||
count_piv_spin = count_piv_spin + 1
|
||||
else
|
||||
piv_pos(iw) = piv(iw) - nrtot
|
||||
piv_spin(iw) = 2
|
||||
end if
|
||||
END DO
|
||||
WRITE(stdout, '(a,I5)') " Number of pivot points with spin up : ", count_piv_spin
|
||||
WRITE(stdout, '(a,I5)') " Number of pivot points with spin down: ", n_wannier - count_piv_spin
|
||||
|
||||
! vv: Compute the points
|
||||
lpt = 0
|
||||
rpos(:,:) = 0.0_DP
|
||||
cpos(:,:) = 0.0_DP
|
||||
DO kpt = 0,dffts%nr3-1
|
||||
DO jpt = 0,dffts%nr2-1
|
||||
DO ipt = 0,dffts%nr1-1
|
||||
lpt = lpt + 1
|
||||
rpos(lpt,1) = DBLE(ipt)/DBLE(dffts%nr1)
|
||||
rpos(lpt,2) = DBLE(jpt)/DBLE(dffts%nr2)
|
||||
rpos(lpt,3) = DBLE(kpt)/DBLE(dffts%nr3)
|
||||
ENDDO
|
||||
ENDDO
|
||||
ENDDO
|
||||
DO iw=1,n_wannier
|
||||
cpos(iw,:) = rpos(piv_pos(iw),:)
|
||||
cpos(iw,:) = cpos(iw,:) - ANINT(cpos(iw,:))
|
||||
ENDDO
|
||||
|
||||
DO ik=1,iknum
|
||||
WRITE (stdout,'(i8)',advance='no') ik
|
||||
IF( MOD(ik,10) == 0 ) WRITE (stdout,*)
|
||||
FLUSH(stdout)
|
||||
ikevc = ik + ikstart - 1
|
||||
|
||||
! vv: SCDM method for generating the Amn matrix
|
||||
! jml: calculate of psi_nk at pivot points using slow FT
|
||||
! This is faster than using invfft because the number of pivot
|
||||
! points is much smaller than the number of FFT grid points.
|
||||
phase(:) = (0.0_DP,0.0_DP)
|
||||
nowfc1(:,:) = (0.0_DP,0.0_DP)
|
||||
nowfc(:,:) = (0.0_DP,0.0_DP)
|
||||
Umat(:,:) = (0.0_DP,0.0_DP)
|
||||
VTmat(:,:) = (0.0_DP,0.0_DP)
|
||||
Amat(:,:) = (0.0_DP,0.0_DP)
|
||||
singval(:) = 0.0_DP
|
||||
rwork2(:) = 0.0_DP
|
||||
|
||||
! jml: calculate phase factors before the loop over bands
|
||||
npw = ngk(ik)
|
||||
ALLOCATE(phase_g(npw, n_wannier))
|
||||
DO iw = 1, n_wannier
|
||||
phase(iw) = cmplx(COS(2.0_DP*pi*(cpos(iw,1)*kpt_latt(1,ik) + &
|
||||
&cpos(iw,2)*kpt_latt(2,ik) + cpos(iw,3)*kpt_latt(3,ik))), & !*ddot(3,cpos(iw,:),1,kpt_latt(:,ik),1)),&
|
||||
&SIN(2.0_DP*pi*(cpos(iw,1)*kpt_latt(1,ik) + &
|
||||
&cpos(iw,2)*kpt_latt(2,ik) + cpos(iw,3)*kpt_latt(3,ik))),kind=DP) !ddot(3,cpos(iw,:),1,kpt_latt(:,ik),1)))
|
||||
|
||||
DO ig_local = 1, npw
|
||||
ig = igk_k(ig_local,ik)
|
||||
tpi_r_dot_g = 2.0_DP * pi * ( cpos(iw,1) * REAL(mill(1,ig), DP) &
|
||||
& + cpos(iw,2) * REAL(mill(2,ig), DP) &
|
||||
& + cpos(iw,3) * REAL(mill(3,ig), DP) )
|
||||
phase_g(ig_local, iw) = cmplx(COS(tpi_r_dot_g), SIN(tpi_r_dot_g), kind=DP)
|
||||
END DO
|
||||
END DO
|
||||
|
||||
locibnd = 0
|
||||
CALL davcio (evc, 2*nwordwfc, iunwfc, ikevc, -1 )
|
||||
DO ibnd=1,nbtot
|
||||
IF (excluded_band(ibnd)) CYCLE
|
||||
locibnd = locibnd + 1
|
||||
! vv: Define the occupation numbers matrix according to scdm_entanglement
|
||||
IF(TRIM(scdm_entanglement) == 'isolated') THEN
|
||||
focc(locibnd) = 1.0_DP
|
||||
ELSEIF (TRIM(scdm_entanglement) == 'erfc') THEN
|
||||
focc(locibnd) = 0.5_DP*ERFC((et(ibnd,ik)*rytoev - scdm_mu)/scdm_sigma)
|
||||
ELSEIF (TRIM(scdm_entanglement) == 'gaussian') THEN
|
||||
focc(locibnd) = EXP(-1.0_DP*((et(ibnd,ik)*rytoev - scdm_mu)**2)/(scdm_sigma**2))
|
||||
ELSE
|
||||
call errore('compute_amn','scdm_entanglement value not recognized.',1)
|
||||
END IF
|
||||
|
||||
norm_psi= REAL(SUM( evc(1:npw,ibnd) * CONJG(evc(1:npw,ibnd)) )) &
|
||||
+ REAL(SUM( evc(1+npwx:npw+npwx,ibnd) * CONJG(evc(1+npwx:npw+npwx,ibnd)) ))
|
||||
CALL mp_sum(norm_psi, intra_pool_comm)
|
||||
norm_psi= sqrt(norm_psi)
|
||||
|
||||
! jml: nowfc = sum_G (psi(G) * exp(i*G*r)) * focc * phase(iw) / norm_psi
|
||||
DO iw = 1, n_wannier
|
||||
if (piv_spin(iw) == 1) then ! spin up
|
||||
nowfc_tmp = sum( evc(1:npw, ibnd) * phase_g(1:npw, iw) )
|
||||
else ! spin down
|
||||
nowfc_tmp = sum( evc(1+npwx:npw+npwx, ibnd) * phase_g(1:npw, iw) )
|
||||
end if
|
||||
|
||||
nowfc(iw, locibnd) = nowfc_tmp * phase(iw) * focc(locibnd) / norm_psi
|
||||
ENDDO
|
||||
|
||||
END DO ! ibnd
|
||||
CALL mp_sum(nowfc, intra_pool_comm) ! jml
|
||||
DEALLOCATE(phase_g) ! jml
|
||||
|
||||
CALL ZGESVD('S','S',numbands,n_wannier,TRANSPOSE(CONJG(nowfc)),numbands,&
|
||||
&singval,Umat,numbands,VTmat,n_wannier,tmp_cwork,-1,rwork2,info)
|
||||
lcwork = AINT(REAL(tmp_cwork(1)))
|
||||
tmp_cwork(:) = (0.0_DP,0.0_DP)
|
||||
ALLOCATE(cwork(lcwork))
|
||||
#if defined(__MPI)
|
||||
IF(ionode) THEN
|
||||
! vv: SVD to generate orthogonal projections
|
||||
CALL ZGESVD('S','S',numbands,n_wannier,TRANSPOSE(CONJG(nowfc)),numbands,&
|
||||
&singval,Umat,numbands,VTmat,n_wannier,cwork,lcwork,rwork2,info)
|
||||
IF(info/=0) CALL errore('compute_amn','Error in computing the SVD of the PSI matrix in the SCDM method',1)
|
||||
ENDIF
|
||||
CALL mp_bcast(Umat,ionode_id,world_comm)
|
||||
CALL mp_bcast(VTmat,ionode_id,world_comm)
|
||||
#else
|
||||
! vv: SVD to generate orthogonal projections
|
||||
CALL ZGESVD('S','S',numbands,n_wannier,TRANSPOSE(CONJG(nowfc)),numbands,&
|
||||
&singval,Umat,numbands,VTmat,n_wannier,cwork,lcwork,rwork2,info)
|
||||
IF(info/=0) CALL errore('compute_amn','Error in computing the SVD of the PSI matrix in the SCDM method',1)
|
||||
#endif
|
||||
DEALLOCATE(cwork)
|
||||
|
||||
Amat = MATMUL(Umat,VTmat)
|
||||
|
||||
CALL start_clock( 'scdm_write' )
|
||||
DO iw = 1,n_wannier
|
||||
locibnd = 0
|
||||
DO ibnd = 1,nbtot
|
||||
IF (excluded_band(ibnd)) CYCLE
|
||||
locibnd = locibnd + 1
|
||||
IF (ionode) WRITE(iun_amn,'(3i5,2f18.12)') locibnd, iw, ik, REAL(Amat(locibnd,iw)), AIMAG(Amat(locibnd,iw))
|
||||
ENDDO
|
||||
ENDDO
|
||||
CALL stop_clock( 'scdm_write' )
|
||||
ENDDO ! k-points
|
||||
|
||||
! vv: Deallocate all the variables for the SCDM method
|
||||
DEALLOCATE(kpt_latt)
|
||||
DEALLOCATE(psi_gamma)
|
||||
DEALLOCATE(nowfc)
|
||||
DEALLOCATE(nowfc1)
|
||||
DEALLOCATE(focc)
|
||||
DEALLOCATE(piv)
|
||||
DEALLOCATE(piv_pos)
|
||||
DEALLOCATE(piv_spin)
|
||||
DEALLOCATE(qr_tau)
|
||||
DEALLOCATE(rwork)
|
||||
DEALLOCATE(rwork2)
|
||||
DEALLOCATE(rpos)
|
||||
DEALLOCATE(cpos)
|
||||
DEALLOCATE(Umat)
|
||||
DEALLOCATE(VTmat)
|
||||
DEALLOCATE(Amat)
|
||||
DEALLOCATE(singval)
|
||||
|
||||
#if defined(__MPI)
|
||||
DEALLOCATE( psic_all )
|
||||
#endif
|
||||
|
||||
IF (ionode .and. wan_mode=='standalone') CLOSE (iun_amn)
|
||||
WRITE(stdout,'(/)')
|
||||
WRITE(stdout,*) ' AMN calculated'
|
||||
CALL stop_clock( 'compute_amn' )
|
||||
|
||||
RETURN
|
||||
END SUBROUTINE compute_amn_with_scdm_spinor
|
||||
|
||||
|
||||
subroutine orient_gf_spinor(npw)
|
||||
use constants, only: eps6
|
||||
use noncollin_module, only: npol
|
||||
|
|
|
|||
Loading…
Reference in New Issue