Skip to content

Commit

Permalink
clean up 3c(in) and 3c(out) output
Browse files Browse the repository at this point in the history
  • Loading branch information
@aoterodelaroza
aoterodelaroza committed Jan 29, 2024
1 parent 7fd3f94 commit 445534b
Showing 1 changed file with 104 additions and 72 deletions.
176 changes: 104 additions & 72 deletions src/[email protected]
View file
Original file line number Diff line number Diff line change
Expand Up @@ -675,12 +675,12 @@ module subroutine int_output_header(bas,res,nomol0,usesym0)
itaux = string(fid)
elseif (sy%propi(i)%itype == itype_deloc_wnr .or. sy%propi(i)%itype == itype_deloc_psink) then
if (res(i)%done .and. sy%propi(i)%di3) then
if (sy%propi(i)%di3_atom1 < 0) then
if (sy%propi(i)%di3_atom1 < 0 .or. sy%propi(i)%di3_atom1 > bas%nattr) then
saux = " 3-center indices calculated "
elseif (sy%propi(i)%di3_atom2(1) < 0) then
saux = " 3-center indices calculated for atom " // string(sy%propi(i)%di3_atom1) // " "
elseif (sy%propi(i)%di3_atom2(1) < 0 .or. sy%propi(i)%di3_atom2(1) > bas%nattr) then
saux = " 3-center indices calculated for attractor " // string(sy%propi(i)%di3_atom1) // " "
else
saux = " 3-center indices calculated for atom pair " // string(sy%propi(i)%di3_atom1) //&
saux = " 3-center indices calculated for attractor pair " // string(sy%propi(i)%di3_atom1) //&
"," // string(sy%propi(i)%di3_atom2(1)) // "+(" //&
string(sy%propi(i)%di3_atom2(2)) // "," //&
string(sy%propi(i)%di3_atom2(3)) // "," //&
Expand Down Expand Up @@ -2584,7 +2584,7 @@ subroutine calc_di3_wannier(res,nmo,nbnd,nlat,nattr,nspin,atom1,atom2)
integer, intent(in) :: nmo, nbnd, nlat(3), nattr, nspin
integer, intent(in) :: atom1, atom2(4)

integer :: i, j, k
integer :: i, j, k, n
integer :: iat_a, iat_b, iat_c, rat_b, rat_c
integer :: iat_a_ini, iat_a_end, iat_b_ini, iat_b_end, rat_b_ini, rat_b_end
integer :: is, imo, jmo, kmo, ia, ja, ka, iba, ib, jb, kb, ibb
Expand All @@ -2596,39 +2596,39 @@ subroutine calc_di3_wannier(res,nmo,nbnd,nlat,nattr,nspin,atom1,atom2)
logical, allocatable :: useovrlp(:,:,:)
real*8, parameter :: ovrlp_thr = 1d-50

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

! set limits for the sums
iat_a_ini = 1
iat_a_end = nattr
if (atom1 >= 0) then
if (atom1 > 0 .and. atom1 <= nattr) then
iat_a_ini = atom1
iat_a_end = atom1
end if
iat_b_ini = 1
iat_b_end = nattr
rat_b_ini = 1
rat_b_end = nlattot
if (atom2(1) >= 0) then
if (atom2(1) > 0 .and. atom2(1) <= nattr) then
iat_b_ini = atom2(1)
iat_b_end = atom2(1)
nlattot = 0
n = 0
loopi: do i = 0, nlat(1)-1
do j = 0, nlat(2)-1
do k = 0, nlat(3)-1
nlattot = nlattot + 1
n = n + 1
if (i == modulo(atom2(2),nlat(1)) .and. j == modulo(atom2(3),nlat(2)) .and.&
k == modulo(atom2(4),nlat(3))) then
rat_b_ini = nlattot
rat_b_end = nlattot
rat_b_ini = n
rat_b_end = n
exit loopi
end if
end do
end do
end do loopi
end if

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

