TABLE OF CONTENTS


ABINIT/m_paw_slater [ Modules ]

[ Top ] [ Modules ]

NAME

  m_paw_slater

FUNCTION

  This module defines objects and procedures to evaluate Slater-like integrals
  using real spherical Harmonics.

COPYRIGHT

 Copyright (C) 2008-2018 ABINIT group (MG)
 This file is distributed under the terms of the
 GNU General Public License, see ~abinit/COPYING
 or http://www.gnu.org/copyleft/gpl.txt .

NOTES

  * Routines tagged with "@type_name" are tightly connected to the definition of the data type.
    Tightly connected means that the proper functioning of the implementation relies on the
    assumption that the tagged procedure is consistent with the type declaration.
    Every time a developer changes the structure "type_name" adding new entries, he/she has to make sure
    that all the tightly connected routines are changed accordingly to accommodate the modification of the data type.
    Typical examples of tightly connected routines are creation, destruction or reset methods.

SOURCE

26 #if defined HAVE_CONFIG_H
27 #include "config.h"
28 #endif
29 
30 #include "abi_common.h"
31 
32 MODULE m_paw_slater
33 
34  use defs_basis
35  use m_abicore
36  use m_errors
37  use m_splines
38 
39  use m_fstrings,     only : basename
40  use m_paw_atomorb,  only : atomorb_type, init_atomorb, print_atomorb, destroy_atomorb, get_overlap
41  use m_crystal,      only : crystal_t
42  use m_paw_io,       only : pawio_print_ij
43  use m_pawang,       only : pawang_type, realgaunt
44  use m_pawrad,       only : pawrad_type, pawrad_free, pawrad_isame, &
45 &                           pawrad_deducer0, simp_gen, calc_slatradl
46  use m_pawtab,       only : pawtab_type
47  use m_pawrhoij,     only : pawrhoij_type
48  use m_paw_lmn,      only : make_kln2ln, make_klm2lm, make_indln, klmn2ijlmn
49 
50  implicit none
51 
52  private
53 
54  public :: paw_sigxcore            ! The onsite matrix elements of the Fock operator generated by (closed) core shells.
55  public :: paw_mkdijexc_core       ! Calculate the onsite matrix element of the Fock operator generated by the core.
56  public :: paw_dijhf               ! Compute the onsite D_{ij} strengths of the exchange parth of the self energy.

m_paw_slater/paw_dijhf [ Functions ]

[ Top ] [ m_paw_slater ] [ Functions ]

NAME

  paw_dihf

FUNCTION

  This routine calculates the onsite D_{ij} strengths of the exchange parth of the self energy.

INPUTS

  ndij=Usually ndij=nspden, except for spin-orbit (where ndij=nspinor**2)
  cplex_dij=1 if sigx_dij is real, 2 if they are complex
  lmn2_size_max=Max Number of (klmn) channels over type of atoms.
  my_natom=number of atoms treated by current process
  ntypat=number of atom types
  Pawtab(ntypat)<pawtab_type>=paw tabulated starting data
  Pawrad(ntypat)<pawrad_type>=paw radial mesh and related data
  Pawang<type(pawang_type)>=paw angular mesh and related data
  pawprtvol=Flags governing the verbosity of the output.

OUTPUT

  sigx_dij(cplex_dij*lmn2_size_max,ndij,my_natom)=
    For each atom, the Pseudopotential strengths of the on-site operator Sigma_x

NOTES

  The on-site contribution to the matrix elements of the exchange part of the self-energy is given by:
  <\tpsi_a| [\sum_{ij} |tprj_i\> D_{ij} \<tprj_j|] |\tpsi_b\>.

  When compensation charges are used one obtains:

  D_{ij} = - sum_{kl} \rho_lk [ \Phi_{ikjl} - \Phihat_{ijkl} =

         = - sum_{kl} \rho_lk \sum_{LM} \Gaunt_{ik}^{LM} \Gaunt_{jl}^{LM} [S_{ikjl}^L - tS{ikjl}^L}]

  where S and tS are Slater-like integrals given by

  1)  S_{ijkl}^L = dfrac{4\pi}{2L+1} \iint u_i(1)u_j(1) u_k(2) u_l(2) \dfrac{r_<^L/}{r_>^{L+1}} d1d2.
  1) tS_{ijkl}^L = dfrac{4\pi}{2L+1} \iint [u_i(1)u_j(1)+ tq_{ij}^L g^L(1)]
                                           [u_k(2)u_l(2)+ tq_{kl}^L g^L(2)] \dfrac{r_<^L/}{r_>^{L+1}} d1d2.

  tq_{ij}^L is defined in terms of q_{ij}^L via: q_{ij]^{LM} = tq_{ij}^L \Gaunt_{ij}^{LM}

PARENTS

CHILDREN

SOURCE

1393 subroutine paw_dijhf(ndij,cplex_dij,lmn2_size_max,my_natom,ntypat,Pawtab,Pawrad,Pawang,Pawrhoij,&
1394 &                    sigx_dij,pawprtvol)
1395 
1396 
1397 !This section has been created automatically by the script Abilint (TD).
1398 !Do not modify the following lines by hand.
1399 #undef ABI_FUNC
1400 #define ABI_FUNC 'paw_dijhf'
1401 !End of the abilint section
1402 
1403  implicit none
1404 
1405 !Arguments ------------------------------------
1406 !scalars
1407  integer,intent(in) :: pawprtvol,ndij,cplex_dij,lmn2_size_max,my_natom,ntypat
1408  type(pawang_type),intent(in) :: Pawang
1409 !arrays
1410  real(dp),target,intent(out) :: sigx_dij(cplex_dij*lmn2_size_max,ndij,my_natom) !TODO use ragged arrays pawij?
1411  type(pawtab_type),intent(in) :: Pawtab(ntypat)
1412  type(pawrad_type),intent(in) :: Pawrad(ntypat)
1413  type(pawrhoij_type),intent(in) :: Pawrhoij(my_natom)
1414 
1415 !Local variables ---------------------------------------
1416 !scalars
1417  integer,parameter :: cplex=1    ! FIXME preliminary implementation
1418  integer :: iatom,itypat,lmn_size,lmn2_size,ispden,nspden,ln2_size
1419  integer :: lm2_size !,isppol ln_size,
1420  integer :: irhoij,jrhoij,dplex
1421  integer :: rho_lmn !,rho_klm,rho_kln,rho_lmin,rho_lmax,rho_iln,rho_jln
1422  integer :: klmn
1423  integer :: i_lmn,j_lmn,k_lmn,l_lmn
1424  integer :: which_intg,l_max,opt_l
1425  real(dp) :: slt_ikjl,slt_iljk
1426  !character(len=500) :: msg
1427 !arrays
1428  integer :: opt_l_index(0,0),pack2ij(0)
1429  real(dp) :: ro(cplex)
1430  real(dp), ABI_CONTIGUOUS pointer :: sigx_atm(:,:)
1431  type(slatrad_t),allocatable :: Slatrad4(:)
1432 
1433 ! *************************************************************************
1434 
1435  DBG_ENTER("COLL")
1436 
1437  ABI_CHECK(ndij/=4,"ndij=4 not coded")
1438  ABI_CHECK(cplex_dij==1,"cplex_dij/=2 not coded")
1439  ABI_CHECK(lmn2_size_max==MAXVAL(Pawtab(:)%lmn2_size),"Wrong lmn2_size_max")
1440 
1441  if (my_natom>0) then
1442    if (pawrhoij(1)%cplex<cplex) then
1443      MSG_BUG('Must have pawrhoij()%cplex >= cplex !')
1444    end if
1445  end if
1446 
1447  dplex=cplex-1 ! used to select the elements of rho_ij.
1448  sigx_dij=zero
1449 
1450  do iatom=1,my_natom
1451    itypat   =Pawrhoij(iatom)%itypat
1452    lmn_size =Pawtab(itypat)%lmn_size
1453    lmn2_size=Pawtab(itypat)%lmn2_size
1454    l_max    =(Pawtab(itypat)%l_size+1)/2
1455    lm2_size = (l_max**2)*(l_max**2+1)/2
1456    !write(std_out,*)"in atom ",iatom,"lm2_size=",lm2_size
1457 
1458    ! Calculate Slater integral for this atom type.
1459    ! TODO obviously these tables should be stored in Pawtab!
1460    ln2_size = Pawtab(itypat)%ij_size
1461    ABI_DATATYPE_ALLOCATE(Slatrad4,(ln2_size*(ln2_size+1)/2))
1462    which_intg=3
1463    call slatrad_init(Slatrad4,which_intg,ln2_size,Pawrad(itypat),Pawtab(itypat))
1464 
1465    sigx_atm => sigx_dij(:,:,iatom)
1466    !
1467    ! * Loop over spin components.
1468    nspden=ndij
1469    do ispden=1,ndij
1470      !
1471      ! ============================================================
1472      ! ==== Summing over the non-zero lk channels of \rho_{lk} ====
1473      ! ============================================================
1474      jrhoij=1
1475      do irhoij=1,pawrhoij(iatom)%nrhoijsel
1476        rho_lmn=pawrhoij(iatom)%rhoijselect(irhoij)
1477 
1478        ! check wheter rho_lmin is consistent with the Indexing used in slatrad
1479        !rho_klm =pawtab(itypat)%indklmn(1,rho_lmn)
1480        !rho_kln =pawtab(itypat)%indklmn(2,rho_lmn)
1481        !rho_lmin=pawtab(itypat)%indklmn(3,rho_lmn)
1482        !rho_lmax=pawtab(itypat)%indklmn(4,rho_lmn)
1483 
1484        ! Retrieve rhoij for this ispden.
1485        if (nspden/=2) then
1486          ro(1:cplex)=pawrhoij(iatom)%rhoijp(jrhoij:jrhoij+dplex,ispden)
1487        else
1488          MSG_ERROR("Recheck this part")
1489          if (ispden==1) then
1490            ro(1:cplex)=pawrhoij(iatom)%rhoijp(jrhoij:jrhoij+dplex,1) + pawrhoij(iatom)%rhoijp(jrhoij:jrhoij+dplex,2)
1491          else if (ispden==2) then
1492            ro(1:cplex)=pawrhoij(iatom)%rhoijp(jrhoij:jrhoij+dplex,1)
1493          end if
1494        end if
1495        !
1496        ! Avoid double-counting the diagonal of rho.
1497        ro(1:cplex)=ro(1:cplex)*pawtab(itypat)%dltij(rho_lmn)*half
1498 
1499        call klmn2ijlmn(rho_lmn,lmn_size,k_lmn,l_lmn)
1500 
1501        ! Loop over the upper triangle of the D_{ij) matrix and accumulate:
1502        ! sum_\lk rho_\kl [ \Phi_{ikjl} + \Phi_{iljk} - \Phihat_{ikjl} - \Phihat_{iljk} ]
1503        do klmn=1,lmn2_size
1504          ! Calculate the indeces in the Slatrad4 structure.
1505          call klmn2ijlmn(klmn,lmn_size,i_lmn,j_lmn)
1506 
1507          ! My formula
1508          slt_ikjl = slat_intg(Slatrad4,Pawtab(itypat),Pawang,i_lmn,k_lmn,j_lmn,l_lmn)
1509          slt_iljk = slat_intg(Slatrad4,Pawtab(itypat),Pawang,i_lmn,l_lmn,j_lmn,k_lmn)
1510 
1511          !slt_ikjl = slat_intg(Slatrad4,Pawtab(itypat),Pawang,i_lmn,k_lmn,l_lmn,j_lmn)
1512          !slt_iljk = slat_intg(Slatrad4,Pawtab(itypat),Pawang,i_lmn,l_lmn,k_lmn,j_lmn)
1513          !slt_iljk = slt_ikjl
1514 
1515          sigx_atm(klmn,ispden) = sigx_atm(klmn,ispden) + ro(1) * (slt_ikjl + slt_iljk)
1516        end do ! klmn
1517 
1518        jrhoij=jrhoij+pawrhoij(iatom)%cplex
1519      end do ! irhoij
1520    end do ! ispden
1521 
1522    if (ABS(pawprtvol)>=1) then ! * Print values
1523      call wrtout(std_out,"   ************** Dij Fock ************ ",'COLL')
1524      opt_l=-1
1525      call pawio_print_ij(std_out,sigx_atm(:,1),lmn2_size,cplex_dij,lmn_size,opt_l,opt_l_index,0,pawprtvol,pack2ij,-one,1)
1526    end if
1527 
1528    call slatrad_free(Slatrad4)
1529    ABI_DATATYPE_DEALLOCATE(Slatrad4)
1530  end do ! iatom
1531 
1532  ! Factor half cancels in the derivation wrt rho_ij.
1533  sigx_dij = - sigx_dij
1534 
1535  DBG_EXIT("COLL")
1536 
1537 end subroutine paw_dijhf

