di3_atom1 and di3_atom2 fields in integrable

src/integration@proc.f90

@@ -674,7 +674,19 @@ module subroutine int_output_header(bas,res,nomol0,usesym0)
           if (res(i)%done) saux = "Lmax = " // string(sy%propi(i)%lmax)
           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) saux = " 3-center indices calculated "
+          if (res(i)%done .and. sy%propi(i)%di3) then
+             if (sy%propi(i)%di3_atom1 < 0) 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) // " "
+             else
+                saux = " 3-center indices calculated for atom 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)) // "," //&
+                   string(sy%propi(i)%di3_atom2(4)) // ") "
+             end if
+          end if
           itaux = string(fid)
        else
           if (res(i)%done) saux = ""
@@ -1598,6 +1610,7 @@ subroutine intgrid_deloc(bas,res)
           end if
        else if ((sy%propi(l)%itype == itype_deloc_wnr .or. sy%propi(l)%itype == itype_deloc_psink).and.&
           sy%propi(l)%fachk .and..not.sy%propi(l)%di3) then
+          ! If DI3 calculation is requested, we need to read/calculate the Sij anyway
           if (read_chk_header(fafname,nbnd,nbndw,nlat,nmo,nlattot,nspin,natt1,isijtype)) then
              fid = sy%propi(l)%fid
 
@@ -1848,6 +1861,8 @@ subroutine intgrid_deloc(bas,res)
 
        ! finished the Sij and Fa successfully
 999    continue
+
+       ! deallocate
        if (allocated(res(l)%sijc)) deallocate(res(l)%sijc)
        if (allocated(res(l)%sij_wnr_imap)) deallocate(res(l)%sij_wnr_imap)
 

src/systemmod@proc.f90

@@ -854,7 +854,7 @@ module subroutine new_integrable_string(s,line,errmsg)
     character(len=:), allocatable, intent(out) :: errmsg
 
     logical :: ok
-    integer :: id, lp, lpold, idum, lp2
+    integer :: id, lp, lpold, idum, idum3(3), lp2
     character(len=:), allocatable :: word, expr, str
     logical :: useexpr, inpsijchk, inpfachk
 
@@ -986,6 +986,20 @@ module subroutine new_integrable_string(s,line,errmsg)
                    end if
                 else if (equal(word,"di3")) then
                    s%propi(s%npropi)%di3 = .true.
+                   s%propi(s%npropi)%di3_atom1 = -1
+                   s%propi(s%npropi)%di3_atom2 = (/-1,0,0,0/)
+                   ok = isinteger(idum,line,lp)
+                   if (ok) then
+                      s%propi(s%npropi)%di3_atom1 = idum
+                      ok = isinteger(idum,line,lp)
+                      if (ok) then
+                         s%propi(s%npropi)%di3_atom2(1) = idum
+                         ok = isinteger(idum3(1),line,lp)
+                         ok = ok .and. isinteger(idum3(2),line,lp)
+                         ok = ok .and. isinteger(idum3(3),line,lp)
+                         if (ok) s%propi(s%npropi)%di3_atom2(2:4) = idum3
+                      end if
+                   end if
                 else
                    lp = lp2
                    exit

src/types.f90

@@ -223,6 +223,8 @@ module types
      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
+     integer :: di3_atom1 = -1 ! first atom for 3-body DI calculation
+     integer :: di3_atom2(4) = (/-1,0,0,0/) ! second atom for 3-body DI calculation
      character(len=mlen) :: sijchkfile = "" ! name of sijchk file
      character(len=mlen) :: fachkfile = ""  ! name of fachk file
   end type integrable