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-2024 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

25 #if defined HAVE_CONFIG_H
26 #include "config.h"
27 #endif
28 
29 #include "abi_common.h"
30 
31 MODULE m_paw_slater
32 
33  use defs_basis
34  use m_abicore
35  use m_errors
36  use m_splines
37 
38  use m_fstrings,     only : basename
39  use m_paw_atomorb,  only : atomorb_type, init_atomorb, print_atomorb, destroy_atomorb, get_overlap
40  use m_crystal,      only : crystal_t
41  use m_paw_io,       only : pawio_print_ij
42  use m_pawang,       only : pawang_type
43  use m_paw_sphharm,   only : 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 part 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
  qphase=2 if dij contains a exp(-i.q.r) phase (as in the q<>0 RF case), 1 if not
  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}

SOURCE

1274 subroutine paw_dijhf(ndij,cplex_dij,qphase,lmn2_size_max,my_natom,ntypat,Pawtab,Pawrad,Pawang,Pawrhoij,&
1275 &                    sigx_dij,pawprtvol)
1276 
1277  implicit none
1278 
1279 !Arguments ------------------------------------
1280 !scalars
1281  integer,intent(in) :: pawprtvol,ndij,cplex_dij,lmn2_size_max,my_natom,ntypat,qphase
1282  type(pawang_type),intent(in) :: Pawang
1283 !arrays
1284  real(dp),target,intent(out) :: sigx_dij(cplex_dij*lmn2_size_max,ndij,my_natom) !TODO use ragged arrays pawij?
1285  type(pawtab_type),intent(in) :: Pawtab(ntypat)
1286  type(pawrad_type),intent(in) :: Pawrad(ntypat)
1287  type(pawrhoij_type),intent(in) :: Pawrhoij(my_natom)
1288 
1289 !Local variables ---------------------------------------
1290 !scalars
1291  integer,parameter :: cplex=1    ! FIXME preliminary implementation
1292  integer :: cplex_rhoij,iatom,iq,iq0_dij,iq0_rhoij,itypat,lmn_size,lmn2_size,ispden,nspden,ln2_size
1293  integer :: lm2_size !,isppol ln_size,
1294  integer :: irhoij,jrhoij
1295  integer :: rho_lmn !,rho_klm,rho_kln,rho_lmin,rho_lmax,rho_iln,rho_jln
1296  integer :: klmn
1297  integer :: i_lmn,j_lmn,k_lmn,l_lmn
1298  integer :: which_intg,l_max,opt_l
1299  real(dp) :: ro,slt_ikjl,slt_iljk
1300  !character(len=500) :: msg
1301 !arrays
1302  integer :: opt_l_index(0,0),pack2ij(0)
1303  real(dp), ABI_CONTIGUOUS pointer :: sigx_atm(:,:)
1304  type(slatrad_t),allocatable :: Slatrad4(:)
1305 
1306 ! *************************************************************************
1307 
1308  DBG_ENTER("COLL")
1309 
1310  ABI_CHECK(ndij/=4,"ndij=4 not coded")
1311  ABI_CHECK(cplex_dij==1,"cplex_dij/=2 not coded")
1312  ABI_CHECK(lmn2_size_max==MAXVAL(Pawtab(:)%lmn2_size),"Wrong lmn2_size_max")
1313 
1314  if (my_natom>0) then
1315    if (pawrhoij(1)%qphase<cplex) then
1316      ABI_BUG('Must have pawrhoij()%qphase >= cplex !')
1317    end if
1318  end if
1319 
1320  sigx_dij=zero
1321 
1322  do iatom=1,my_natom
1323    itypat   =Pawrhoij(iatom)%itypat
1324    lmn_size =Pawtab(itypat)%lmn_size
1325    lmn2_size=Pawtab(itypat)%lmn2_size
1326    l_max    =(Pawtab(itypat)%l_size+1)/2
1327    lm2_size = (l_max**2)*(l_max**2+1)/2
1328    cplex_rhoij=Pawrhoij(iatom)%cplex_rhoij
1329    !write(std_out,*)"in atom ",iatom,"lm2_size=",lm2_size
1330 
1331    ! Calculate Slater integral for this atom type.
1332    ! TODO obviously these tables should be stored in Pawtab!
1333    ln2_size = Pawtab(itypat)%ij_size
1334    ABI_MALLOC(Slatrad4,(ln2_size*(ln2_size+1)/2))
1335    which_intg=3
1336    call slatrad_init(Slatrad4,which_intg,ln2_size,Pawrad(itypat),Pawtab(itypat))
1337 
1338    sigx_atm => sigx_dij(:,:,iatom)
1339 
1340 !  Loop over phase exp(iqr) phase real/imaginary part, if any
1341    do iq=1,qphase
1342      !First loop: we store the real part in dij(1 -> lmn2_size)
1343      !2nd loop: we store the imaginary part in dij(lmn2_size+1 -> 2*lmn2_size)
1344      iq0_dij=merge(0,cplex_dij*lmn2_size,iq==1)
1345      iq0_rhoij=merge(0,cplex_rhoij*lmn2_size,iq==1)
1346 
1347      ! * Loop over spin components.
1348      nspden=ndij
1349      do ispden=1,ndij
1350        !
1351        ! ============================================================
1352        ! ==== Summing over the non-zero lk channels of \rho_{lk} ====
1353        ! ============================================================
1354        jrhoij=1+iq0_rhoij
1355        do irhoij=1,pawrhoij(iatom)%nrhoijsel
1356          rho_lmn=pawrhoij(iatom)%rhoijselect(irhoij)
1357 
1358          ! check wheter rho_lmin is consistent with the Indexing used in slatrad
1359          !rho_klm =pawtab(itypat)%indklmn(1,rho_lmn)
1360          !rho_kln =pawtab(itypat)%indklmn(2,rho_lmn)
1361          !rho_lmin=pawtab(itypat)%indklmn(3,rho_lmn)
1362          !rho_lmax=pawtab(itypat)%indklmn(4,rho_lmn)
1363 
1364          ! Retrieve rhoij for this ispden.
1365          if (nspden/=2) then
1366            ro=pawrhoij(iatom)%rhoijp(jrhoij,ispden)
1367          else
1368            ABI_ERROR("Recheck this part")
1369            if (ispden==1) then
1370              ro=pawrhoij(iatom)%rhoijp(jrhoij,1) + pawrhoij(iatom)%rhoijp(jrhoij,2)
1371            else if (ispden==2) then
1372              ro=pawrhoij(iatom)%rhoijp(jrhoij,1)
1373            end if
1374          end if
1375          !
1376          ! Avoid double-counting the diagonal of rho.
1377          ro=ro*pawtab(itypat)%dltij(rho_lmn)*half
1378 
1379          call klmn2ijlmn(rho_lmn,lmn_size,k_lmn,l_lmn)
1380 
1381          ! Loop over the upper triangle of the D_{ij) matrix and accumulate:
1382          ! sum_\lk rho_\kl [ \Phi_{ikjl} + \Phi_{iljk} - \Phihat_{ikjl} - \Phihat_{iljk} ]
1383          do klmn=1,lmn2_size
1384            ! Calculate the indices in the Slatrad4 structure.
1385            call klmn2ijlmn(klmn,lmn_size,i_lmn,j_lmn)
1386 
1387            ! My formula
1388            slt_ikjl = slat_intg(Slatrad4,Pawtab(itypat),Pawang,i_lmn,k_lmn,j_lmn,l_lmn)
1389            slt_iljk = slat_intg(Slatrad4,Pawtab(itypat),Pawang,i_lmn,l_lmn,j_lmn,k_lmn)
1390 
1391            !slt_ikjl = slat_intg(Slatrad4,Pawtab(itypat),Pawang,i_lmn,k_lmn,l_lmn,j_lmn)
1392            !slt_iljk = slat_intg(Slatrad4,Pawtab(itypat),Pawang,i_lmn,l_lmn,k_lmn,j_lmn)
1393            !slt_iljk = slt_ikjl
1394 
1395            sigx_atm(klmn+iq0_dij,ispden) = sigx_atm(klmn,ispden) + ro * (slt_ikjl + slt_iljk)
1396          end do ! klmn
1397 
1398          jrhoij=jrhoij+cplex_rhoij
1399        end do ! irhoij
1400      end do ! iq
1401    end do ! ispden
1402 
1403    if (ABS(pawprtvol)>=1) then ! * Print values
1404      call wrtout(std_out,"   ************** Dij Fock ************ ",'COLL')
1405      opt_l=-1
1406      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)
1407    end if
1408 
1409    call slatrad_free(Slatrad4)
1410    ABI_FREE(Slatrad4)
1411  end do ! iatom
1412 
1413  ! Factor half cancels in the derivation wrt rho_ij.
1414  sigx_dij = - sigx_dij
1415 
1416  DBG_EXIT("COLL")
1417 
1418 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.

SOURCE

831 subroutine paw_mkdijexc_core(ndij,cplex_dij,lmn2_size_max,Cryst,Pawtab,Pawrad,dijexc_core,pawprtvol,filpsp)
832 
833  implicit none
834 
835 !Arguments ------------------------------------
836 !scalars
837  integer,intent(in) :: pawprtvol,ndij,cplex_dij,lmn2_size_max
838  type(crystal_t),intent(in) :: Cryst
839 !arrays
840  real(dp),intent(out) :: dijexc_core(cplex_dij*lmn2_size_max,ndij,Cryst%ntypat) !TODO use ragged arrays pawij?
841  character(len=fnlen) :: filpsp(Cryst%ntypat)
842  type(Pawtab_type),target,intent(in) :: Pawtab(Cryst%ntypat)
843  type(Pawrad_type),intent(in) :: Pawrad(Cryst%ntypat)
844 
845 !Local variables ---------------------------------------
846 !scalars
847  integer :: itypat,ic,ierr,lmn_size,lmn2_size,ln_size,isppol
848  real(dp) :: rcut
849  character(len=500) :: header,msg
850  character(len=fnlen) :: fcore,string
851 !arrays
852  integer,allocatable :: phi_indln(:,:)
853  real(dp),ABI_CONTIGUOUS pointer :: phi(:,:)
854  real(dp),allocatable :: overlap(:,:)
855  type(atomorb_type),allocatable :: Atm(:)
856  type(Pawrad_type),allocatable :: Radatm(:)
857 
858 ! *************************************************************************
859 
860  ABI_MALLOC(Atm,(Cryst%ntypat))
861  ABI_MALLOC(Radatm,(Cryst%ntypat))
862 
863  ABI_CHECK(ndij==1     ,"spinor+HF not available")
864  ABI_CHECK(cplex_dij==1,"spinor+HF not available")
865  ABI_CHECK(lmn2_size_max==MAXVAL(Pawtab(:)%lmn2_size),"Wrong lmn2_size_max")
866 
867  !allocate(dijexc_core(cplex_dij*lmn2_size_max,ndij,Cryst%ntypat)) !TODO use ragged arrays pawij?
868  dijexc_core=zero
869 
870  do itypat=1,Cryst%ntypat
871 
872    ! Read core orbitals for this atom type.
873    string = filpsp(itypat)
874    fcore = "CORE_"//TRIM(basename(string))
875    ic = INDEX (TRIM(string), "/" , back=.TRUE.) ! if string is a path, prepend path to fcore.
876    if (ic>0 .and. ic<LEN_TRIM(string)) fcore = filpsp(itypat)(1:ic)//TRIM(fcore)
877 
878    rcut=Pawtab(itypat)%rpaw
879    call init_atomorb(Atm(itypat),Radatm(itypat),rcut,fcore,pawprtvol,ierr)
880 
881    if (ierr/=0) then
882      msg = " Error reading core orbitals from file: "//TRIM(fcore)
883      ABI_ERROR(msg)
884    end if
885    write(header,'(a,i4,a)')" === Atom type = ",itypat," === "
886    call print_atomorb(Atm(itypat),header,unit=std_out,prtvol=pawprtvol)
887    !
888    ! * Calculate $ \<\phi_i|Sigma_x^\core|\phi_j\> $ for this atom type.
889    lmn_size  = Pawtab(itypat)%lmn_size
890    lmn2_size = Pawtab(itypat)%lmn2_size
891 
892    call paw_sigxcore(cplex_dij,lmn2_size,ndij,&
893 &    Pawrad(itypat),Pawtab(itypat),Atm(itypat),Radatm(itypat),dijexc_core(1:lmn2_size,:,itypat))
894 
895    ln_size =  Pawtab(itypat)%basis_size
896    phi     => Pawtab(itypat)%phi
897 
898    ABI_MALLOC(phi_indln,(2,ln_size))
899    call make_indln(lmn_size,ln_size,Pawtab(itypat)%indlmn(:,:),phi_indln)
900 
901    ABI_MALLOC(overlap,(Atm(itypat)%ln_size,ln_size))
902    isppol=1 ! hardcoded
903    call get_overlap(Atm(itypat),Radatm(itypat),Pawrad(itypat),isppol,ln_size,phi,phi_indln,overlap)
904 
905    ABI_FREE(phi_indln)
906    ABI_FREE(overlap)
907  end do ! ntypat
908 
909  ! Free memory
910  call pawrad_free(Radatm)
911  do itypat=1,Cryst%ntypat
912    call destroy_atomorb(Atm(itypat))
913  end do
914 
915  ABI_FREE(Atm)
916  ABI_FREE(Radatm)
917 
918 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)

SOURCE

693 subroutine paw_sigxcore(cplex_dij,lmn2_size,ndij,Pawrad,Pawtab,Atm,Atmrad,dijexc_core)
694 
695  implicit none
696 
697 !Arguments ------------------------------------
698 !scalars
699  integer,intent(in) :: lmn2_size,cplex_dij,ndij
700 !arrays
701  real(dp),intent(out) :: dijexc_core(cplex_dij*lmn2_size,ndij)
702  type(atomorb_type),intent(in) :: Atm
703  type(pawrad_type),intent(in) :: Atmrad
704  type(pawrad_type),intent(in) :: Pawrad
705  type(pawtab_type),target,intent(in) :: Pawtab
706 
707 !Local variables ---------------------------------------
708 !scalars
709  integer :: ilnc,ilc,lnc_size,l_max
710  integer :: lm2_size,ln_size,ln2_size,lmn_size
711  integer :: lm_size,klmn,kln,klm
712  integer :: lc_max,ilsl,isgg,israd,opt_l,pawprtvol
713  real(dp) :: tmp,sgg,intgrl
714 !character(len=500) :: msg
715 !arrays
716  integer :: opt_l_index(0,0),pack2ij(0)
717  integer,allocatable :: kln2ln(:,:),klm2lm(:,:)
718  integer, pointer :: indklmn(:,:),indlmn(:,:)
719  type(slatang_cshell_t),allocatable :: Slatang3l(:)
720  type(slatrad_cshell_t),allocatable :: Slatrad3l(:)
721 
722 ! *************************************************************************
723 
724  ! * Consistency check
725  ABI_CHECK(cplex_dij==1,"cplex_dij must be 1")
726 
727  ABI_CHECK(ndij==1,"ndij must be 1")
728 
729  ABI_CHECK(lmn2_size==Pawtab%lmn2_size,"Wrong lmn2_size")
730 
731  lmn_size  = Pawtab%lmn_size
732  ln_size   = Pawtab%basis_size
733  ln2_size  = Pawtab%ij_size
734  l_max     = (Pawtab%l_size-1)/2 +1
735  lm_size   = l_max**2
736  lm2_size  = lm_size*(lm_size+1)/2
737 
738  indlmn  => Pawtab%indlmn(1:6,1:lmn_size)
739  indklmn => Pawtab%indklmn(1:8,1:lmn2_size)
740 
741  ! * Setup of useful tables.
742  ABI_MALLOC(kln2ln,(6,ln2_size))
743  call make_kln2ln(lmn_size,lmn2_size,ln2_size,indlmn,indklmn,kln2ln)
744 
745  ABI_MALLOC(klm2lm,(6,lm2_size))
746  call make_klm2lm(lmn_size,lmn2_size,lm2_size,indlmn,indklmn,klm2lm)
747 
748  ! * Integrate angular part.
749  lnc_size  = Atm%ln_size
750  lc_max    = Atm%l_max
751 
752  ABI_MALLOC(Slatang3l,(lm2_size))
753  call slatang_cshell_init(Slatang3l,l_max,lm2_size,lc_max,klm2lm)
754 
755  ABI_FREE(klm2lm)
756 
757  ! * Integrate radial part.
758  ABI_MALLOC(Slatrad3l,(ln2_size))
759 
760  call slatrad_cshell_init(Slatrad3l,ln2_size,Pawrad,Pawtab,Atm,Atmrad)
761 
762  ! === Calculate matrix elements of Sigma_x^core ===
763  ! * $<\phi_i|\Sigma_x^\core|\phi_j>$
764  dijexc_core = zero
765  do klmn=1,lmn2_size
766    klm = Pawtab%indklmn(1,klmn)
767    kln = Pawtab%indklmn(2,klmn)
768    !
769    ! * Summing over (lc,nc) and lslat
770    tmp = zero
771    if (Slatang3l(klm)%nsggsel >0) then
772      do ilnc=1,Atm%ln_size
773        ilc = 1+Atm%indln(1,ilnc)
774        do ilsl=1,Slatang3l(klm)%lslat_max !FIXME check this
775          !do ilsl=Slatang3l(klm)%lslat_min,Slatang3l(klm)%lslat_max
776          isgg  = Slatang3l(klm)%sggselect(ilsl,ilc)
777          israd = Slatrad3l(kln)%rlphic_select(ilsl,ilnc)
778          if (isgg>0 .and. israd>0) then
779            sgg    = Slatang3l(klm)%sgg(isgg)
780            intgrl = Slatrad3l(kln)%rlphic_int(israd)
781            tmp = tmp + intgrl*sgg
782          end if
783        end do
784      end do
785    end if
786 
787    dijexc_core(klmn,1) = -tmp  ! Store results.
788  end do
789 
790  ! * Print values
791  call wrtout(std_out,"   ************** Dij Fock_core ************ ",'COLL')
792  opt_l=-1; pawprtvol=-1
793  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)
794 
795  ! * Free memory.
796  ABI_FREE(kln2ln)
797  call slatang_cshell_free(Slatang3l)
798  ABI_FREE(Slatang3l)
799  call slatrad_cshell_free(Slatrad3l)
800  ABI_FREE(Slatrad3l)
801 
802 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

SOURCE

1496 function slat_intg(Slatrad4,Pawtab,Pawang,i_lmn,j_lmn,k_lmn,l_lmn)
1497 
1498  implicit none
1499 
1500 !Arguments ------------------------------------
1501 !scalars
1502  integer,intent(in) :: i_lmn,j_lmn,k_lmn,l_lmn
1503  real(dp) :: slat_intg
1504  type(pawtab_type),intent(in) :: Pawtab
1505  type(pawang_type),intent(in) :: Pawang
1506 !arrays
1507  type(slatrad_t),intent(in) :: Slatrad4(:)
1508 
1509 !Local variables-------------------------------
1510 !scalars
1511  integer :: ij_lmn,kl_lmn,kl_ln,ij_lm,kl_lm,ilsum,ij_ln
1512  integer :: isel,slt_idx
1513  integer :: iln,jln,kln,lln,ii
1514  real(dp) :: sltL_ijkl,angintL_ijkl
1515  !character(len=500) :: msg
1516 
1517 !************************************************************************
1518 
1519  ! The lmn packed indices for (ij) and (kl).
1520  if (j_lmn>=i_lmn) then
1521    ij_lmn = i_lmn + j_lmn*(j_lmn-1)/2
1522  else
1523    ij_lmn = j_lmn + i_lmn*(i_lmn-1)/2
1524  end if
1525 
1526  if (l_lmn>=k_lmn) then
1527    kl_lmn = k_lmn + l_lmn*(l_lmn-1)/2
1528  else
1529    kl_lmn = l_lmn + k_lmn*(k_lmn-1)/2
1530  end if
1531  !
1532  ! The lm indices for (ij) and (kl) in packed storage.
1533  ij_lm = pawtab%indklmn(1,ij_lmn)
1534  ij_ln = pawtab%indklmn(2,ij_lmn)
1535 
1536  kl_lm = pawtab%indklmn(1,kl_lmn)
1537  kl_ln = pawtab%indklmn(2,kl_lmn)
1538  !
1539  ! The index of (ijkl) in the Slatrad4 database.
1540  if (kl_ln>=ij_ln) then
1541    slt_idx = ij_ln +kl_ln*(kl_ln-1)/2
1542  else
1543    slt_idx = kl_ln +ij_ln*(ij_ln-1)/2
1544  end if
1545 
1546 !BEGIN DEBUG
1547  iln = Slatrad4(slt_idx)%iln
1548  jln = Slatrad4(slt_idx)%jln
1549  kln = Slatrad4(slt_idx)%kln
1550  lln = Slatrad4(slt_idx)%lln
1551 
1552  ii = kln + lln*(lln-1)/2
1553  if (slt_idx /=  (iln + jln*(jln-1)/2 + ii*(ii-1)/2 )) then
1554    write(std_out,*)"slt_idx, iln, jln, kln, lln",slt_idx, iln, jln, kln, lln
1555    ABI_BUG("Check indices")
1556  end if
1557 !END DEBUG
1558  !
1559  ! Calculate the integral by summing over ilsum.
1560  slat_intg=zero
1561  if (Slatrad4(slt_idx)%nintgl>0) then
1562    do ilsum=Slatrad4(slt_idx)%lslat_min,Slatrad4(slt_idx)%lslat_max
1563    !% do ilsum=Slatrad4(slt_idx)%lslat_min,Slatrad4(slt_idx)%lslat_max,2
1564      isel = Slatrad4(slt_idx)%intgl_select(ilsum)
1565      if (isel/=0) then
1566        sltL_ijkl = Slatrad4(slt_idx)%intgl(isel)
1567        angintL_ijkl = summ_2gaunt(Pawang,ij_lm,kl_lm,ilsum)
1568        slat_intg = slat_intg + sltL_ijkl * angintL_ijkl
1569      end if
1570    end do
1571  end if
1572 
1573 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

SOURCE

406 subroutine slatang_cshell_free(Slatang3l)
407 
408  implicit none
409 
410 !Arguments ------------------------------------
411 !scalars
412  type(slatang_cshell_t),intent(inout) :: Slatang3l(:)
413 
414 !Local variables-------------------------------
415  integer :: ii
416 ! *********************************************************************
417 
418  !@slatang_cshell_t
419  do ii=1,SIZE(Slatang3l)
420    if (allocated(Slatang3l(ii)%sggselect)) then
421      ABI_FREE(Slatang3l(ii)%sggselect)
422    end if
423    if (allocated(Slatang3l(ii)%sgg)) then
424      ABI_FREE(Slatang3l(ii)%sgg)
425    end if
426  end do
427 
428 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)} } $