m_paw_slater/paw_mkdijexc_core [ Functions ]

[ Top ] [ m_paw_slater ] [ Functions ]

NAME

  paw_mkdijexc_core

FUNCTION

  Driver routine to calculate the onsite matrix element of the Fock operator between two
  all-electron partial waves.

INPUTS

  ndij=Usually ndij=nspden, except for spin-orbit (where ndij=nspinor**2)
  cplex_dij=1 if dijexc_core is real, 2 if they are complex
  lmn2_size_max=Max Number of (klmn) channels over type of atoms.
  Cryst<crystal_t>=Structure describing the crystal structure and its symmmetries.
  Pawtab(ntypat)<pawtab_type>=paw tabulated starting data
  Pawrad(ntypat)<pawrad_type>=paw radial mesh and related data
  pawprtvol=Flags governing the verbosity of the output.
  filpsp(ntypat)=names of the files containing the all-electron core WF

OUTPUT

  dijexc_core(cplex_dij*lmn2_size_max,ndij,ntypat)= On-site matrix elements $ \<\phi_i|Sigma_x^\core|\phi_j\>
    for each type of atom.

PARENTS

      sigma

CHILDREN

      klmn2ijlmn,pawio_print_ij,slatrad_free,slatrad_init,wrtout

SOURCE

902 subroutine paw_mkdijexc_core(ndij,cplex_dij,lmn2_size_max,Cryst,Pawtab,Pawrad,dijexc_core,pawprtvol,filpsp)
903 
904 
905 !This section has been created automatically by the script Abilint (TD).
906 !Do not modify the following lines by hand.
907 #undef ABI_FUNC
908 #define ABI_FUNC 'paw_mkdijexc_core'
909 !End of the abilint section
910 
911  implicit none
912 
913 !Arguments ------------------------------------
914 !scalars
915  integer,intent(in) :: pawprtvol,ndij,cplex_dij,lmn2_size_max
916  type(crystal_t),intent(in) :: Cryst
917 !arrays
918  real(dp),intent(out) :: dijexc_core(cplex_dij*lmn2_size_max,ndij,Cryst%ntypat) !TODO use ragged arrays pawij?
919  character(len=fnlen) :: filpsp(Cryst%ntypat)
920  type(Pawtab_type),target,intent(in) :: Pawtab(Cryst%ntypat)
921  type(Pawrad_type),intent(in) :: Pawrad(Cryst%ntypat)
922 
923 !Local variables ---------------------------------------
924 !scalars
925  integer :: itypat,ic,ierr,lmn_size,lmn2_size,ln_size,isppol
926  real(dp) :: rcut
927  character(len=500) :: header,msg
928  character(len=fnlen) :: fcore,string
929 !arrays
930  integer,allocatable :: phi_indln(:,:)
931  real(dp),ABI_CONTIGUOUS pointer :: phi(:,:)
932  real(dp),allocatable :: overlap(:,:)
933  type(atomorb_type),allocatable :: Atm(:)
934  type(Pawrad_type),allocatable :: Radatm(:)
935 
936 ! *************************************************************************
937 
938  ABI_DATATYPE_ALLOCATE(Atm,(Cryst%ntypat))
939  ABI_DATATYPE_ALLOCATE(Radatm,(Cryst%ntypat))
940 
941  ABI_CHECK(ndij==1     ,"spinor+HF not available")
942  ABI_CHECK(cplex_dij==1,"spinor+HF not available")
943  ABI_CHECK(lmn2_size_max==MAXVAL(Pawtab(:)%lmn2_size),"Wrong lmn2_size_max")
944 
945  !allocate(dijexc_core(cplex_dij*lmn2_size_max,ndij,Cryst%ntypat)) !TODO use ragged arrays pawij?
946  dijexc_core=zero
947 
948  do itypat=1,Cryst%ntypat
949 
950    ! Read core orbitals for this atom type.
951    string = filpsp(itypat)
952    fcore = "CORE_"//TRIM(basename(string))
953    ic = INDEX (TRIM(string), "/" , back=.TRUE.) ! if string is a path, prepend path to fcore.
954    if (ic>0 .and. ic<LEN_TRIM(string)) fcore = filpsp(itypat)(1:ic)//TRIM(fcore)
955 
956    rcut=Pawtab(itypat)%rpaw
957    call init_atomorb(Atm(itypat),Radatm(itypat),rcut,fcore,pawprtvol,ierr)
958 
959    if (ierr/=0) then
960      msg = " Error reading core orbitals from file: "//TRIM(fcore)
961      MSG_ERROR(msg)
962    end if
963    write(header,'(a,i4,a)')" === Atom type = ",itypat," === "
964    call print_atomorb(Atm(itypat),header,unit=std_out,prtvol=pawprtvol)
965    !
966    ! * Calculate $ \<\phi_i|Sigma_x^\core|\phi_j\> $ for this atom type.
967    lmn_size  = Pawtab(itypat)%lmn_size
968    lmn2_size = Pawtab(itypat)%lmn2_size
969 
970    call paw_sigxcore(cplex_dij,lmn2_size,ndij,&
971 &    Pawrad(itypat),Pawtab(itypat),Atm(itypat),Radatm(itypat),dijexc_core(1:lmn2_size,:,itypat))
972 
973    ln_size =  Pawtab(itypat)%basis_size
974    phi     => Pawtab(itypat)%phi
975 
976    ABI_MALLOC(phi_indln,(2,ln_size))
977    call make_indln(lmn_size,ln_size,Pawtab(itypat)%indlmn(:,:),phi_indln)
978 
979    ABI_MALLOC(overlap,(Atm(itypat)%ln_size,ln_size))
980    isppol=1 ! hardcoded
981    call get_overlap(Atm(itypat),Radatm(itypat),Pawrad(itypat),isppol,ln_size,phi,phi_indln,overlap)
982 
983    ABI_FREE(phi_indln)
984    ABI_FREE(overlap)
985  end do ! ntypat
986 
987  ! Free memory
988  call pawrad_free(Radatm)
989  do itypat=1,Cryst%ntypat
990    call destroy_atomorb(Atm(itypat))
991  end do
992 
993  ABI_DATATYPE_DEALLOCATE(Atm)
994  ABI_DATATYPE_DEALLOCATE(Radatm)
995 
996 end subroutine paw_mkdijexc_core

m_paw_slater/paw_sigxcore [ Functions ]

[ Top ] [ m_paw_slater ] [ Functions ]

NAME

  paw_sigxcore

FUNCTION

  Calculate the integrals:
  \dfrac{4\pi}{2L+1} \int \phi_{\ni\li}(1) orb_{nl}(1) \dfrac{r_<^L}{r_>^{L+1}} orb_{nl}(2) \phi_{\nj\lj} d1d2!!
  for given (in,il) and (jn,jl) as a function of (nc,lc) and L = |il-jl|, |il-jl|+2, ..., |il+il|

INPUTS

  cplex_dij=1 if dijexc_core is real, 2 if they are complex
  lmn2_size=Number of (klmn) channels
  ndij=Usually ndij=nspden, except for spin-orbit (where ndij=nspinor**2)
  Pawtab<pawtab_type>=paw tabulated starting data
  Atm<atomorb_type>=Structure containing core orbitals
  Atmrad<pawrad_type>=The radial mesh for core orbitals

OUTPUT

  dijexc_core(cplex_dij*lmn2_size,ndij)

PARENTS

      m_paw_slater

CHILDREN

      klmn2ijlmn,pawio_print_ij,slatrad_free,slatrad_init,wrtout

SOURCE

751 subroutine paw_sigxcore(cplex_dij,lmn2_size,ndij,Pawrad,Pawtab,Atm,Atmrad,dijexc_core)
752 
753 
754 !This section has been created automatically by the script Abilint (TD).
755 !Do not modify the following lines by hand.
756 #undef ABI_FUNC
757 #define ABI_FUNC 'paw_sigxcore'
758 !End of the abilint section
759 
760  implicit none
761 
762 !Arguments ------------------------------------
763 !scalars
764  integer,intent(in) :: lmn2_size,cplex_dij,ndij
765 !arrays
766  real(dp),intent(out) :: dijexc_core(cplex_dij*lmn2_size,ndij)
767  type(atomorb_type),intent(in) :: Atm
768  type(pawrad_type),intent(in) :: Atmrad
769  type(pawrad_type),intent(in) :: Pawrad
770  type(pawtab_type),target,intent(in) :: Pawtab
771 
772 !Local variables ---------------------------------------
773 !scalars
774  integer :: ilnc,ilc,lnc_size,l_max
775  integer :: lm2_size,ln_size,ln2_size,lmn_size
776  integer :: lm_size,klmn,kln,klm
777  integer :: lc_max,ilsl,isgg,israd,opt_l,pawprtvol
778  real(dp) :: tmp,sgg,intgrl
779 !character(len=500) :: msg
780 !arrays
781  integer :: opt_l_index(0,0),pack2ij(0)
782  integer,allocatable :: kln2ln(:,:),klm2lm(:,:)
783  integer, pointer :: indklmn(:,:),indlmn(:,:)
784  type(slatang_cshell_t),allocatable :: Slatang3l(:)
785  type(slatrad_cshell_t),allocatable :: Slatrad3l(:)
786 
787 ! *************************************************************************
788 
789  ! * Consistency check
790  ABI_CHECK(cplex_dij==1,"cplex_dij must be 1")
791 
792  ABI_CHECK(ndij==1,"ndij must be 1")
793 
794  ABI_CHECK(lmn2_size==Pawtab%lmn2_size,"Wrong lmn2_size")
795 
796  lmn_size  = Pawtab%lmn_size
797  ln_size   = Pawtab%basis_size
798  ln2_size  = Pawtab%ij_size
799  l_max     = (Pawtab%l_size-1)/2 +1
800  lm_size   = l_max**2
801  lm2_size  = lm_size*(lm_size+1)/2
802 
803  indlmn  => Pawtab%indlmn(1:6,1:lmn_size)
804  indklmn => Pawtab%indklmn(1:8,1:lmn2_size)
805 
806  ! * Setup of useful tables.
807  ABI_MALLOC(kln2ln,(6,ln2_size))
808  call make_kln2ln(lmn_size,lmn2_size,ln2_size,indlmn,indklmn,kln2ln)
809 
810  ABI_MALLOC(klm2lm,(6,lm2_size))
811  call make_klm2lm(lmn_size,lmn2_size,lm2_size,indlmn,indklmn,klm2lm)
812 
813  ! * Integrate angular part.
814  lnc_size  = Atm%ln_size
815  lc_max    = Atm%l_max
816 
817  ABI_DATATYPE_ALLOCATE(Slatang3l,(lm2_size))
818  call slatang_cshell_init(Slatang3l,l_max,lm2_size,lc_max,klm2lm)
819 
820  ABI_FREE(klm2lm)
821 
822  ! * Integrate radial part.
823  ABI_DATATYPE_ALLOCATE(Slatrad3l,(ln2_size))
824 
825  call slatrad_cshell_init(Slatrad3l,ln2_size,Pawrad,Pawtab,Atm,Atmrad)
826 
827  ! === Calculate matrix elements of Sigma_x^core ===
828  ! * $<\phi_i|\Sigma_x^\core|\phi_j>$
829  dijexc_core = zero
830  do klmn=1,lmn2_size
831    klm = Pawtab%indklmn(1,klmn)
832    kln = Pawtab%indklmn(2,klmn)
833    !
834    ! * Summing over (lc,nc) and lslat
835    tmp = zero
836    if (Slatang3l(klm)%nsggsel >0) then
837      do ilnc=1,Atm%ln_size
838        ilc = 1+Atm%indln(1,ilnc)
839        do ilsl=1,Slatang3l(klm)%lslat_max !FIXME check this
840          !do ilsl=Slatang3l(klm)%lslat_min,Slatang3l(klm)%lslat_max
841          isgg  = Slatang3l(klm)%sggselect(ilsl,ilc)
842          israd = Slatrad3l(kln)%rlphic_select(ilsl,ilnc)
843          if (isgg>0 .and. israd>0) then
844            sgg    = Slatang3l(klm)%sgg(isgg)
845            intgrl = Slatrad3l(kln)%rlphic_int(israd)
846            tmp = tmp + intgrl*sgg
847          end if
848        end do
849      end do
850    end if
851 
852    dijexc_core(klmn,1) = -tmp  ! Store results.
853  end do
854 
855  ! * Print values
856  call wrtout(std_out,"   ************** Dij Fock_core ************ ",'COLL')
857  opt_l=-1; pawprtvol=-1
858  call pawio_print_ij(std_out,dijexc_core(:,1),lmn2_size,cplex_dij,lmn_size,opt_l,opt_l_index,0,pawprtvol,pack2ij,-one,1)
859 
860  ! * Free memory.
861  ABI_FREE(kln2ln)
862  call slatang_cshell_free(Slatang3l)
863  ABI_DT_FREE(Slatang3l)
864  call slatrad_cshell_free(Slatrad3l)
865  ABI_DT_FREE(Slatrad3l)
866 
867 end subroutine paw_sigxcore

m_paw_slater/slat_intg [ Functions ]

[ Top ] [ m_paw_slater ] [ Functions ]

NAME

   slat_intg

FUNCTION

  Helper function returning the slater integral
    \int_\Omega \phi_i(1)\phi_j(1) \dfrac{1}{|1-2|} \phi_k(2)\phi_l(2) d1d2

INPUTS

PARENTS

CHILDREN

SOURCE

1630 function slat_intg(Slatrad4,Pawtab,Pawang,i_lmn,j_lmn,k_lmn,l_lmn)
1631 
1632 
1633 !This section has been created automatically by the script Abilint (TD).
1634 !Do not modify the following lines by hand.
1635 #undef ABI_FUNC
1636 #define ABI_FUNC 'slat_intg'
1637 !End of the abilint section
1638 
1639  implicit none
1640 
1641 !Arguments ------------------------------------
1642 !scalars
1643  integer,intent(in) :: i_lmn,j_lmn,k_lmn,l_lmn
1644  real(dp) :: slat_intg
1645  type(pawtab_type),intent(in) :: Pawtab
1646  type(pawang_type),intent(in) :: Pawang
1647 !arrays
1648  type(slatrad_t),intent(in) :: Slatrad4(:)
1649 
1650 !Local variables-------------------------------
1651 !scalars
1652  integer :: ij_lmn,kl_lmn,kl_ln,ij_lm,kl_lm,ilsum,ij_ln
1653  integer :: isel,slt_idx
1654  integer :: iln,jln,kln,lln,ii
1655  real(dp) :: sltL_ijkl,angintL_ijkl
1656  !character(len=500) :: msg
1657 
1658 !************************************************************************
1659 
1660  ! The lmn packed indeces for (ij) and (kl).
1661  if (j_lmn>=i_lmn) then
1662    ij_lmn = i_lmn + j_lmn*(j_lmn-1)/2
1663  else
1664    ij_lmn = j_lmn + i_lmn*(i_lmn-1)/2
1665  end if
1666 
1667  if (l_lmn>=k_lmn) then
1668    kl_lmn = k_lmn + l_lmn*(l_lmn-1)/2
1669  else
1670    kl_lmn = l_lmn + k_lmn*(k_lmn-1)/2
1671  end if
1672  !
1673  ! The lm indeces for (ij) and (kl) in packed storage.
1674  ij_lm = pawtab%indklmn(1,ij_lmn)
1675  ij_ln = pawtab%indklmn(2,ij_lmn)
1676 
1677  kl_lm = pawtab%indklmn(1,kl_lmn)
1678  kl_ln = pawtab%indklmn(2,kl_lmn)
1679  !
1680  ! The index of (ijkl) in the Slatrad4 database.
1681  if (kl_ln>=ij_ln) then
1682    slt_idx = ij_ln +kl_ln*(kl_ln-1)/2
1683  else
1684    slt_idx = kl_ln +ij_ln*(ij_ln-1)/2
1685  end if
1686 
1687 !BEGIN DEBUG
1688  iln = Slatrad4(slt_idx)%iln
1689  jln = Slatrad4(slt_idx)%jln
1690  kln = Slatrad4(slt_idx)%kln
1691  lln = Slatrad4(slt_idx)%lln
1692 
1693  ii = kln + lln*(lln-1)/2
1694  if (slt_idx /=  (iln + jln*(jln-1)/2 + ii*(ii-1)/2 )) then
1695    write(std_out,*)"slt_idx, iln, jln, kln, lln",slt_idx, iln, jln, kln, lln
1696    MSG_BUG("Check indeces")
1697  end if
1698 !END DEBUG
1699  !
1700  ! Calculate the integral by summing over ilsum.
1701  slat_intg=zero
1702  if (Slatrad4(slt_idx)%nintgl>0) then
1703    do ilsum=Slatrad4(slt_idx)%lslat_min,Slatrad4(slt_idx)%lslat_max
1704    !% do ilsum=Slatrad4(slt_idx)%lslat_min,Slatrad4(slt_idx)%lslat_max,2
1705      isel = Slatrad4(slt_idx)%intgl_select(ilsum)
1706      if (isel/=0) then
1707        sltL_ijkl = Slatrad4(slt_idx)%intgl(isel)
1708        angintL_ijkl = summ_2gaunt(Pawang,ij_lm,kl_lm,ilsum)
1709        slat_intg = slat_intg + sltL_ijkl * angintL_ijkl
1710      end if
1711    end do
1712  end if
1713 
1714 end function slat_intg