! flag which overlaps will be used
allocate(useovrlp(nmo,nmo,nspin))
useovrlp = .true.
Expand Down Expand Up @@ -2711,20 +2711,20 @@ subroutine calc_di3_wannier(res,nmo,nbnd,nlat,nattr,nspin,atom1,atom2)
end do
write (*,*) "fin"

write (*,*) "xx checking consistency between fa and fa3"
do is = 1, nspin
do iat_a = iat_a_ini, iat_a_end
do iat_b = iat_b_ini, iat_b_end
do rat_b = rat_b_ini, rat_b_end
write (*,'(I3,X,I3,X,I3,X,4(F20.12,X))') iat_a, iat_b, rat_b,&
abs(res%fa(iat_a,iat_b,rat_b,is)-sum(res%fa3(iat_a,iat_b,rat_b,:,:,is))),&
abs(res%fa(iat_a,iat_b,rat_b,is)/sum(res%fa3(iat_a,iat_b,rat_b,:,:,is))),&
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
! write (*,*) "xx checking consistency between fa and fa3"
! do is = 1, nspin
! do iat_a = iat_a_ini, iat_a_end
! do iat_b = iat_b_ini, iat_b_end
! do rat_b = rat_b_ini, rat_b_end
! write (*,'(I3,X,I3,X,I3,X,4(F20.12,X))') iat_a, iat_b, rat_b,&
! abs(res%fa(iat_a,iat_b,rat_b,is)-sum(res%fa3(iat_a,iat_b,rat_b,:,:,is))),&
! abs(res%fa(iat_a,iat_b,rat_b,is)/sum(res%fa3(iat_a,iat_b,rat_b,:,:,is))),&
! 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

Expand Down Expand Up @@ -2812,38 +2812,40 @@ subroutine calc_di3_psink(res,nmo,nbnd,nlat,nattr,nspin,kpt,occ,atom1,atom2)

real*8, allocatable :: lvec(:,:)

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

! set limits for the sums
iat_a_ini = 1
iat_a_end = nattr
if (atom1 >= 0) then
if (atom1 > 0 .and. atom1 <= nattr) then
iat_a_ini = atom1
iat_a_end = atom1
end if
iat_b_ini = 1
iat_b_end = nattr
rat_b_ini = 1
rat_b_end = nlattot
if (atom2(1) >= 0) then
if (atom2(1) > 0 .and. atom2(1) <= nattr) then
iat_b_ini = atom2(1)
iat_b_end = atom2(1)
nlattot = 0
n = 0
loopi: do i = 0, nlat(1)-1
do j = 0, nlat(2)-1
do k = 0, nlat(3)-1
nlattot = nlattot + 1
n = n + 1
if (i == modulo(atom2(2),nlat(1)) .and. j == modulo(atom2(3),nlat(2)) .and.&
k == modulo(atom2(4),nlat(3))) then
rat_b_ini = nlattot
rat_b_end = nlattot
rat_b_ini = n
rat_b_end = n
exit loopi
end if
end do
end do
end do loopi
end if

! number and list of lattice vectors
nlattot = nlat(1)*nlat(2)*nlat(3)
! list of lattice vectors
allocate(lvec(3,nlattot))
nlattot = 0
do i = 0, nlat(1)-1
Expand Down Expand Up @@ -2932,36 +2934,36 @@ subroutine calc_di3_psink(res,nmo,nbnd,nlat,nattr,nspin,kpt,occ,atom1,atom2)
end if
res%fa3 = 0.5d0 * res%fa3 / (fspin*fspin*fspin)

do is = 1, nspin
do iat_a = iat_a_ini, iat_a_end
do iat_b = iat_b_ini, iat_b_end
do rat_b = rat_b_ini, rat_b_end
do iat_c = 1, nattr
do rat_c = 1, nlattot
write (*,'("is=",I1," A:",I2," B:",I2,"/",I2," C:",I2,"/",I2," = ",F15.8)') &
is, iat_a, iat_b, rat_b, iat_c, rat_c, res%fa3(iat_a,iat_b,rat_b,iat_c,rat_c,is)
end do
end do
end do
end do
end do
end do
write (*,*) "fin"

write (*,*) "xx checking consistency between fa and fa3"
do is = 1, nspin
do iat_a = iat_a_ini, iat_a_end
do iat_b = iat_b_ini, iat_b_end
do rat_b = rat_b_ini, rat_b_end
write (*,'(I3,X,I3,X,I3,X,4(F20.12,X))') iat_a, iat_b, rat_b,&
abs(res%fa(iat_a,iat_b,rat_b,is)-sum(res%fa3(iat_a,iat_b,rat_b,:,:,is))),&
abs(res%fa(iat_a,iat_b,rat_b,is)/sum(res%fa3(iat_a,iat_b,rat_b,:,:,is))),&
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
! do is = 1, nspin
! do iat_a = iat_a_ini, iat_a_end
! do iat_b = iat_b_ini, iat_b_end
! do rat_b = rat_b_ini, rat_b_end
! do iat_c = 1, nattr
! do rat_c = 1, nlattot
! write (*,'("is=",I1," A:",I2," B:",I2,"/",I2," C:",I2,"/",I2," = ",F15.8)') &
! is, iat_a, iat_b, rat_b, iat_c, rat_c, res%fa3(iat_a,iat_b,rat_b,iat_c,rat_c,is)
! end do
! end do
! end do
! end do
! end do
! end do
! write (*,*) "fin"

! write (*,*) "xx checking consistency between fa and fa3"
! do is = 1, nspin
! do iat_a = iat_a_ini, iat_a_end
! do iat_b = iat_b_ini, iat_b_end
! do rat_b = rat_b_ini, rat_b_end
! write (*,'(I3,X,I3,X,I3,X,4(F20.12,X))') iat_a, iat_b, rat_b,&
! abs(res%fa(iat_a,iat_b,rat_b,is)-sum(res%fa3(iat_a,iat_b,rat_b,:,:,is))),&
! abs(res%fa(iat_a,iat_b,rat_b,is)/sum(res%fa3(iat_a,iat_b,rat_b,:,:,is))),&
! 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_psink

Expand Down Expand Up @@ -3231,13 +3233,14 @@ subroutine int_output_deloc(bas,res)
real*8 :: fspin, xli, xnn, r1(3), asum, d2, raux, sfac
real*8, allocatable :: dimol(:,:,:,:,:), limol(:), namol(:)
real*8, allocatable :: dist(:), diout(:), xcm(:,:), diout_i(:), diout_o(:)
logical, allocatable :: diout_ioprint1(:), diout_ioprint2(:)
integer, allocatable :: io(:), ilvec(:,:), idat(:), idxmol(:,:)
integer :: ia, ja
integer :: ic, jc, kc, lvec1(3), lvec2(3), lvec3(3)
character(len=:), allocatable :: sncp, scp, sname, sz, smult
type(crystal) :: cr1
type(crystalseed) :: ncseed
logical :: di3
logical :: di3, ok

if (bas%imtype == imtype_isosurface) return

Expand Down Expand Up @@ -3296,13 +3299,40 @@ subroutine int_output_deloc(bas,res)
write (uout,'(" with all atoms in the environment. Last line: sum of LI + 0.5 * DIs,")')
write (uout,'(" equal to the atomic population. Distances are in ",A,".")') iunitname0(iunit)
allocate(dist(natt*nlattot),io(natt*nlattot),diout(natt*nlattot),ilvec(3,natt*nlattot),idat(natt*nlattot))
if (di3) allocate(diout_i(natt*nlattot),diout_o(natt*nlattot))
if (di3) then
allocate(diout_i(natt*nlattot),diout_o(natt*nlattot),diout_ioprint1(natt),diout_ioprint2(natt*nlattot))
k = 0
m = 0
do ic = 0, nlat(1)-1
do jc = 0, nlat(2)-1
do kc = 0, nlat(3)-1
k = k + 1
do j = 1, natt
m = m + 1
ok = sy%propi(l)%di3_atom2(1) < 0 .or. sy%propi(l)%di3_atom2(1) > natt .or.&
(sy%propi(l)%di3_atom2(1) == j) .and.&
ic == modulo(sy%propi(l)%di3_atom2(2),nlat(1)) .and.&
jc == modulo(sy%propi(l)%di3_atom2(3),nlat(2)) .and.&
kc == modulo(sy%propi(l)%di3_atom2(4),nlat(3))
diout_ioprint2(m) = ok
end do
end do
end do
end do
if (sy%propi(l)%di3_atom1 < 0 .or. sy%propi(l)%di3_atom1 > natt) then
diout_ioprint1 = .true.
else
diout_ioprint1 = .false.
diout_ioprint1(sy%propi(l)%di3_atom1) = .true.
end if
end if

do i = 1, natt
if (.not.bas%docelatom(bas%icp(i))) cycle
call assign_strings(i,bas%icp(i),.false.,scp,sncp,sname,smult,sz)
write (uout,'("# Attractor ",A," (cp=",A,", ncp=",A,", name=",A,", Z=",A,") at: ",3(A," "))') &
string(i), trim(scp), trim(sncp), trim(adjustl(sname)), trim(sz), (trim(string(bas%xattr(j,i),'f',12,7)),j=1,3)
if (di3) then
if (di3.and.diout_ioprint1(i)) then
write (uout,'("# Id cp ncp Name Z Latt. vec. &
&---- Cryst. coordinates ---- Distance LI/DI 3c(in) 3c(out)")')
else
Expand All @@ -3329,7 +3359,7 @@ subroutine int_output_deloc(bas,res)
if (dist(m) < 1d-5) sfac = 0.5d0

diout(m) = 2d0 * sum(abs(res(l)%fa(i,j,k,:))) * fspin * sfac
if (di3) then
if (di3.and.diout_ioprint1(i).and.diout_ioprint2(m)) then
diout_i(m) = sum(abs(res(l)%fa3(i,j,k,i,1,:))) + sum(abs(res(l)%fa3(i,j,k,j,k,:)))
diout_o(m) = sum(abs(res(l)%fa3(i,j,k,:,:,:))) - diout_i(m)
diout_i(m) = diout_i(m) * 2d0 * fspin * sfac
Expand All @@ -3349,7 +3379,7 @@ subroutine int_output_deloc(bas,res)
j = io(m)
if (.not.bas%docelatom(bas%icp(idat(j)))) cycle
if (dist(j) < 1d-5) then
if (di3) then
if (di3.and.diout_ioprint1(i).and.diout_ioprint2(j)) then
write (uout,'(" Localization index",71("."),3(A," "))') string(diout(j),'f',12,8,4),&
string(diout_i(j),'f',12,8,4), string(diout_o(j),'f',12,8,4)
else
Expand All @@ -3359,7 +3389,7 @@ subroutine int_output_deloc(bas,res)
else
call assign_strings(j,bas%icp(idat(j)),.false.,scp,sncp,sname,smult,sz)
r1 = bas%xattr(:,idat(j)) + ilvec(:,j)
if (di3) then
if (di3.and.diout_ioprint1(i).and.diout_ioprint2(j)) then
write (uout,'(" ",99(A," "))') string(j,4,ioj_left), scp, sncp, sname, sz,&
(string(ilvec(k,j),3,ioj_right),k=1,3), (string(r1(k),'f',12,7,4),k=1,3),&
string(dist(j),'f',12,7,4), string(diout(j),'f',12,8,4),&
Expand Down Expand Up @@ -3503,6 +3533,8 @@ subroutine int_output_deloc(bas,res)
deallocate(dist,io,diout,ilvec,idat)
if (allocated(diout_i)) deallocate(diout_i)
if (allocated(diout_o)) deallocate(diout_o)
if (allocated(diout_ioprint1)) deallocate(diout_ioprint1)
if (allocated(diout_ioprint2)) deallocate(diout_ioprint2)
end do

end subroutine int_output_deloc
Expand Down

0 comments on commit 445534b

Please sign in to comment.