Commit 3f84f373 authored by Alexander Neukirchen's avatar Alexander Neukirchen
Browse files

Merge branch 'juPhon/DFPT' into 'develop'

Juphon/dfpt to dev

See merge request fleur/fleur!202
parents e27f5060 3b6d67b5
......@@ -133,7 +133,7 @@ CONTAINS
CALL setabc1lo(atoms,iType,usdus,jspin,alo1,blo1,clo1)
! generate the spinors (chi)
IF(noco%l_noco) ccchi=nococonv%chi(itype)
IF(noco%l_noco) ccchi=conjg(nococonv%umat(itype))
nintsp = 1
......@@ -212,7 +212,7 @@ CONTAINS
!$acc host_data use_device(work_c,abCoeffs,abTemp)
CALL zgemm_acc("T","C",ne,2*abSize,nvmax,CMPLX(1.0,0.0),work_c,MAXVAL(lapw%nv),abCoeffs,2*atoms%lmaxd*(atoms%lmaxd+2)+2,CMPLX(0.0,0.0),abTemp,acof_size)
CALL zgemm_acc("T","T",ne,2*abSize,nvmax,CMPLX(1.0,0.0),work_c,MAXVAL(lapw%nv),abCoeffs,2*atoms%lmaxd*(atoms%lmaxd+2)+2,CMPLX(0.0,0.0),abTemp,acof_size)
!$acc end host_data
!stop "DEBUG"
......@@ -302,8 +302,8 @@ CONTAINS
lm = ll1 + m
lmp = ll1 - m
inv_f = (-1)**(l-m)
acof(:,lmp,jatom)=acof(:,lmp,jatom)+inv_f*matmul(abCoeffs(lm+1,:),work_c(:nvmax,:))
bcof(:,lmp,jatom)=bcof(:,lmp,jatom)+inv_f*matmul(abCoeffs(lm+1+abSize,:),work_c(:nvmax,:))
acof(:,lmp,jatom)=acof(:,lmp,jatom)+inv_f*matmul(CONJG(abCoeffs(lm+1,:)),work_c(:nvmax,:)) !TODO: Is this conjugation costly?
bcof(:,lmp,jatom)=bcof(:,lmp,jatom)+inv_f*matmul(CONJG(abCoeffs(lm+1+abSize,:)),work_c(:nvmax,:)) !TODO: Is this conjugation costly?
!CALL zaxpy(ne,c_1,workTrans_c(:,iLAPW),1, acof(:,lmp,jatom),1)
!CALL zaxpy(ne,c_2,workTrans_c(:,iLAPW),1, bcof(:,lmp,jatom),1)
END DO
......@@ -336,10 +336,10 @@ CONTAINS
fgpl(:,iLAPW) = MATMUL(fgr,cell%bmat)
ENDDO
helpMat_c = abCoeffs(1+abSize:,:)
helpMat_c = CONJG(abCoeffs(1+abSize:,:)) !TODO: Is this conjugation costly?
workTrans_cf = 0.0
CALL zgemm("N","C",ne,atoms%lmaxd*(atoms%lmaxd+2)+1,nvmax,CMPLX(1.0,0.0),s2h_e,ne,abCoeffs,size(abcoeffs,1),CMPLX(1.0,0.0),force%e1cof(:,:,iAtom),ne)
CALL zgemm("N","T",ne,atoms%lmaxd*(atoms%lmaxd+2)+1,nvmax,CMPLX(1.0,0.0),s2h_e,ne,abCoeffs,size(abcoeffs,1),CMPLX(1.0,0.0),force%e1cof(:,:,iAtom),ne)
CALL zgemm("N","C",ne,atoms%lmaxd*(atoms%lmaxd+2)+1,nvmax,CMPLX(1.0,0.0),s2h_e,ne,helpMat_c,size(helpMat_c,1),CMPLX(1.0,0.0),force%e2cof(:,:,iAtom),ne)
DO i =1,3
IF (zmat%l_real) THEN
......@@ -351,7 +351,7 @@ CONTAINS
workTrans_cf(:,iLAPW) = workTrans_c(:,iLAPW) * fgpl(i,iLAPW)
ENDDO
ENDIF
CALL zgemm("N","C",ne,atoms%lmaxd*(atoms%lmaxd+2)+1,nvmax,CMPLX(1.0,0.0),workTrans_cf,ne,abCoeffs,size(abCoeffs,1),CMPLX(0.0,0.0),helpMat_force,ne)
CALL zgemm("N","T",ne,atoms%lmaxd*(atoms%lmaxd+2)+1,nvmax,CMPLX(1.0,0.0),workTrans_cf,ne,abCoeffs,size(abCoeffs,1),CMPLX(0.0,0.0),helpMat_force,ne)
force%aveccof(i,:,:,iAtom) = force%aveccof(i,:,:,iAtom) + helpMat_force(:,:)
CALL zgemm("N","C",ne,atoms%lmaxd*(atoms%lmaxd+2)+1,nvmax,CMPLX(1.0,0.0),workTrans_cf,ne,helpMat_c,size(helpMat_c,1),CMPLX(0.0,0.0),helpMat_force,ne)
force%bveccof(i,:,:,iAtom) = force%bveccof(i,:,:,iAtom) + helpMat_force(:,:)
......@@ -363,8 +363,8 @@ CONTAINS
ll1 = l* (l+1)
DO m = -l,l
lm = ll1 + m
c_1 = CONJG(abCoeffs(lm+1,iLAPW))
c_2 = CONJG(abCoeffs(lm+1+abSize,iLAPW))
c_1 = abCoeffs(lm+1,iLAPW)
c_2 = abCoeffs(lm+1+abSize,iLAPW)
jatom = sym%invsatnr(iAtom)
lmp = ll1 - m
inv_f = (-1)**(l-m)
......
......@@ -30,7 +30,6 @@ eigen/tlmplm.F90
eigen/tlo.f90
eigen/vacfun.F90
eigen/eigen_redist_matrix.f90
eigen/hs_int_onespin.F90
eigen/hs_int_direct.F90
)
#if (FLEUR_USE_GPU)
......
......@@ -100,13 +100,14 @@ CONTAINS
REAL, ALLOCATABLE :: bkpt(:)
REAL, ALLOCATABLE :: eig(:), eigBuffer(:,:,:)
TYPE(t_tlmplm) :: td, tdV1
TYPE(t_usdus) :: ud
TYPE(t_tlmplm) :: td, tdV1, tdmod
TYPE(t_usdus) :: ud, uddummy
TYPE(t_lapw) :: lapw
CLASS(t_mat), ALLOCATABLE :: zMat
CLASS(t_mat), ALLOCATABLE :: hmat,smat
CLASS(t_mat), ALLOCATABLE :: smat_unfold !used for unfolding bandstructure
TYPE(t_kpts) :: kqpts ! basically kpts, but with q added onto each one.
TYPE(t_hub1data) :: hub1datadummy
! Variables for HF or fi%hybinp functional calculation
INTEGER :: comm(fi%kpts%nkpt),irank2(fi%kpts%nkpt),isize2(fi%kpts%nkpt), dealloc_stat
......@@ -126,6 +127,7 @@ CONTAINS
END IF
call ud%init(fi%atoms,fi%input%jspins)
call uddummy%init(fi%atoms,fi%input%jspins)
ALLOCATE(eig(fi%input%neig))
ALLOCATE(bkpt(3))
ALLOCATE(eigBuffer(fi%input%neig,fi%kpts%nkpt,fi%input%jspins))
......@@ -144,9 +146,11 @@ CONTAINS
alpha_hybrid = MERGE(xcpot%get_exchange_weight(),0.0,hybdat%l_subvxc)
CALL mt_setup(fi%atoms,fi%sym,sphhar,fi%input,fi%noco,nococonv,enpara,fi%hub1inp,hub1data,inden,v,vx,fmpi,td,ud,alpha_hybrid)
! Get matrix elements of perturbed potential in DFPT case.
! Get matrix elements of perturbed potential and modified H/S in DFPT case.
IF (l_dfpteigen) THEN
hub1datadummy = hub1data
CALL dfpt_tlmplm(fi%atoms,fi%sym,sphhar,fi%input,fi%noco,enpara,fi%hub1inp,hub1data,v,fmpi,tdV1,v1real,v1imag)
CALL mt_setup(fi%atoms,fi%sym,sphhar,fi%input,fi%noco,nococonv,enpara,fi%hub1inp,hub1datadummy,inden,v,vx,fmpi,tdmod,uddummy,alpha_hybrid,.TRUE.)
END IF
neigBuffer = 0
......@@ -199,12 +203,14 @@ CONTAINS
select type(smat_unfold)
type is (t_mat)
smat_unfold=smat
smat_unfold%data_c=CONJG(smat%data_c)
end select
type is (t_mpimat)
allocate(t_mpimat::smat_unfold)
select type(smat_unfold)
type is (t_mpimat)
smat_unfold=smat
smat_unfold%data_c=CONJG(smat%data_c)
end select
end select
END IF
......@@ -235,9 +241,6 @@ CONTAINS
CALL MPI_ALLREDUCE(ne_found,ne_all,1,MPI_INTEGER,MPI_SUM, fmpi%diag_sub_comm,ierr)
ne_all=MIN(fi%input%neig,ne_all)
#endif
IF (.NOT.zMat%l_real) THEN
zMat%data_c(:lapw%nmat,:ne_found) = CONJG(zMat%data_c(:lapw%nmat,:ne_found))
END IF
endif
IF (fmpi%n_rank == 0) THEN
......
......@@ -70,7 +70,7 @@ CONTAINS
!Generate interstitial part of Hamiltonian
ALLOCATE(vpw_wTemp(SIZE(v%pw_w,1),SIZE(v%pw_w,2)))
vpw_wTemp = merge(v%pw_w - xcpot%get_exchange_weight() * vx%pw_w, v%pw_w, hybdat%l_subvxc)
CALL hs_int(fi%input, fi%noco, nococonv, stars, lapw, fmpi, fi%cell, isp, vpw_wTemp, smat, hmat)
CALL hs_int(fi%input, fi%noco, nococonv, stars, lapw, fmpi, fi%cell%bbmat, isp, vpw_wTemp, smat, hmat)
DEALLOCATE(vpw_wTemp)
CALL timestop("Interstitial part")
......@@ -103,8 +103,9 @@ CONTAINS
if (any(shape(smat) /= 1)) then
call judft_error("Hybrid doesn't do noco.")
end if
smat(1,1)%data_c = CONJG(smat(1,1)%data_c)
CALL write_eig(hybdat%eig_id, nk, isp, smat=smat(1, 1), n_start=fmpi%n_size, n_end=fmpi%n_rank)
smat(1,1)%data_c = CONJG(smat(1,1)%data_c)
END IF
IF (fi%hybinp%l_hybrid) THEN
......
This diff is collapsed.
......@@ -7,8 +7,18 @@
MODULE m_hs_int
CONTAINS
!Subroutine to construct the interstitial Hamiltonian and overlap matrix
SUBROUTINE hs_int(input, noco, nococonv, stars, lapw, fmpi, cell, isp, vpw, &
& smat,hmat)
SUBROUTINE hs_int(input, noco, nococonv, stars, lapw, fmpi, bbmat, isp, vpw, &
& smat, hmat)
! Control subroutine for the calculation of Hamiltonian/overlap matrix
! elements in the interstitial. The spin logic and case selections are
! found here while the actual calculation loop is one layer deeper in
! hs_int_direct.
!
! <\Phi'|H/S|\Phi>
!
! All primed indices are to be understood as corresponding to the Bra
! basis function, e.g. iSpinPr is the left and iSpin the right hand spin.
USE m_types
USE m_hs_int_direct
......@@ -18,14 +28,14 @@ CONTAINS
TYPE(t_noco), INTENT(IN) :: noco
TYPE(t_nococonv), INTENT(IN) :: nococonv
TYPE(t_stars), INTENT(IN) :: stars
TYPE(t_cell), INTENT(IN) :: cell
REAL, INTENT(IN) :: bbmat(3, 3)
TYPE(t_lapw), INTENT(IN) :: lapw
TYPE(t_mpi), INTENT(IN) :: fmpi
INTEGER, INTENT(IN) :: isp
COMPLEX, INTENT(IN) :: vpw(:,:)
CLASS(t_mat), INTENT(INOUT) :: smat(:,:),hmat(:,:)
INTEGER :: ispin, jspin, iispin, jjspin
INTEGER :: iSpinPr, iSpin, igSpin, igSpinPr
INTEGER :: iTkin, fact, iQss
LOGICAL :: l_smat
......@@ -35,24 +45,24 @@ CONTAINS
IF (noco%l_noco.AND.isp==2) RETURN !was done already
DO ispin=MERGE(1,isp,noco%l_noco),MERGE(2,isp,noco%l_noco)
iispin=MIN(ispin,SIZE(smat,1))
DO jspin=MERGE(1,isp,noco%l_noco),MERGE(2,isp,noco%l_noco)
jjspin=MIN(jspin,SIZE(smat,1))
IF (ispin.EQ.1.AND.jspin.EQ.2) THEN
DO iSpinPr=MERGE(1,isp,noco%l_noco),MERGE(2,isp,noco%l_noco)
igSpinPr=MIN(iSpinPr,SIZE(smat,1))
DO iSpin=MERGE(1,isp,noco%l_noco),MERGE(2,isp,noco%l_noco)
igSpin=MIN(iSpin,SIZE(smat,1))
IF (iSpinPr.EQ.1.AND.iSpin.EQ.2) THEN
vpw_temp = conjg(vpw(:, 3))
l_smat = .FALSE. ! Offdiagonal part --> No step function part.
iTkin = 0 ! Offdiagonal part --> No T part.
fact = -1 ! (12)-element --> (-1) prefactor
iQss = 0 ! No spin-spiral considered (no T).
ELSE IF (ispin.EQ.2.AND.jspin.EQ.1) THEN
ELSE IF (iSpinPr.EQ.2.AND.iSpin.EQ.1) THEN
vpw_temp = vpw(:, 3)
l_smat = .FALSE.
iTkin = 0
fact = 1
iQss = 0
ELSE
vpw_temp = vpw(:, ispin)
vpw_temp = vpw(:, iSpin)
l_smat = .TRUE.
IF (input%l_useapw) THEN
iTkin = 1 ! Dirac form.
......@@ -63,9 +73,9 @@ CONTAINS
END IF
fact = 1
END IF
CALL hs_int_direct(fmpi, lapw%gvec(:,:,ispin), lapw%gvec(:,:,jspin), &
& lapw%bkpt+iQss*(2*ispin - 3)/2.0*nococonv%qss, lapw%bkpt+iQss*(2*jspin - 3)/2.0*nococonv%qss, &
& lapw%nv(ispin), lapw%nv(jspin), stars, cell, vpw_temp, hmat(jjspin,iispin), smat(jjspin,iispin), l_smat, .FALSE., iTkin, fact)
CALL hs_int_direct(fmpi, stars, bbmat, lapw%gvec(:,:,iSpinPr), lapw%gvec(:,:,iSpin), &
& lapw%bkpt+iQss*(2*iSpinPr - 3)/2.0*nococonv%qss, lapw%bkpt+iQss*(2*iSpin - 3)/2.0*nococonv%qss, &
& lapw%nv(iSpinPr), lapw%nv(iSpin), iTkin, fact, l_smat, .FALSE., vpw_temp, hmat(igSpinPr,igSpin), smat(igSpinPr,igSpin))
END DO
END DO
END SUBROUTINE hs_int
......
......@@ -6,100 +6,99 @@
MODULE m_hs_int_direct
CONTAINS
SUBROUTINE hs_int_direct(fmpi, gvec, gvecPr, kvec, kvecPr, nv, nvPr, stars, cell, vpw, hmat, smat, l_smat, l_fullj, iTkin, fact)
! Calculates matrix elements of the form
! <\phi_{k'G'}|M|\phi_{kG}>
! for different use cases in the DFT/DFPT scf loop and operators M.
!
! The spin loop and distinction is found on a higher level and translates
! into new switches.
!
! DFT:
! M = \Theta_{IR} * (T + V) goes into hmat, M = \Theta_{IR} into smat.
! [vpw = \Theta_{IR} * V]
! [stars%ustep = \Theta_{IR}]
! [l_smat = F for offdiags, l_fullj = F]
! [iTkin = 0 for offdiags, 1 for l_useapw, 2 else]
!
! DFPT:
! M = \Theta_{IR}^{(1)} * (T + V) + \Theta_{IR} * V^{(1)} goes into hmat,
! M = \Theta_{IR}^{(1)} into smat.
! [vpw = \Theta_{IR}^{(1)} * V + \Theta_{IR} * V^{(1)}]
! [stars%ustep = \Theta_{IR}^{(1)}]
! [l_smat = F for offdiags, l_fullj = T]
! [iTkin = 0 for offdiags, 1 else]
USE m_types
IMPLICIT NONE
TYPE(t_stars),INTENT(IN) :: stars
TYPE(t_cell),INTENT(IN) :: cell
INTEGER ,INTENT(IN) :: gvecPr(:, :), gvec(:, :)
INTEGER, INTENT(IN) :: nvPr, nv, iTkin, fact
REAL, INTENT(IN) :: kvecPr(3), kvec(3)
TYPE(t_mpi),INTENT(IN) :: fmpi
COMPLEX,INTENT(IN) :: vpw(:)
CLASS(t_mat),INTENT(INOUT) :: hmat, smat
LOGICAL, INTENT(IN) :: l_smat, l_fullj
INTEGER :: i, j, i0, jmax, gPrG(3)
INTEGER :: gInd
COMPLEX :: th, ts, phase
REAL :: bvecPr(3), bvec(3), r2
!$OMP PARALLEL DO SCHEDULE(dynamic) DEFAULT(none) &
!$OMP SHARED(fmpi, stars, cell, vpw, gvecPr, gvec, kvecPr, kvec) &
!$OMP SHARED(nvPr, nv, l_smat, l_fullj, iTkin, fact)&
!$OMP SHARED(hmat, smat)&
!$OMP PRIVATE(gPrG, i0, i, j, jmax, gInd, phase, bvecPr, bvec, r2, th, ts)
DO i = fmpi%n_rank + 1, nv, fmpi%n_size
i0 = (i-1) / fmpi%n_size + 1
jmax = MERGE(nvPr, MIN(i, nvPr), l_fullj)
DO j = 1, jmax
gPrG = fact * (gvec(:,i) - gvecPr(:,j))
gInd = stars%ig(gPrG(1), gPrG(2), gPrG(3))
IF (gInd.EQ.0) CYCLE
phase = stars%rgphs(gPrG(1), gPrG(2), gPrG(3))
th = phase * vpw(gInd)
IF (iTkin.GT.0) THEN
bvecPr = kvecPr + gvecPr(:,j)
bvec = kvec + gvec(:,i)
IF (iTkin.EQ.1) THEN ! Symmetric Dirac form
r2 = 0.5 * DOT_PRODUCT(MATMUL(bvecPr,cell%bbmat),bvec)
ELSE IF (iTkin.EQ.2) THEN ! Symmetrized Laplace form
r2 = 0.25 * DOT_PRODUCT(MATMUL(bvecPr,cell%bbmat),bvecPr)
r2 = r2 + 0.25 * DOT_PRODUCT(MATMUL(bvec,cell%bbmat),bvec)
! Old form:
! 0.25* (rk(i)**2+rkPr(j)**2); rk(Pr)=lapw(Pr)%rk
ELSE ! Pure Laplace form
r2 = 0.5 * DOT_PRODUCT(MATMUL(bvec,cell%bbmat),bvec)
END IF
th = th + phase * r2 * stars%ustep(gInd)
END IF
IF (l_smat) THEN
ts = phase*stars%ustep(gInd)
ELSE
ts = 0.0
END IF
IF (hmat%l_real) THEN
hmat%data_r(j,i0) = REAL(th)
smat%data_r(j,i0) = REAL(ts)
ELSE
hmat%data_c(j,i0) = th
smat%data_c(j,i0) = ts
END IF
END DO
END DO
!$OMP END PARALLEL DO
END SUBROUTINE hs_int_direct
SUBROUTINE hs_int_direct(fmpi, stars, bbmat, gvecPr, gvec, kvecPr, kvec, nvPr, nv, &
& iTkin, fact, l_smat, l_fullj, vpw, hmat, smat)
! Calculates matrix elements of the form
! <\phi_{k'G'}|M|\phi_{kG}>
! for different use cases in the DFT/DFPT scf loop and operators M.
!
! The spin loop and distinction is found on a higher level and translates
! into new switches.
!
! DFT:
! M = \Theta_{IR} * (T + V) goes into hmat, M = \Theta_{IR} into smat.
! [vpw = \Theta_{IR} * V]
! [stars%ustep = \Theta_{IR}]
! [l_smat = F for offdiags, l_fullj = F]
! [iTkin = 0 for offdiags, 1 for l_useapw, 2 else]
!
! DFPT:
! M = \Theta_{IR}^{(1)} * (T + V) + \Theta_{IR} * V^{(1)} goes into hmat,
! M = \Theta_{IR}^{(1)} into smat.
! [vpw = \Theta_{IR}^{(1)} * V + \Theta_{IR} * V^{(1)}]
! [stars%ustep = \Theta_{IR}^{(1)}]
! [l_smat = F for offdiags, l_fullj = T]
! [iTkin = 0 for offdiags, 1 else]
USE m_types
IMPLICIT NONE
TYPE(t_mpi), INTENT(IN) :: fmpi
TYPE(t_stars), INTENT(IN) :: stars
REAL, INTENT(IN) :: bbmat(3, 3)
INTEGER, INTENT(IN) :: gvecPr(:, :), gvec(:, :)
REAL, INTENT(IN) :: kvecPr(3), kvec(3)
INTEGER, INTENT(IN) :: nvPr, nv, iTkin, fact
LOGICAL, INTENT(IN) :: l_smat, l_fullj
COMPLEX, INTENT(IN) :: vpw(:)
CLASS(t_mat), INTENT(INOUT) :: hmat, smat
INTEGER :: ikGPr, ikG, ikG0, gPrG(3), gInd
COMPLEX :: th, ts, phase
REAL :: bvecPr(3), bvec(3), r2
!$OMP PARALLEL DO SCHEDULE(dynamic) DEFAULT(none) &
!$OMP SHARED(fmpi, stars, bbmat, gvecPr, gvec, kvecPr, kvec) &
!$OMP SHARED(nvPr, nv, iTkin, fact, l_smat, l_fullj, vpw, hmat, smat) &
!$OMP PRIVATE(ikGPr, ikG, ikG0, gPrG, gInd, th, ts, phase, bvecPr, bvec, r2)
DO ikG = fmpi%n_rank + 1, nv, fmpi%n_size
ikG0 = (ikG-1) / fmpi%n_size + 1
DO ikGPr = 1, MERGE(nvPr, MIN(ikG, nvPr), l_fullj)
gPrG = fact * (gvecPr(:, ikGPr) - gvec(:, ikG))
gInd = stars%ig(gPrG(1), gPrG(2), gPrG(3))
IF (gInd.EQ.0) CYCLE
phase = stars%rgphs(gPrG(1), gPrG(2), gPrG(3))
th = phase * vpw(gInd)
IF (iTkin.GT.0) THEN
bvecPr = kvecPr + gvecPr(:, ikGPr)
bvec = kvec + gvec(:, ikG)
IF (iTkin.EQ.1) THEN ! Symmetric Dirac form
r2 = 0.5 * DOT_PRODUCT(MATMUL(bvecPr, bbmat), bvec)
ELSE IF (iTkin.EQ.2) THEN ! Symmetrized Laplace form
r2 = 0.25 * DOT_PRODUCT(MATMUL(bvecPr, bbmat), bvecPr)
r2 = r2 + 0.25 * DOT_PRODUCT(MATMUL(bvec, bbmat), bvec)
! Old form:
! 0.25* (rk(i)**2+rkPr(j)**2); rk(Pr)=lapw(Pr)%rk
ELSE ! Pure Laplace form
r2 = 0.5 * DOT_PRODUCT(MATMUL(bvec, bbmat), bvec)
END IF
th = th + phase * r2 * stars%ustep(gInd)
END IF
IF (l_smat) THEN
ts = phase * stars%ustep(gInd)
ELSE
ts = 0.0
END IF
IF (hmat%l_real) THEN
hmat%data_r(ikGPr, ikG0) = REAL(th)
smat%data_r(ikGPr, ikG0) = REAL(ts)
ELSE
hmat%data_c(ikGPr, ikG0) = th ! This corresponds to the
smat%data_c(ikGPr, ikG0) = ts ! conjugation in eigen.
END IF
END DO
END DO
!$OMP END PARALLEL DO
END SUBROUTINE hs_int_direct
END MODULE m_hs_int_direct
!--------------------------------------------------------------------------------
! Copyright (c) 2022 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
! This file is part of FLEUR and available as free software under the conditions
! of the MIT license as expressed in the LICENSE file in more detail.
!--------------------------------------------------------------------------------
MODULE m_hs_int_onespin
CONTAINS
SUBROUTINE hs_int_onespin(input, fmpi, gvec1, gvec2, bkpt1, bkpt2, nv1, nv2, rk1, rk2, stars, ispin, jspin, cell, vpw, hmat, smat, l_dfpths)
USE m_types
IMPLICIT NONE
TYPE(t_input),INTENT(IN) :: input
TYPE(t_stars),INTENT(IN) :: stars
TYPE(t_cell),INTENT(IN) :: cell
INTEGER ,INTENT(IN) :: gvec1(:, :), gvec2(:, :)
INTEGER, INTENT(IN) :: nv1, nv2
REAL, INTENT(IN) :: bkpt1(3), bkpt2(3), rk1(:), rk2(:)
TYPE(t_mpi),INTENT(IN) :: fmpi
INTEGER,INTENT(IN) :: ispin, jspin
COMPLEX,INTENT(IN) :: vpw(:,:)
CLASS(t_mat),INTENT(INOUT) :: hmat, smat
LOGICAL, INTENT(IN) :: l_dfpths
INTEGER :: i, j, i0, jmax, ii(3)
INTEGER :: in
COMPLEX :: th,ts,phase
REAL :: b1(3),b2(3),r2
!$OMP PARALLEL DO SCHEDULE(dynamic) DEFAULT(none) &
!$OMP SHARED(fmpi,stars,input,cell,vpw,gvec1,gvec2,bkpt1,bkpt2,rk1,rk2) &
!$OMP SHARED(nv1,nv2,ispin,jspin,l_dfpths)&
!$OMP SHARED(hmat,smat)&
!$OMP PRIVATE(ii,i0,i,j,jmax,in,phase,b1,b2,r2,th,ts)
DO i = fmpi%n_rank+1,nv1,fmpi%n_size
i0=(i-1)/fmpi%n_size+1
!---> loop over (k+g)
jmax = MERGE(nv2,MIN(i,nv2),l_dfpths)
DO j = 1, jmax
ii = gvec1(:,i) - gvec2(:,j)
IF (ispin==1.AND.jspin==2.AND..NOT.l_dfpths) THEN
ii=-1*ii
in = stars%ig(ii(1),ii(2),ii(3))
IF (in.EQ.0) CYCLE
th = stars%rgphs(ii(1),ii(2),ii(3))*conjg(vpw(in,3))
ts=0.0
ELSEIF(ispin==2.and.jspin==1.AND..NOT.l_dfpths) THEN
! ii = -1*ii
in = stars%ig(ii(1),ii(2),ii(3))
IF (in.EQ.0) CYCLE
th = stars%rgphs(ii(1),ii(2),ii(3))*vpw(in,3)
ts=0.0
ELSE
!--> determine index and phase factor
in = stars%ig(ii(1),ii(2),ii(3))
IF (in.EQ.0) CYCLE
phase = stars%rgphs(ii(1),ii(2),ii(3))
ts = phase*stars%ustep(in)
IF (input%l_useapw.OR.l_dfpths) THEN
b1=bkpt1+gvec1(:,i)
b2=bkpt2+gvec2(:,j)
r2 = DOT_PRODUCT(MATMUL(b2,cell%bbmat),b1)
th = phase*(0.5*r2*stars%ustep(in)+vpw(in,ispin))
ELSE
th = phase* (0.25* (rk1(i)**2+rk2(j)**2)*stars%ustep(in) + vpw(in,ispin))
ENDIF
ENDIF
!---> determine matrix element and store
IF (hmat%l_real) THEN
hmat%data_r(j,i0) = REAL(th)
smat%data_r(j,i0) = REAL(ts)
else
hmat%data_c(j,i0) = th
smat%data_c(j,i0) = ts
endif
ENDDO
ENDDO
!$OMP END PARALLEL DO
END SUBROUTINE hs_int_onespin
END MODULE m_hs_int_onespin
......@@ -52,8 +52,8 @@ CONTAINS
!locals
TYPE(t_fjgj)::fjgj
INTEGER :: ispin,jspin !local spin in atom
INTEGER :: iintsp,jintsp,n
INTEGER :: ilSpinPr,ilSpin !local spin in atom
INTEGER :: igSpinPr,igSpin,n
COMPLEX :: chi(2,2),chi_one
CLASS(t_mat),ALLOCATABLE::smat_tmp,hmat_tmp
......@@ -74,61 +74,61 @@ CONTAINS
CALL fjgj%alloc(MAXVAL(lapw%nv),atoms%lmaxd,isp,noco)
!$acc data copyin(fjgj) create(fjgj%fj,fjgj%gj)
iintsp=1;jintsp=1;chi_one=1.0 !Defaults in non-noco case
igSpinPr=1;igSpin=1;chi_one=1.0 !Defaults in non-noco case
DO n=1,atoms%ntype
DO ispin=MERGE(1,isp,noco%l_noco),MERGE(2,isp,noco%l_noco)
DO ilSpinPr=MERGE(1,isp,noco%l_noco),MERGE(2,isp,noco%l_noco)
CALL timestart("fjgj coefficients")
CALL fjgj%calculate(input,atoms,cell,lapw,noco,usdus,n,ispin)
CALL fjgj%calculate(input,atoms,cell,lapw,noco,usdus,n,ilSpinPr)
!$acc update device(fjgj%fj,fjgj%gj)
CALL timestop("fjgj coefficients")
DO jspin=ispin,MERGE(2,isp,noco%l_noco)
DO ilSpin=ilSpinPr,MERGE(2,isp,noco%l_noco)
IF (.NOT.noco%l_noco) THEN
!This is for collinear calculations: the (1,1) element of the matrices is all
!that is needed and allocated
CALL hsmt_sph(n,atoms,fmpi,ispin,input,nococonv,1,1,chi_one,lapw,enpara%el0,td%e_shift(n,ispin),usdus,fjgj,smat(1,1),hmat(1,1),.FALSE.)
CALL hsmt_nonsph(n,fmpi,sym,atoms,ispin,jspin,1,1,chi_one,noco,nococonv,cell,lapw,td,fjgj,hmat(1,1),.FALSE.)
CALL hsmt_lo(input,atoms,sym,cell,fmpi,noco,nococonv,lapw,usdus,td,fjgj,n,chi_one,ispin,jspin,iintsp,jintsp,hmat(1,1),.FALSE.,smat(1,1))
CALL hsmt_sph(n,atoms,fmpi,ilSpinPr,input,nococonv,1,1,chi_one,lapw,enpara%el0,td%e_shift(n,ilSpinPr),usdus,fjgj,smat(1,1),hmat(1,1),.FALSE.,.FALSE.,cell%bmat,1)
CALL hsmt_nonsph(n,fmpi,sym,atoms,ilSpinPr,ilSpin,1,1,chi_one,noco,nococonv,cell,lapw,td,fjgj,hmat(1,1),.FALSE.)
CALL hsmt_lo(input,atoms,sym,cell,fmpi,noco,nococonv,lapw,usdus,td,fjgj,n,chi_one,ilSpinPr,ilSpin,igSpinPr,igSpin,hmat(1,1),.FALSE.,smat(1,1))
ELSEIF(noco%l_noco.AND..NOT.noco%l_ss) THEN
!The NOCO but non-spinspiral setup follows:
!The Matrix-elements are first calculated in the local frame of the atom and
!stored in tmp-variables. Then these are distributed (rotated) into the 2x2
!global spin-matrices.
IF (ispin==jspin) THEN !local spin-diagonal contribution