m_paw_slater/slatang_cshell_free [ Functions ]

[ Top ] [ m_paw_slater ] [ Functions ]

NAME

  slatang_cshell_free

FUNCTION

  Free the dynamic memory allocated in a structure of type slatang_cshell_t

SIDE EFFECTS

  Slatang3l(lm2_size) <type(slatang_cshell_t)> = Object containing radial integrals

PARENTS

      m_paw_slater

CHILDREN

      klmn2ijlmn,pawio_print_ij,slatrad_free,slatrad_init,wrtout

SOURCE

425 subroutine slatang_cshell_free(Slatang3l)
426 
427 
428 !This section has been created automatically by the script Abilint (TD).
429 !Do not modify the following lines by hand.
430 #undef ABI_FUNC
431 #define ABI_FUNC 'slatang_cshell_free'
432 !End of the abilint section
433 
434  implicit none
435 
436 !Arguments ------------------------------------
437 !scalars
438  type(slatang_cshell_t),intent(inout) :: Slatang3l(:)
439 
440 !Local variables-------------------------------
441  integer :: ii
442 ! *********************************************************************
443 
444  !@slatang_cshell_t
445  do ii=1,SIZE(Slatang3l)
446    if (allocated(Slatang3l(ii)%sggselect)) then
447      ABI_FREE(Slatang3l(ii)%sggselect)
448    end if
449    if (allocated(Slatang3l(ii)%sgg)) then
450      ABI_FREE(Slatang3l(ii)%sgg)
451    end if
452  end do
453 
454 end subroutine slatang_cshell_free

m_paw_slater/slatang_cshell_init [ Functions ]

[ Top ] [ m_paw_slater ] [ Functions ]

NAME

  slatang_cshell_init

FUNCTION

  Initialize the structure slatrad_cshell_t containing radial integrals, see below.

INPUTS

  l_max= max. value of ang. momentum l+1;
   Gaunt coeffs up to [(2*l_max-1,m),(l_max,m),(l_max,m)] are computed
  lc_max=Max Lc+1 for core states used to contrusct \Sigma_x^\core.
  lm2_size=Number of symmetrix elements in the (l,m) basis set.
   NB: lm2_size = (l_max**2)*(l_max**2+1)/2.
  klm2lm(6,lm2_size)=Table giving il, jl ,im, jm, ilm and jlm for each klm=(ilm,jlm)
  where ilm=(il,im) and ilm<=jlm. NB: klm2lm is an application and not a bijection.

OUTPUT

  Slatang3l(lm2_size) <type(slatang_cshell_t)> = Object storing :

   $ F^{ilsl,ilc}_{klm} = sum_{msl,mc}  <li mi|lsl msl;lc mc> <lsl msl;lc mc|lj mj> $

  where klm = runs over the upper triangle of the ((il,im),(jl,jm)) matrix.
  ilc runs from 1 up to lc_max and |li-lc| <= lsl <= |li+lc|

NOTES

  Selection rules for F
   1) mi = mj
   2) In the case of closed shells, one sums for all possible mc"s values from -lc up to +lc.
      In this particular case, one can use the symmetry properties of Clebsch-Gordan
      coefficients to show that F is non null only if li==lj. In particular, F can be rewritten as:

      $ F^{ilsl,ilc}_{klm} = \delta{li,lj}\delta{mi,mj} \times
         \Gaunt^{lsl,0}_{lc,0;li,0} \sqrt{ \dfrac{(2*lc+1) (2*lsl+1)}{4\pi*(2*li+1)} } $

PARENTS

      m_paw_slater

CHILDREN

      klmn2ijlmn,pawio_print_ij,slatrad_free,slatrad_init,wrtout

SOURCE

259 subroutine slatang_cshell_init(Slatang3l,l_max,lm2_size,lc_max,klm2lm)
260 
261 
262 !This section has been created automatically by the script Abilint (TD).
263 !Do not modify the following lines by hand.
264 #undef ABI_FUNC
265 #define ABI_FUNC 'slatang_cshell_init'
266 !End of the abilint section
267 
268  implicit none
269 
270 !Arguments ------------------------------------
271 !scalars
272  integer,intent(in) :: l_max,lc_max,lm2_size
273 !arrays
274  integer,intent(in) :: klm2lm(6,lm2_size)
275  type(slatang_cshell_t),intent(out) :: Slatang3l(lm2_size)
276 
277 !Local variables-------------------------------
278 !scalars
279  integer :: ilm,ilm0,jlm,lgnt_max,ngnt,ilsl,ilc,lc,ilm0c,lsl
280  integer :: klm_ci,ilm0sl,li,il,jl,im,jm,ig000
281  integer :: klm,k0lm_i,k0lm_j,k0lm_c,nsggsel,lslat_max,lslat_min
282  real(dp) :: dum
283 !arrays
284  integer,allocatable :: gntselect(:,:)
285  real(dp),allocatable :: realgnt(:),tmp_sgg(:)
286  real(dp),allocatable :: g000(:,:,:)
287 
288 ! *************************************************************************
289 
290  !@slatang_cshell_t
291 
292  ! * Calculate $\Gaunt^{lsl,msl}_{lc,mc;li,mi}$
293  lgnt_max = MAX(l_max,lc_max)
294  ABI_MALLOC(  realgnt,((2*lgnt_max-1)**2*(lgnt_max)**4))
295  ABI_MALLOC(gntselect,((2*lgnt_max-1)**2, lgnt_max**2*(lgnt_max**2+1)/2))
296 
297  call realgaunt(lgnt_max,ngnt,gntselect,realgnt)
298 
299  ! Below we need $\Gaunt_{lsl,0}_{lc,0;li,0}$
300  ABI_MALLOC(g000,(2*lgnt_max-1,lc_max,l_max))
301  g000 = zero
302 
303  do ilsl=1,2*lgnt_max-1
304    lsl    = ilsl-1
305    ilm0sl = 1+lsl**2+lsl
306    do il=1,l_max
307      li    = il-1
308      ilm0  = 1+li**2+li
309      k0lm_i = ilm0 *(ilm0-1)/2
310      do ilc=1,lc_max
311        lc    = ilc-1
312        ilm0c = 1+lc**2+lc
313        k0lm_c= ilm0c * (ilm0c-1)/2
314        if (ilm0c > ilm0) then
315          klm_ci = k0lm_c + ilm0
316        else
317          klm_ci = k0lm_i + ilm0c
318        end if
319        ig000 = gntselect(ilm0sl,klm_ci) ! Index of $\Gaunt_{lsl,0}_{lc,0;li,0}$
320        if (ig000 > 0) g000(ilsl,ilc,il)=realgnt(ig000)
321      end do
322    end do
323  end do
324 
325  ABI_FREE(realgnt)
326  ABI_FREE(gntselect)
327 
328  ! === Loop over klm channels in packed form ===
329  do klm=1,lm2_size
330    il = klm2lm(1,klm); im = klm2lm(3,klm)
331    jl = klm2lm(2,klm); jm = klm2lm(4,klm)
332 
333    nsggsel=0
334    lslat_min = 1 !FIXME find better way
335    lslat_max = il+lc_max-1
336 
337    Slatang3l(klm)%lslat_min = lslat_min
338    Slatang3l(klm)%lslat_max = lslat_max
339    Slatang3l(klm)%lc_max    = lc_max
340 
341    ABI_MALLOC(tmp_sgg,(lslat_max*lc_max))
342    tmp_sgg = zero
343 
344    ! === Calculate F^{lsl,lc}_{li,mi;lj,mj} ===
345    ! * Selection rule: mi = mj and li==lj
346    if (im == jm .and. il==jl) then
347      li  = il-1
348      ilm = klm2lm(5,klm); k0lm_i = ilm *(ilm -1)/2
349      jlm = klm2lm(6,klm); k0lm_j = jlm *(jlm -1)/2
350 
351      ABI_MALLOC(Slatang3l(klm)%sggselect,(lslat_max,lc_max))
352      Slatang3l(klm)%sggselect = 0
353 
354      do ilsl=lslat_min,lslat_max
355      !% do ilsl=lslat_min,lslat_max,2
356        lsl = ilsl-1
357        do ilc=1,lc_max
358          lc = ilc-1
359          dum = SQRT( (two*lc+1)*(two*lsl+1) / (four_pi*(two*li+1)) ) * g000(ilsl,ilc,il)
360          if (ABS(dum)>=tol12) then ! * Store results and progressive index if non null.
361            nsggsel = nsggsel + 1
362            tmp_sgg(nsggsel) = dum
363            Slatang3l(klm)%sggselect(ilsl,ilc) = nsggsel
364          end if
365        end do !ilc
366      end do !ilsl
367    end if ! Selection rule li=lj and mi == mj
368    !
369    ! * Finalize the object.
370    Slatang3l(klm)%nsggsel = nsggsel
371    if (nsggsel > 0) then
372      ABI_MALLOC(Slatang3l(klm)%sgg,(nsggsel))
373      Slatang3l(klm)%sgg = tmp_sgg(1:nsggsel)
374    end if
375    ABI_FREE(tmp_sgg)
376  end do !klm
377 
378  ABI_FREE(g000)
379 
380 #if 0
381 ! Debugging code
382  do klm=1,lm2_size
383    if (Slatang3l(klm)%nsggsel>0) then
384      il  = klm2lm(1,klm)
385      jl  = klm2lm(2,klm)
386      im  = klm2lm(3,klm)
387      jm  = klm2lm(4,klm)
388      write(std_out,*)"--for li, mi",il-1,im-il
389      lslat_min = Slatang3l(klm)%lslat_min
390      lslat_max = Slatang3l(klm)%lslat_max
391 
392      do ilc=1,lc_max
393        do ilsl=lslat_min,lslat_max
394          ii = Slatang3l(klm)%sggselect(ilsl,ilc)
395          if (ii>0) write(std_out,*)"   lc, lslat, sgg",ilc-1,ilsl-1,Slatang3l(klm)%sgg(ii)
396        end do
397      end do
398    end if
399  end do
400 #endif
401 
402 end subroutine slatang_cshell_init

m_paw_slater/slatang_cshell_t [ Types ]

[ Top ] [ m_paw_slater ] [ Types ]

NAME

  slatang_cshell_t

FUNCTION

  Object used to store:
   $ F^{lsl,lc}_{li,lj,mi,mj} = sum_{msl mc} \<li mi|lsl msl;lc mc\> \<lsl msl;lc mc| lj mj\> $
  This (less general) type of radial integral is needed to evaluate the Exchange term generated
  by a closed-shell atom. In the equation, (lc,mc) are the set of angular quantum number associated
  to (closed) core electrons while  (lsl,msl) comes from the expansion of 1/|r1-r2|.
  Since the F is invariant under exchange of i and j we use an array of structures indexed
  by kln = (iln,jln) in packed form.

SOURCE

137  type, public :: slatang_cshell_t
138 
139   integer :: nsggsel
140   ! Number of non null matrix elements
141 
142   integer :: lslat_max
143   ! Max l+1 in the expansion of the Coulomb potential
144 
145   integer :: lslat_min
146   ! Min l+1 in the expansion of the Coulomb potential
147 
148   integer :: lc_max
149   ! Max l+1 for orbitals summed over (usually core orbitals)
150 
151   integer,allocatable :: sggselect(:,:)
152   ! sggselect(lslat_max,lc_max)
153   ! Index of non null sgg, 0 if sgg is zero by symmetry.
154 
155   real(dp),allocatable :: sgg(:)
156   ! sgg(nsggsel)
157   ! Non null matrix elements in packed form. The index is given by sggselect.
158 
159  end type slatang_cshell_t
160 
161  public :: slatang_cshell_init  ! Creation method for slatang_cshell_t.
162  public :: slatang_cshell_free  ! Destruction method for the slatang_cshell_t.

m_paw_slater/slatrad_cshell_free [ Functions ]

[ Top ] [ m_paw_slater ] [ Functions ]

NAME

  slatrad_cshell_free

FUNCTION

  Free the dynamic memory allocated in a structure of type slatrad_cshell_t

SIDE EFFECTS

  Slatrad3l(ln2_size) <type(slarad3l_type)> = Object containing radial integrals

PARENTS

      m_paw_slater

CHILDREN

      klmn2ijlmn,pawio_print_ij,slatrad_free,slatrad_init,wrtout

SOURCE

477 subroutine slatrad_cshell_free(Slatrad3l)
478 
479 
480 !This section has been created automatically by the script Abilint (TD).
481 !Do not modify the following lines by hand.
482 #undef ABI_FUNC
483 #define ABI_FUNC 'slatrad_cshell_free'
484 !End of the abilint section
485 
486  implicit none
487 
488 !Arguments ------------------------------------
489 !scalars
490  type(slatrad_cshell_t),intent(inout) :: Slatrad3l(:)
491 
492 !Local variables-------------------------------
493  integer :: ii
494 ! *********************************************************************
495 
496  !@slatrad_cshell_t
497  do ii=1,SIZE(Slatrad3l)
498    if (allocated(Slatrad3l(ii)%rlphic_select)) then
499      ABI_FREE(Slatrad3l(ii)%rlphic_select)
500    end if
501    if (allocated(Slatrad3l(ii)%rlphic_int)) then
502      ABI_FREE(Slatrad3l(ii)%rlphic_int)
503    end if
504  end do
505 
506 end subroutine slatrad_cshell_free

m_paw_slater/slatrad_cshell_init [ Functions ]

[ Top ] [ m_paw_slater ] [ Functions ]

NAME

  slatrad_cshell_init

FUNCTION

  Initialize the structure storing the radial part of Slater"s integrals.

INPUTS

  ln2_size=Number of symmetrical (l,n) channels
  Pawrad<pawrad_type>=paw radial mesh and related data
  Pawtab<pawtab_type>=paw tabulated starting data
  Atm<atomorb_type>=Object containing core orbitals.
  Atmrad<pawrad_type>=paw radial mesh and related data for the atom.
  kln_mask

OUTPUT

  Slatrad3l<slatrad_cshell_t>=The object completely initialized.

PARENTS

      m_paw_slater

CHILDREN

      klmn2ijlmn,pawio_print_ij,slatrad_free,slatrad_init,wrtout

SOURCE

537 subroutine slatrad_cshell_init(Slatrad3l,ln2_size,Pawrad,Pawtab,Atm,Atmrad,kln_mask)
538 
539 
540 !This section has been created automatically by the script Abilint (TD).
541 !Do not modify the following lines by hand.
542 #undef ABI_FUNC
543 #define ABI_FUNC 'slatrad_cshell_init'
544 !End of the abilint section
545 
546  implicit none
547 
548 !Arguments ------------------------------------
549 !scalars
550  integer,intent(in) :: ln2_size
551 !arrays
552  integer,optional,intent(in) :: kln_mask(ln2_size)
553  type(atomorb_type),intent(in) :: Atm
554  type(pawrad_type),target,intent(in) :: Atmrad,Pawrad
555  type(pawtab_type),target,intent(in) :: Pawtab
556  type(slatrad_cshell_t),intent(out) :: Slatrad3l(ln2_size)
557 
558 !Local variables ---------------------------------------
559 !scalars
560  integer :: cmesh_size,dmesh_size
561  integer :: il,iln,ilnc,isl,in,jl,jln,jn,kln,ll,lnc_size
562  integer :: lslat_max,lslat_min,lc_max,nintg
563  integer :: lmn_size,lmn2_size,do_spline,ln_size,whichdenser,isppol
564  real(dp) :: intg,intg1,ybcbeg,ybcend
565  logical :: hasameq
566 !arrays
567  integer,allocatable :: kln2ln(:,:)
568  integer, pointer :: indklmn(:,:),indlmn(:,:)
569  real(dp),allocatable :: ff1(:),ff2(:),tmp_integrals(:)
570  real(dp),ABI_CONTIGUOUS pointer :: phi_i(:),phi_j(:)
571  real(dp),allocatable,target :: phi_spl(:,:)
572  real(dp),allocatable :: der(:),ypp(:)
573  real(dp),ABI_CONTIGUOUS pointer :: crad(:),drad(:),phi_in(:)
574 
575 ! *************************************************************************
576 
577  ABI_CHECK(ln2_size==Pawtab%ij_size,"Wrong ln2_size")
578  if (PRESENT(kln_mask)) then
579    MSG_ERROR("kln_mask is present")
580  end if
581 
582  !@slatrad_cshell_t
583  lmn_size  = Pawtab%lmn_size
584  lmn2_size = Pawtab%lmn2_size
585  ln_size   = Pawtab%basis_size
586 
587  lnc_size = Atm%ln_size
588  lc_max   = Atm%l_max
589 
590  call pawrad_isame(Atmrad,Pawrad,hasameq,whichdenser)
591 
592  do_spline=0
593  if (.not.hasameq) then
594    do_spline=1
595    if (whichdenser/=1) &
596 &    MSG_COMMENT("Pawrad is denser than Atmrad!")
597  else
598    ABI_CHECK(whichdenser==1,"Pawrad is denser than Atmrad!")
599  end if
600 
601  dmesh_size = Atmrad%mesh_size
602  cmesh_size = Pawtab%mesh_size
603 
604  drad  => Atmrad%rad(1:dmesh_size)
605  crad  => Pawrad%rad(1:cmesh_size)
606 
607  ! === Spline valence basis set onto core mesh (natural spline) ===
608  if (do_spline==1) then
609    MSG_COMMENT("Splining in init_slatrad3l")
610    ABI_MALLOC(phi_spl,(dmesh_size,ln_size))
611    ABI_MALLOC(der,(cmesh_size))
612    ABI_MALLOC(ypp,(cmesh_size))
613 
614    do iln=1,ln_size
615      phi_in => Pawtab%phi(:,iln)
616      ypp(:) = zero; ybcbeg = zero; ybcend = zero
617      call spline(crad,phi_in,cmesh_size,ybcbeg,ybcend,ypp)
618      call splint(cmesh_size,crad,phi_in,ypp,dmesh_size,drad,phi_spl(:,iln))
619    end do
620 
621    ABI_FREE(der)
622    ABI_FREE(ypp)
623  end if
624 
625  indlmn  => Pawtab%indlmn(1:6,1:lmn_size)
626  indklmn => Pawtab%indklmn(1:8,1:lmn2_size)
627 
628  ABI_MALLOC(kln2ln,(6,ln2_size))
629 
630  call make_kln2ln(lmn_size,lmn2_size,ln2_size,indlmn,indklmn,kln2ln)
631 
632  ABI_MALLOC(ff1,(dmesh_size))
633  ABI_MALLOC(ff2,(dmesh_size))
634 
635  ! * Loop over the upper triangle of the [(in,il), (jn,il)] matrix.
636  ABI_CHECK(Atm%nsppol==1,"nsppol==2 not tested")
637 
638  do isppol=1,Atm%nsppol
639    do kln=1,ln2_size
640      il  = kln2ln(1,kln)
641      jl  = kln2ln(2,kln)
642      in  = kln2ln(3,kln)
643      jn  = kln2ln(4,kln)
644      iln = kln2ln(5,kln)
645      jln = kln2ln(6,kln)
646 
647      lslat_max = MAX((il+lc_max),(jl+lc_max))       - 1   ! These are indeces, not l-values.
648      !lslat_min = MIN(ABS(il-lc_max),ABS(jl-lc_max)) + 1
649      lslat_min = 1 ! FIXME find better way
650 
651      Slatrad3l(kln)%lnc_size    = lnc_size
652      Slatrad3l(kln)%lslat_min   = lslat_min
653      Slatrad3l(kln)%lslat_max   = lslat_max
654 
655      Slatrad3l(kln)%nrlphic_int = 0
656 
657      ABI_MALLOC(Slatrad3l(kln)%rlphic_select,(lslat_max,lnc_size))
658      Slatrad3l(kln)%rlphic_select(:,:) = 0
659 
660      !if (PRESENT(kln_mask)) then  !FIXME THIS IS WRONG, move it below in case
661      ! if (kln_mask(kln)==0) CYCLE
662      !end if
663 
664      if (do_spline==1) then
665        MSG_COMMENT("Performing spline of valence phi")
666        phi_i => phi_spl(:,iln)
667        phi_j => phi_spl(:,jln)
668      else
669        phi_i => Pawtab%phi(:,iln)
670        phi_j => Pawtab%phi(:,jln)
671      end if
672 
673      ! * Loop over (n,l) channels for Atom orbitals
674      ABI_MALLOC(tmp_integrals,(lslat_max*lnc_size))
675      tmp_integrals(:) = zero
676      nintg=0
677 
678      do ilnc=1,lnc_size
679        ! phicore => Atm%phi(:,ilnc,isppol)
680        ff1 = phi_i * Atm%phi(:,ilnc,isppol)
681        ff2 = phi_j * Atm%phi(:,ilnc,isppol)
682        do isl=lslat_min,lslat_max ! L coming from Coulomb expansion
683          ll = isl-1
684          call calc_slatradl(ll,dmesh_size,ff2,ff1,Atmrad,intg1)
685          call calc_slatradl(ll,dmesh_size,ff1,ff2,Atmrad,intg)
686 
687          !call calc_slatradl(ll,cmesh_size,ff2,ff1,Pawrad,intg1)
688          !call calc_slatradl(ll,cmesh_size,ff1,ff2,Pawrad,intg)
689 
690          if (ABS(intg1-intg)>tol6) write(std_out,*)"DEBUG ",ll,il,in,jl,jn,intg1,intg
691 
692          ! * Store results
693          if (ABS(intg)>=tol12) then
694            nintg = nintg +1
695            Slatrad3l(kln)%rlphic_select(isl,ilnc) = nintg
696            tmp_integrals(nintg) = intg
697          end if
698        end do !ll
699      end do ! ilnc
700 
701      ! Finalize the object
702      Slatrad3l(kln)%nrlphic_int = nintg
703      ABI_MALLOC(Slatrad3l(kln)%rlphic_int,(nintg))
704      if (nintg>0) Slatrad3l(kln)%rlphic_int(:) = tmp_integrals(1:nintg)
705 
706      ABI_FREE(tmp_integrals)
707    end do !kln
708  end do !isppol
709 
710  ABI_FREE(ff1)
711  ABI_FREE(ff2)
712  ABI_FREE(kln2ln)
713 
714  if (do_spline==1)  then
715    ABI_FREE(phi_spl)
716  end if
717 
718 end subroutine slatrad_cshell_init

m_paw_slater/slatrad_cshell_t [ Types ]

[ Top ] [ m_paw_slater ] [ Types ]

NAME

  slatrad_cshell_t

FUNCTION

  Object used to store the set of radial integrals:

  $ \dfrac{4\pi}{2L+1} \times
    \int \phi_{\ni\li}(1) \phic_{\nc\lc}(1) \dfrac{r_<^L}{r_>^{L+1}} \phic_{\nc\lc}(2) \phi_{\nj\lj} d1d2 $

  for given (in,il) and (jn,jl) as a function L = |il-jl|, |il-jl|+2, ..., |il+il| and ilnc = (lc,nc).
  This (less general) type of radial integral is needed to evaluate the Exchange term generated
  by a closed-shell atom. In the equation, (\nc,\lc) are the set of angular quantum number associated
  to core electrons while  (lsl,msl) comes from the expansion of 1/|r1-r2|.
  Since the F is invariant under exchange of i and j we use an array of structures indexed
  by kln = (iln,jln) in packed form.

SOURCE

186  type,public :: slatrad_cshell_t
187 
188   integer :: lnc_size
189   ! Number of (n,l) channel for core orbitals.
190 
191   integer :: lslat_max
192   ! Max l+1 in the expansion of the Coulomb potential
193 
194   integer :: lslat_min
195   ! Min l+1 in the expansion of the Coulomb potential
196 
197   integer :: nrlphic_int
198   ! The number of non-zero integrals stored in rlphic_int.
199 
200   integer,allocatable :: rlphic_select(:,:)
201   ! rlphic_select(lslat_max,lnc_size)  TODO should be allocated with lslat_min:lslat_max
202   ! Index of the non-zero integrals in rlphic_int, 0 if rlphic_int has not been
203   ! calculated thanks to selection rules coming from the angular integration.
204 
205   real(dp),allocatable :: rlphic_int(:)
206   ! rlphic_int(1:nrlphic_int)
207   ! The integrals:
208   ! \dfrac{4\pi}{2L+1} \int \phi_{\ni\li}(1) \phic_{nc\lc}(1) \dfrac{r_<^L}{r_>^{L+1}} \phic_{nc\lc}(2) \phi_{\nj\lj} d1d2
209   ! for given (in,il) and (jn,jl) as a function L = |il-jl|, |il-jl|+2, ..., |il+il| and ilnc = (lc,nc).
210 
211  end type slatrad_cshell_t
212 
213 
214 CONTAINS  !========================================================================================

m_paw_slater/slatrad_free_0D [ Functions ]

[ Top ] [ m_paw_slater ] [ Functions ]

NAME

  slatrad_free_0D

FUNCTION

  Free the dynamic memory allocated in a structure of type slatrad_t

PARENTS

      m_paw_slater

CHILDREN

      klmn2ijlmn,pawio_print_ij,slatrad_free,slatrad_init,wrtout

SOURCE

1016 subroutine slatrad_free_0D(Slatrad)
1017 
1018 
1019 !This section has been created automatically by the script Abilint (TD).
1020 !Do not modify the following lines by hand.
1021 #undef ABI_FUNC
1022 #define ABI_FUNC 'slatrad_free_0D'
1023 !End of the abilint section
1024 
1025  implicit none
1026 
1027 !Arguments ------------------------------------
1028 !scalars
1029  type(slatrad_t),intent(inout) :: Slatrad
1030 
1031 ! *********************************************************************
1032 
1033  !@slatrad_t
1034  if (allocated(Slatrad%intgl_select)) then
1035    ABI_FREE(Slatrad%intgl_select)
1036  end if
1037  if (allocated(Slatrad%intgl)) then
1038    ABI_FREE(Slatrad%intgl)
1039  end if
1040 
1041 end subroutine slatrad_free_0D

m_paw_slater/slatrad_free_1D [ Functions ]

[ Top ] [ m_paw_slater ] [ Functions ]

NAME

  slatrad_free_1D

FUNCTION

  Free the dynamic memory allocated in a structure of type slatrad_t

PARENTS

CHILDREN

      klmn2ijlmn,pawio_print_ij,slatrad_free,slatrad_init,wrtout

SOURCE

1060 subroutine slatrad_free_1D(Slatrad)
1061 
1062 
1063 !This section has been created automatically by the script Abilint (TD).
1064 !Do not modify the following lines by hand.
1065 #undef ABI_FUNC
1066 #define ABI_FUNC 'slatrad_free_1D'
1067 !End of the abilint section
1068 
1069  implicit none
1070 
1071 !Arguments ------------------------------------
1072 !scalars
1073  type(slatrad_t),intent(inout) :: Slatrad(:)
1074 
1075 !Local variables-------------------------------
1076  integer :: ii
1077 ! *********************************************************************
1078 
1079  do ii=1,SIZE(Slatrad)
1080    call slatrad_free_0D(Slatrad(ii))
1081  end do
1082 
1083 end subroutine slatrad_free_1D

