Commit 88960a9a authored by Alexander Neukirchen's avatar Alexander Neukirchen
Browse files

DFPT eigenstuff

parent e3de0257
......@@ -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
......
......@@ -137,7 +137,7 @@ CONTAINS
& abCoeffs, SIZE(abCoeffs, 1), h_loc, size(td%h_loc_nonsph, 1), &
& cmplx(0.0, 0.0), ab1, size_ab)
!$acc end host_data
ELSE ! Needed, because t^H .NE. t! TODO: Do we possibly need to do this for fully fully noco as well?
ELSE ! Needed, because t^H .NE. t!
!!$acc update device(ab)
!$acc host_data use_device(abCoeffs,ab1,h_loc)
CALL CPP_zgemm("C", "C", lapw%nv(igSpin), ab_size, ab_size, cmplx(1.0, 0.0), &
......@@ -274,7 +274,7 @@ CONTAINS
! data_c += chi * ab2 * ab_select^H
! = chi * aPr'^H * H * a
ELSE
!$acc host_data use_device(abCoeffs,ab1,data_c,ab_select)
!$acc host_data use_device(abCoeffs,data_c,ab_select)
IF (set0 .AND. nn == 1) THEN
!CPP_data_c = CMPLX(0.0,0.0)
CALL CPP_zgemm("C", "C", lapwPr%nv(igSpinPr), lapwPr%num_local_cols(igSpin), ab_size, cchi, &
......
......@@ -125,7 +125,7 @@ CONTAINS
ddnln = usdus%ddn(l,n,isp)
elall = el(l,n,isp)
IF (l<=atoms%lnonsph(n)) elall=elall-e_shift!(isp)
IF (l<=atoms%lnonsph(n).AND..NOT.l_fullj) elall=elall-e_shift!(isp)
! Legendre polynomials
l3 = modulo(l, 3)
......@@ -139,10 +139,14 @@ CONTAINS
& - fleg2(l-1)*plegend(modulo(l-2,3))
END IF ! l
fct = plegend(l3)*fl2p1(l) * ( fjkiln*fjgj%fj(ikGPr,l,isp,igSpinPr) &
& + gjkiln*fjgj%gj(ikGPr,l,isp,igSpinPr)*ddnln )
fct2 = plegend(l3)*fl2p1(l) * 0.5 * ( gjkiln*fjgj%fj(ikGPr,l,isp,igSpinPr) &
& + fjkiln*fjgj%gj(ikGPr,l,isp,igSpinPr) )
fct = plegend(l3) * fl2p1(l) * ( fjkiln*fjgj%fj(ikGPr,l,isp,igSpinPr) &
& + gjkiln*fjgj%gj(ikGPr,l,isp,igSpinPr)*ddnln )
IF (.NOT.l_fullj) THEN
fct2 = plegend(l3)*fl2p1(l) * 0.5 * ( gjkiln*fjgj%fj(ikGPr,l,isp,igSpinPr) &
& + fjkiln*fjgj%gj(ikGPr,l,isp,igSpinPr) )
ELSE
fct2 = plegend(l3)*fl2p1(l) * gjkiln*fjgj%fj(ikGPr,l,isp,igSpinPr)
END IF
VecHelpS = VecHelpS + fct
VecHelpH = VecHelpH + fct*elall + fct2
......@@ -325,7 +329,7 @@ CONTAINS
ddnln = usdus%ddn(l,n,isp)
elall = el(l,n,isp)
IF (l<=atoms%lnonsph(n)) elall=elall-e_shift!(isp)
IF (l<=atoms%lnonsph(n).AND..NOT.l_fullj) elall=elall-e_shift!(isp)
! Legendre polynomials
l3 = modulo(l, 3)
......@@ -338,10 +342,14 @@ CONTAINS
& - fleg2(l-1)*plegend(:NVEC_REM,modulo(l-2,3))
END IF ! l
fct(:NVEC_REM) = plegend(:NVEC_REM,l3)*fl2p1(l) * ( fjkiln*fjgj%fj(kj_off:kj_vec,l,isp,igSpinPr) &
& + gjkiln*fjgj%gj(kj_off:kj_vec,l,isp,igSpinPr)*ddnln )
fct2(:NVEC_REM) = plegend(:NVEC_REM,l3)*fl2p1(l) * 0.5 * ( gjkiln*fjgj%fj(kj_off:kj_vec,l,isp,igSpinPr) &
& + fjkiln*fjgj%gj(kj_off:kj_vec,l,isp,igSpinPr) )
fct(:NVEC_REM) = plegend(:NVEC_REM,l3) * fl2p1(l) * ( fjkiln*fjgj%fj(kj_off:kj_vec,l,isp,igSpinPr) &
& + gjkiln*fjgj%gj(kj_off:kj_vec,l,isp,igSpinPr)*ddnln )
IF (.NOT.l_fullj) THEN
fct2(:NVEC_REM) = plegend(:NVEC_REM,l3) * fl2p1(l) * 0.5 * ( gjkiln*fjgj%fj(kj_off:kj_vec,l,isp,igSpinPr) &
& + fjkiln*fjgj%gj(kj_off:kj_vec,l,isp,igSpinPr) )
ELSE
fct2(:NVEC_REM) = plegend(:NVEC_REM,l3) * fl2p1(l) * gjkiln*fjgj%fj(kj_off:kj_vec,l,isp,igSpinPr)
END IF
VecHelpS(:NVEC_REM) = VecHelpS(:NVEC_REM) + fct(:NVEC_REM)
VecHelpH(:NVEC_REM) = VecHelpH(:NVEC_REM) + fct(:NVEC_REM)*elall + fct2(:NVEC_REM)
......
......@@ -7,7 +7,7 @@
MODULE m_mt_setup
CONTAINS
SUBROUTINE mt_setup(atoms,sym,sphhar,input,noco,nococonv,enpara,hub1inp,hub1data,inden,vTot,vx,fmpi,td,ud,alpha_hybrid)
SUBROUTINE mt_setup(atoms,sym,sphhar,input,noco,nococonv,enpara,hub1inp,hub1data,inden,vTot,vx,fmpi,td,ud,alpha_hybrid,l_dfptmod)
USE m_types
USE m_tlmplm_cholesky
USE m_spnorb
......@@ -29,6 +29,7 @@ CONTAINS
TYPE(t_hub1inp),INTENT(IN) :: hub1inp
TYPE(t_hub1data),INTENT(INOUT)::hub1data
REAl,INTENT(IN) :: alpha_hybrid
LOGICAL,INTENT(IN),OPTIONAL :: l_dfptmod
INTEGER:: jsp
......@@ -38,7 +39,11 @@ CONTAINS
DO jsp=1,MERGE(4,input%jspins,any(noco%l_unrestrictMT).OR.any(noco%l_spinoffd_ldau))
!CALL tlmplm_cholesky(sphhar,atoms,DIMENSION,enpara, jsp,1,fmpi,vTot%mt(:,0,1,jsp),input,vTot%mmpMat, td,ud)
CALL tlmplm_cholesky(sphhar,atoms,sym,noco,nococonv,enpara,jsp,fmpi,vTot,vx,inDen,input,hub1inp,hub1data,td,ud,alpha_hybrid)
IF (PRESENT(l_dfptmod)) THEN
CALL tlmplm_cholesky(sphhar,atoms,sym,noco,nococonv,enpara,jsp,fmpi,vTot,vx,inDen,input,hub1inp,hub1data,td,ud,alpha_hybrid,l_dfptmod)
ELSE
CALL tlmplm_cholesky(sphhar,atoms,sym,noco,nococonv,enpara,jsp,fmpi,vTot,vx,inDen,input,hub1inp,hub1data,td,ud,alpha_hybrid)
END IF
END DO
CALL timestop("tlmplm")
......
......@@ -9,7 +9,7 @@ MODULE m_tlmplm_cholesky
!*********************************************************************
CONTAINS
SUBROUTINE tlmplm_cholesky(sphhar,atoms,sym,noco,nococonv,enpara,&
jspin,fmpi,v,vx,inden,input,hub1inp,hub1data,td,ud,alpha_hybrid)
jspin,fmpi,v,vx,inden,input,hub1inp,hub1data,td,ud,alpha_hybrid,l_dfptmod)
USE m_tlmplm
USE m_types
USE m_radovlp
......@@ -29,6 +29,7 @@ CONTAINS
! .. Scalar Arguments ..
INTEGER, INTENT (IN) :: jspin!physical spin&spin index for data
REAL,INTENT(IN) :: alpha_hybrid
LOGICAL,INTENT(IN),OPTIONAL :: l_dfptmod
! ..
TYPE(t_potden),INTENT(IN) :: v,vx,inden
TYPE(t_tlmplm),INTENT(INOUT) :: td
......@@ -39,7 +40,7 @@ CONTAINS
REAL temp
INTEGER i,l,lm,lmin,lmin0,lmp,lmplm,lp,info,in,jsp,j1,j2
INTEGER lpl ,mp,n,m,s,i_u,jmin,jmax,i_opc,lh0
LOGICAL OK, isRoot, l_call_tlmplm
LOGICAL OK, isRoot, l_call_tlmplm, l_dfpt
COMPLEX :: one
! ..
! .. Local Arrays ..
......@@ -49,6 +50,8 @@ CONTAINS
REAL,PARAMETER:: e_shift_min=0.5
REAL,PARAMETER:: e_shift_max=65.0
l_dfpt = PRESENT(l_dfptmod)
! ..e_shift
jsp=jspin
IF (jsp<3) THEN
......@@ -85,7 +88,7 @@ CONTAINS
!$OMP PRIVATE(lmplm,lp,m,mp,n)&
!$OMP PRIVATE(OK,s,in,info,i_u,i_opc)&
!$OMP SHARED(one,nococonv,atoms,jspin,jsp,sym,sphhar,enpara,td,ud,v,vx,alpha_hybrid,isRoot,l_call_tlmplm)&
!$OMP SHARED(fmpi,input,hub1inp,hub1data,uun21,udn21,dun21,ddn21,opc_corrections,j1,j2)
!$OMP SHARED(fmpi,input,hub1inp,hub1data,uun21,udn21,dun21,ddn21,opc_corrections,j1,j2,l_dfpt)
DO n = 1,atoms%ntype
IF (l_call_tlmplm) THEN
lh0 = MERGE(1,0,jsp<3.and.alpha_hybrid==0)
......@@ -143,8 +146,9 @@ CONTAINS
td%h_loc_nonsph(lm ,lm ,n,j1,j2) = td%h_loc_nonsph(lm ,lm ,n,j1,j2) + opc_corrections(i_opc) * m
td%h_loc_nonsph(lm+s,lm+s,n,j1,j2) = td%h_loc_nonsph(lm+s,lm+s,n,j1,j2) + opc_corrections(i_opc) * m * ud%ddn(l,n,jsp)
END DO
END DO
IF (jsp<3) THEN
END DO
! For DFPT, do not decompose
IF (jsp<3.AND..NOT.l_dfpt) THEN
!Create Cholesky decomposition of local hamiltonian
!---> Add diagonal terms to make matrix positive definite
......@@ -189,6 +193,11 @@ CONTAINS
td%h_loc(lm+s,lm+s,n,jsp,jsp)=td%h_loc(lm+s,lm+s,n,jsp,jsp) + enpara%el0(l,n,jsp)*ud%ddn(l,n,jsp)
td%h_loc(lm+s,lm,n,jsp,jsp)=td%h_loc(lm+s,lm,n,jsp,jsp) + 0.5
td%h_loc(lm,lm+s,n,jsp,jsp)=td%h_loc(lm,lm+s,n,jsp,jsp) + 0.5
! For DFPT we need a non-symmetrized local Hamiltonian
IF (l_dfpt) THEN
td%h_loc(lm+s,lm,n,jsp,jsp)=td%h_loc(lm+s,lm,n,jsp,jsp) - 0.5 !+0.0
td%h_loc(lm,lm+s,n,jsp,jsp)=td%h_loc(lm,lm+s,n,jsp,jsp) + 0.5 !+1.0
END IF
ENDDO
ENDDO
ENDIF
......
......@@ -3,7 +3,7 @@
! 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_hsmt
MODULE m_dfpt_hsmt
USE m_juDFT
IMPLICIT NONE
CONTAINS
......@@ -19,9 +19,11 @@ CONTAINS
!! - In the spin-spiral case, a loop over the global spin is performed and the four parts of the matrix are calculated one-by-one
!! @todo
!! The off-diagonal contribution in first-variation soc and constraint calculations is still missing
! DFPT: Handle H/S elements with shifted k+q on the lhs i(k+G-k'-G'-q) prefactor through sph,
! V1 elements fully with nonsph. LO?
SUBROUTINE hsmt(atoms,sym,enpara,&
isp,input,fmpi,noco,nococonv,cell,lapw,usdus,td,smat,hmat)
SUBROUTINE dfpt_hsmt(atoms,sym,enpara,&
isp,input,fmpi,noco,nococonv,cell,lapw,usdus,td,smat,hmat,tdV1)
USE m_types
USE m_types_mpimat
USE m_hsmt_nonsph
......@@ -43,7 +45,7 @@ CONTAINS
TYPE(t_atoms),INTENT(IN) :: atoms
TYPE(t_enpara),INTENT(IN) :: enpara
TYPE(t_lapw),INTENT(IN) :: lapw
TYPE(t_tlmplm),INTENT(IN) :: td
TYPE(t_tlmplm),INTENT(IN) :: td, tdV1
TYPE(t_usdus),INTENT(IN) :: usdus
CLASS(t_mat),INTENT(INOUT) :: smat(:,:),hmat(:,:)
! ..
......@@ -83,13 +85,10 @@ CONTAINS
CALL timestop("fjgj coefficients")
DO jspin=ispin,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))
ELSEIF(noco%l_noco.AND..NOT.noco%l_ss) THEN
ELSE IF(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
......@@ -118,10 +117,10 @@ CONTAINS
ENDDO
END DO
!$acc end data
if (noco%l_noco.AND..NOT.noco%l_ss) then
if (noco%l_noco) then
!$acc exit data delete(smat_tmp%data_c,smat_tmp%data_r,hmat_tmp%data_c,hmat_tmp%data_r)
!$acc exit data delete(smat_tmp,hmat_tmp)
endif
RETURN
END SUBROUTINE hsmt
END MODULE m_hsmt
END SUBROUTINE dfpt_hsmt
END MODULE m_dfpt_hsmt
......@@ -7,7 +7,7 @@
MODULE m_dfpt_tlmplm
CONTAINS
SUBROUTINE dfpt_tlmplm(atoms,sym,sphhar,input,noco,enpara,hub1inp,hub1data,vTot,fmpi,td,v1real,v1imag)
SUBROUTINE dfpt_tlmplm(atoms,sym,sphhar,input,noco,enpara,hub1inp,hub1data,vTot,fmpi,tdV1,v1real,v1imag)
USE m_types
USE m_tlmplm
......@@ -17,33 +17,33 @@ CONTAINS
TYPE(t_enpara),INTENT(IN) :: enpara
TYPE(t_input),INTENT(IN) :: input
TYPE(t_noco),INTENT(IN) :: noco
TYPE(t_sym),INTENT(IN) :: sym
TYPE(t_sphhar),INTENT(IN) :: sphhar
TYPE(t_atoms),INTENT(IN) :: atoms
TYPE(t_potden),INTENT(IN) :: vTot
TYPE(t_tlmplm),INTENT(INOUT) :: td
TYPE(t_hub1inp),INTENT(IN) :: hub1inp
TYPE(t_hub1data),INTENT(INOUT)::hub1data
TYPE(t_sym),INTENT(IN) :: sym
TYPE(t_sphhar),INTENT(IN) :: sphhar
TYPE(t_atoms),INTENT(IN) :: atoms
TYPE(t_potden),INTENT(IN) :: vTot
TYPE(t_tlmplm),INTENT(INOUT) :: tdV1
TYPE(t_hub1inp),INTENT(IN) :: hub1inp
TYPE(t_hub1data),INTENT(INOUT)::hub1data
TYPE(t_potden), INTENT(IN) :: v1real, v1imag
TYPE(t_potden), INTENT(IN) :: v1real, v1imag
INTEGER :: iSpinV1, iSpinPr, iSpin, iPart, n
COMPLEX :: one
INTEGER :: iSpinV1, iSpinPr, iSpin, iPart, n, offs
COMPLEX :: one
REAL, ALLOCATABLE :: vr1(:, :)
REAL, ALLOCATABLE :: vr1(:, :)
TYPE(t_usdus) :: uddummy
TYPE(t_potden) :: vxdummy
TYPE(t_nococonv) :: nococonvdummy
TYPE(t_usdus) :: uddummy
TYPE(t_potden) :: vxdummy
TYPE(t_nococonv) :: nococonvdummy
ALLOCATE( vr1(SIZE(v1real%mt,1),0:SIZE(v1real%mt,2)-1))
call uddummy%init(atoms,input%jspins)
CALL timestart("tlmplm")
CALL td%init(atoms,input%jspins,.FALSE.)
CALL tdV1%init(atoms,input%jspins,.FALSE.)
!$OMP PARALLEL DO DEFAULT(NONE)&
!$OMP PRIVATE(n,one,iSpinV1,iSpinPr,iSpin,vr1)&
!$OMP SHARED(noco,nococonvdummy,atoms,sym,sphhar,enpara,td,uddummy,vTot,vxdummy,v1real,v1imag)&
!$OMP PRIVATE(n,one,iSpinV1,iSpinPr,iSpin,vr1,offs)&
!$OMP SHARED(noco,nococonvdummy,atoms,sym,sphhar,enpara,tdV1,uddummy,vTot,vxdummy,v1real,v1imag)&
!$OMP SHARED(fmpi,input,hub1inp,hub1data)
DO n = 1,atoms%ntype
DO iSpinV1 = 1, MERGE(4, input%jspins, any(noco%l_unrestrictMT))
......@@ -55,9 +55,15 @@ CONTAINS
IF (iPart.EQ.1) one = CMPLX(1.0, 0.0); vr1 = v1real%mt(:, :, n, iSpinV1)
IF (iPart.EQ.2) one = CMPLX(0.0, 1.0); vr1 = v1imag%mt(:, :, n, iSpinV1)
CALL tlmplm(n, sphhar, atoms, sym, enpara, nococonvdummy, iSpinPr, iSpin, iSpinV1, fmpi, &
& vTot, vxdummy, input, hub1inp, hub1data, td, uddummy, 0.0, 0, one, vr1)
& vTot, vxdummy, input, hub1inp, hub1data, tdV1, uddummy, 0.0, 0, one, vr1)
END DO
END DO
offs = tdV1%h_loc2_nonsph(n)
tdV1%h_loc_nonsph(0:offs-1,0:offs-1,n,:,:) = tdV1%h_loc(0:offs-1,0:offs-1,n,:,:)
tdV1%h_loc_nonsph(offs:offs+offs-1,0:offs-1,n,:,:) = tdV1%h_loc(tdV1%h_loc2(n):offs+tdV1%h_loc2(n)-1,0:offs-1,n,:,:)
tdV1%h_loc_nonsph(0:offs-1,offs:offs+offs-1,n,:,:) = tdV1%h_loc(0:offs-1,tdV1%h_loc2(n):offs+tdV1%h_loc2(n)-1,n,:,:)
tdV1%h_loc_nonsph(offs:offs+offs-1,offs:offs+offs-1,n,:,:)= tdV1%h_loc(tdV1%h_loc2(n):offs+tdV1%h_loc2(n)-1,tdV1%h_loc2(n):offs+tdV1%h_loc2(n)-1,n,:,:)
END DO
!$OMP END PARALLEL DO
CALL timestop("tlmplm")
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment