Skip to content

Commit

Permalink
new DI3 keyword
Browse files Browse the repository at this point in the history
  • Loading branch information
@aoterodelaroza
aoterodelaroza committed Dec 18, 2023
1 parent e6ef468 commit cb1e143
Show file tree
Hide file tree
Showing 3 changed files with 114 additions and 70 deletions.
179 changes: 109 additions & 70 deletions src/[email protected]
View file
Original file line number Diff line number Diff line change
Expand Up @@ -36,7 +36,7 @@
! subroutine read_fachk_body(fafname,fa)
! subroutine calc_sij_wannier(fid,wancut,useu,restart,imtype,natt1,iatt,ilvec,idg1,xattr,dat,luevc,luevc_ibnd,sij)
! subroutine calc_sij_psink(fid,imtype,natt1,iatt,ilvec,idg1,xattr,dat,sij)
! subroutine calc_fa_wannier(res,nmo,nbnd,nlat,nattr,nspin)
! subroutine calc_fa_wannier(res,nmo,nbnd,nlat,nattr,nspin,di3)
! subroutine calc_fa_psink(res,nmo,nbnd,nlat,nattr,nspin,kpt,occ)
! subroutine find_sij_translations(res,nmo,nbnd,nlat,nlattot)
! subroutine check_sij_sanity(res,nspin,nmo,nbnd,nlat,nlattot)
Expand Down Expand Up @@ -1591,7 +1591,7 @@ subroutine intgrid_deloc(bas,res)
cycle
end if
else if ((sy%propi(l)%itype == itype_deloc_wnr .or. sy%propi(l)%itype == itype_deloc_psink).and.&
sy%propi(l)%fachk) then
sy%propi(l)%fachk .and..not.sy%propi(l)%di3) then
if (read_chk_header(fafname,nbnd,nbndw,nlat,nmo,nlattot,nspin,natt1,isijtype)) then
fid = sy%propi(l)%fid

Expand Down Expand Up @@ -1749,6 +1749,9 @@ subroutine intgrid_deloc(bas,res)
elseif (sy%propi(l)%itype == itype_deloc_psink) then
write (uout,'(99(A," "))') " Calculating overlaps using Bloch states (psink)"
end if
if (sy%propi(l)%di3) then
write (uout,'(99(A," "))') " Calculating three-center delocalization indices"
end if
if (sy%propi(l)%sijrestart) &
write (uout,'(99(A," "))') " Restart file read/written: ", trim(sy%f(fid)%file) // "-sijrestart"

Expand Down Expand Up @@ -1832,6 +1835,12 @@ subroutine intgrid_deloc(bas,res)
call write_fachk(fafname,nbnd,nbndw,nlat,nmo,nlattot,nspin,bas%nattr,res(l)%sijtype,res(l)%fa)
end if

! calculate three-body DIs, if requested
if (sy%propi(l)%di3) then
write (uout,'("# Calculating three-atom delocalization indices")')
call calc_di3_wannier(res(l),nmo,nbnd,nlat,bas%nattr,nspin)
end if

! finished the Sij and Fa successfully
999 continue
if (allocated(res(l)%sijc)) deallocate(res(l)%sijc)
Expand Down Expand Up @@ -2486,11 +2495,9 @@ subroutine calc_fa_wannier(res,nmo,nbnd,nlat,nattr,nspin)
type(int_result), intent(inout) :: res
integer, intent(in) :: nmo, nbnd, nlat(3), nattr, nspin

integer :: iat_a, iat_b, iat_c, rat_b, rat_c
integer :: is, imo, jmo, kmo, ia, ja, ka, iba, ib, jb, kb, ibb
integer :: ic, jc, kc, icc, nlattot
integer :: iat_a, iat_b, rat_b, nlattot
integer :: is, imo, jmo, ia, ja, ka, iba, ib, jb, kb, ibb
real*8 :: fac
complex*16 :: f3temp

real*8, allocatable :: ftemp(:,:,:), f3temp(:,:,:,:,:)
logical, allocatable :: useovrlp(:,:,:)
Expand Down Expand Up @@ -2546,72 +2553,104 @@ subroutine calc_fa_wannier(res,nmo,nbnd,nlat,nattr,nspin)
!$omp end parallel do
end do

! ! three-atom indices
! allocate(res%fa3(nattr,nattr,nlattot,nattr,nlattot,nspin))
! allocate(f3temp(nattr,nattr,nlattot,nattr,nlattot))
! res%fa3 = 0d0
! f3temp = 0d0
! do is = 1, nspin

! !$omp parallel do private(ia,ja,ka,iba,ib,jb,kb,ibb,ic,jc,kc,icc,fac) firstprivate(f3temp)
! do imo = 1, nmo
! call unpackidx(imo,ia,ja,ka,iba,nmo,nbnd,nlat)
! do jmo = 1, nmo
! call unpackidx(jmo,ib,jb,kb,ibb,nmo,nbnd,nlat)
! if (.not.useovrlp(imo,jmo,is)) cycle
! do kmo = imo, nmo
! call unpackidx(jmo,ic,jc,kc,icc,nmo,nbnd,nlat)
! if (.not.useovrlp(imo,kmo,is)) cycle
! if (.not.useovrlp(jmo,kmo,is)) cycle

! f3temp = 0d0
! do iat_a = 1, nattr
! do iat_b = 1, nattr
! do iat_c = 1, nattr
! do rat_b = 1, nlattot
! do rat_c = 1, nlattot
! f3temp(iat_a,iat_b,rat_b,iat_c,rat_c) = &
! real(res%sijc(imo,kmo,iat_a,is) * &
! res%sijc(res%sij_wnr_imap(jmo,rat_b),res%sij_wnr_imap(imo,rat_b),iat_b,is) * &
! res%sijc(res%sij_wnr_imap(kmo,rat_c),res%sij_wnr_imap(jmo,rat_c),iat_c,is),8) + &
! real(res%sijc(imo,kmo,iat_a,is) * &
! res%sijc(res%sij_wnr_imap(kmo,rat_b),res%sij_wnr_imap(jmo,rat_b),iat_b,is) * &
! res%sijc(res%sij_wnr_imap(jmo,rat_c),res%sij_wnr_imap(imo,rat_c),iat_c,is),8)
! end do ! rat_c
! end do ! rat_b
! end do ! iat_c
! end do ! iat_b
! end do ! iat_a

! if (imo == kmo) then
! fac = 1d0
! else
! fac = 2d0
! end if

! !$omp critical (addfa)
! res%fa3(:,:,:,:,:,is) = res%fa3(:,:,:,:,:,is) + fac * f3temp
! !$omp end critical (addfa)
! end do ! kmo
! end do ! jmo
! end do ! imo
! !$omp end parallel do
! end do ! is
! res%fa3 = 0.5d0 * res%fa3

! do is = 1, nspin
! do iat_a = 1, nattr
! do iat_b = 1, nattr
! do rat_b = 1, nlattot
! write (*,*) "xxxx ", abs(res%fa(iat_a,iat_b,rat_b,is)-sum(res%fa3(iat_a,iat_b,rat_b,:,:,is)))
! end do
! end do
! end do
! end do
! stop 1

end subroutine calc_fa_wannier

!> Calculate the three-body DIs, wannier version
subroutine calc_di3_wannier(res,nmo,nbnd,nlat,nattr,nspin)
use types, only: int_result
type(int_result), intent(inout) :: res
integer, intent(in) :: nmo, nbnd, nlat(3), nattr, nspin

integer :: iat_a, iat_b, iat_c, rat_b, rat_c
integer :: is, imo, jmo, kmo, ia, ja, ka, iba, ib, jb, kb, ibb
integer :: ic, jc, kc, icc, nlattot
real*8 :: fac
complex*16 :: f3temp

real*8, allocatable :: f3temp(:,:,:,:,:)
logical, allocatable :: useovrlp(:,:,:)
real*8, parameter :: ovrlp_thr = 1d-50

! number of lattice vectors
nlattot = nlat(1)*nlat(2)*nlat(3)

! flag which overlaps will be used
allocate(useovrlp(nmo,nmo,nspin))
useovrlp = .true.
do is = 1, nspin
do imo = 1, nmo
do jmo = imo, nmo
useovrlp(imo,jmo,is) = (maxval(abs(res%sijc(imo,jmo,:,is))) >= ovrlp_thr)
end do
end do
end do

! three-atom indices
allocate(res%fa3(nattr,nattr,nlattot,nattr,nlattot,nspin))
allocate(f3temp(nattr,nattr,nlattot,nattr,nlattot))
res%fa3 = 0d0
f3temp = 0d0
do is = 1, nspin

!$omp parallel do private(ia,ja,ka,iba,ib,jb,kb,ibb,ic,jc,kc,icc,fac) firstprivate(f3temp)
do imo = 1, nmo
call unpackidx(imo,ia,ja,ka,iba,nmo,nbnd,nlat)
do jmo = 1, nmo
call unpackidx(jmo,ib,jb,kb,ibb,nmo,nbnd,nlat)
if (.not.useovrlp(imo,jmo,is)) cycle
do kmo = imo, nmo
call unpackidx(jmo,ic,jc,kc,icc,nmo,nbnd,nlat)
if (.not.useovrlp(imo,kmo,is)) cycle
if (.not.useovrlp(jmo,kmo,is)) cycle

f3temp = 0d0
do iat_a = 1, nattr
do iat_b = 1, nattr
do iat_c = 1, nattr
do rat_b = 1, nlattot
do rat_c = 1, nlattot
f3temp(iat_a,iat_b,rat_b,iat_c,rat_c) = &
real(res%sijc(imo,kmo,iat_a,is) * &
res%sijc(res%sij_wnr_imap(jmo,rat_b),res%sij_wnr_imap(imo,rat_b),iat_b,is) * &
res%sijc(res%sij_wnr_imap(kmo,rat_c),res%sij_wnr_imap(jmo,rat_c),iat_c,is),8) + &
real(res%sijc(imo,kmo,iat_a,is) * &
res%sijc(res%sij_wnr_imap(kmo,rat_b),res%sij_wnr_imap(jmo,rat_b),iat_b,is) * &
res%sijc(res%sij_wnr_imap(jmo,rat_c),res%sij_wnr_imap(imo,rat_c),iat_c,is),8)
end do ! rat_c
end do ! rat_b
end do ! iat_c
end do ! iat_b
end do ! iat_a

if (imo == kmo) then
fac = 1d0
else
fac = 2d0
end if

!$omp critical (addfa)
res%fa3(:,:,:,:,:,is) = res%fa3(:,:,:,:,:,is) + fac * f3temp
!$omp end critical (addfa)
end do ! kmo
end do ! jmo
end do ! imo
!$omp end parallel do
end do ! is
res%fa3 = 0.5d0 * res%fa3

do is = 1, nspin
do iat_a = 1, nattr
do iat_b = 1, nattr
do rat_b = 1, nlattot
write (*,*) "xxxx ", abs(res%fa(iat_a,iat_b,rat_b,is)-sum(res%fa3(iat_a,iat_b,rat_b,:,:,is)))
end do
end do
end do
end do
stop 1

end subroutine calc_di3_wannier

!> Calculate the Fa matrix from the Sij matrix, psink version
subroutine calc_fa_psink(res,nmo,nbnd,nlat,nattr,nspin,kpt,occ)
use types, only: int_result
Expand Down
4 changes: 4 additions & 0 deletions src/[email protected]
View file
Original file line number Diff line number Diff line change
Expand Up @@ -960,6 +960,7 @@ module subroutine new_integrable_string(s,line,errmsg)
s%propi(s%npropi)%wancut = 4d0
s%propi(s%npropi)%sijchkfile = ""
s%propi(s%npropi)%fachkfile = ""
s%propi(s%npropi)%di3 = .false.

do while (.true.)
lp2 = lp
Expand All @@ -983,11 +984,14 @@ module subroutine new_integrable_string(s,line,errmsg)
s%npropi = s%npropi - 1
return
end if
else if (equal(word,"di3")) then
s%propi(s%npropi)%di3 = .true.
else
lp = lp2
exit
end if
end do

elseif (equal(word,"name")) then
word = getword(line,lp)
s%propi(s%npropi)%prop_name = string(word)
Expand Down
1 change: 1 addition & 0 deletions src/types.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -222,6 +222,7 @@ module types
logical :: fachk = .true. ! read and write the fachk file
logical :: sijrestart = .false. ! read and write the sij restart file
real*8 :: wancut = 4d0 ! Wannier center distance cutoff
logical :: di3 ! calculate the 3-body DIs
character(len=mlen) :: sijchkfile = "" ! name of sijchk file
character(len=mlen) :: fachkfile = "" ! name of fachk file
end type integrable
Expand Down

0 comments on commit cb1e143

Please sign in to comment.