m_paw_slater/slatrad_init [ Functions ]

[ Top ] [ m_paw_slater ] [ Functions ]

NAME

  slatrad_init

FUNCTION

  Initialize the structure storing the radial part of Slater"s integrals.

INPUTS

  which_intg= Option defining what kind of integrals have to be calculated:
   -- 1 for Slater integral of AE partial waves only.
          A = \frac{4\pi}{2L+1} \int u_i(1) u_j(1) \frac{r_<^L}{r_^{L+1}} u_k(2)u_l(2) d1d2
   -- 2 for Slater integral of (AE-PS) partial waves
          B = \frac{4\pi}{2L+1} \int  u_i(1)  u_j(1) \frac{r_<^L}{r_^{L+1}}  u_k(2)  u_l(2) d1d2 -
              \frac{4\pi}{2L+1} \int tu_i(1) tu_j(1) \frac{r_<^L}{r_^{L+1}} tu_k(2) tu_l(2) d1d2

   -- 3 for Slater integral of (AE-PS-compensation charges)
          C = A -
          \frac{4\pi}{2L+1} \int [tu_i(1) tu_j(1) + qhat^L_\ij r_1^2 g^L(1) ] \frac{r_<^L}{r_^{L+1}}*
                                 [tu_k(2) tu_l(2) + qhat^L_\kl r_2^2 g^L(2) ] d1d2

  where u = \phi/r; tu = \tphi/r; and qhat^L_\ij are related to q^\LM_\ij via
    q^\LM_\ij = \Gaunt_\ij^\LM qhat^L_\ij => qhat^L\ij = \int (u_i*u_j - tu_i*tu_j) r^L

  ln2_size=Number of symmetrical (l,n) channels for this atom type type.
  Pawrad<pawrad_type>=paw radial mesh and related data
  Pawtab<pawtab_type>=paw tabulated starting data

OUTPUT

  Slatrad4<slatrad_t>=The object completely initialized.

NOTES

  Slater integrals S_ij are invariant under exchage of the indeces,
  but the results reported by calc_slatradl are not due to numerical roundoff errors (err < 10^-9).
  However this does not cause any problem since only the upper triangle of the S_ij matrix
  is stored and used in the other routines.

PARENTS

      m_paw_slater

CHILDREN

      klmn2ijlmn,pawio_print_ij,slatrad_free,slatrad_init,wrtout

SOURCE

1133 subroutine slatrad_init(Slatrad4,which_intg,ln2_size,Pawrad,Pawtab)
1134 
1135 
1136 !This section has been created automatically by the script Abilint (TD).
1137 !Do not modify the following lines by hand.
1138 #undef ABI_FUNC
1139 #define ABI_FUNC 'slatrad_init'
1140 !End of the abilint section
1141 
1142  implicit none
1143 
1144 !Arguments ------------------------------------
1145 !scalars
1146  integer,intent(in) :: ln2_size,which_intg
1147 !arrays
1148  type(pawrad_type),target,intent(in) :: Pawrad
1149  type(pawtab_type),target,intent(in) :: Pawtab
1150  type(slatrad_t),intent(out) :: Slatrad4(ln2_size*(ln2_size+1)/2)
1151 
1152 !Local variables ---------------------------------------
1153 !scalars
1154  integer :: mesh_size,il,iln,isl,in,jl,jln,jn,sln1,sln2,l_slat
1155  integer :: kn,kl,ln,ll,kln,lln,lslat_max,lslat_min,nintgl
1156  integer :: lmn_size,lmn2_size,ln_size,slt_idx
1157  real(dp) :: ae_intg,ps_intg,pshat_intg,intg,tqij_L,tqkl_L !intg1
1158  character(len=500) :: msg
1159 !arrays
1160  integer,allocatable :: kln2ln(:,:)
1161  integer, pointer :: indklmn(:,:),indlmn(:,:)
1162  real(dp),allocatable :: uiuj(:),ukul(:),tuituj(:),tuktul(:),tuituj_tqgl(:),tuktul_tqgl(:)
1163  real(dp),allocatable :: tmp_integrals(:),ff(:)
1164  real(dp),ABI_CONTIGUOUS pointer :: phi_i(:),phi_j(:),phi_k(:),phi_l(:)
1165  real(dp),ABI_CONTIGUOUS pointer :: tphi_i(:),tphi_j(:),tphi_k(:),tphi_l(:)
1166  real(dp),ABI_CONTIGUOUS pointer :: shapefunc(:),rad(:)
1167 
1168 ! *************************************************************************
1169 
1170  DBG_ENTER("COLL")
1171 
1172  ABI_CHECK(ln2_size==Pawtab%ij_size,"Wrong ln2_size")
1173 
1174  if ( ALL(which_intg /= (/1,2,3/)) ) then
1175    write(msg,'(a,i0)')"Wrong value for which_intg: ",which_intg
1176    MSG_ERROR(msg)
1177  end if
1178 
1179  !@slatrad_t
1180  lmn_size   = Pawtab%lmn_size
1181  lmn2_size  = Pawtab%lmn2_size
1182  ln_size    = Pawtab%basis_size
1183  mesh_size  = Pawtab%mesh_size
1184  !
1185  ! Useful table for looping.
1186  indlmn  => Pawtab%indlmn(1:6,1:lmn_size)
1187  indklmn => Pawtab%indklmn(1:8,1:lmn2_size)
1188 
1189  ABI_MALLOC(kln2ln,(6,ln2_size))
1190  call make_kln2ln(lmn_size,lmn2_size,ln2_size,indlmn,indklmn,kln2ln)
1191 
1192  ABI_MALLOC(uiuj,(mesh_size))
1193  ABI_MALLOC(ukul,(mesh_size))
1194  ABI_MALLOC(ff,(mesh_size))
1195  ABI_MALLOC(tuituj,(mesh_size))
1196  ABI_MALLOC(tuktul,(mesh_size))
1197  ABI_MALLOC(tuituj_tqgl,(mesh_size))
1198  ABI_MALLOC(tuktul_tqgl,(mesh_size))
1199  rad => Pawrad%rad
1200  !
1201  ! * Loop over (k,l) channels in packed form.
1202  do sln2=1,ln2_size
1203    kl  = kln2ln(1,sln2)
1204    ll  = kln2ln(2,sln2)
1205    kn  = kln2ln(3,sln2)
1206    ln  = kln2ln(4,sln2)
1207    kln = kln2ln(5,sln2)
1208    lln = kln2ln(6,sln2)
1209    !write(std_out,*)"sln2, kln, lln",sln2,kln,lln
1210 
1211    phi_k  => Pawtab%phi (:,kln)
1212    tphi_k => Pawtab%tphi(:,kln)
1213 
1214    phi_l  => Pawtab%phi (:,lln)
1215    tphi_l => Pawtab%tphi(:,lln)
1216    !
1217    ! * Loop over (i,j) channels in packed form AND only for the upper triangle of sln2, sln1
1218    do sln1=1,sln2
1219      il  = kln2ln(1,sln1)
1220      jl  = kln2ln(2,sln1)
1221      in  = kln2ln(3,sln1)
1222      jn  = kln2ln(4,sln1)
1223      iln = kln2ln(5,sln1)
1224      jln = kln2ln(6,sln1)
1225      !write(std_out,*)"sln1, iln, jln",sln1,iln,jln
1226 
1227      slt_idx = sln1 + sln2*(sln2-1)/2 ! index for packed storage.
1228 
1229      phi_i  => Pawtab%phi (:,iln)
1230      tphi_i => Pawtab%tphi(:,iln)
1231 
1232      phi_j  => Pawtab%phi (:,jln)
1233      tphi_j => Pawtab%tphi(:,jln)
1234 
1235      lslat_min = MAX(ABS(il-jl),ABS(kl-ll)) + 1  ! We use indeces not l-values.
1236      lslat_max = MIN((il+jl),(kl+ll)) - 1
1237 
1238      !lslat_min = MIN(ABS(il-jl),ABS(kl-ll)) + 1
1239      !lslat_max = MAX((il+jl),(kl+ll)) - 1
1240 
1241      Slatrad4(slt_idx)%lslat_min = lslat_min
1242      Slatrad4(slt_idx)%lslat_max = lslat_max
1243 
1244      Slatrad4(slt_idx)%iln = iln
1245      Slatrad4(slt_idx)%jln = jln
1246      Slatrad4(slt_idx)%kln = kln
1247      Slatrad4(slt_idx)%lln = lln
1248 
1249      ABI_MALLOC(Slatrad4(slt_idx)%intgl_select,(lslat_min:lslat_max))
1250      Slatrad4(slt_idx)%intgl_select(:) = 0
1251      Slatrad4(slt_idx)%nintgl          = 0
1252 
1253      if (lslat_min > lslat_max) then
1254        ! e.g. (1 2) (1 1). Due to angular selection rules, this integral do not contribue
1255        !write(std_out,*)"lslat_min, lslat_max",lslat_min,lslat_max
1256        !write(std_out,*)"il,jl,kl,ll",il,jl,kl,ll
1257        !MSG_ERROR("")
1258        ABI_MALLOC(Slatrad4(slt_idx)%intgl,(0))
1259        CYCLE
1260      end if
1261 
1262      uiuj   =  phi_i *  phi_j  ! The AE part.
1263      ukul   =  phi_k *  phi_l
1264      tuituj = tphi_i * tphi_j  ! The pseudized part.
1265      tuktul = tphi_k * tphi_l
1266      !
1267      ! Calculate L-depedent integrals where L come from the expansion the Coulomb interaction.
1268      ABI_MALLOC(tmp_integrals,(MAX(lslat_min,lslat_max)))
1269      tmp_integrals=zero
1270      nintgl=0
1271 
1272      do isl=lslat_min,lslat_max
1273      !do isl=lslat_min,lslat_max,2  ! TODO Here I can reduce the number of iterations using a step of 2.
1274        l_slat = isl-1
1275        call calc_slatradl(l_slat,mesh_size,uiuj,ukul,Pawrad,ae_intg)
1276        intg = ae_intg
1277 
1278 #if 0
1279        call calc_slatradl(l_slat,mesh_size,ukul,uiuj,Pawrad,intg1)
1280        if (ABS(intg1-ae_intg)>tol12) then
1281          write(msg,'(a,es16.8)')"s_ij and s_ij differ by ",intg1-ae_intg
1282          MSG_WARNING(msg)
1283        end if
1284 #endif
1285        if (which_intg == 2) then ! Subtract the pseudo part.
1286          call calc_slatradl(l_slat,mesh_size,tuituj,tuktul,Pawrad,ps_intg)
1287          intg = intg - ps_intg
1288 
1289        else if (which_intg == 3) then ! Subtract (pseudo + compensation charges)
1290          !
1291          ! Evaluate tqij_L and tqkl_L (without M-dependent part).
1292          ff(1)=zero
1293          ff(2:mesh_size)=(pawtab%phiphj(2:mesh_size,sln1)-pawtab%tphitphj(2:mesh_size,sln1))*rad(2:mesh_size)**l_slat
1294          if (l_slat==0.and.kl==1.and.ll==1) then
1295            call pawrad_deducer0(ff,mesh_size,pawrad)
1296          end if
1297          call simp_gen(tqij_L,ff,pawrad)
1298 
1299          ff(1)=zero
1300          ff(2:mesh_size)=(pawtab%phiphj(2:mesh_size,sln2)-pawtab%tphitphj(2:mesh_size,sln2))*rad(2:mesh_size)**l_slat
1301          if (l_slat==0.and.il==1.and.jl==1) then
1302            call pawrad_deducer0(ff,mesh_size,pawrad)
1303          end if
1304          call simp_gen(tqkl_L,ff,pawrad)
1305 
1306          shapefunc   => Pawtab%shapefunc(:,isl)  ! Recheck this part, in particular the convention
1307          tuituj_tqgl = tuituj + tqij_L * shapefunc * rad**2
1308          tuktul_tqgl = tuktul + tqkl_L * shapefunc * rad**2
1309 
1310          call calc_slatradl(l_slat,mesh_size,tuituj_tqgl,tuktul_tqgl,Pawrad,pshat_intg)
1311          intg = intg - pshat_intg
1312        end if
1313        !
1314        ! * Store results
1315        if (ABS(intg)>=tol12) then
1316          nintgl = nintgl +1
1317          Slatrad4(slt_idx)%intgl_select(isl) = nintgl
1318          tmp_integrals(nintgl) = intg
1319        end if
1320      end do !isl
1321      !
1322      ! Finalize the object.
1323      Slatrad4(slt_idx)%nintgl = nintgl
1324      ABI_MALLOC(Slatrad4(slt_idx)%intgl,(nintgl))
1325      if (nintgl>0) Slatrad4(slt_idx)%intgl(:) = tmp_integrals(1:nintgl)
1326      ABI_FREE(tmp_integrals)
1327    end do !sln1
1328  end do !sln2
1329  !
1330  ! Free memory
1331  ABI_FREE(kln2ln)
1332  ABI_FREE(uiuj)
1333  ABI_FREE(ukul)
1334  ABI_FREE(tuituj)
1335  ABI_FREE(tuktul)
1336  ABI_FREE(ff)
1337  ABI_FREE(tuituj_tqgl)
1338  ABI_FREE(tuktul_tqgl)
1339 
1340  DBG_EXIT("COLL")
1341 
1342 end subroutine slatrad_init