SOURCE

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

SOURCE

445 subroutine slatrad_cshell_free(Slatrad3l)
446 
447  implicit none
448 
449 !Arguments ------------------------------------
450 !scalars
451  type(slatrad_cshell_t),intent(inout) :: Slatrad3l(:)
452 
453 !Local variables-------------------------------
454  integer :: ii
455 ! *********************************************************************
456 
457  !@slatrad_cshell_t
458  do ii=1,SIZE(Slatrad3l)
459    if (allocated(Slatrad3l(ii)%rlphic_select)) then
460      ABI_FREE(Slatrad3l(ii)%rlphic_select)
461    end if
462    if (allocated(Slatrad3l(ii)%rlphic_int)) then
463      ABI_FREE(Slatrad3l(ii)%rlphic_int)
464    end if
465  end do
466 
467 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.

SOURCE

492 subroutine slatrad_cshell_init(Slatrad3l,ln2_size,Pawrad,Pawtab,Atm,Atmrad,kln_mask)
493 
494  implicit none
495 
496 !Arguments ------------------------------------
497 !scalars
498  integer,intent(in) :: ln2_size
499 !arrays
500  integer,optional,intent(in) :: kln_mask(ln2_size)
501  type(atomorb_type),intent(in) :: Atm
502  type(pawrad_type),target,intent(in) :: Atmrad,Pawrad
503  type(pawtab_type),target,intent(in) :: Pawtab
504  type(slatrad_cshell_t),intent(out) :: Slatrad3l(ln2_size)
505 
506 !Local variables ---------------------------------------
507 !scalars
508  integer :: cmesh_size,dmesh_size
509  integer :: il,iln,ilnc,isl,in,jl,jln,jn,kln,ll,lnc_size
510  integer :: lslat_max,lslat_min,lc_max,nintg
511  integer :: lmn_size,lmn2_size,do_spline,ln_size,whichdenser,isppol
512  real(dp) :: intg,intg1,ybcbeg,ybcend
513  logical :: hasameq
514 !arrays
515  integer,allocatable :: kln2ln(:,:)
516  integer, pointer :: indklmn(:,:),indlmn(:,:)
517  real(dp),allocatable :: ff1(:),ff2(:),tmp_integrals(:)
518  real(dp),ABI_CONTIGUOUS pointer :: phi_i(:),phi_j(:)
519  real(dp),allocatable,target :: phi_spl(:,:)
520  real(dp),allocatable :: der(:),ypp(:)
521  real(dp),ABI_CONTIGUOUS pointer :: crad(:),drad(:),phi_in(:)
522 
523 ! *************************************************************************
524 
525  ABI_CHECK(ln2_size==Pawtab%ij_size,"Wrong ln2_size")
526  if (PRESENT(kln_mask)) then
527    ABI_ERROR("kln_mask is present")
528  end if
529 
530  !@slatrad_cshell_t
531  lmn_size  = Pawtab%lmn_size
532  lmn2_size = Pawtab%lmn2_size
533  ln_size   = Pawtab%basis_size
534 
535  lnc_size = Atm%ln_size
536  lc_max   = Atm%l_max
537 
538  call pawrad_isame(Atmrad,Pawrad,hasameq,whichdenser)
539 
540  do_spline=0
541  if (.not.hasameq) then
542    do_spline=1
543    if (whichdenser/=1) &
544 &    ABI_COMMENT("Pawrad is denser than Atmrad!")
545  else
546    ABI_CHECK(whichdenser==1,"Pawrad is denser than Atmrad!")
547  end if
548 
549  dmesh_size = Atmrad%mesh_size
550  cmesh_size = Pawtab%mesh_size
551 
552  drad  => Atmrad%rad(1:dmesh_size)
553  crad  => Pawrad%rad(1:cmesh_size)
554 
555  ! === Spline valence basis set onto core mesh (natural spline) ===
556  if (do_spline==1) then
557    ABI_COMMENT("Splining in init_slatrad3l")
558    ABI_MALLOC(phi_spl,(dmesh_size,ln_size))
559    ABI_MALLOC(der,(cmesh_size))
560    ABI_MALLOC(ypp,(cmesh_size))
561 
562    do iln=1,ln_size
563      phi_in => Pawtab%phi(:,iln)
564      ypp(:) = zero; ybcbeg = zero; ybcend = zero
565      call spline(crad,phi_in,cmesh_size,ybcbeg,ybcend,ypp)
566      call splint(cmesh_size,crad,phi_in,ypp,dmesh_size,drad,phi_spl(:,iln))
567    end do
568 
569    ABI_FREE(der)
570    ABI_FREE(ypp)
571  end if
572 
573  indlmn  => Pawtab%indlmn(1:6,1:lmn_size)
574  indklmn => Pawtab%indklmn(1:8,1:lmn2_size)
575 
576  ABI_MALLOC(kln2ln,(6,ln2_size))
577 
578  call make_kln2ln(lmn_size,lmn2_size,ln2_size,indlmn,indklmn,kln2ln)
579 
580  ABI_MALLOC(ff1,(dmesh_size))
581  ABI_MALLOC(ff2,(dmesh_size))
582 
583  ! * Loop over the upper triangle of the [(in,il), (jn,il)] matrix.
584  ABI_CHECK(Atm%nsppol==1,"nsppol==2 not tested")
585 
586  do isppol=1,Atm%nsppol
587    do kln=1,ln2_size
588      il  = kln2ln(1,kln)
589      jl  = kln2ln(2,kln)
590      in  = kln2ln(3,kln)
591      jn  = kln2ln(4,kln)
592      iln = kln2ln(5,kln)
593      jln = kln2ln(6,kln)
594 
595      lslat_max = MAX((il+lc_max),(jl+lc_max))       - 1   ! These are indices, not l-values.
596      !lslat_min = MIN(ABS(il-lc_max),ABS(jl-lc_max)) + 1
597      lslat_min = 1 ! FIXME find better way
598 
599      Slatrad3l(kln)%lnc_size    = lnc_size
600      Slatrad3l(kln)%lslat_min   = lslat_min
601      Slatrad3l(kln)%lslat_max   = lslat_max
602 
603      Slatrad3l(kln)%nrlphic_int = 0
604 
605      ABI_MALLOC(Slatrad3l(kln)%rlphic_select,(lslat_max,lnc_size))
606      Slatrad3l(kln)%rlphic_select(:,:) = 0
607 
608      !if (PRESENT(kln_mask)) then  !FIXME THIS IS WRONG, move it below in case
609      ! if (kln_mask(kln)==0) CYCLE
610      !end if
611 
612      if (do_spline==1) then
613        ABI_COMMENT("Performing spline of valence phi")
614        phi_i => phi_spl(:,iln)
615        phi_j => phi_spl(:,jln)
616      else
617        phi_i => Pawtab%phi(:,iln)
618        phi_j => Pawtab%phi(:,jln)
619      end if
620 
621      ! * Loop over (n,l) channels for Atom orbitals
622      ABI_MALLOC(tmp_integrals,(lslat_max*lnc_size))
623      tmp_integrals(:) = zero
624      nintg=0
625 
626      do ilnc=1,lnc_size
627        ! phicore => Atm%phi(:,ilnc,isppol)
628        ff1 = phi_i * Atm%phi(:,ilnc,isppol)
629        ff2 = phi_j * Atm%phi(:,ilnc,isppol)
630        do isl=lslat_min,lslat_max ! L coming from Coulomb expansion
631          ll = isl-1
632          call calc_slatradl(ll,dmesh_size,ff2,ff1,Atmrad,intg1)
633          call calc_slatradl(ll,dmesh_size,ff1,ff2,Atmrad,intg)
634 
635          !call calc_slatradl(ll,cmesh_size,ff2,ff1,Pawrad,intg1)
636          !call calc_slatradl(ll,cmesh_size,ff1,ff2,Pawrad,intg)
637 
638          if (ABS(intg1-intg)>tol6) write(std_out,*)"DEBUG ",ll,il,in,jl,jn,intg1,intg
639 
640          ! * Store results
641          if (ABS(intg)>=tol12) then
642            nintg = nintg +1
643            Slatrad3l(kln)%rlphic_select(isl,ilnc) = nintg
644            tmp_integrals(nintg) = intg
645          end if
646        end do !ll
647      end do ! ilnc
648 
649      ! Finalize the object
650      Slatrad3l(kln)%nrlphic_int = nintg
651      ABI_MALLOC(Slatrad3l(kln)%rlphic_int,(nintg))
652      if (nintg>0) Slatrad3l(kln)%rlphic_int(:) = tmp_integrals(1:nintg)
653 
654      ABI_FREE(tmp_integrals)
655    end do !kln
656  end do !isppol
657 
658  ABI_FREE(ff1)
659  ABI_FREE(ff2)
660  ABI_FREE(kln2ln)
661 
662  if (do_spline==1)  then
663    ABI_FREE(phi_spl)
664  end if
665 
666 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

SOURCE

932 subroutine slatrad_free_0D(Slatrad)
933 
934  implicit none
935 
936 !Arguments ------------------------------------
937 !scalars
938  type(slatrad_t),intent(inout) :: Slatrad
939 
940 ! *********************************************************************
941 
942  !@slatrad_t
943  if (allocated(Slatrad%intgl_select)) then
944    ABI_FREE(Slatrad%intgl_select)
945  end if
946  if (allocated(Slatrad%intgl)) then
947    ABI_FREE(Slatrad%intgl)
948  end if
949 
950 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

SOURCE

964 subroutine slatrad_free_1D(Slatrad)
965 
966  implicit none
967 
968 !Arguments ------------------------------------
969 !scalars
970  type(slatrad_t),intent(inout) :: Slatrad(:)
971 
972 !Local variables-------------------------------
973  integer :: ii
974 ! *********************************************************************
975 
976  do ii=1,SIZE(Slatrad)
977    call slatrad_free_0D(Slatrad(ii))
978  end do
979 
980 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 exchange of the indices,
  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.

SOURCE

1024 subroutine slatrad_init(Slatrad4,which_intg,ln2_size,Pawrad,Pawtab)
1025 
1026  implicit none
1027 
1028 !Arguments ------------------------------------
1029 !scalars
1030  integer,intent(in) :: ln2_size,which_intg
1031 !arrays
1032  type(pawrad_type),target,intent(in) :: Pawrad
1033  type(pawtab_type),target,intent(in) :: Pawtab
1034  type(slatrad_t),intent(out) :: Slatrad4(ln2_size*(ln2_size+1)/2)
1035 
1036 !Local variables ---------------------------------------
1037 !scalars
1038  integer :: mesh_size,il,iln,isl,in,jl,jln,jn,sln1,sln2,l_slat
1039  integer :: kn,kl,ln,ll,kln,lln,lslat_max,lslat_min,nintgl
1040  integer :: lmn_size,lmn2_size,ln_size,slt_idx
1041  real(dp) :: ae_intg,ps_intg,pshat_intg,intg,tqij_L,tqkl_L !intg1
1042  character(len=500) :: msg
1043 !arrays
1044  integer,allocatable :: kln2ln(:,:)
1045  integer, pointer :: indklmn(:,:),indlmn(:,:)
1046  real(dp),allocatable :: uiuj(:),ukul(:),tuituj(:),tuktul(:),tuituj_tqgl(:),tuktul_tqgl(:)
1047  real(dp),allocatable :: tmp_integrals(:),ff(:)
1048  real(dp),ABI_CONTIGUOUS pointer :: phi_i(:),phi_j(:),phi_k(:),phi_l(:)
1049  real(dp),ABI_CONTIGUOUS pointer :: tphi_i(:),tphi_j(:),tphi_k(:),tphi_l(:)
1050  real(dp),ABI_CONTIGUOUS pointer :: shapefunc(:),rad(:)
1051 
1052 ! *************************************************************************
1053 
1054  DBG_ENTER("COLL")
1055 
1056  ABI_CHECK(ln2_size==Pawtab%ij_size,"Wrong ln2_size")
1057 
1058  if ( ALL(which_intg /= (/1,2,3/)) ) then
1059    write(msg,'(a,i0)')"Wrong value for which_intg: ",which_intg
1060    ABI_ERROR(msg)
1061  end if
1062 
1063  !@slatrad_t
1064  lmn_size   = Pawtab%lmn_size
1065  lmn2_size  = Pawtab%lmn2_size
1066  ln_size    = Pawtab%basis_size
1067  mesh_size  = Pawtab%mesh_size
1068  !
1069  ! Useful table for looping.
1070  indlmn  => Pawtab%indlmn(1:6,1:lmn_size)
1071  indklmn => Pawtab%indklmn(1:8,1:lmn2_size)
1072 
1073  ABI_MALLOC(kln2ln,(6,ln2_size))
1074  call make_kln2ln(lmn_size,lmn2_size,ln2_size,indlmn,indklmn,kln2ln)
1075 
1076  ABI_MALLOC(uiuj,(mesh_size))
1077  ABI_MALLOC(ukul,(mesh_size))
1078  ABI_MALLOC(ff,(mesh_size))
1079  ABI_MALLOC(tuituj,(mesh_size))
1080  ABI_MALLOC(tuktul,(mesh_size))
1081  ABI_MALLOC(tuituj_tqgl,(mesh_size))
1082  ABI_MALLOC(tuktul_tqgl,(mesh_size))
1083  rad => Pawrad%rad
1084  !
1085  ! * Loop over (k,l) channels in packed form.
1086  do sln2=1,ln2_size
1087    kl  = kln2ln(1,sln2)
1088    ll  = kln2ln(2,sln2)
1089    kn  = kln2ln(3,sln2)
1090    ln  = kln2ln(4,sln2)
1091    kln = kln2ln(5,sln2)
1092    lln = kln2ln(6,sln2)
1093    !write(std_out,*)"sln2, kln, lln",sln2,kln,lln
1094 
1095    phi_k  => Pawtab%phi (:,kln)
1096    tphi_k => Pawtab%tphi(:,kln)
1097 
1098    phi_l  => Pawtab%phi (:,lln)
1099    tphi_l => Pawtab%tphi(:,lln)
1100    !
1101    ! * Loop over (i,j) channels in packed form AND only for the upper triangle of sln2, sln1
1102    do sln1=1,sln2
1103      il  = kln2ln(1,sln1)
1104      jl  = kln2ln(2,sln1)
1105      in  = kln2ln(3,sln1)
1106      jn  = kln2ln(4,sln1)
1107      iln = kln2ln(5,sln1)
1108      jln = kln2ln(6,sln1)
1109      !write(std_out,*)"sln1, iln, jln",sln1,iln,jln
1110 
1111      slt_idx = sln1 + sln2*(sln2-1)/2 ! index for packed storage.
1112 
1113      phi_i  => Pawtab%phi (:,iln)
1114      tphi_i => Pawtab%tphi(:,iln)
1115 
1116      phi_j  => Pawtab%phi (:,jln)
1117      tphi_j => Pawtab%tphi(:,jln)
1118 
1119      lslat_min = MAX(ABS(il-jl),ABS(kl-ll)) + 1  ! We use indices not l-values.
1120      lslat_max = MIN((il+jl),(kl+ll)) - 1
1121 
1122      !lslat_min = MIN(ABS(il-jl),ABS(kl-ll)) + 1
1123      !lslat_max = MAX((il+jl),(kl+ll)) - 1
1124 
1125      Slatrad4(slt_idx)%lslat_min = lslat_min
1126      Slatrad4(slt_idx)%lslat_max = lslat_max
1127 
1128      Slatrad4(slt_idx)%iln = iln
1129      Slatrad4(slt_idx)%jln = jln
1130      Slatrad4(slt_idx)%kln = kln
1131      Slatrad4(slt_idx)%lln = lln
1132 
1133      ABI_MALLOC(Slatrad4(slt_idx)%intgl_select,(lslat_min:lslat_max))
1134      Slatrad4(slt_idx)%intgl_select(:) = 0
1135      Slatrad4(slt_idx)%nintgl          = 0
1136 
1137      if (lslat_min > lslat_max) then
1138        ! e.g. (1 2) (1 1). Due to angular selection rules, this integral do not contribue
1139        !write(std_out,*)"lslat_min, lslat_max",lslat_min,lslat_max
1140        !write(std_out,*)"il,jl,kl,ll",il,jl,kl,ll
1141        !ABI_ERROR("")
1142        ABI_MALLOC(Slatrad4(slt_idx)%intgl,(0))
1143        CYCLE
1144      end if
1145 
1146      uiuj   =  phi_i *  phi_j  ! The AE part.
1147      ukul   =  phi_k *  phi_l
1148      tuituj = tphi_i * tphi_j  ! The pseudized part.
1149      tuktul = tphi_k * tphi_l
1150      !
1151      ! Calculate L-depedent integrals where L come from the expansion the Coulomb interaction.
1152      ABI_MALLOC(tmp_integrals,(MAX(lslat_min,lslat_max)))
1153      tmp_integrals=zero
1154      nintgl=0
1155 
1156      do isl=lslat_min,lslat_max
1157      !do isl=lslat_min,lslat_max,2  ! TODO Here I can reduce the number of iterations using a step of 2.
1158        l_slat = isl-1
1159        call calc_slatradl(l_slat,mesh_size,uiuj,ukul,Pawrad,ae_intg)
1160        intg = ae_intg
1161 
1162 #if 0
1163        call calc_slatradl(l_slat,mesh_size,ukul,uiuj,Pawrad,intg1)
1164        if (ABS(intg1-ae_intg)>tol12) then
1165          write(msg,'(a,es16.8)')"s_ij and s_ij differ by ",intg1-ae_intg
1166          ABI_WARNING(msg)
1167        end if
1168 #endif
1169        if (which_intg == 2) then ! Subtract the pseudo part.
1170          call calc_slatradl(l_slat,mesh_size,tuituj,tuktul,Pawrad,ps_intg)
1171          intg = intg - ps_intg
1172 
1173        else if (which_intg == 3) then ! Subtract (pseudo + compensation charges)
1174          !
1175          ! Evaluate tqij_L and tqkl_L (without M-dependent part).
1176          ff(1)=zero
1177          ff(2:mesh_size)=(pawtab%phiphj(2:mesh_size,sln1)-pawtab%tphitphj(2:mesh_size,sln1))*rad(2:mesh_size)**l_slat
1178          if (l_slat==0.and.kl==1.and.ll==1) then
1179            call pawrad_deducer0(ff,mesh_size,pawrad)
1180          end if
1181          call simp_gen(tqij_L,ff,pawrad)
1182 
1183          ff(1)=zero
1184          ff(2:mesh_size)=(pawtab%phiphj(2:mesh_size,sln2)-pawtab%tphitphj(2:mesh_size,sln2))*rad(2:mesh_size)**l_slat
1185          if (l_slat==0.and.il==1.and.jl==1) then
1186            call pawrad_deducer0(ff,mesh_size,pawrad)
1187          end if
1188          call simp_gen(tqkl_L,ff,pawrad)
1189 
1190          shapefunc   => Pawtab%shapefunc(:,isl)  ! Recheck this part, in particular the convention
1191          tuituj_tqgl = tuituj + tqij_L * shapefunc * rad**2
1192          tuktul_tqgl = tuktul + tqkl_L * shapefunc * rad**2
1193 
1194          call calc_slatradl(l_slat,mesh_size,tuituj_tqgl,tuktul_tqgl,Pawrad,pshat_intg)
1195          intg = intg - pshat_intg
1196        end if
1197        !
1198        ! * Store results
1199        if (ABS(intg)>=tol12) then
1200          nintgl = nintgl +1
1201          Slatrad4(slt_idx)%intgl_select(isl) = nintgl
1202          tmp_integrals(nintgl) = intg
1203        end if
1204      end do !isl
1205      !
1206      ! Finalize the object.
1207      Slatrad4(slt_idx)%nintgl = nintgl
1208      ABI_MALLOC(Slatrad4(slt_idx)%intgl,(nintgl))
1209      if (nintgl>0) Slatrad4(slt_idx)%intgl(:) = tmp_integrals(1:nintgl)
1210      ABI_FREE(tmp_integrals)
1211    end do !sln1
1212  end do !sln2
1213  !
1214  ! Free memory
1215  ABI_FREE(kln2ln)
1216  ABI_FREE(uiuj)
1217  ABI_FREE(ukul)
1218  ABI_FREE(tuituj)
1219  ABI_FREE(tuktul)
1220  ABI_FREE(ff)
1221  ABI_FREE(tuituj_tqgl)
1222  ABI_FREE(tuktul_tqgl)
1223 
1224  DBG_EXIT("COLL")
1225 
1226 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 indices 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) indices 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

SOURCE

1438 function summ_2gaunt(Pawang,ij_lm,kl_lm,ll_idx)
1439 
1440  implicit none
1441 
1442 !Arguments ------------------------------------
1443 !scalars
1444  integer,intent(in) :: ij_lm,kl_lm,ll_idx
1445  real(dp) :: summ_2gaunt
1446  type(pawang_type),intent(in) :: Pawang
1447 !arrays
1448 
1449 !Local variables-------------------------------
1450 !scalars
1451  integer :: ignt1,ignt2,idx_LM,max_klm,mm,ii,ll
1452  character(len=500) :: msg
1453 
1454 !************************************************************************
1455 
1456  ! FIXME: size of gntselect depends on pawxcdev!
1457  ! Consistency check on input arguments.
1458  max_klm = pawang%l_max**2*(pawang%l_max**2+1)/2
1459  if (ij_lm>max_klm.or.kl_lm>max_klm.or.ij_lm<1.or.kl_lm<1.or.&
1460 &    ll_idx>pawang%l_size_max.or.ll_idx<1) then
1461    write(msg,'(a,3i0)')"Wrong indices, check pawxcdev ",ij_lm,kl_lm,ll_idx
1462    ABI_ERROR(msg)
1463  end if
1464 
1465  ll = ll_idx-1
1466  summ_2gaunt=zero; ii=0
1467  do mm=-ll,ll
1468    idx_LM = 1 + ll**2 + ll + mm
1469    ignt1 = Pawang%gntselect(idx_LM,ij_lm)
1470    ignt2 = Pawang%gntselect(idx_LM,kl_lm)
1471    if (ignt1>0 .and. ignt2>0) then
1472      summ_2gaunt = summ_2gaunt + Pawang%realgnt(ignt1)*Pawang%realgnt(ignt2)
1473      ii=ii+1
1474      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
1475      if (ii/=1) ABI_WARNING("ii>1")
1476    end if
1477  end do
1478 
1479 end function summ_2gaunt