m_paw_slater/slatrad_t [ Types ]

[ Top ] [ m_paw_slater ] [ Types ]

NAME

  slatrad_t

FUNCTION

  Object used to store radial integrals of the form.

  $ F_{ijkl}^L = \dfrac{4\pi}{2L+1} \int u_i(1) u_j(1) \dfrac{r_<^L}{r_>^{L+1}} u_k(2) u_l(2) d1d2 $

  for a given quadruple (i,j,k,l) as a function L \in [L_min, L_max].
  i,j,k,l are shorthand indeces for (nn,ll) quantum numbers.

NOTES

   Basic symmetry properties:
   1) invariant under the exchange (i<-->j) and (k<-->l).
   2) invariant under the exchange (i,j) <--> (k,l).

  Memory saving is achieved by storing the upper triangle of the (ij) (kl) matrix
  and, for each dimension, only the upper triangle of the two matrices (iln,jln) (kln,lln).

  Some matrix elements will never contribute to <ij|1/|1-2||kl> due to selection rules
  introduced by the integration of the angular part.

SOURCE

 84  type,public :: slatrad_t
 85 
 86   integer :: iln,jln,kln,lln
 87   ! The (l,n) indeces associated to the partial waves.
 88 
 89   integer :: lslat_min
 90   ! Min l+1 in the expansion of the Coulomb potential.
 91 
 92   integer :: lslat_max
 93   ! Max l+1 in the expansion of the Coulomb potential.
 94 
 95   integer :: nintgl
 96   ! The number of non-zero integrals stored in intgl.
 97 
 98   integer,allocatable :: intgl_select(:)
 99   ! intgl_select(lslat_min:lslat_max)
100   ! Index of the non-zero integrals in intgl, 0 if intgl has not been
101   ! calculated thanks to selection rules coming from the angular integration.
102 
103   real(dp),allocatable :: intgl(:)
104   ! intgl(1:nintgl)
105   ! The integrals:
106   ! \dfrac{4\pi}{2L+1} \int \phi_{\ni\li}(1) \phi_{\nj\lj}(1) \dfrac{r_<^L}{r_>^{L+1}} \phi_{\nk\lk}(2) \phi_{\nl\ll} d1d2
107   ! for given (i,j,k,l) as a function L = |il-jl|, |il-jl|+2, ..., |il+il| and ilnc = (lc,nc).
108 
109  end type slatrad_t
110 
111  public :: slatrad_init      ! Creation method
112  public :: slatrad_free      ! Free memory

m_paw_slater/summ_2gaunt [ Functions ]

[ Top ] [ m_paw_slater ] [ Functions ]

NAME

   summ_2gaunt

FUNCTION

  Helper function returning \sum_M G_{ij}^{LM} G_{kl}^{LM}

INPUTS

  ij_lm=index of (i_lm,j_lm) element in packed form.
  kl_lm=index of (k_lm,l_lm) element in packed form.
  ll_idx=Index for L (thus L+1).
  Pawang<type(pawang_type)>=paw angular mesh and related data

PARENTS

CHILDREN

SOURCE

1561 function summ_2gaunt(Pawang,ij_lm,kl_lm,ll_idx)
1562 
1563 
1564 !This section has been created automatically by the script Abilint (TD).
1565 !Do not modify the following lines by hand.
1566 #undef ABI_FUNC
1567 #define ABI_FUNC 'summ_2gaunt'
1568 !End of the abilint section
1569 
1570  implicit none
1571 
1572 !Arguments ------------------------------------
1573 !scalars
1574  integer,intent(in) :: ij_lm,kl_lm,ll_idx
1575  real(dp) :: summ_2gaunt
1576  type(pawang_type),intent(in) :: Pawang
1577 !arrays
1578 
1579 !Local variables-------------------------------
1580 !scalars
1581  integer :: ignt1,ignt2,idx_LM,max_klm,mm,ii,ll
1582  character(len=500) :: msg
1583 
1584 !************************************************************************
1585 
1586  ! FIXME: size of gntselect depends on pawxcdev!
1587  ! Consistency check on input arguments.
1588  max_klm = pawang%l_max**2*(pawang%l_max**2+1)/2
1589  if (ij_lm>max_klm.or.kl_lm>max_klm.or.ij_lm<1.or.kl_lm<1.or.&
1590 &    ll_idx>pawang%l_size_max.or.ll_idx<1) then
1591    write(msg,'(a,3i0)')"Wrong indeces, check pawxcdev ",ij_lm,kl_lm,ll_idx
1592    MSG_ERROR(msg)
1593  end if
1594 
1595  ll = ll_idx-1
1596  summ_2gaunt=zero; ii=0
1597  do mm=-ll,ll
1598    idx_LM = 1 + ll**2 + ll + mm
1599    ignt1 = Pawang%gntselect(idx_LM,ij_lm)
1600    ignt2 = Pawang%gntselect(idx_LM,kl_lm)
1601    if (ignt1>0 .and. ignt2>0) then
1602      summ_2gaunt = summ_2gaunt + Pawang%realgnt(ignt1)*Pawang%realgnt(ignt2)
1603      ii=ii+1
1604      write(std_out,'(a,4(i2,1x),f8.5,i2)')"ll, mm, ij_lm, kl_lm: ",ll,mm,ij_lm,kl_lm,summ_2gaunt,ii
1605      if (ii/=1) MSG_WARNING("ii>1")
1606    end if
1607  end do
1608 
1609 end function summ_2gaunt