TABLE OF CONTENTS


ABINIT/m_paw_denpot [ Modules ]

[ Top ] [ Modules ]

NAME

  m_paw_denpot

FUNCTION

  This module contains routines related to PAW on-site densities and on-site potentials.

COPYRIGHT

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

SOURCE

16 #if defined HAVE_CONFIG_H
17 #include "config.h"
18 #endif
19 
20 #include "abi_common.h"
21 
22 MODULE m_paw_denpot
23 
24  use defs_basis
25  use m_abicore
26  use m_errors
27  use m_xmpi
28  use m_time, only : timab
29 
30  use m_pawang,           only : pawang_type
31  use m_pawrad,           only : pawrad_type,pawrad_deducer0,poisson,simp_gen
32  use m_pawtab,           only : pawtab_type
33  use m_paw_an,           only : paw_an_type
34  use m_paw_ij,           only : paw_ij_type
35  use m_pawfgrtab,        only : pawfgrtab_type
36  use m_pawrhoij,         only : pawrhoij_type
37  use m_pawdij,           only : pawdijhartree,pawdiju_euijkl,pawdijnd,pawdijso,pawxpot,pawdijfock,symdij,symdij_all
38  use m_pawxc,            only : pawxc,pawxc_dfpt,pawxcm,pawxcm_dfpt,pawxcpositron,pawxcmpositron,pawxc_get_usekden
39  use m_paw_finegrid,     only : pawgylm
40  use m_paral_atom,       only : get_my_atmtab,free_my_atmtab
41  use m_paw_correlations, only : pawuenergy,pawxenergy,setnoccmmp
42  use m_paral_atom,       only : get_my_atmtab,free_my_atmtab
43 
44  use m_crystal,          only : crystal_t
45  use m_electronpositron, only : electronpositron_type,electronpositron_calctype
46 
47 #ifdef HAVE_FC_ISO_C_BINDING
48  use, intrinsic :: iso_c_binding, only : c_ptr,c_loc,c_f_pointer
49 #endif
50 
51  implicit none
52 
53  private
54 
55 !public procedures.
56  public :: pawdenpot           ! Compute different (PAW) energies, densities and potentials inside PAW spheres
57  public :: pawdensities        ! Compute PAW on-site densities (all-electron, pseudo and compensation)
58  public :: pawkindensities     ! Compute PAW on-site kinetic energy densities (all-electron, pseudo)
59  public :: pawaccenergy        ! Accumulate the atomic contribution of a PAW on-site energy
60  public :: pawaccenergy_nospin ! As pawaccenergy, but with no spin polarization
61  public :: paw_mknewh0         ! Compute bare PAW on-site Hamiltonian (-> GW calculations)
62 
63 CONTAINS  !========================================================================================

m_paw_denpot/paw_mknewh0 [ Functions ]

[ Top ] [ m_paw_denpot ] [ Functions ]

NAME

 paw_mknewh0

FUNCTION

 Calculates the new bare PAW Hamiltonian in the case of quasi-particle self-consistent GW calculations.

INPUTS

  mpi_atmtab(:)=--optional-- indexes of the atoms treated by current proc
  comm_atom=--optional-- MPI communicator over atoms
  my_natom=number of atoms treated by current processor
  nsppol=1 for unpolarized, 2 for spin-polarized
  nspden=number of spin-density components
  nfftf=(effective) number of FFT grid points (for this proc) for the "fine" grid
  pawspnorb=flag: 1 if spin-orbit coupling is activated
  pawprtvol=control print volume and debugging output for PAW
  Cryst<crystal_t>=Info on unit cell and its symmetries
  Pawtab(ntypat*usepaw)<type(pawtab_type)>=paw tabulated starting data
  Paw_an(natom) <type(paw_an_type)>=paw arrays given on angular mesh
  Pawang<type(pawang_type)>=paw angular mesh and related data
  Pawfgrtab(natom) <type(pawfgrtab_type)>=atomic data given on fine rectangular grid
  vxc(nfftf,nspden)=exchange-correlation potential
  vxc_val(nfftf,nspden)=valence only exchange-correlation potential
  vtrial(nfftf,nspden)=potential (Hartree+XC+loc)

SIDE EFFECTS

  Paw_ij(natom*usepaw)<Paw_ij_type)>=paw arrays given on (i,j) channels
     At output: new value for Paw_ij()%dij

SOURCE

2071 subroutine paw_mknewh0(my_natom,nsppol,nspden,nfftf,pawspnorb,pawprtvol,Cryst,&
2072 &          Pawtab,Paw_an,Paw_ij,Pawang,Pawfgrtab,vxc,vxc_val,vtrial,&
2073 &          mpi_atmtab,comm_atom) ! optional arguments (parallelism)
2074 
2075 !Arguments ------------------------------------
2076 !scalars
2077  integer,intent(in) :: my_natom,nsppol,nspden,nfftf,pawprtvol,pawspnorb
2078  integer,optional,intent(in) :: comm_atom
2079 !arrays
2080  integer,optional,target,intent(in) :: mpi_atmtab(:)
2081  real(dp),intent(in) :: vxc(nfftf,nspden),vxc_val(nfftf,nspden),vtrial(nfftf,nspden)
2082  type(crystal_t),intent(in) :: Cryst
2083  type(Pawang_type),intent(in) :: Pawang
2084  type(Pawtab_type),target,intent(in) :: Pawtab(Cryst%ntypat)
2085  type(Paw_an_type),intent(in) :: Paw_an(my_natom)
2086  type(Paw_ij_type),intent(inout) :: Paw_ij(my_natom)
2087  type(Pawfgrtab_type),intent(inout) :: Pawfgrtab(my_natom)
2088 
2089 !Local variables-------------------------------
2090 !scalars
2091  integer,parameter :: ipert0=0
2092  integer :: iat,iat_tot,idij,ndij,option_dij
2093  integer :: itypat,lmn_size,j0lmn,jlmn,ilmn,klmn,klmn1,klm
2094  integer :: lmin,lmax,mm,isel,lm_size,lmn2_size,my_comm_atom,cplex_dij
2095  integer :: ils,ilslm,ic,lm0
2096  integer :: nsploop,is2fft,qphase
2097  real(dp) :: gylm,qijl
2098  logical :: ltest,my_atmtab_allocated,paral_atom
2099  character(len=500) :: msg
2100 !arrays
2101  integer, pointer :: indklmn_(:,:)
2102  integer,pointer :: my_atmtab(:)
2103  real(dp) :: rdum(1),rdum2(1)
2104  real(dp),allocatable :: prod_hloc(:,:),prodhxc_core(:,:)
2105  real(dp),allocatable :: dijhl_hat(:,:),dijhmxc_val(:,:)
2106 
2107 ! *************************************************************************
2108 
2109  DBG_ENTER("COLL")
2110 
2111  call wrtout(std_out,'Assembling PAW strengths for the bare Hamiltonian','COLL')
2112 
2113 !== Set up parallelism over atoms ===
2114  paral_atom=(present(comm_atom).and.(my_natom/=Cryst%natom))
2115  nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab
2116  my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom
2117  call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,Cryst%natom,my_natom_ref=my_natom)
2118 
2119  if (my_natom>0) then
2120 
2121 !  === Test if required pointers in paw_ij are allocated ===
2122    ltest = (allocated(Paw_ij(1)%dijxc).and.allocated(Paw_ij(1)%dijxc_val) )
2123    ABI_CHECK(ltest,'dijxc or dijxc_val not calculated')
2124 
2125    ltest=(allocated(Paw_ij(1)%dijhat)) !.and.Paw_ij(1)%has_dijhat==2)
2126    ABI_CHECK(ltest,'dijhat not calculated')
2127 
2128    ltest=(allocated(Paw_ij(1)%dijhartree)) !.and.Paw_ij(1)%has_dijhartree==2)
2129    ABI_CHECK(ltest,'dijhartree not calculated')
2130 
2131    if (ANY(Pawtab(:)%usepawu/=0)) then
2132      do iat=1,my_natom
2133        iat_tot=iat;if (paral_atom) iat_tot=my_atmtab(iat)
2134        itypat=Cryst%typat(iat_tot)
2135        if (Pawtab(itypat)%usepawu/=0) then
2136          ltest=(allocated(Paw_ij(iat)%dijU) ) !.and.Paw_ij(iat)%has_dijU==2)
2137          write(msg,'(a,i3,a)')" For atom no. ",iat," %dijU(iat) has not been calculated."
2138          ABI_CHECK(ltest,msg)
2139        end if
2140      end do
2141    end if
2142 
2143    if (pawspnorb>0) then
2144      do iat=1,my_natom
2145        ltest=(allocated(Paw_ij(iat)%dijso) ) !.and.Paw_ij(iat)%has_dijso==2)
2146        write(msg,'(a,i3,a)')" For atom no. ",iat," %dijso(iat) has not been calculated."
2147        ABI_CHECK(ltest,msg)
2148      end do
2149    end if
2150  end if ! my_natom>0
2151 
2152 !== Construct the new PAW H0 Hamiltonian ===
2153  do iat=1,my_natom
2154    iat_tot=iat;if (paral_atom) iat_tot=my_atmtab(iat)
2155 
2156    itypat    = Cryst%typat(iat_tot)
2157    lmn_size  = Pawtab(itypat)%lmn_size
2158    lmn2_size = Pawtab(itypat)%lmn2_size
2159    lm_size   = Paw_an(iat)%lm_size
2160    cplex_dij = Paw_ij(iat)%cplex_dij
2161    qphase    = Paw_ij(iat)%qphase
2162    ndij      = Paw_ij(iat)%ndij
2163 
2164    ABI_CHECK(cplex_dij==1,'cplex_dij/=1 not implemented')
2165    ABI_CHECK(qphase==1,'qphase/=1 not implemented')
2166 !
2167 !  Eventually compute g_l(r).Y_lm(r) factors for the current atom (if not already done)
2168    if (Pawfgrtab(iat)%gylm_allocated==0) then
2169      if (allocated(Pawfgrtab(iat)%gylm))  then
2170        ABI_FREE(Pawfgrtab(iat)%gylm)
2171      end if
2172      ABI_MALLOC(Pawfgrtab(iat)%gylm,(Pawfgrtab(iat)%nfgd,lm_size))
2173      Pawfgrtab(iat)%gylm_allocated=2
2174 
2175      call pawgylm(Pawfgrtab(iat)%gylm,rdum,rdum2,lm_size,&
2176 &     Pawfgrtab(iat)%nfgd,1,0,0,Pawtab(itypat),Pawfgrtab(iat)%rfgd)
2177    end if
2178 
2179 !  === Calculate LM contribution to dijhmxc_val for this atom ===
2180 !  * Dijxc contains also the Hat term on the FFT mesh while Dijxc_val does not
2181 !  contain neither the hat term nor the LM sum of onsite terms (they should cancel each other)
2182 !  FIXME change paw_dij,  otherwise I miss tnc in vxc
2183 !  * prodhxc_core is used to assemble $\int g_l Ylm (vtrial - vxc_val[tn+nhat] dr$ on the FFT mesh ===
2184 !  * The following quantities do not depend on ij
2185    ABI_MALLOC(prod_hloc   ,(lm_size,ndij))
2186    ABI_MALLOC(prodhxc_core,(lm_size,ndij))
2187    prod_hloc   =zero
2188    prodhxc_core=zero
2189    do idij=1,ndij
2190      do ilslm=1,lm_size
2191        do ic=1,Pawfgrtab(iat)%nfgd
2192          is2fft=Pawfgrtab(iat)%ifftsph(ic)
2193          gylm=Pawfgrtab(iat)%gylm(ic,ilslm)
2194          prod_hloc (ilslm,idij)=prod_hloc (ilslm,idij) + (vtrial(is2fft,idij)-vxc(is2fft,idij))*gylm
2195 !        prodhxc_core(ilslm,idij)=prodhxc_core(ilslm,idij) + (vxc_val(is2fft,idij))*gylm
2196          prodhxc_core(ilslm,idij)=prodhxc_core(ilslm,idij) + (vtrial(is2fft,idij)-vxc_val(is2fft,idij))*gylm
2197        end do
2198      end do
2199    end do !idij
2200 
2201 !  === Assembly the "Hat" contribution for this atom ====
2202    ABI_MALLOC(dijhl_hat  ,(cplex_dij*lmn2_size,ndij))
2203    ABI_MALLOC(dijhmxc_val,(cplex_dij*lmn2_size,ndij))
2204    dijhl_hat  =zero
2205    dijhmxc_val=zero
2206    indklmn_ => Pawtab(itypat)%indklmn(1:6,1:lmn2_size)
2207 
2208    do idij=1,ndij
2209      do klmn=1,lmn2_size
2210        klm =indklmn_(1,klmn)
2211        lmin=indklmn_(3,klmn)
2212        lmax=indklmn_(4,klmn)
2213 
2214 !      === $\sum_lm q_ij^l prod* for each idij$ ===
2215        do ils=lmin,lmax,2
2216          lm0=ils**2+ils+1
2217          do mm=-ils,ils
2218            ilslm=lm0+mm
2219            isel=Pawang%gntselect(lm0+mm,klm)
2220            if (isel>0) then
2221              qijl=Pawtab(itypat)%qijl(ilslm,klmn)
2222              dijhl_hat  (klmn,idij)=dijhl_hat  (klmn,idij) +  prod_hloc (ilslm,idij)*qijl
2223              dijhmxc_val(klmn,idij)=dijhmxc_val(klmn,idij) +prodhxc_core(ilslm,idij)*qijl
2224            end if
2225          end do
2226        end do
2227      end do
2228    end do
2229 
2230    ABI_FREE(prod_hloc)
2231    ABI_FREE(prodhxc_core)
2232 
2233 !  * Normalization factor due to integration on the FFT mesh
2234    dijhl_hat  = dijhl_hat  *Cryst%ucvol/DBLE(nfftf)
2235    dijhmxc_val= dijhmxc_val*Cryst%ucvol/DBLE(nfftf)
2236 
2237 !  === Now assembly the bare Hamiltonian ===
2238 !  * Loop over density components overwriting %dij
2239    nsploop=nsppol; if (Paw_ij(iat)%ndij==4) nsploop=4
2240 
2241    do idij=1,nsploop
2242      klmn1=1
2243 
2244      do jlmn=1,lmn_size
2245        j0lmn=jlmn*(jlmn-1)/2
2246        do ilmn=1,jlmn
2247          klmn=j0lmn+ilmn
2248 
2249 !        The following gives back the input dij.
2250 !        since dijxc contains the hat term done on the FFT mesh
2251          if (.FALSE.) then
2252            Paw_ij(iat)%dij(klmn,idij) =        &
2253 &           Pawtab(itypat)%dij0    (klmn)      &
2254 &           +Paw_ij(iat)%dijhartree(klmn)      &
2255 &           +Paw_ij(iat)%dijxc     (klmn,idij) &
2256 &           +dijhl_hat   (klmn,idij)
2257 
2258          else
2259 !          === Make nonlocal part of h0 removing the valence contribution ===
2260 !          Remeber that XC contains already the Hat contribution
2261            Paw_ij(iat)%dij(klmn,idij) =        &
2262 &           Pawtab(itypat)%dij0      (klmn)    &
2263 &           +Paw_ij(iat)%dijhartree(klmn)      &
2264 &           +Paw_ij(iat)%dijxc     (klmn,idij) &  ! 2 lines to get the d1-dt1 XC core contribution + XC hat (core+val)
2265 &          -Paw_ij(iat)%dijxc_val (klmn,idij) &  ! I suppose that the "hat" term on the FFT mesh in included in both.
2266 &          +dijhmxc_val(klmn,idij)               ! Local + Hartree - XC val contribution to the "hat" term.
2267 
2268 !          Add the U contribution to the
2269 !          if (.FALSE. .and. Pawtab(itypat)%usepawu/=0) then
2270            if (.TRUE. .and. Pawtab(itypat)%usepawu/=0) then
2271              Paw_ij(iat)%dij(klmn,idij) = Paw_ij(iat)%dij(klmn,idij) + Paw_ij(iat)%dijU(klmn,idij)
2272            end if
2273          end if
2274 !        TODO dijso, dijU, vpawx?
2275 !        Just to be consistent, update some values.
2276 !$Paw_ij(iat)%dijhat(klmn,idij)=Paw_ij(iat)%dijhat(klmn,idij)-dijhmxc_val(klmn,idij)
2277 
2278        end do !ilmn
2279      end do !jlmn
2280    end do !idij
2281 
2282 !  this is to be consistent?
2283 !  deallocate(Paw_ij(iat)%dijvxc_val)
2284    ABI_FREE(dijhl_hat)
2285    ABI_FREE(dijhmxc_val)
2286  end do !iat
2287 
2288 !=== Symmetrize total Dij ===
2289  option_dij=0 ! For total Dij.
2290 #if 0
2291  if (paral_atom) then
2292    call symdij(Cryst%gprimd,Cryst%indsym,ipert0,my_natom,Cryst%natom,Cryst%nsym,Cryst%ntypat,option_dij,&
2293 &   Paw_ij,Pawang,pawprtvol,Pawtab,Cryst%rprimd,Cryst%symafm,Cryst%symrec,&
2294 &   comm_atom=my_comm_atom,mpi_atmtab=my_atmtab)
2295  else
2296    call symdij(Cryst%gprimd,,Cryst%indsym,ipert0,my_natom,Cryst%natom,Cryst%nsym,Cryst%ntypat,option_dij,&
2297 &   Paw_ij,Pawang,pawprtvol,Pawtab,Cryst%rprimd,Cryst%symafm,Cryst%symrec)
2298  end if
2299 #else
2300  if (paral_atom) then
2301    call symdij_all(Cryst%gprimd,Cryst%indsym,ipert0,my_natom,Cryst%natom,Cryst%nsym,Cryst%ntypat,&
2302 &   Paw_ij,Pawang,pawprtvol,Pawtab,Cryst%rprimd,Cryst%symafm,Cryst%symrec,&
2303 &   comm_atom=my_comm_atom,mpi_atmtab=my_atmtab)
2304  else
2305    call symdij_all(Cryst%gprimd,Cryst%indsym,ipert0,my_natom,Cryst%natom,Cryst%nsym,Cryst%ntypat,&
2306 &   Paw_ij,Pawang,pawprtvol,Pawtab,Cryst%rprimd,Cryst%symafm,Cryst%symrec)
2307  end if
2308 #endif
2309 
2310 !Destroy atom table used for parallelism
2311  call free_my_atmtab(my_atmtab,my_atmtab_allocated)
2312 
2313  DBG_EXIT("COLL")
2314 
2315 end subroutine paw_mknewh0

m_paw_denpot/pawaccenergy [ Functions ]

[ Top ] [ m_paw_denpot ] [ Functions ]

NAME

 pawaccenergy

FUNCTION

 Accumulate an on-site PAW energy by adding the contribution of the current atom.
 This contribution has the form: Sum_ij[Rhoij.Dij]

INPUTS

  pawrhoij<type(pawrhoij_type)>= datastructure containing Rho_ij values
  dij(cplex_dij*qphase_dij*lmn2_size,nspden_dij)= array containing D_ij values
  cplex_dij= 2 if dij is COMPLEX (as in the spin-orbit case), 1 if dij is REAL
  qphase_dij= 2 if dij has a exp(iqR) phase, 1 if not
  nspden_dij= number of spin components for dij
  pawtab<type(pawtab_type)>=paw tabulated starting data

OUTPUT

SIDE EFFECTS

  epaw= PAW on-site energy. At output, the contribution of the current atom
        has been added to epaw.
  [epaw_im]= imaginary part of PAW on-site energy. At output, the contribution
             of the current atom has been added to epaw.
             This imaginary p    rt only exists in a few cases (f.i. non-stationnary
             expression of 2nd-order energy)

NOTES

 * The general form for Dij is:
   D^{s1,s2}_ij = D1^{s1,s2}_ij.cos(qr) + i.D2^{s1,s2}_ij.sin(qr)
       =   [D1re^{s1,s2}_ij + i.D1im^{s1,s2}_ij).cos(qr)]
       + i.[D2re^{s1,s2}_ij + i.D2im^{s1,s2}_ij).sin(qr)]
    where
      ij are the partial waves channels
      s1,s2 are spin/spinor components
      q is the wave vector of the phase
   D1^{s1,s2}_ij.cos(qr) is stored in the the first half of paw_ij%dij and corresponds to iq=1
   D2^{s1,s2}_ij.sin(qr) is stored in the the 2nd half of paw_ij%dij and corresponds to iq=2
   D1^{s1,s2}_ij.cos(qr) and D2^{s1,s2}_ij.sin(qr) are complex if cplex_dij=2

 * The same for Rho_ij

 * The contribution to the PAW on-site energy is:
     Sum_ij_s1s2[Rho^{s2,s1}_ij * D^{s1,s2}_ij]
   Note the order of s1/s2 indices, especially for Rho_ij.
   The present implementation follows eq(15) in Hobbs et al, PRB 62, 11556(2000)
     rho^{s1,s2}^_ij = Sum[<Psi^s2|pi><pj|Psi^s1]  (s1 and s2 exponents inverted)

SOURCE

1824 subroutine pawaccenergy(epaw,pawrhoij,dij,cplex_dij,qphase_dij,nspden_dij,pawtab,epaw_im)
1825 
1826 !Arguments ---------------------------------------------
1827 !scalars
1828  integer,intent(in) :: cplex_dij,qphase_dij,nspden_dij
1829  real(dp),intent(inout) :: epaw
1830  real(dp),intent(inout),optional :: epaw_im
1831  type(pawrhoij_type),intent(in),target :: pawrhoij
1832  type(pawtab_type),intent(in) :: pawtab
1833 !arrays
1834  real(dp),intent(in) :: dij(cplex_dij*qphase_dij*pawtab%lmn2_size,nspden_dij)
1835 
1836 !Local variables ---------------------------------------
1837 !scalars
1838  integer :: cplex_rhoij,iq,iq0_dij,iq0_rhoij,irhoij,isp_dij,isp_rhoij,jrhoij
1839  integer :: klmn,kklmn,krhoij,lmn2_size,nspden_rhoij,nsploop
1840  logical :: add_imaginary
1841  real(dp) :: etmp
1842  character(len=500) :: msg
1843 !arrays
1844  real(dp),pointer :: rhoij(:,:)
1845 
1846 ! *************************************************************************
1847 
1848  DBG_ENTER("COLL")
1849 
1850 !Compatibility tests
1851  if (pawrhoij%qphase/=qphase_dij) then
1852    msg='pawaccenergy: pawrhoij%qphase/=qphase_dij!'
1853    ABI_BUG(msg)
1854  end if
1855  if (pawrhoij%nspden>nspden_dij.and.nspden_dij/=1) then
1856    msg='pawaccenergy: pawrhoij%nspden>nspden_dij!'
1857    ABI_BUG(msg)
1858  end if
1859 
1860 !Useful data
1861  nspden_rhoij=pawrhoij%nspden
1862  lmn2_size=pawtab%lmn2_size
1863 
1864 !Special treatment for nspden
1865  nsploop=nspden_rhoij
1866  if (nspden_dij==1.and.nspden_rhoij==4) nsploop=1
1867 
1868 !Non-collinear case: need a temporary rhoij
1869  if (nspden_rhoij==4.and.nspden_dij==4) then
1870    cplex_rhoij=2
1871    ABI_MALLOC(rhoij,(2*lmn2_size,4))
1872  else
1873    cplex_rhoij=pawrhoij%cplex_rhoij
1874    rhoij => pawrhoij%rhoijp
1875  end if
1876 
1877  add_imaginary=(cplex_dij==2.and.cplex_rhoij==2)
1878 
1879 !Loop over qphase components
1880  do iq=1,qphase_dij
1881    iq0_rhoij=(iq-1)*lmn2_size*cplex_rhoij
1882    iq0_dij  =(iq-1)*lmn2_size*cplex_dij
1883 
1884 !  Non-collinear case
1885    if (nspden_rhoij==4.and.nspden_dij==4) then
1886      rhoij(:,:)=zero
1887      jrhoij=(iq-1)*lmn2_size*pawrhoij%cplex_rhoij+1 ; krhoij=1
1888      do irhoij=1,pawrhoij%nrhoijsel
1889        klmn=pawrhoij%rhoijselect(irhoij)
1890        rhoij(krhoij  ,1)= half*(pawrhoij%rhoijp(jrhoij,1)+pawrhoij%rhoijp(jrhoij,4))
1891        rhoij(krhoij  ,2)= half*(pawrhoij%rhoijp(jrhoij,1)-pawrhoij%rhoijp(jrhoij,4))
1892        !Be careful we store rhoij^21 in rhoij(:,3) and rhoij^12 in rhoij(:,4)
1893        !because of the inversion of spins in rhoij definition
1894        rhoij(krhoij  ,3)= half*pawrhoij%rhoijp(jrhoij,2)
1895        rhoij(krhoij+1,3)= half*pawrhoij%rhoijp(jrhoij,3)
1896        rhoij(krhoij  ,4)= half*pawrhoij%rhoijp(jrhoij,2)
1897        rhoij(krhoij+1,4)=-half*pawrhoij%rhoijp(jrhoij,3)
1898        if (pawrhoij%cplex_rhoij==2) then
1899          rhoij(krhoij+1,1)= half*(pawrhoij%rhoijp(jrhoij+1,1)+pawrhoij%rhoijp(jrhoij+1,4))
1900          rhoij(krhoij+1,2)= half*(pawrhoij%rhoijp(jrhoij+1,1)-pawrhoij%rhoijp(jrhoij+1,4))
1901          !Be careful we store rhoij^21 in rhoij(:,3) and rhoij^12 in rhoij(:,4)
1902          !because of the inversion of spins in rhoij definition
1903          rhoij(krhoij  ,3)= rhoij(krhoij  ,3)-half*pawrhoij%rhoijp(jrhoij+1,3)
1904          rhoij(krhoij+1,3)= rhoij(krhoij+1,3)+half*pawrhoij%rhoijp(jrhoij+1,2)
1905          rhoij(krhoij  ,4)= rhoij(krhoij  ,4)+half*pawrhoij%rhoijp(jrhoij+1,3)
1906          rhoij(krhoij+1,4)= rhoij(krhoij+1,4)+half*pawrhoij%rhoijp(jrhoij+1,2)
1907        end if
1908        jrhoij=jrhoij+pawrhoij%cplex_rhoij ; krhoij=krhoij+2
1909      end do
1910      iq0_rhoij=0
1911    end if
1912 
1913 !  Contribution to on-site energy (real part)
1914    do isp_rhoij=1,nsploop
1915      isp_dij=min(isp_rhoij,nspden_dij)
1916      jrhoij=iq0_rhoij+1
1917      do irhoij=1,pawrhoij%nrhoijsel
1918        klmn=pawrhoij%rhoijselect(irhoij)
1919        kklmn=iq0_dij+cplex_dij*(klmn-1)+1
1920        etmp=rhoij(jrhoij,isp_rhoij)*dij(kklmn,isp_dij)
1921        if (add_imaginary) etmp=etmp-rhoij(jrhoij+1,isp_rhoij)*dij(kklmn+1,isp_dij)
1922        epaw=epaw+etmp*pawtab%dltij(klmn)
1923        jrhoij=jrhoij+cplex_rhoij
1924      end do
1925    end do ! nsploop
1926 
1927 !  Contribution to on-site energy (imaginary part)
1928    if (present(epaw_im).and.qphase_dij==2) then
1929      do isp_rhoij=1,nsploop
1930        isp_dij=min(isp_rhoij,nspden_dij)
1931        jrhoij=iq0_rhoij+1
1932        do irhoij=1,pawrhoij%nrhoijsel
1933          klmn=pawrhoij%rhoijselect(irhoij)
1934          if (iq==1) then
1935            kklmn=lmn2_size*cplex_dij+cplex_dij*(klmn-1)+1
1936            etmp=-rhoij(jrhoij,isp_rhoij)*dij(kklmn,isp_dij)
1937            if (add_imaginary) etmp=etmp+rhoij(jrhoij+1,isp_rhoij)*dij(kklmn+1,isp_dij)
1938          end if
1939          if (iq==2) then
1940            kklmn=cplex_dij*(klmn-1)+1
1941            etmp=rhoij(jrhoij,isp_rhoij)*dij(kklmn,isp_dij)
1942            if (add_imaginary) etmp=etmp-rhoij(jrhoij+1,isp_rhoij)*dij(kklmn+1,isp_dij)
1943          end if
1944          epaw_im=epaw_im+etmp*pawtab%dltij(klmn)
1945          jrhoij=jrhoij+cplex_rhoij
1946        end do
1947      end do ! nsploop
1948    end if
1949 
1950  end do ! qphase
1951 
1952  if (nspden_rhoij==4.and.nspden_dij==4) then
1953    ABI_FREE(rhoij)
1954  end if
1955 
1956  DBG_EXIT("COLL")
1957 
1958 end subroutine pawaccenergy

m_paw_denpot/pawaccenergy_nospin [ Functions ]

[ Top ] [ m_paw_denpot ] [ Functions ]

NAME

 pawaccenergy_nospin

FUNCTION

 Accumulate an on-site PAW energy by adding the contribution of the current atom.
 This contribution has the form: Sum_ij[Rhoij.Dij]
 Applies only for Dij without spin components (as f.i. Dij^Hartree).
 This routine is a wrapper to pawaccenergy.

INPUTS

  pawrhoij<type(pawrhoij_type)>= datastructure containing Rho_ij values
  dij(cplex_dij*qphase_dij*lmn2_size)= array containing D_ij values
  cplex_dij= 2 if dij is COMPLEX (as in the spin-orbit case), 1 if dij is REAL
  qphase_dij= 2 if dij has a exp(iqR) phase, 1 if not
  pawtab<type(pawtab_type)>=paw tabulated starting data

OUTPUT

SIDE EFFECTS

  epaw= PAW on-site energy. At output, the contribution of the current atom
        has been added to epaw.
  [epaw_im]= imaginary part of PAW on-site energy. At output, the contribution
             of the current atom has been added to epaw.
             This imaginary part only exists in a few cases (f.i. non-stationnary
             expression of 2nd-order energy)

SOURCE

1992 subroutine pawaccenergy_nospin(epaw,pawrhoij,dij,cplex_dij,qphase_dij,pawtab,epaw_im)
1993 
1994 !Arguments ---------------------------------------------
1995 !scalars
1996  integer,intent(in) :: cplex_dij,qphase_dij
1997  real(dp),intent(inout) :: epaw
1998  real(dp),intent(inout),optional :: epaw_im
1999  type(pawrhoij_type),intent(in),target :: pawrhoij
2000  type(pawtab_type),intent(in) :: pawtab
2001 !arrays
2002  real(dp),intent(in),target :: dij(cplex_dij*qphase_dij*pawtab%lmn2_size)
2003 
2004 !Local variables ---------------------------------------
2005 !scalars
2006  integer :: size_dij
2007 #ifdef HAVE_FC_ISO_C_BINDING
2008  type(C_PTR) :: cptr
2009 #endif
2010 !arrays
2011 real(dp), ABI_CONTIGUOUS pointer :: dij_2D(:,:)
2012 
2013 ! *************************************************************************
2014 
2015  size_dij=size(dij)
2016 
2017 #ifdef HAVE_FC_ISO_C_BINDING
2018  cptr=c_loc(dij(1))
2019  call c_f_pointer(cptr,dij_2D,shape=[size_dij,1])
2020 #else
2021  ABI_MALLOC(dij_2D,(size_dij,1))
2022  dij_2D=reshape(dij,[size_dij,1])
2023 #endif
2024 
2025  if (present(epaw_im)) then
2026    call pawaccenergy(epaw,pawrhoij,dij_2D,cplex_dij,qphase_dij,1,pawtab,epaw_im=epaw_im)
2027  else
2028    call pawaccenergy(epaw,pawrhoij,dij_2D,cplex_dij,qphase_dij,1,pawtab)
2029  end if
2030 
2031 #ifndef HAVE_FC_ISO_C_BINDING
2032  ABI_FREE(dij_2D)
2033 #endif
2034 
2035 end subroutine pawaccenergy_nospin

m_paw_denpot/pawdenpot [ Functions ]

[ Top ] [ m_paw_denpot ] [ Functions ]

NAME

 pawdenpot

FUNCTION

 Compute different (PAW) energies, densities and potentials (or potential-like quantities)
 inside PAW spheres
 Can also compute first-order densities potentials and second-order energies (RF calculations).

INPUTS

  electronpositron <type(electronpositron_type)>=quantities for the electron-positron annihilation (optional argument)
  [hyb_mixing, hyb_mixing_sr]= -- optional-- mixing factors for the global (resp. screened) XC hybrid functional
  ipert=index of perturbation (used only for RF calculation ; set ipert<=0 for GS calculations.
  ixc= choice of exchange-correlation scheme (see above, and below)
  mpi_atmtab(:)=--optional-- indexes of the atoms treated by current proc
  comm_atom=--optional-- MPI communicator over atoms
  my_natom=number of atoms treated by current processor
  natom=total number of atoms in cell
  nspden=number of spin-density components
  ntypat=number of types of atoms in unit cell.
  nucdipmom(3,natom) nuclear dipole moments
  nzlmopt= if -1, compute all LM-moments of densities
                  initialize "lmselect" (index of non-zero LM-moments of densities)
           if  0, compute all LM-moments of densities
                  force "lmselect" to .true. (index of non-zero LM-moments of densities)
           if  1, compute only non-zero LM-moments of densities (stored before)
  option=0: compute both energies and potentials
         1: compute only potentials
         2: compute only energies
  paw_an(my_natom) <type(paw_an_type)>=paw arrays given on angular mesh
  paw_an0(my_natom) <type(paw_an_type)>=paw arrays given on angular mesh for Ground-State
  paw_ij(my_natom) <type(paw_ij_type)>=paw arrays given on (i,j) channels
  pawang <type(pawang_type)>=paw angular mesh and related data
  pawprtvol=control print volume and debugging output for PAW
  pawrad(ntypat) <type(pawrad_type)>=paw radial mesh and related data
  pawrhoij(my_natom) <type(pawrhoij_type)>= paw rhoij occupancies and related data
  pawspnorb=flag: 1 if spin-orbit coupling is activated
  pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data
  pawxcdev=Choice of XC development (0=no dev. (use of angular mesh) ; 1 or 2=dev. on moments)
  ucvol=unit cell volume (bohr^3)
  xclevel= XC functional level
  xc_denpos= lowest allowe density (usually for the computation of the XC functionals)
  znucl(ntypat)=gives the nuclear charge for all types of atoms

OUTPUT

  paw_ij(my_natom)%dijhartree(qphase*lmn2_size)=Hartree contribution to dij;
                                      Enters into calculation of hartree energy
  ==== if option=0 or 2
    epaw=contribution to total energy from the PAW "on-site" part
    epawdc=contribution to total double-counting energy from the PAW "on-site" part
  ==== if option=0 or 2 and ipert<=0
    compch_sph=compensation charge integral inside spheres computed over spherical meshes
  ==== if (option=0 or 1) and paw_an(:)%has_vxc=1
    paw_an(my_natom)%vxc1(cplex*mesh_size,:,nspden)=XC potential calculated from "on-site" density
    paw_an(my_natom)%vxct1(cplex*mesh_size,:,nspden)=XC potential calculated from "on-site" pseudo density
  ==== if (option=0 or 1) and paw_an(:)%has_vxctau=1
    paw_an(my_natom)%vxctau1(cplex*mesh_size,:,nspden)=1st deriv. of XC energy wrt to kinetic energy density (all electron)
    paw_an(my_natom)%vxcttau1(cplex*mesh_size,:,nspden)=1st deriv. of XC energy wrt to kinetic energy density (pseudo)
  ==== if paw_an(iatom_tot)%has_vxcval==1 compute also XC potentials neglecting core charge
      paw_an(my_natom)%vxc1_val(cplex*mesh_size,:nspden)=XC potential calculated from spherical valence density
      paw_an(my_natom)%vxct1_val(cplex*mesh_size,:nspden)=XC potential calculated from spherical valence pseudo density
  ==== if nzlmopt==-1,
    paw_an(iatom_tot)%lnmselect(lm_size,nspden)=select the non-zero LM-moments of rho1 and trho1
  ==== if paw_an(:)%has_vhartree=1
    paw_an(my_natom)%vh1(cplex*mesh_size,1,1)=Hartree total potential calculated from "on-site" density
  ==== if pawspnorb>0
    paw_ij(my_natom)%dijso(qphase*cplex_dij*lmn2_size,nspden)=spin-orbit contribution to dij

NOTES

  Response function calculations:
    In order to compute first- or second-order quantities, paw_an (resp. paw_ij) datastructures
    must contain first-order quantities, namely paw_an1 (resp. paw_ij1).

SOURCE

 144 subroutine pawdenpot(compch_sph,epaw,epawdc,ipert,ixc,&
 145 & my_natom,natom,nspden,ntypat,nucdipmom,nzlmopt,option,paw_an,paw_an0,&
 146 & paw_ij,pawang,pawprtvol,pawrad,pawrhoij,pawspnorb,pawtab,pawxcdev,spnorbscl,xclevel,xc_denpos,ucvol,znucl,&
 147 & electronpositron,mpi_atmtab,comm_atom,vpotzero,hyb_mixing,hyb_mixing_sr) ! optional arguments
 148 
 149 !Arguments ---------------------------------------------
 150 !scalars
 151  integer,intent(in) :: ipert,ixc,my_natom,natom,nspden,ntypat,nzlmopt,option,pawprtvol
 152  integer,intent(in) :: pawspnorb,pawxcdev,xclevel
 153  integer,optional,intent(in) :: comm_atom
 154  real(dp), intent(in) :: spnorbscl,xc_denpos,ucvol
 155  real(dp),intent(in),optional :: hyb_mixing,hyb_mixing_sr
 156  real(dp),intent(out) :: compch_sph,epaw,epawdc
 157  type(electronpositron_type),pointer,optional :: electronpositron
 158  type(pawang_type),intent(in) :: pawang
 159 !arrays
 160  integer,optional,target,intent(in) :: mpi_atmtab(:)
 161  real(dp),intent(in) :: nucdipmom(3,natom),znucl(ntypat)
 162  real(dp),intent(out),optional :: vpotzero(2)
 163  type(paw_an_type),intent(inout) :: paw_an(my_natom)
 164  type(paw_an_type), intent(in) :: paw_an0(my_natom)
 165  type(paw_ij_type),intent(inout) :: paw_ij(my_natom)
 166  type(pawrad_type),intent(in) :: pawrad(ntypat)
 167  type(pawrhoij_type),intent(inout) :: pawrhoij(my_natom)
 168  type(pawtab_type),intent(in) :: pawtab(ntypat)
 169 
 170 !Local variables ---------------------------------------
 171 !scalars
 172  integer, parameter :: PAWU_ALGO_1=1,PAWU_ALGO_2=2
 173  integer, parameter :: PAWU_FLL=1,PAWU_AMF=2
 174  integer :: cplex,cplex_dij,cplex_rhoij,has_kxc,has_k3xc,has_vxctau
 175  integer :: iatom,iatom_tot,idum,ierr,ii,ipositron,iq,iq0_dij,iq0_rhoij
 176  integer :: itypat,itypat0,lm_size,lmn2_size,mesh_size
 177  integer :: my_comm_atom,ndij,nkxc1,nk3xc1,nsppol,opt_compch,pawu_algo,pawu_dblec
 178  integer :: qphase,usecore,usekden,usetcore,usepawu,usexcnhat,usenhat,usefock
 179  logical :: keep_vhartree,my_atmtab_allocated,need_kxc,need_k3xc,need_vxctau
 180  logical :: non_magnetic_xc,paral_atom,temp_vxc
 181  real(dp) :: e1t10,e1xc,e1xcdc,efock,efockdc,eexc,eexcdc,eexdctemp
 182  real(dp) :: eexc_val,eexcdc_val,eexex,eexexdc,eextemp,eh2
 183  real(dp) :: edftumdc,edftumdcdc,edftufll,enucdip,etmp,espnorb,etild1xc,etild1xcdc
 184  real(dp) :: exccore,exchmix,hyb_mixing_,hyb_mixing_sr_,rdum
 185  character(len=3) :: pertstrg
 186  character(len=500) :: msg
 187 !arrays
 188  integer :: idum1(0),idum3(0,0,0)
 189  integer,pointer :: my_atmtab(:)
 190  logical,allocatable :: lmselect_cur(:),lmselect_cur_ep(:),lmselect_ep(:),lmselect_tmp(:)
 191  real(dp) :: mpiarr(7),tsec(2)
 192  real(dp),allocatable :: dij_ep(:),dijfock_vv(:,:),dijfock_cv(:,:)
 193  real(dp),allocatable :: one_over_rad2(:),kxc_tmp(:,:,:),k3xc_tmp(:,:,:)
 194  real(dp),allocatable :: nhat1(:,:,:),nhat1_ep(:,:,:)
 195  real(dp) :: rdum2(0,0),rdum3(0,0,0),rdum3a(0,0,0),rdum4(0,0,0,0)
 196  real(dp),allocatable :: rho(:),rho1(:,:,:),rho1_ep(:,:,:),rho1xx(:,:,:)
 197  real(dp),allocatable :: tau1(:,:,:),ttau1(:,:,:), trho1(:,:,:),trho1_ep(:,:,:)
 198  real(dp),allocatable :: vh(:),vxc_tmp(:,:,:),vxctau_tmp(:,:,:)
 199 
 200 ! *************************************************************************
 201 
 202  DBG_ENTER("COLL")
 203 
 204  call timab(560,1,tsec)
 205 
 206 !Various inits
 207  hyb_mixing_   =zero ; if(present(hyb_mixing))    hyb_mixing_   =hyb_mixing
 208  hyb_mixing_sr_=zero ; if(present(hyb_mixing_sr)) hyb_mixing_sr_=hyb_mixing_sr
 209  usefock=0;if (abs(hyb_mixing_)>tol8.or.abs(hyb_mixing_sr_)>tol8) usefock=1
 210  usexcnhat=maxval(pawtab(1:ntypat)%usexcnhat)
 211  usekden=pawxc_get_usekden(ixc)
 212  usenhat = usexcnhat
 213  keep_vhartree=(maxval(paw_an(:)%has_vhartree)>0)
 214  if (keep_vhartree) usenhat = 1
 215  compch_sph=-1.d5
 216  opt_compch=0;if (option/=1.and.ipert<=0) opt_compch=1
 217  if (opt_compch==1) compch_sph=zero
 218  nsppol=1;if (my_natom>0) nsppol=pawrhoij(1)%nsppol
 219  pertstrg=" ";if (ipert>0) pertstrg="(1)"
 220 
 221 !Various checks
 222  if(nzlmopt/=0.and.nzlmopt/=1.and.nzlmopt/=-1) then
 223    msg='invalid value for variable "nzlmopt"!'
 224    ABI_BUG(msg)
 225  end if
 226  if (my_natom>0) then
 227    if(paw_ij(1)%has_dijhartree==0.and.ipert/=natom+1.and.ipert/=natom+10) then
 228      msg='dijhartree must be allocated!'
 229      ABI_BUG(msg)
 230    end if
 231    if(paw_ij(1)%has_dijU==0.and.pawtab(1)%usepawu/=0.and. &
 232 &    ((ipert>0.and.ipert/=natom+1.and.ipert/=natom+10).or.pawtab(1)%usepawu<0)) then
 233      msg='dijU must be allocated!'
 234      ABI_BUG(msg)
 235    end if
 236    if (pawrhoij(1)%qphase<paw_an(1)%cplex) then
 237      msg='pawrhoij()%qphase must be >=paw_an()%cplex!'
 238      ABI_BUG(msg)
 239    end if
 240    if (ipert>0.and.(ipert<=natom.or.ipert==natom+2).and.paw_an0(1)%has_kxc/=2) then
 241      msg='XC kernels for ground state must be in memory!'
 242      ABI_BUG(msg)
 243    end if
 244    if(paw_an(1)%has_vxc==0.and.(option==0.or.option==1).and. &
 245 &   .not.(ipert==natom+1.or.ipert==natom+10)) then
 246      msg='vxc1 and vxct1 must be allocated!'
 247      ABI_BUG(msg)
 248    end if
 249    if(paw_an(1)%has_vxctau==0.and.(option==0.or.option==1).and.usekden==1) then
 250      msg='vxctau1 and vxcttau1 must be allocated!'
 251      ABI_BUG(msg)
 252    end if
 253    if (ipert>0.and.paw_an(1)%has_vxctau==1.and.usekden==1) then
 254      msg='computation of vxctau not compatible with RF (ipert>0)!'
 255      ABI_BUG(msg)
 256    end if
 257    if (ipert>0.and.paw_an(1)%has_vhartree==1) then
 258      msg='computation of vhartree not compatible with RF (ipert>0)!'
 259      ABI_BUG(msg)
 260    end if
 261    if (ipert>0.and.paw_an(1)%has_vxcval==1.and.(option==0.or.option==1)) then
 262      msg='computation of vxc_val not compatible with RF (ipert>0)!'
 263      ABI_BUG(msg)
 264    end if
 265  end if
 266 
 267  ipositron=0
 268  if (present(electronpositron)) then
 269    ipositron=electronpositron_calctype(electronpositron)
 270    if (ipositron==1.and.pawtab(1)%has_kij/=2) then
 271      msg='kij must be in memory for electronpositron%calctype=1!'
 272      ABI_BUG(msg)
 273    end if
 274    if (ipert>0) then
 275      msg='electron-positron calculation not available for ipert>0!'
 276      ABI_ERROR(msg)
 277    end if
 278  end if
 279 
 280 !Set up parallelism over atoms
 281  paral_atom=(present(comm_atom).and.(my_natom/=natom))
 282  nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab
 283  my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom
 284  call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,&
 285 & my_natom_ref=my_natom)
 286 
 287 !For some perturbations, nothing to do
 288  if (ipert==natom+1.or.ipert==natom+10) then
 289    if (option/=1) then
 290      epaw=zero;epawdc=zero
 291    end if
 292    return
 293  end if
 294 
 295 !Init energies
 296  if (option/=1) then
 297    e1xc=zero     ; e1xcdc=zero
 298    etild1xc=zero ; etild1xcdc=zero
 299    exccore=zero  ; eh2=zero ; e1t10=zero
 300    edftumdc=zero ; edftumdcdc=zero ; edftufll=zero
 301    eexex=zero    ; eexexdc=zero
 302    eextemp=zero  ; eexdctemp=zero
 303    espnorb=zero  ; enucdip=zero
 304    efock=zero    ; efockdc=zero
 305    if (ipositron/=0) then
 306      electronpositron%e_paw  =zero
 307      electronpositron%e_pawdc=zero
 308    end if
 309  end if
 310 !vpotzero is needed for both the energy and the potential
 311  if (present(vpotzero)) vpotzero(:)=zero
 312 
 313 !Select PAW+U algo, different for DFT and DFPT
 314  usepawu=maxval(pawtab(1:ntypat)%usepawu)
 315  ii=minval(pawtab(1:ntypat)%usepawu);if (ii<0) usepawu=ii
 316  non_magnetic_xc=(mod(abs(usepawu),10)==4)
 317 
 318 !if PAW+U, compute noccmmp^{\sigma}_{m,m'} occupation matrix
 319  if (usepawu/=0.and.ipert<=0.and.ipositron/=1) then
 320    if (paral_atom) then
 321      call setnoccmmp(1,0,rdum4,0,0,idum3,my_natom,natom,0,1,nsppol,0,ntypat,&
 322 &     paw_ij,pawang,pawprtvol,pawrhoij,pawtab,rdum2,idum1,idum1,0,usepawu,&
 323 &     comm_atom=my_comm_atom,mpi_atmtab=mpi_atmtab)
 324    else
 325      call setnoccmmp(1,0,rdum4,0,0,idum3,my_natom,natom,0,1,nsppol,0,ntypat,&
 326 &     paw_ij,pawang,pawprtvol,pawrhoij,pawtab,rdum2,idum1,idum1,0,usepawu)
 327    end if
 328  end if
 329 
 330 !Print some titles
 331  if (abs(pawprtvol)>=2) then
 332    if (nzlmopt<1) write(msg, '(6a)') ch10,' PAW TEST:',ch10,&
 333    ' ====== Moments of (n1-tn1)',trim(pertstrg),' ========='
 334    if (nzlmopt==1) write(msg, '(6a)') ch10,' PAW TEST:',ch10,&
 335    ' ==== Non-zero Moments of (n1-tn1)',trim(pertstrg),' ===='
 336    call wrtout(std_out,msg,'COLL')
 337    if (usexcnhat/=0) then
 338      write(msg, '(6a)')' The moments of (n1-tn1-nhat1)',trim(pertstrg),' must be very small...'
 339      call wrtout(std_out,msg,'COLL')
 340    end if
 341  end if
 342 
 343 !================ Big loop on atoms =======================
 344 !==========================================================
 345 
 346  do iatom=1,my_natom
 347    iatom_tot=iatom;if (paral_atom) iatom_tot=my_atmtab(iatom)
 348    itypat=pawrhoij(iatom)%itypat
 349    exchmix=pawtab(itypat)%exchmix
 350    lmn2_size=paw_ij(iatom)%lmn2_size
 351    lm_size=paw_an(iatom)%lm_size
 352    mesh_size=pawtab(itypat)%mesh_size
 353    usecore=1;usetcore =pawtab(itypat)%usetcore
 354    if (ipert/=0) usecore=0  ! This is true for phonons and Efield pert.
 355    if (ipert/=0) usetcore=0 ! This is true for phonons and Efield pert.
 356    has_kxc =paw_an(iatom)%has_kxc ;need_kxc =(has_kxc ==1)
 357    has_k3xc=paw_an(iatom)%has_k3xc;need_k3xc=(has_k3xc==1)
 358    has_vxctau=paw_an(iatom)%has_vxctau ;need_vxctau =(has_vxctau>=1.and.usekden==1)
 359    cplex=paw_an(iatom)%cplex
 360    cplex_dij=paw_ij(iatom)%cplex_dij
 361    cplex_rhoij=pawrhoij(iatom)%cplex_rhoij
 362    qphase=pawrhoij(iatom)%qphase
 363    ndij=paw_ij(iatom)%ndij
 364    iq0_rhoij=cplex_rhoij*lmn2_size
 365    iq0_dij=cplex_dij*lmn2_size
 366    usepawu=pawtab(itypat)%usepawu
 367    pawu_algo=merge(PAWU_ALGO_1,PAWU_ALGO_2,ipert<=0.and.usepawu>=0)
 368    pawu_dblec=merge(PAWU_FLL,PAWU_AMF,abs(usepawu)==1.or.abs(usepawu)==4)
 369 
 370 !  Allocations of "on-site" densities
 371    ABI_MALLOC(rho1 ,(cplex*mesh_size,lm_size,nspden))
 372    ABI_MALLOC(trho1,(cplex*mesh_size,lm_size,nspden))
 373    ABI_MALLOC(nhat1,(cplex*mesh_size,lm_size,nspden*usenhat))
 374    rho1(:,:,:)=zero;trho1(:,:,:)=zero;nhat1(:,:,:)=zero
 375    if (usekden==1) then
 376      ABI_MALLOC(tau1 ,(cplex*mesh_size,lm_size,nspden))
 377      ABI_MALLOC(ttau1,(cplex*mesh_size,lm_size,nspden))
 378    end if
 379    if (ipositron/=0) then ! Additional allocation for the electron-positron case
 380      ABI_MALLOC(rho1_ep ,(cplex*mesh_size,lm_size,nspden))
 381      ABI_MALLOC(trho1_ep,(cplex*mesh_size,lm_size,nspden))
 382      ABI_MALLOC(nhat1_ep,(cplex*mesh_size,lm_size,nspden*usenhat))
 383    end if
 384    ABI_MALLOC(lmselect_cur,(lm_size))
 385    lmselect_cur(:)=.true.
 386    if (nzlmopt==1) lmselect_cur(:)=paw_an(iatom)%lmselect(:)
 387 
 388 !  Store some usefull quantities
 389    itypat0=0;if (iatom>1) itypat0=pawrhoij(iatom-1)%itypat
 390    if (itypat/=itypat0) then
 391      ABI_MALLOC(one_over_rad2,(mesh_size))
 392      one_over_rad2(1)=zero
 393      one_over_rad2(2:mesh_size)=one/pawrad(itypat)%rad(2:mesh_size)**2
 394    end if
 395 
 396 !  Need to allocate vxc1 in particular cases
 397    if (pawspnorb>0.and.ipert==0.and.option==2.and.ipositron/=1.and. &
 398 &      cplex_rhoij==2.and.paw_an(iatom)%has_vxc==0) then
 399 !    These should already be allocated in paw_an_init!
 400      if (allocated(paw_an(iatom)%vxc1))  then
 401        ABI_FREE(paw_an(iatom)%vxc1)
 402      end if
 403      if (pawxcdev==0)then
 404        ABI_MALLOC(paw_an(iatom)%vxc1,(cplex*mesh_size,paw_an(iatom)%angl_size,nspden))
 405      else
 406        ABI_MALLOC(paw_an(iatom)%vxc1,(cplex*mesh_size,lm_size,nspden))
 407      end if
 408      paw_an(iatom)%has_vxc=1
 409      temp_vxc=.true.
 410    else
 411      temp_vxc=.false.
 412    end if
 413 
 414 !  ===== Compute "on-site" densities (n1, ntild1, nhat1) =====
 415 !  ==========================================================
 416 
 417    call pawdensities(compch_sph,cplex,iatom_tot,lmselect_cur,paw_an(iatom)%lmselect,lm_size,&
 418 &   nhat1,nspden,nzlmopt,opt_compch,1-usenhat,-1,1,pawang,pawprtvol,pawrad(itypat),&
 419 &   pawrhoij(iatom),pawtab(itypat),rho1,trho1,one_over_rad2=one_over_rad2)
 420 
 421    if (usekden==1) then
 422      ABI_MALLOC(lmselect_tmp,(lm_size))
 423      lmselect_tmp(:)=.true.
 424      call pawkindensities(cplex,lmselect_tmp,lm_size,nspden,-1,1,-1,&
 425 &     pawang,pawrad(itypat),pawrhoij(iatom),pawtab(itypat),tau1,ttau1,&
 426 &     one_over_rad2=one_over_rad2)
 427      ABI_FREE(lmselect_tmp)
 428    end if
 429 
 430    if (ipositron/=0) then
 431 !    Electron-positron calculation: need additional on-site densities:
 432 !    if ipositron==1, need electronic on-site densities
 433 !    if ipositron==2, need positronic on-site densities
 434      ABI_MALLOC(lmselect_ep,(lm_size))
 435      ABI_MALLOC(lmselect_cur_ep,(lm_size))
 436      lmselect_cur_ep(:)=.true.
 437      if (nzlmopt==1) lmselect_cur_ep(:)=electronpositron%lmselect_ep(1:lm_size,iatom)
 438      call pawdensities(rdum,cplex,iatom_tot,lmselect_cur_ep,lmselect_ep,&
 439 &     lm_size,nhat1_ep,nspden,nzlmopt,0,1-usenhat,-1,0,pawang,0,pawrad(itypat),&
 440 &     electronpositron%pawrhoij_ep(iatom),pawtab(itypat),&
 441 &     rho1_ep,trho1_ep,one_over_rad2=one_over_rad2)
 442      if (nzlmopt<1) electronpositron%lmselect_ep(1:lm_size,iatom)=lmselect_ep(1:lm_size)
 443      ABI_FREE(lmselect_ep)
 444      ABI_FREE(lmselect_cur_ep)
 445    end if
 446 
 447 !  =========== Compute XC potentials and energies ===========
 448 !  ==========================================================
 449 
 450 !  Temporary storage
 451    nkxc1 =0;if (paw_an(iatom)%has_kxc /=0) nkxc1 =paw_an(iatom)%nkxc1
 452    nk3xc1=0;if (paw_an(iatom)%has_k3xc/=0.and.pawxcdev==0) nk3xc1=paw_an(iatom)%nk3xc1
 453    if (pawxcdev/=0) then
 454      ABI_MALLOC(vxc_tmp,(cplex*mesh_size,lm_size,nspden))
 455      if (need_kxc) then
 456        ABI_MALLOC(kxc_tmp,(mesh_size,lm_size,nkxc1))
 457      end if
 458      if (need_k3xc) then
 459        msg = 'Computation of k3xc with pawxcdev/=0 is not implemented yet!'
 460        ABI_BUG(msg)
 461      end if
 462    end if
 463    if (pawxcdev==0) then
 464      ABI_MALLOC(vxc_tmp,(cplex*mesh_size,pawang%angl_size,nspden))
 465      vxc_tmp(:,:,:)=zero
 466      if (need_kxc) then
 467        ABI_MALLOC(kxc_tmp,(mesh_size,pawang%angl_size,nkxc1))
 468      end if
 469      if (need_k3xc) then
 470        ABI_MALLOC(k3xc_tmp,(mesh_size,pawang%angl_size,nk3xc1))
 471      end if
 472      if (need_vxctau) then
 473        ABI_MALLOC(vxctau_tmp,(cplex*mesh_size,pawang%angl_size,nspden))
 474        vxctau_tmp(:,:,:)=zero
 475      end if
 476    end if
 477    idum=0
 478    if (.not.allocated(vxc_tmp))  then
 479      ABI_MALLOC(vxc_tmp,(0,0,0))
 480    end if
 481    if (.not.allocated(kxc_tmp))  then
 482      ABI_MALLOC(kxc_tmp,(0,0,0))
 483    end if
 484    if (.not.allocated(k3xc_tmp))  then
 485      ABI_MALLOC(k3xc_tmp,(0,0,0))
 486    end if
 487    if (.not.allocated(vxctau_tmp))  then
 488      ABI_MALLOC(vxctau_tmp,(0,0,0))
 489    end if
 490 
 491 !  ===== Vxc1 term =====
 492    if (ipositron/=1) then
 493      if (pawxcdev/=0) then
 494        if (ipert==0) then
 495          call pawxcm(pawtab(itypat)%coredens,eexc,eexcdc,idum,hyb_mixing_,ixc,kxc_tmp,lm_size,&
 496 &         paw_an(iatom)%lmselect,nhat1,nkxc1,non_magnetic_xc,mesh_size,nspden,option,&
 497 &         pawang,pawrad(itypat),pawxcdev,rho1,usecore,0,vxc_tmp,xclevel,xc_denpos)
 498        else
 499          call pawxcm_dfpt(pawtab(itypat)%coredens,cplex,cplex,eexc,ixc,paw_an0(iatom)%kxc1,lm_size,&
 500 &         paw_an(iatom)%lmselect,nhat1,paw_an0(iatom)%nkxc1,non_magnetic_xc,mesh_size,nspden,option,&
 501 &         pawang,pawrad(itypat),rho1,usecore,0,vxc_tmp,xclevel)
 502          eexcdc=zero
 503        end if
 504      else
 505        if (ipert==0) then
 506          call pawxc(pawtab(itypat)%coredens,eexc,eexcdc,hyb_mixing_,ixc,kxc_tmp,k3xc_tmp,lm_size,&
 507 &         paw_an(iatom)%lmselect,nhat1,nkxc1,nk3xc1,non_magnetic_xc,mesh_size,nspden,option,&
 508 &         pawang,pawrad(itypat),rho1,usecore,0,vxc_tmp,xclevel,xc_denpos,&
 509 &         coretau=pawtab(itypat)%coretau,taur=tau1,vxctau=vxctau_tmp)
 510        else
 511          call pawxc_dfpt(pawtab(itypat)%coredens,cplex,cplex,eexc,ixc,paw_an0(iatom)%kxc1,lm_size,&
 512 &         paw_an(iatom)%lmselect,nhat1,paw_an0(iatom)%nkxc1,non_magnetic_xc,mesh_size,nspden,option,&
 513 &         pawang,pawrad(itypat),rho1,usecore,0,paw_an0(iatom)%vxc1,vxc_tmp,xclevel)
 514          eexcdc=zero
 515        end if
 516      end if
 517      if (option/=1) then
 518        e1xc=e1xc+eexc
 519        e1xcdc=e1xcdc+eexcdc
 520      end if
 521      if (option<2.or.temp_vxc) paw_an(iatom)%vxc1(:,:,:)=vxc_tmp(:,:,:)
 522      if (option<2.and.need_vxctau) paw_an(iatom)%vxctau1(:,:,:)=vxctau_tmp(:,:,:)
 523      if (need_kxc .and.nkxc1>0 ) paw_an(iatom)%kxc1(:,:,:) =kxc_tmp(:,:,:)
 524      if (need_k3xc.and.nk3xc1>0) paw_an(iatom)%k3xc1(:,:,:)=k3xc_tmp(:,:,:)
 525    else ! ipositron==1
 526      if (option<2.or.temp_vxc) paw_an(iatom)%vxc1(:,:,:)=zero
 527      if (need_kxc.and.nkxc1>0) paw_an(iatom)%kxc1(:,:,:)=zero
 528    end if
 529 
 530 
 531 !  Additional electron-positron XC term (if ipositron/=0)
 532    if (ipositron/=0) then
 533      if (pawxcdev/=0) then
 534        call pawxcmpositron(ipositron,pawtab(itypat)%coredens,eexc,eexcdc,electronpositron%ixcpositron,&
 535 &       lm_size,paw_an(iatom)%lmselect,electronpositron%lmselect_ep(1:lm_size,iatom),&
 536 &       nhat1,nhat1_ep,mesh_size,nspden,option,pawang,pawrad(itypat),pawxcdev,&
 537 &       electronpositron%posdensity0_limit,rho1,rho1_ep,usecore,0,vxc_tmp,xc_denpos)
 538      else
 539        call pawxcpositron(ipositron,pawtab(itypat)%coredens,eexc,eexcdc,electronpositron%ixcpositron,&
 540 &       lm_size,paw_an(iatom)%lmselect,electronpositron%lmselect_ep(1:lm_size,iatom),&
 541 &       nhat1,nhat1_ep,mesh_size,nspden,option,pawang,pawrad(itypat),&
 542 &       electronpositron%posdensity0_limit,rho1,rho1_ep,usecore,0,vxc_tmp,xc_denpos)
 543      end if
 544      if (option/=1) then
 545        electronpositron%e_paw  =electronpositron%e_paw  +eexc
 546        electronpositron%e_pawdc=electronpositron%e_pawdc+eexcdc
 547      end if
 548      if (option<2.or.temp_vxc) paw_an(iatom)%vxc1(:,:,:)=paw_an(iatom)%vxc1(:,:,:)+vxc_tmp(:,:,:)
 549      if (need_kxc.and.nkxc1>0) paw_an(iatom)%kxc1(:,:,:)=paw_an(iatom)%kxc1(:,:,:)+kxc_tmp(:,:,:)
 550    end if
 551 
 552 !  ===== tVxc1 term =====
 553    if (ipositron/=1) then
 554      if (pawxcdev/=0) then
 555        if (ipert==0) then
 556          call pawxcm(pawtab(itypat)%tcoredens(:,1),&
 557 &         eexc,eexcdc,idum,hyb_mixing_,ixc,kxc_tmp,lm_size,&
 558 &         paw_an(iatom)%lmselect,nhat1,nkxc1,non_magnetic_xc,mesh_size,nspden,option,&
 559 &         pawang,pawrad(itypat),pawxcdev,trho1,usetcore,2*usexcnhat,vxc_tmp,xclevel,xc_denpos)
 560        else
 561          call pawxcm_dfpt(pawtab(itypat)%tcoredens(:,1),&
 562 &         cplex,cplex,eexc,ixc,paw_an0(iatom)%kxct1,lm_size,&
 563 &         paw_an(iatom)%lmselect,nhat1,paw_an0(iatom)%nkxc1,non_magnetic_xc,mesh_size,nspden,option,&
 564 &         pawang,pawrad(itypat),trho1,usetcore,2*usexcnhat,vxc_tmp,xclevel)
 565          eexcdc=zero
 566        end if
 567      else
 568        if (ipert==0) then
 569          call pawxc(pawtab(itypat)%tcoredens(:,1),&
 570 &         eexc,eexcdc,hyb_mixing_,ixc,kxc_tmp,k3xc_tmp,lm_size,&
 571 &         paw_an(iatom)%lmselect,nhat1,nkxc1,nk3xc1,non_magnetic_xc,mesh_size,nspden,option,&
 572 &         pawang,pawrad(itypat),trho1,usetcore,2*usexcnhat,vxc_tmp,xclevel,xc_denpos,&
 573 &         coretau=pawtab(itypat)%tcoretau,taur=ttau1,vxctau=vxctau_tmp)
 574        else
 575          call pawxc_dfpt(pawtab(itypat)%tcoredens(:,1),&
 576 &         cplex,cplex,eexc,ixc,paw_an0(iatom)%kxct1,lm_size,&
 577 &         paw_an(iatom)%lmselect,nhat1,paw_an0(iatom)%nkxc1,non_magnetic_xc,mesh_size,nspden,option,&
 578 &         pawang,pawrad(itypat),trho1,usetcore,2*usexcnhat,paw_an0(iatom)%vxct1,vxc_tmp,xclevel)
 579          eexcdc=zero
 580        end if
 581      end if
 582      if (option/=1) then
 583        etild1xc=etild1xc+eexc
 584        etild1xcdc=etild1xcdc+eexcdc
 585      end if
 586      if (option<2) paw_an(iatom)%vxct1(:,:,:)=vxc_tmp(:,:,:)
 587      if (option<2.and.need_vxctau) paw_an(iatom)%vxcttau1(:,:,:)=vxctau_tmp(:,:,:)
 588      if (need_kxc.and. nkxc1>0 ) paw_an(iatom)%kxct1(:,:,:) =kxc_tmp(:,:,:)
 589      if (need_k3xc.and.nk3xc1>0) paw_an(iatom)%k3xct1(:,:,:)=k3xc_tmp(:,:,:)
 590    else ! ipositron==1
 591      if (option<2) paw_an(iatom)%vxct1(:,:,:)=zero
 592      if (need_kxc.and.nkxc1>0) paw_an(iatom)%kxct1(:,:,:)=zero
 593    end if
 594 
 595 !  Additional electron-positron XC term (if ipositron/=0)
 596    if (ipositron/=0) then
 597      if (pawxcdev/=0) then
 598        call pawxcmpositron(ipositron,pawtab(itypat)%tcoredens(:,1),&
 599 &       eexc,eexcdc,electronpositron%ixcpositron,&
 600 &       lm_size,paw_an(iatom)%lmselect,electronpositron%lmselect_ep(1:lm_size,iatom),&
 601 &       nhat1,nhat1_ep,mesh_size,nspden,option,pawang,pawrad(itypat),pawxcdev,&
 602 &       electronpositron%posdensity0_limit,trho1,trho1_ep,usetcore,2*usexcnhat,vxc_tmp,xc_denpos)
 603      else
 604        call pawxcpositron(ipositron,pawtab(itypat)%tcoredens(:,1),&
 605 &       eexc,eexcdc,electronpositron%ixcpositron,&
 606 &       lm_size,paw_an(iatom)%lmselect,electronpositron%lmselect_ep(1:lm_size,iatom),&
 607 &       nhat1,nhat1_ep,mesh_size,nspden,option,pawang,pawrad(itypat),&
 608 &       electronpositron%posdensity0_limit,trho1,trho1_ep,usetcore,2*usexcnhat,vxc_tmp,xc_denpos)
 609      end if
 610      if (option/=1) then
 611        electronpositron%e_paw  =electronpositron%e_paw  -eexc
 612        electronpositron%e_pawdc=electronpositron%e_pawdc-eexcdc
 613      end if
 614      if (option<2) paw_an(iatom)%vxct1(:,:,:)=paw_an(iatom)%vxct1(:,:,:)+vxc_tmp(:,:,:)
 615      if (need_kxc.and.nkxc1>0) paw_an(iatom)%kxct1(:,:,:)=paw_an(iatom)%kxct1(:,:,:)+kxc_tmp(:,:,:)
 616    end if
 617 
 618 !  Update flags defining the state of vxc and kxc
 619    if (option<2) paw_an(iatom)%has_vxc=2
 620    if (option<2.and.need_vxctau) paw_an(iatom)%has_vxctau=2
 621    if (need_kxc.and.nkxc1>0) paw_an(iatom)%has_kxc=2
 622 
 623 !  Update core XC contribution to energy
 624    if (option/=1.and.ipositron/=1) exccore=exccore+pawtab(itypat)%exccore
 625 
 626 !  =========== Compute valence-only XC potentials ===========
 627 !  ==========================================================
 628    if (ipert==0.and.paw_an(iatom)%has_vxcval==1.and.(option==0.or.option==1)) then
 629      if (.not.allocated(paw_an(iatom)%vxc1_val).or..not.allocated(paw_an(iatom)%vxct1_val)) then
 630        msg=' vxc1_val and vxct1_val must be associated'
 631        ABI_BUG(msg)
 632      end if
 633 !    ===== Vxc1_val term, vxc[n1] =====
 634      if (pawxcdev/=0) then
 635        write(msg,'(4a,es16.6)')ch10,&
 636 &       ' pawdenpot : Computing valence-only v_xc[n1] using moments ',ch10,&
 637 &       '             Min density rho1 = ',MINVAL(rho1)
 638        call wrtout(std_out,msg,'COLL')
 639        call pawxcm(pawtab(itypat)%coredens,eexc_val,eexcdc_val,idum,hyb_mixing_,ixc,kxc_tmp,lm_size,&
 640 &       paw_an(iatom)%lmselect,nhat1,nkxc1,non_magnetic_xc,mesh_size,nspden,option,&
 641 &       pawang,pawrad(itypat),pawxcdev,rho1,0,0,vxc_tmp,xclevel,xc_denpos)
 642      else
 643        write(msg,'(2a)')ch10,' pawdenpot : Computing valence-only v_xc[n1] using angular mesh '
 644        call wrtout(std_out,msg,'COLL')
 645 
 646        call pawxc(pawtab(itypat)%coredens,eexc_val,eexcdc_val,hyb_mixing_,ixc,kxc_tmp,k3xc_tmp,lm_size,&
 647 &       paw_an(iatom)%lmselect,nhat1,nkxc1,nk3xc1,non_magnetic_xc,mesh_size,nspden,option,&
 648 &       pawang,pawrad(itypat),rho1,0,0,vxc_tmp,xclevel,xc_denpos)
 649      end if
 650      if (option<2) paw_an(iatom)%vxc1_val(:,:,:)=vxc_tmp(:,:,:)
 651 
 652 !    ===== tVxc1_val term =====
 653      if (pawxcdev/=0) then
 654        if (usexcnhat/=0) then
 655          write(msg,'(4a,e16.6,2a,es16.6)')ch10,&
 656 &         ' pawdenpot : Computing valence-only v_xc[tn1+nhat] using moments ',ch10,&
 657 &         '             Min density trho1        = ',MINVAL(trho1),ch10,&
 658 &         '             Min density trho1 + nhat = ',MINVAL(trho1+nhat1)
 659        else
 660          write(msg,'(4a,e16.6)')ch10,&
 661 &         ' pawdenpot : Computing valence-only v_xc[tn1] using moments ',ch10,&
 662 &         '             Min density trho1        = ',MINVAL(trho1)
 663        end if
 664        call wrtout(std_out,msg,'COLL')
 665        call pawxcm(pawtab(itypat)%tcoredens(:,1),&
 666 &       eexc_val,eexcdc_val,idum,hyb_mixing_,ixc,kxc_tmp,lm_size,&
 667 &       paw_an(iatom)%lmselect,nhat1,nkxc1,non_magnetic_xc,mesh_size,nspden,option,&
 668 &       pawang,pawrad(itypat),pawxcdev,trho1,0,2*usexcnhat,vxc_tmp,xclevel,xc_denpos)
 669      else
 670        write(msg,'(2a)')ch10,' pawdenpot : Computing valence-only v_xc[tn1+nhat] using angular mesh'
 671        call wrtout(std_out,msg,'COLL')
 672        call pawxc(pawtab(itypat)%tcoredens(:,1),&
 673 &       eexc_val,eexcdc_val,hyb_mixing_,ixc,kxc_tmp,k3xc_tmp,lm_size,&
 674 &       paw_an(iatom)%lmselect,nhat1,nkxc1,nk3xc1,non_magnetic_xc,mesh_size,nspden,option,&
 675 &       pawang,pawrad(itypat),trho1,0,2*usexcnhat,vxc_tmp,xclevel,xc_denpos)
 676      end if
 677      if (option<2) then
 678        paw_an(iatom)%vxct1_val(:,:,:)=vxc_tmp(:,:,:)
 679        paw_an(iatom)%has_vxcval=2
 680      end if
 681    end if ! valence-only XC potentials
 682 
 683    ABI_FREE(vxc_tmp)
 684    if (usekden==1) then
 685      ABI_FREE(vxctau_tmp)
 686    end if
 687    ABI_FREE(kxc_tmp)
 688    ABI_FREE(k3xc_tmp)
 689 
 690 !  ===== Compute first part of local exact-exchange energy term =====
 691 !  ===== Also compute corresponding potential                   =====
 692 !  ==================================================================
 693 
 694    if (pawtab(itypat)%useexexch/=0.and.ipert==0.and.ipositron/=1) then
 695 
 696 !    ===== Re-compute a partial "on-site" density n1 (only l=lexexch contrib.)
 697      ABI_MALLOC(rho1xx,(mesh_size,lm_size,nspden))
 698      ABI_MALLOC(lmselect_tmp,(lm_size))
 699      lmselect_tmp(:)=lmselect_cur(:)
 700      call pawdensities(rdum,cplex,iatom_tot,lmselect_cur,lmselect_tmp,lm_size,rdum3,nspden,&
 701 &     1,0,2,pawtab(itypat)%lexexch,0,pawang,pawprtvol,pawrad(itypat),&
 702 &     pawrhoij(iatom),pawtab(itypat),rho1xx,rdum3a,one_over_rad2=one_over_rad2)
 703      ABI_FREE(lmselect_tmp)
 704 !    ===== Re-compute Exc1 and Vxc1; for local exact-exchange, this is done in GGA only
 705      ABI_MALLOC(vxc_tmp,(mesh_size,lm_size,nspden))
 706      ABI_MALLOC(kxc_tmp,(mesh_size,lm_size,nkxc1))
 707      call pawxcm(pawtab(itypat)%coredens,eextemp,eexdctemp,pawtab(itypat)%useexexch,hyb_mixing_,ixc,kxc_tmp,lm_size,&
 708 &     paw_an(iatom)%lmselect,nhat1,nkxc1,non_magnetic_xc,mesh_size,nspden,option,pawang,pawrad(itypat),pawxcdev,&
 709 &     rho1xx,0,0,vxc_tmp,xclevel,xc_denpos)
 710      if (option/=1) then
 711        e1xc=e1xc-eextemp*exchmix
 712        e1xcdc=e1xcdc-eexdctemp*exchmix
 713      end if
 714      if (option<2) then
 715        paw_an(iatom)%vxc_ex(:,:,:)=vxc_tmp(:,:,:)
 716        paw_an(iatom)%has_vxc_ex=2
 717      end if
 718      ABI_FREE(rho1xx)
 719      ABI_FREE(vxc_tmp)
 720      ABI_FREE(kxc_tmp)
 721 
 722    end if ! useexexch
 723 
 724    itypat0=0;if (iatom<my_natom) itypat0=pawrhoij(iatom+1)%itypat
 725    if (itypat/=itypat0) then
 726      ABI_FREE(one_over_rad2)
 727    end if
 728 
 729    ABI_FREE(lmselect_cur)
 730 
 731 !  ==== Compute Hartree potential terms and some energy terms ====
 732 !  ===============================================================
 733 
 734 !  Hartree Dij computation
 735    if (ipositron/=1) then
 736      call pawdijhartree(paw_ij(iatom)%dijhartree,cplex,nspden,pawrhoij(iatom),pawtab(itypat))
 737    else
 738      paw_ij(iatom)%dijhartree(:)=zero
 739    end if
 740    paw_ij(iatom)%has_dijhartree=2
 741 
 742 !  Hartree energy computation
 743    if (option/=1) then
 744      call pawaccenergy_nospin(eh2,pawrhoij(iatom),paw_ij(iatom)%dijhartree,1,qphase,pawtab(itypat))
 745    end if
 746 
 747 !  Electron-positron calculation:
 748 !    - Compute Dij due to fixed particles (elec. or pos. depending on calctype)
 749 !    - Compute contribution to energy
 750 !    - Add electron and positron
 751    if (ipositron/=0) then
 752      ABI_CHECK(qphase==1,'qphase should be 1 for electron-positron!')
 753      ABI_MALLOC(dij_ep,(qphase*lmn2_size))
 754      call pawdijhartree(dij_ep,qphase,nspden,electronpositron%pawrhoij_ep(iatom),pawtab(itypat))
 755      if (option/=1) then
 756        etmp=zero
 757        call pawaccenergy_nospin(etmp,pawrhoij(iatom),dij_ep,1,1,pawtab(itypat))
 758        electronpositron%e_paw  =electronpositron%e_paw  -etmp
 759        electronpositron%e_pawdc=electronpositron%e_pawdc-etmp
 760      end if
 761      paw_ij(iatom)%dijhartree(:)=paw_ij(iatom)%dijhartree(:)-dij_ep(:)
 762      ABI_FREE(dij_ep)
 763    end if
 764 
 765 !  Compute 1st moment of total Hartree potential VH(n_Z+n_core+n1)
 766 !  Equation 10 (density) and up to 43 (Hartree potential of density)
 767 !    of Kresse and Joubert PRB 59 1758 (1999) [[cite:Kresse1999]]
 768    keep_vhartree=(paw_an(iatom)%has_vhartree>0)
 769    if ((pawspnorb>0.and.ipert==0.and.ipositron/=1).or.keep_vhartree) then
 770 
 771      !In the first clause case, would it not be simpler just to turn on has_vhartree?
 772      if (.not. allocated(paw_an(iatom)%vh1)) then
 773        ABI_MALLOC(paw_an(iatom)%vh1,(cplex*mesh_size,1,1))
 774      end if
 775      if (.not. allocated(paw_an(iatom)%vht1)) then
 776        ABI_MALLOC(paw_an(iatom)%vht1,(cplex*mesh_size,1,1))
 777      end if
 778      ABI_MALLOC(rho,(mesh_size))
 779      ABI_MALLOC(vh,(mesh_size))
 780 
 781 !    Construct vh1 and tvh1
 782      do iq=1,cplex
 783        !Construct vh1
 784        !  The sqrt(4pi) factor comes from the fact we are calculating the spherical moments,
 785        !   and for the 00 channel the prefactor of Y_00 = 2 sqrt(pi)
 786        rho(1:mesh_size)=rho1(iq:cplex*mesh_size:cplex,1,1)*four_pi*pawrad(itypat)%rad(1:mesh_size)**2
 787        if (usecore==1) then
 788          rho(1:mesh_size)=rho(1:mesh_size)+sqrt(four_pi)*pawtab(itypat)%coredens(1:mesh_size) &
 789 &                                         *four_pi*pawrad(itypat)%rad(1:mesh_size)**2
 790        end if
 791        call poisson(rho,0,pawrad(itypat),vh)
 792        vh(2:mesh_size)=(vh(2:mesh_size)-sqrt(four_pi)*znucl(itypat))/pawrad(itypat)%rad(2:mesh_size)
 793        call pawrad_deducer0(vh,mesh_size,pawrad(itypat))
 794        paw_an(iatom)%vh1(iq:cplex*mesh_size:cplex,1,1)=vh(1:mesh_size)
 795 !      TODO: check this is equivalent to the previous version (commented) which explicitly recalculated VH(coredens)
 796 !      DONE: numerically there are residual differences on abiref (7th digit).
 797 !            paw_an(iatom)%vh1(2:mesh_size,1,1)=paw_an(iatom)%vh1(2:mesh_size,1,1)/pawrad(itypat)%rad(2:mesh_size) &
 798 !&                                          +sqrt(four_pi) * pawtab(itypat)%VHnZC(2:mesh_size)
 799 
 800        !Same for vht1
 801        rho(1:mesh_size)=trho1(iq:cplex*mesh_size:cplex,1,1)*four_pi*pawrad(itypat)%rad(1:mesh_size)**2
 802        if (usenhat/=0) then
 803          rho(1:mesh_size)=rho(1:mesh_size)+nhat1(iq:cplex*mesh_size:cplex,1,1) &
 804 &                                         *four_pi*pawrad(itypat)%rad(1:mesh_size)**2
 805        end if
 806        if (usetcore==1) then
 807          rho(1:mesh_size)=rho(1:mesh_size)+sqrt(four_pi)*pawtab(itypat)%tcoredens(1:mesh_size,1) &
 808 &                                         *four_pi*pawrad(itypat)%rad(1:mesh_size)**2
 809        end if
 810        call poisson(rho,0,pawrad(itypat),vh)
 811        vh(2:mesh_size)=(vh(2:mesh_size)-sqrt(four_pi)*znucl(itypat))/pawrad(itypat)%rad(2:mesh_size)
 812        call pawrad_deducer0(vh,mesh_size,pawrad(itypat))
 813        paw_an(iatom)%vht1(iq:cplex*mesh_size:cplex,1,1)=vh(1:mesh_size)
 814 
 815      end do ! cplex phase
 816 
 817      paw_an(iatom)%has_vhartree=2
 818      ABI_FREE(rho)
 819      ABI_FREE(vh)
 820    end if
 821 
 822 !  ========= Compute PAW+U and energy contribution  =========
 823 !  ==========================================================
 824 
 825    if (usepawu/=0.and.usepawu<10.and.ipositron/=1.and.option/=1) then
 826 
 827      if (pawu_algo==PAWU_ALGO_1) then
 828 
 829 !      PAW+U energy computation from nocc_m_mp
 830        call pawuenergy(iatom_tot,edftumdc,edftumdcdc,paw_ij(iatom)%noccmmp, &
 831 &                      paw_ij(iatom)%nocctot,pawprtvol,pawtab(itypat))
 832      else
 833 
 834 !      PAW+U energy computation from eU_ijkl
 835        !First, compute DijU
 836        call pawdiju_euijkl(paw_ij(iatom)%dijU,cplex_dij,cplex,ndij, &
 837 &                          pawrhoij(iatom),pawtab(itypat))
 838        paw_ij(iatom)%has_dijU=2
 839        !Then, compute energy
 840        if (option/=1) then
 841          etmp=zero
 842          call pawaccenergy(etmp,pawrhoij(iatom),paw_ij(iatom)%dijU,cplex_dij,qphase,ndij,pawtab(itypat))
 843          edftumdc=edftumdc+half*etmp ; edftumdcdc=edftumdcdc-half*etmp
 844          !Add FLL double-counting part
 845          if (pawu_dblec==PAWU_FLL.and.ipert==0) then
 846            ABI_CHECK(qphase==1,'BUG in pawdenpot: qphase should be 1 for Dble-C FLL term!')
 847            call pawaccenergy_nospin(edftufll,pawrhoij(iatom),pawtab(itypat)%euij_fll,1,1,pawtab(itypat))
 848          end if
 849        end if
 850 
 851      end if ! DFT+U algo
 852    end if ! Dij Hartree
 853 
 854 !  ========= Compute nuclear dipole moment energy contribution  ========
 855 !  =====================================================================
 856 
 857    if (any(abs(nucdipmom(:,iatom))>tol8).and.ipert==0.and.ipositron/=1) then
 858 
 859      ABI_CHECK(cplex_rhoij==2,'BUG in pawdenpot: rhoij must be complex for ND moments!')
 860      ABI_CHECK(qphase==1,'BUG in pawdenpot: qphase should be 1 for ND moments!')
 861 
 862 !    Compute nuclear dipole contribution to Dij if necessary
 863      if (paw_ij(iatom)%has_dijnd/=2) then
 864        call pawdijnd(paw_ij(iatom)%dijnd,cplex_dij,ndij,nucdipmom(:,iatom),pawrad(itypat),pawtab(itypat))
 865        paw_ij(iatom)%has_dijnd=2
 866      end if
 867 
 868 !    Compute nuclear dipole contribution to energy
 869      if (option/=1) then
 870        call pawaccenergy_nospin(enucdip,pawrhoij(iatom),paw_ij(iatom)%dijnd,cplex_dij,1,pawtab(itypat))
 871      end if
 872 
 873    end if
 874 
 875 !  ========= Compute spin-orbit energy contribution  ========
 876 !  ==========================================================
 877 
 878    if (pawspnorb>0.and.ipert==0.and.ipositron/=1) then
 879 
 880 !    Compute spin-orbit contribution to Dij
 881      if (option/=2.or.cplex_rhoij==2) then
 882        call pawdijso(paw_ij(iatom)%dijso,cplex_dij,cplex,ndij,nspden,pawang,pawrad(itypat),pawtab(itypat), &
 883 &                    pawxcdev,spnorbscl,paw_an(iatom)%vh1,paw_an(iatom)%vxc1)
 884        paw_ij(iatom)%has_dijso=2
 885      end if
 886 
 887 !    Compute spin-orbit contribution to on-site energy
 888      if (option/=1.and.cplex_rhoij==2) then
 889        call pawaccenergy(espnorb,pawrhoij(iatom),paw_ij(iatom)%dijso,cplex_dij,qphase,ndij,pawtab(itypat))
 890      end if
 891 
 892    end if
 893 
 894 !  === Compute 2nd part of local exact-exchange energy and potential  ===
 895 !  ======================================================================
 896 
 897    if (pawtab(itypat)%useexexch/=0.and.ipert==0.and.ipositron/=1) then
 898 
 899      ABI_CHECK(paw_ij(iatom)%nspden/=4,'BUG in pawdenpot: Local ex-exch. not implemented for nspden=4!')
 900      if (option<2) then
 901        call pawxpot(ndij,pawprtvol,pawrhoij(iatom),pawtab(itypat),paw_ij(iatom)%vpawx)
 902        paw_ij(iatom)%has_exexch_pot=2
 903      end if
 904      if (option/=1) then
 905        if (abs(pawprtvol)>=2) then
 906          write(msg, '(2a)' )ch10,'======= PAW local exact exchange terms (in Hartree) ===='
 907          call wrtout(std_out,  msg,'COLL')
 908          write(msg, '(2a,i4)' )ch10,' For Atom',iatom_tot
 909          call wrtout(std_out,  msg,'COLL')
 910        end if
 911        call pawxenergy(eexex,pawprtvol,pawrhoij(iatom),pawtab(itypat))
 912      end if
 913 
 914    end if ! useexexch
 915 
 916 !  ==== Compute Fock Dij term and Fock energy terms ====
 917 !  =====================================================
 918 
 919    if (usefock==1) then
 920 
 921      if (ipositron/=1) then
 922 
 923 !      Fock contribution to Dij
 924        ABI_MALLOC(dijfock_vv,(cplex_dij*qphase*lmn2_size,ndij))
 925        ABI_MALLOC(dijfock_cv,(cplex_dij*qphase*lmn2_size,ndij))
 926        call pawdijfock(dijfock_vv,dijfock_cv,cplex_dij,cplex,hyb_mixing_,hyb_mixing_sr_, &
 927 &                      ndij,pawrhoij(iatom),pawtab(itypat))
 928        paw_ij(iatom)%dijfock(:,:)=dijfock_vv(:,:)+dijfock_cv(:,:)
 929        paw_ij(iatom)%has_dijfock=2
 930 
 931 !      Fock contribution to energy
 932        if (option/=1) then
 933          dijfock_vv(:,:)=half*dijfock_vv(:,:) ; dijfock_cv(:,:)=dijfock_vv(:,:)+dijfock_cv(:,:)
 934          call pawaccenergy(efock  ,pawrhoij(iatom),dijfock_cv,cplex_dij,qphase,ndij,pawtab(itypat))
 935          call pawaccenergy(efockdc,pawrhoij(iatom),dijfock_vv,cplex_dij,qphase,ndij,pawtab(itypat))
 936        end if
 937 
 938        ABI_FREE(dijfock_vv)
 939        ABI_FREE(dijfock_cv)
 940      end if
 941 
 942 !    Special case for positron
 943      if (ipositron==1) then
 944        paw_ij(iatom)%dijfock(:,:)=zero
 945        paw_ij(iatom)%has_dijfock=2
 946      end if
 947 
 948    end if
 949 
 950 !  === Compute the zero of the potentials if requested ==================
 951 !  ======================================================================
 952 
 953    if (pawtab(itypat)%usepotzero==1.and.present(vpotzero).and.ipert<=0) then
 954 
 955      !Term 1 : beta
 956      vpotzero(1)=vpotzero(1)-pawtab(itypat)%beta/ucvol
 957 
 958      !Term 2 : \sum_ij rho_ij gamma_ij
 959      etmp=zero
 960      call pawaccenergy_nospin(etmp,pawrhoij(iatom),pawtab(itypat)%gammaij,1,1,pawtab(itypat))
 961      vpotzero(2)=vpotzero(2)-etmp/ucvol
 962 
 963    end if
 964 
 965 !  ======= Compute atomic contribution to the energy (Dij0)   ===========
 966 !  ======================================================================
 967 
 968    if (option/=1.and.ipert<=0) then
 969      call pawaccenergy_nospin(e1t10,pawrhoij(iatom),pawtab(itypat)%dij0,1,1,pawtab(itypat))
 970 
 971 !    Positron special case (dij0 is opposite, except for kinetic term)
 972      if (ipositron==1) then
 973        ABI_MALLOC(dij_ep,(lmn2_size))
 974        dij_ep(:)=two*(pawtab(itypat)%kij(:)-pawtab(itypat)%dij0(:))
 975        call pawaccenergy_nospin(e1t10,pawrhoij(iatom),dij_ep,1,1,pawtab(itypat))
 976        ABI_FREE(dij_ep)
 977      end if
 978 
 979    end if
 980 
 981 !  ==========================================================
 982 !  No more need of some densities/potentials
 983 
 984 !  Deallocate densities
 985    ABI_FREE(rho1)
 986    ABI_FREE(trho1)
 987    ABI_FREE(nhat1)
 988    if (usekden==1) then
 989      ABI_FREE(tau1)
 990      ABI_FREE(ttau1)
 991    end if
 992    if (ipositron/=0)  then
 993      ABI_FREE(rho1_ep)
 994      ABI_FREE(trho1_ep)
 995      ABI_FREE(nhat1_ep)
 996    end if
 997 
 998 !  Deallocate potentials
 999    if (.not.keep_vhartree) then
1000      paw_an(iatom)%has_vhartree=0
1001      if (allocated(paw_an(iatom)%vh1)) then
1002        ABI_FREE(paw_an(iatom)%vh1)
1003      end if
1004    end if
1005    if (temp_vxc) then
1006      paw_an(iatom)%has_vxc=0
1007      if (allocated(paw_an(iatom)%vxc1)) then
1008        ABI_FREE(paw_an(iatom)%vxc1)
1009      end if
1010    end if
1011 
1012 !  =========== End loop on atoms ============================
1013 !  ==========================================================
1014 
1015  end do
1016 
1017 !========== Assemble "on-site" energy terms ===============
1018 !==========================================================
1019 
1020  if (option/=1) then
1021    if (ipert==0) then
1022      epaw=e1xc+half*eh2+e1t10-exccore-etild1xc+edftumdc+edftufll+eexex+espnorb+efock+enucdip
1023      epawdc=e1xc-e1xcdc-half*eh2-exccore-etild1xc+etild1xcdc+edftumdcdc-eexex-efockdc
1024    else
1025      epaw=e1xc-etild1xc+eh2+two*edftumdc
1026      epawdc=zero
1027    end if
1028  end if
1029 
1030 !========== Reduction in case of parallelism ==============
1031 !==========================================================
1032 
1033  if (paral_atom) then
1034    if (option/=1)  then
1035      call timab(48,1,tsec)
1036      mpiarr=zero
1037      mpiarr(1)=compch_sph;mpiarr(2)=epaw;mpiarr(3)=epawdc
1038      if (ipositron/=0) then
1039        mpiarr(4)=electronpositron%e_paw
1040        mpiarr(5)=electronpositron%e_pawdc
1041      end if
1042      if (present(vpotzero)) then
1043        mpiarr(6)=vpotzero(1)
1044        mpiarr(7)=vpotzero(2)
1045      end if
1046      call xmpi_sum(mpiarr,my_comm_atom,ierr)
1047      compch_sph=mpiarr(1);epaw=mpiarr(2);epawdc=mpiarr(3)
1048      if (ipositron/=0) then
1049        electronpositron%e_paw=mpiarr(4)
1050        electronpositron%e_pawdc=mpiarr(5)
1051      end if
1052      if (present(vpotzero)) then
1053        vpotzero(1)=mpiarr(6)
1054        vpotzero(2)=mpiarr(7)
1055      end if
1056      call timab(48,2,tsec)
1057    end if
1058  end if
1059 
1060 !Destroy atom table used for parallelism
1061  call free_my_atmtab(my_atmtab,my_atmtab_allocated)
1062 
1063  call timab(560,2,tsec)
1064 
1065  DBG_EXIT("COLL")
1066 
1067 end subroutine pawdenpot

m_paw_denpot/pawdensities [ Functions ]

[ Top ] [ m_paw_denpot ] [ Functions ]

NAME

 pawdensities

FUNCTION

 Compute PAW on-site densities (all-electron, pseudo and compensation) for a given atom

INPUTS

  cplex: if 1, on-site densities are REAL, if 2, COMPLEX (response function only)
  iatom=index of current atom (note: this is the absolute index, not the index on current proc)
  lm_size=number of (l,m) moments
  lmselectin(lm_size)=flags selecting the non-zero LM-moments of on-site densities
                      (value of these flags at input; must be .TRUE. for nzlmopt/=1)
  nspden=number of spin-density components
  nzlmopt=if -1, compute all LM-moments of densities (lmselectin=.true. forced)
                 initialize "lmselectout" (index of non-zero LM-moments of densities)
          if  0, compute all LM-moments of densities (lmselectin=.true. forced)
                 force "lmselectout" to .true. (index of non-zero LM-moments of densities)
          if  1, compute only non-zero LM-moments of densities (stored before in "lmselectin")
  one_over_rad2(mesh_size)= contains 1/r**2 for each point of the radial grid -optional argument-
  opt_compch=flag controlling the accumulation of compensation charge density moments
             inside PAW spheres (compch_sph)
  opt_dens=flag controlling which on-site density(ies) is (are) computed
           0: all on-site densities (all-electron, pseudo and compensation)
           1: all-electron and pseudo densities (no compensation)
           2: only all-electron density
  opt_l=controls which l-moment(s) contribute to the density:
        <0 : all l contribute
        >=0: only l=opt_l contributes
        Note: opt_l>=0 is only compatible with opt_dens=2
  opt_print=1 if the densities moments have to be printed out (only if pawprtvol>=2)
  pawang <type(pawang_type)>=paw angular mesh and related data
  pawprtvol=control print volume and debugging output for PAW
  pawrad <type(pawrad_type)>=paw radial mesh and related data (for the current atom type)
  pawrhoij <type(pawrhoij_type)>= paw rhoij occupancies and related data (for the current atom)
  pawtab <type(pawtab_type)>=paw tabulated starting data (for the current atom type)

OUTPUT

  nhat1(cplex*mesh_size,lm_size,nspden)= compensation charge on-site density for current atom
  rho1(cplex*mesh_size,lm_size,nspden)= all electron on-site density for current atom
  trho1(cplex*mesh_size,lm_size,nspden)= pseudo on-site density for current atom
  ==== if nzlmopt/=1
    lmselectout(lm_size)=flags selecting the non-zero LM-moments of on-site densities
                         (value of these flags at output if updated, i.e. if nzlmopt<1)

SIDE EFFECTS

  ==== if opt_compch==1
    compch_sph=compensation charge integral inside spheres computed over spherical meshes
               updated with the contribution of current atom

SOURCE

1124 subroutine pawdensities(compch_sph,cplex,iatom,lmselectin,lmselectout,lm_size,nhat1,nspden,nzlmopt,&
1125 &          opt_compch,opt_dens,opt_l,opt_print,pawang,pawprtvol,pawrad,pawrhoij,pawtab,rho1,trho1,&
1126 &          one_over_rad2) ! optional
1127 
1128 !Arguments ---------------------------------------------
1129 !scalars
1130  integer,intent(in) :: cplex,iatom,lm_size,nspden,nzlmopt,opt_compch,opt_dens,opt_l,opt_print,pawprtvol
1131  real(dp),intent(inout) :: compch_sph
1132  type(pawang_type),intent(in) :: pawang
1133  type(pawrad_type),intent(in) :: pawrad
1134  type(pawrhoij_type),intent(in) :: pawrhoij
1135  type(pawtab_type),intent(in) :: pawtab
1136 !arrays
1137  logical,intent(in) :: lmselectin(lm_size)
1138  logical,intent(inout) :: lmselectout(lm_size)
1139  real(dp),intent(in),target,optional :: one_over_rad2(pawtab%mesh_size)
1140  real(dp),intent(out) :: nhat1(cplex*pawtab%mesh_size,lm_size,nspden*(1-((opt_dens+1)/2)))
1141  real(dp),intent(out) ::  rho1(cplex*pawtab%mesh_size,lm_size,nspden)
1142  real(dp),intent(out) :: trho1(cplex*pawtab%mesh_size,lm_size,nspden*(1-(opt_dens/2)))
1143 !Local variables ---------------------------------------
1144 !scalars
1145  integer :: dplex,ii,ilm,iplex,iq0,ir,irhoij,isel,ispden,jrhoij
1146  integer :: klm,klmn,kln,ll,lmax,lmin,mesh_size
1147  real(dp) :: m1,mt1,rdum
1148  character(len=500) :: msg
1149 !arrays
1150  real(dp) :: compchspha(cplex),compchsphb(cplex),ro(cplex),ro_ql(cplex),ro_rg(cplex)
1151  real(dp),allocatable :: aa(:),bb(:)
1152  real(dp),pointer :: one_over_rad2_(:)
1153 
1154 ! *************************************************************************
1155 
1156  DBG_ENTER("COLL")
1157 
1158 !Compatibility tests
1159  if (opt_dens/=2.and.opt_l>=0) then
1160    msg='opt_dens/=2 incompatible with opt_l>=0!'
1161    ABI_BUG(msg)
1162  end if
1163  if(nzlmopt/=0.and.nzlmopt/=1.and.nzlmopt/=-1) then
1164    msg='invalid value for variable "nzlmopt"!'
1165    ABI_BUG(msg)
1166  end if
1167  if(nspden>pawrhoij%nspden) then
1168    msg='nspden must be <= pawrhoij%nspden!'
1169    ABI_BUG(msg)
1170  end if
1171  if (cplex>pawrhoij%qphase) then
1172    msg='cplex must be <= pawrhoij%qphase!'
1173    ABI_BUG(msg)
1174  end if
1175  if (nzlmopt/=1) then
1176    if (any(.not.lmselectin(1:lm_size))) then
1177      msg='With nzlmopt/=1, lmselectin must be true!'
1178      ABI_BUG(msg)
1179    end if
1180  end if
1181  if (pawang%gnt_option==0) then
1182    msg='pawang%gnt_option=0!'
1183    ABI_BUG(msg)
1184  end if
1185 
1186 !Various inits
1187  rho1=zero
1188  if (opt_dens<2) trho1=zero
1189  if (opt_dens==0) nhat1=zero
1190  mesh_size=pawtab%mesh_size;dplex=cplex-1
1191  iq0=pawrhoij%cplex_rhoij*pawrhoij%lmn2_size
1192  if (nzlmopt<1) lmselectout(1:lm_size)=.true.
1193  if (present(one_over_rad2)) then
1194    one_over_rad2_ => one_over_rad2
1195  else
1196    ABI_MALLOC(one_over_rad2_,(mesh_size))
1197    one_over_rad2_(1)=zero
1198    one_over_rad2_(2:mesh_size)=one/pawrad%rad(2:mesh_size)**2
1199  end if
1200 
1201 !===== Compute "on-site" densities (n1, ntild1, nhat1) =====
1202 !===========================================================
1203 
1204  do ispden=1,nspden
1205 
1206 !  -- Loop over ij channels (basis components)
1207    jrhoij=1
1208    do irhoij=1,pawrhoij%nrhoijsel
1209      klmn=pawrhoij%rhoijselect(irhoij)
1210      klm =pawtab%indklmn(1,klmn)
1211      kln =pawtab%indklmn(2,klmn)
1212      lmin=pawtab%indklmn(3,klmn)
1213      lmax=pawtab%indklmn(4,klmn)
1214 
1215 !    Retrieve rhoij
1216      if (pawrhoij%nspden/=2) then
1217        ro(1)=pawrhoij%rhoijp(jrhoij,ispden)
1218        if (cplex==2) ro(2)=pawrhoij%rhoijp(iq0+jrhoij,ispden)
1219      else
1220        if (ispden==1) then
1221          ro(1)=pawrhoij%rhoijp(jrhoij,1)+pawrhoij%rhoijp(jrhoij,2)
1222          if (cplex==2) ro(2)=pawrhoij%rhoijp(iq0+jrhoij,1)+pawrhoij%rhoijp(iq0+jrhoij,2)
1223        else if (ispden==2) then
1224          ro(1)=pawrhoij%rhoijp(jrhoij,1)
1225          if (cplex==2) ro(2)=pawrhoij%rhoijp(iq0+jrhoij,1)
1226        end if
1227      end if
1228      ro(1:cplex)=pawtab%dltij(klmn)*ro(1:cplex)
1229 
1230 !    First option: all on-site densities are computed (opt_dens==0)
1231 !    --------------------------------------------------------------
1232      if (opt_dens==0) then
1233        do ll=lmin,lmax,2
1234          do ilm=ll**2+1,(ll+1)**2
1235            if (lmselectin(ilm)) then
1236              isel=pawang%gntselect(ilm,klm)
1237              if (isel>0) then
1238                ro_ql(1:cplex)=ro(1:cplex)*pawtab%qijl(ilm,klmn)
1239                ro_rg(1:cplex)=ro(1:cplex)*pawang%realgnt(isel)
1240 !              == nhat1(r=0)
1241                nhat1(1:cplex,ilm,ispden)=nhat1(1:cplex,ilm,ispden) &
1242 &               +ro_ql(1:cplex)*pawtab%shapefunc(1,ll+1)
1243 !              == rho1(r>0), trho1(r>0), nhat1(r>0)
1244                do ir=2,mesh_size
1245                  rho1(cplex*ir-dplex:ir*cplex,ilm,ispden) =rho1(cplex*ir-dplex:ir*cplex,ilm,ispden)&
1246 &                 +ro_rg(1:cplex)*pawtab%phiphj(ir,kln)*one_over_rad2_(ir)
1247                  trho1(cplex*ir-dplex:ir*cplex,ilm,ispden)=trho1(cplex*ir-dplex:ir*cplex,ilm,ispden)&
1248 &                 +ro_rg(1:cplex)*pawtab%tphitphj(ir,kln)*one_over_rad2_(ir)
1249                  nhat1(cplex*ir-dplex:ir*cplex,ilm,ispden)=nhat1(cplex*ir-dplex:ir*cplex,ilm,ispden)&
1250 &                 +ro_ql(1:cplex)*pawtab%shapefunc(ir,ll+1)
1251                end do
1252              end if
1253            end if
1254          end do  ! End loops over ll,lm
1255        end do
1256 
1257 !      2nd option: AE and pseudo densities are computed (opt_dens==1)
1258 !      --------------------------------------------------------------
1259      else if (opt_dens==1) then
1260        do ll=lmin,lmax,2
1261          do ilm=ll**2+1,(ll+1)**2
1262            if (lmselectin(ilm)) then
1263              isel=pawang%gntselect(ilm,klm)
1264              if (isel>0) then
1265                ro_rg(1:cplex)=ro(1:cplex)*pawang%realgnt(isel)
1266 !              == rho1(r>0), trho1(r>0)
1267                do ir=2,mesh_size
1268                  rho1(cplex*ir-dplex:ir*cplex,ilm,ispden) =rho1(cplex*ir-dplex:ir*cplex,ilm,ispden)&
1269 &                 +ro_rg(1:cplex)*pawtab%phiphj  (ir,kln)*one_over_rad2_(ir)
1270                  trho1(cplex*ir-dplex:ir*cplex,ilm,ispden)=trho1(cplex*ir-dplex:ir*cplex,ilm,ispden)&
1271 &                 +ro_rg(1:cplex)*pawtab%tphitphj(ir,kln)*one_over_rad2_(ir)
1272                end do
1273              end if
1274            end if
1275          end do  ! End loops over ll,lm
1276        end do
1277 
1278 !      3rd option: only all-electron on-site density is computed (opt_dens==2)
1279 !      -----------------------------------------------------------------------
1280      else if (opt_dens==2) then
1281        if (opt_l<0.or.(pawtab%indklmn(3,klmn)==0.and.pawtab%indklmn(4,klmn)==2*opt_l)) then
1282          do ll=lmin,lmax,2
1283            do ilm=ll**2+1,(ll+1)**2
1284              if (lmselectin(ilm)) then
1285                isel=pawang%gntselect(ilm,klm)
1286                if (isel>0) then
1287                  ro_rg(1:cplex)=ro(1:cplex)*pawang%realgnt(isel)
1288 !                == rho1(r>0)
1289                  do ir=2,mesh_size
1290                    rho1(cplex*ir-dplex:ir*cplex,ilm,ispden) =rho1(cplex*ir-dplex:ir*cplex,ilm,ispden)&
1291 &                   +ro_rg(1:cplex)*pawtab%phiphj(ir,kln)*one_over_rad2_(ir)
1292                  end do
1293                end if
1294              end if
1295            end do  ! End loops over ll, lm
1296          end do
1297        end if
1298      end if
1299 
1300 
1301 !    -- End loop over ij channels
1302      jrhoij=jrhoij+pawrhoij%cplex_rhoij
1303    end do
1304 
1305 !  Compute rho1(r=0) and trho1(r=0)
1306    if (cplex==2)  then
1307      ABI_MALLOC(aa,(5))
1308      ABI_MALLOC(bb,(5))
1309    end if
1310    if (opt_dens==0.or.opt_dens==1) then
1311      do ll=0,pawtab%lcut_size-1
1312        do ilm=ll**2+1,(ll+1)**2
1313          if (lmselectin(ilm)) then
1314            if (cplex==1) then
1315              call pawrad_deducer0(rho1 (:,ilm,ispden),mesh_size,pawrad)
1316              call pawrad_deducer0(trho1(:,ilm,ispden),mesh_size,pawrad)
1317            else
1318              do ii=0,1
1319                do ir=2,5
1320                  aa(ir)=rho1 (2*ir-ii,ilm,ispden)
1321                  bb(ir)=trho1(2*ir-ii,ilm,ispden)
1322                end do
1323                call pawrad_deducer0(aa,5,pawrad)
1324                call pawrad_deducer0(bb,5,pawrad)
1325                rho1 (2-ii,ilm,ispden)=aa(1)
1326                trho1(2-ii,ilm,ispden)=bb(1)
1327              end do
1328            end if
1329          end if
1330        end do
1331      end do
1332    else
1333      do ll=0,pawtab%lcut_size-1
1334        do ilm=ll**2+1,(ll+1)**2
1335          if (lmselectin(ilm)) then
1336            if (cplex==1) then
1337              call pawrad_deducer0(rho1(:,ilm,ispden),mesh_size,pawrad)
1338            else
1339              do ii=0,1
1340                do ir=2,5
1341                  aa(ir)=rho1 (2*ir-ii,ilm,ispden)
1342                end do
1343                call pawrad_deducer0(aa,5,pawrad)
1344                rho1(2-ii,ilm,ispden)=aa(1)
1345              end do
1346            end if
1347          end if
1348        end do
1349      end do
1350    end if
1351    if (cplex==2)  then
1352      ABI_FREE(aa)
1353      ABI_FREE(bb)
1354    end if
1355 
1356 !  -- Test moments of densities and store non-zero ones
1357    if (nzlmopt==-1) then
1358      do ll=0,pawtab%lcut_size-1
1359        do ilm=ll**2+1,(ll+1)**2
1360          m1=zero;mt1=zero
1361          if (cplex==1) then
1362            m1=maxval(abs(rho1 (1:mesh_size,ilm,ispden)))
1363            if (opt_dens<2) mt1=maxval(abs(trho1(1:mesh_size,ilm,ispden)))
1364          else
1365            do ir=1,mesh_size
1366              rdum=sqrt(rho1(2*ir-1,ilm,ispden)**2+rho1(2*ir,ilm,ispden)**2)
1367              m1=max(m1,rdum)
1368            end do
1369            if (opt_dens<2) then
1370              do ir=1,mesh_size
1371                rdum=sqrt(trho1(2*ir-1,ilm,ispden)**2+trho1(2*ir,ilm,ispden)**2)
1372                mt1=max(mt1,rdum)
1373              end do
1374            end if
1375          end if
1376          if (ispden==1) then
1377            if ((ilm>1).and.(m1<tol16).and.(mt1<tol16)) then
1378              lmselectout(ilm)=.false.
1379            end if
1380          else if (.not.(lmselectout(ilm))) then
1381            lmselectout(ilm)=((m1>=tol16).or.(mt1>=tol16))
1382          end if
1383        end do
1384      end do
1385    end if
1386 
1387 !  -- Compute integral of (n1-tn1) inside spheres
1388    if (opt_compch==1.and.ispden==1.and.opt_dens<2) then
1389      ABI_MALLOC(aa,(mesh_size))
1390      aa(1:mesh_size)=(rho1(1:mesh_size,1,1)-trho1(1:mesh_size,1,1)) &
1391 &     *pawrad%rad(1:mesh_size)**2
1392      call simp_gen(compchspha(1),aa,pawrad)
1393      compch_sph=compch_sph+compchspha(1)*sqrt(four_pi)
1394      ABI_FREE(aa)
1395    end if
1396 
1397 !  -- Print out moments of densities (if requested)
1398    if (abs(pawprtvol)>=2.and.opt_print==1.and.opt_dens<2) then
1399      ABI_MALLOC(aa,(cplex*mesh_size))
1400      ABI_MALLOC(bb,(cplex*mesh_size))
1401      if (opt_dens==0) then
1402        write(msg,'(2a,i3,a,i1,3a)') ch10, &
1403 &       ' Atom ',iatom,' (ispden=',ispden,'):',ch10,&
1404 &       '  ******* Moment of (n1-tn1)  ** Moment of (n1-tn1-nhat1)'
1405      else
1406        write(msg,'(2a,i3,a,i1,3a)') ch10, &
1407 &       ' Atom ',iatom,' (ispden=',ispden,'):',ch10,&
1408 &       '  ******* Moment of (n1-tn1)'
1409      end if
1410      call wrtout(std_out,msg,'PERS')
1411      do ll=0,pawtab%lcut_size-1
1412        do ilm=ll**2+1,(ll+1)**2
1413          if (lmselectin(ilm)) then
1414            do iplex=1,cplex
1415              if (opt_dens==0) then
1416                do ir=1,mesh_size
1417                  ii=cplex*(ir-1)+iplex
1418                  ro(1)=pawrad%rad(ir)**(2+ll)
1419                  aa(ir)=ro(1)*(rho1(ii,ilm,ispden)-trho1(ii,ilm,ispden))
1420                  bb(ir)=ro(1)*nhat1(ii,ilm,ispden)
1421                end do
1422                call simp_gen(compchspha(iplex),aa,pawrad)
1423                call simp_gen(compchsphb(iplex),bb,pawrad)
1424              else
1425                do ir=1,mesh_size
1426                  ii=cplex*(ir-1)+iplex
1427                  ro(1)=pawrad%rad(ir)**(2+ll)
1428                  aa(ir)=ro(1)*(rho1(ii,ilm,ispden)-trho1(ii,ilm,ispden))
1429                end do
1430                call simp_gen(compchspha(iplex),aa,pawrad)
1431              end if
1432            end do
1433            if (opt_dens==0) then
1434              if (cplex==1) then
1435                write(msg,'(3x,a,2i2,2(a,es14.7))') &
1436 &               'l,m=',ll,ilm-(ll**2+ll+1),': M=',compchspha(1),&
1437 &               ' **    M=',compchspha(1)-compchsphb(1)
1438              else
1439                write(msg,'(3x,a,2i2,2(a,2es14.7))') &
1440 &               'l,m=',ll,ilm-(ll**2+ll+1),': M=',compchspha(1:2),&
1441 &               ' **    M=',compchspha(1:2)-compchsphb(1:2)
1442              end if
1443            else
1444              if (cplex==1) then
1445                write(msg,'(3x,a,2i2,a,es14.7)') &
1446 &               'l,m=',ll,ilm-(ll**2+ll+1),': M=',compchspha(1)
1447              else
1448                write(msg,'(3x,a,2i2,a,2es14.7)') &
1449 &               'l,m=',ll,ilm-(ll**2+ll+1),': M=',compchspha(1:2)
1450              end if
1451            end if
1452            call wrtout(std_out,msg,'PERS')
1453          end if
1454        end do
1455      end do
1456      ABI_FREE(aa)
1457      ABI_FREE(bb)
1458    end if
1459 
1460 !  ----- End loop over spin components
1461  end do
1462 
1463  if (.not.present(one_over_rad2))  then
1464    ABI_FREE(one_over_rad2_)
1465  end if
1466 
1467  DBG_EXIT("COLL")
1468 
1469 end subroutine pawdensities

m_paw_denpot/pawkindensities [ Functions ]

[ Top ] [ m_paw_denpot ] [ Functions ]

NAME

 pawkindensities

FUNCTION

 Compute PAW on-site kinetic energy densities (all-electron, pseudo) for a given atom

INPUTS

  cplex: if 1, on-site densities are REAL, if 2, COMPLEX (response function only)
  lm_size=number of (l,m) moments
  lmselectin(lm_size)=flags selecting the non-zero LM-moments of on-site kinetic energy densities
                      (value of these flags at input; must be .TRUE. for nzlmopt/=1)
  nspden=number of spin-density components
  nzlmopt=if -1, compute all LM-moments of densities (lmselectin=.true. forced)
                 initialize "lmselectout" (index of non-zero LM-moments of densities)
          if  0, compute all LM-moments of densities (lmselectin=.true. forced)
                 force "lmselectout" to .true. (index of non-zero LM-moments of densities)
          if  1, compute only non-zero LM-moments of densities (stored before in "lmselectin")
  one_over_rad2(mesh_size)= contains 1/r**2 for each point of the radial grid -optional argument-
  opt_dens=flag controlling which on-site kinetic energy density(ies) is (are) computed
           0,1: all-electron and pseudo on-site kinetic energy densities
             2: only all-electron density
  opt_l=controls which l-moment(s) contribute to the kinetic energy density:
        <0 : all l contribute
        >=0: only l=opt_l contributes
        Note: opt_l>=0 is only compatible with opt_dens=2
  pawang <type(pawang_type)>=paw angular mesh and related data
  pawrad <type(pawrad_type)>=paw radial mesh and related data (for the current atom type)
  pawrhoij <type(pawrhoij_type)>= paw rhoij occupancies and related data (for the current atom)
  pawtab <type(pawtab_type)>=paw tabulated starting data (for the current atom type)

OUTPUT

  tau1(cplex*mesh_size,lm_size,nspden)= on site kinetic energy density
  ttau1(cplex*mesh_size,lm_size,nspden)]= pseudo on site kinetic energy density

SOURCE

1511 subroutine pawkindensities(cplex,lmselectin,lm_size,nspden,nzlmopt,&
1512 &          opt_dens,opt_l,pawang,pawrad,pawrhoij,pawtab,tau1,ttau1,&
1513 &          one_over_rad2) ! optional
1514 
1515 !Arguments ---------------------------------------------
1516 !scalars
1517  integer,intent(in) :: cplex,lm_size,nspden,nzlmopt,opt_dens,opt_l
1518  type(pawang_type),intent(in) :: pawang
1519  type(pawrad_type),intent(in) :: pawrad
1520  type(pawrhoij_type),intent(in) :: pawrhoij
1521  type(pawtab_type),intent(in) :: pawtab
1522 !arrays
1523  logical,intent(in) :: lmselectin(lm_size)
1524  real(dp),intent(in),target,optional :: one_over_rad2(pawtab%mesh_size)
1525  real(dp),intent(out),optional :: tau1(cplex*pawtab%mesh_size,lm_size,nspden*(1-(opt_dens/2)))
1526  real(dp),intent(out),optional :: ttau1(cplex*pawtab%mesh_size,lm_size,nspden*(1-(opt_dens/2)))
1527 !Local variables ---------------------------------------
1528 !scalars
1529  integer :: dplex,ii,iq0,ir,irhoij,isel,ispden,jrhoij
1530  integer :: ilmn,ilm,ilm1,iln,jlmn,jlm1,jln,klm,klmn,ll,lmax,lmin,mesh_size
1531  real(dp) :: phiphj,tphitphj
1532  character(len=500) :: msg
1533 !arrays
1534  real(dp) :: ro(cplex),ro_rg(cplex)
1535  real(dp),allocatable :: aa(:),bb(:)
1536  real(dp),pointer :: one_over_rad2_(:)
1537 
1538 ! *************************************************************************
1539 
1540  DBG_ENTER("COLL")
1541 
1542 !Compatibility tests
1543  if (nzlmopt/=-1) then
1544    msg='nzlmopt/=-1 has not not been tested (might be wrong)!'
1545    ABI_BUG(msg)
1546  end if
1547  if (opt_dens/=2.and.opt_l>=0) then
1548    msg='opt_dens/=2 incompatible with opt_l>=0!'
1549    ABI_BUG(msg)
1550  end if
1551  if(nzlmopt/=0.and.nzlmopt/=1.and.nzlmopt/=-1) then
1552    msg='invalid value for variable "nzlmopt"!'
1553    ABI_BUG(msg)
1554  end if
1555  if(nspden>pawrhoij%nspden) then
1556    msg='nspden must be <= pawrhoij%nspden!'
1557    ABI_BUG(msg)
1558  end if
1559  if (cplex>pawrhoij%qphase) then
1560    msg='cplex must be <= pawrhoij%qphase!'
1561    ABI_BUG(msg)
1562  end if
1563  if (nzlmopt/=1) then
1564    if (any(.not.lmselectin(1:lm_size))) then
1565      msg='With nzlmopt/=1, lmselectin must be true!'
1566      ABI_BUG(msg)
1567    end if
1568  end if
1569  if (pawang%gnt_option==0) then
1570    msg='pawang%gnt_option=0!'
1571    ABI_BUG(msg)
1572  end if
1573  if (pawang%nabgnt_option==0) then
1574    msg='pawang%nabgnt_option=0!'
1575    ABI_BUG(msg)
1576  end if
1577 
1578 !Various inits
1579  tau1=zero
1580  if (opt_dens<2) ttau1=zero
1581  mesh_size=pawtab%mesh_size;dplex=cplex-1
1582  iq0=pawrhoij%cplex_rhoij*pawrhoij%lmn2_size
1583  if (present(one_over_rad2)) then
1584    one_over_rad2_ => one_over_rad2
1585  else
1586    ABI_MALLOC(one_over_rad2_,(mesh_size))
1587    one_over_rad2_(1)=zero
1588    one_over_rad2_(2:mesh_size)=one/pawrad%rad(2:mesh_size)**2
1589  end if
1590 
1591 !=== Compute "on-site" kin. energy densities (n1, ntild1) =====
1592 !==============================================================
1593 
1594  do ispden=1,nspden
1595 
1596 !  -- Loop over ij channels (basis components)
1597    jrhoij=1
1598    do irhoij=1,pawrhoij%nrhoijsel
1599      klmn=pawrhoij%rhoijselect(irhoij)
1600      klm =pawtab%indklmn(1,klmn)
1601      lmin=pawtab%indklmn(3,klmn)
1602      lmax=pawtab%indklmn(4,klmn)
1603      ilmn=pawtab%indklmn(7,klmn) ! (l,m,n) orbital 1
1604      jlmn=pawtab%indklmn(8,klmn) ! (l,m,n) orbital 2
1605      ilm1=pawtab%indklmn(5,klmn) ! (l,m) orbital 1
1606      jlm1=pawtab%indklmn(6,klmn) ! (l,m) orbital 2
1607      iln=pawtab%indlmn(5,ilmn)   ! (l,n) orbital 1
1608      jln=pawtab%indlmn(5,jlmn)   ! (l,n) orbital 2
1609 
1610 !    Retrieve rhoij
1611      if (pawrhoij%nspden/=2) then
1612        ro(1)=pawrhoij%rhoijp(jrhoij,ispden)
1613        if (cplex==2) ro(2)=pawrhoij%rhoijp(iq0+jrhoij,ispden)
1614      else
1615        if (ispden==1) then
1616          ro(1)=pawrhoij%rhoijp(jrhoij,1)+pawrhoij%rhoijp(jrhoij,2)
1617          if (cplex==2) ro(2)=pawrhoij%rhoijp(iq0+jrhoij,1)+pawrhoij%rhoijp(iq0+jrhoij,2)
1618        else if (ispden==2) then
1619          ro(1)=pawrhoij%rhoijp(jrhoij,1)
1620          if (cplex==2) ro(2)=pawrhoij%rhoijp(iq0+jrhoij,1)
1621        end if
1622      end if
1623 !    Apply factor 1/2 (because tau=1/2 * Sum_ij[rhoij.Nabla_phi_i*Nabla_phi_j])
1624      ro(1:cplex)=half*pawtab%dltij(klmn)*ro(1:cplex)
1625 
1626 !    First option: AE and PS on-site kin. energy densities (opt_dens==0 or 1)
1627 !    ------------------------------------------------------------------------
1628      if (opt_dens==0.or.opt_dens==1) then
1629 
1630 !      Compute part of tau_lm depending on gaunt coefficients
1631        do ll=lmin,lmax,2
1632          do ilm=ll**2+1,(ll+1)**2
1633            if (lmselectin(ilm)) then
1634              isel=pawang%gntselect(ilm,klm)
1635              if (isel>0) then
1636                ro_rg(1:cplex)=ro(1:cplex)*pawang%realgnt(isel)
1637                do ir=2,mesh_size
1638                  phiphj=pawtab%nablaphi(ir,iln)*pawtab%nablaphi(ir,jln)
1639                  tphitphj=pawtab%tnablaphi(ir,iln)*pawtab%tnablaphi(ir,jln)
1640                  tau1(cplex*ir-dplex:ir*cplex,ilm,ispden)=tau1(cplex*ir-dplex:ir*cplex,ilm,ispden)&
1641 &                 +ro_rg(1:cplex)*phiphj*one_over_rad2_(ir)
1642                  ttau1(cplex*ir-dplex:ir*cplex,ilm,ispden)=ttau1(cplex*ir-dplex:ir*cplex,ilm,ispden)&
1643 &                 +ro_rg(1:cplex)*tphitphj*one_over_rad2_(ir)
1644                end do
1645              end if
1646            end if
1647         end do  ! End loops over ll,lm
1648        end do
1649 
1650 !      Compute the part of tau_lm depending on nablagaunt coefficients
1651        do ilm=1,(1+lmax)**2
1652          if (lmselectin(ilm)) then
1653            isel=pawang%nablagntselect(ilm,ilm1,jlm1)
1654            if (isel>0) then
1655              ro_rg(1:cplex)=ro(1:cplex)*pawang%nablarealgnt(isel)
1656              do ir=2,mesh_size
1657                phiphj=pawtab%phi(ir,iln)*pawtab%phi(ir,jln)
1658                tphitphj=pawtab%tphi(ir,iln)*pawtab%tphi(ir,jln)
1659                tau1(cplex*ir-dplex:ir*cplex,ilm,ispden)=tau1(cplex*ir-dplex:ir*cplex,ilm,ispden)&
1660 &               +ro_rg(1:cplex)*phiphj*one_over_rad2_(ir)**2
1661                ttau1(cplex*ir-dplex:ir*cplex,ilm,ispden)=ttau1(cplex*ir-dplex:ir*cplex,ilm,ispden)&
1662 &               +ro_rg(1:cplex)*tphitphj*one_over_rad2_(ir)**2
1663               end do
1664            end if
1665          end if
1666        end do
1667 
1668 !    2nd option: AE on-site kinetic energy density only (opt_dens==2)
1669 !    ----------------------------------------------------------------
1670      else if (opt_dens==2) then
1671 
1672 !      Compute part of tau_lm depending on gaunt coefficients
1673        do ll=lmin,lmax,2
1674          do ilm=ll**2+1,(ll+1)**2
1675            if (lmselectin(ilm)) then
1676              isel=pawang%gntselect(ilm,klm)
1677              if (isel>0) then
1678                ro_rg(1:cplex)=ro(1:cplex)*pawang%realgnt(isel)
1679                do ir=2,mesh_size
1680                  phiphj=pawtab%nablaphi(ir,iln)*pawtab%nablaphi(ir,jln)
1681                  tau1(cplex*ir-dplex:ir*cplex,ilm,ispden)=tau1(cplex*ir-dplex:ir*cplex,ilm,ispden)&
1682     &             +ro_rg(1:cplex)*phiphj*one_over_rad2_(ir)
1683                end do
1684              end if
1685            end if
1686         end do  ! End loops over ll,lm
1687        end do
1688 
1689 !      Compute the part of tau_lm depending on nablagaunt coefficients
1690        do ilm=1,(1+lmax)**2
1691          if (lmselectin(ilm)) then
1692            isel=pawang%nablagntselect(ilm,ilm1,jlm1)
1693            if (isel>0) then
1694              ro_rg(1:cplex)=ro(1:cplex)*pawang%nablarealgnt(isel)
1695              do ir=2,mesh_size
1696                phiphj=pawtab%phi(ir,iln)*pawtab%phi(ir,jln)
1697                tau1(cplex*ir-dplex:ir*cplex,ilm,ispden)=tau1(cplex*ir-dplex:ir*cplex,ilm,ispden)&
1698 &               +ro_rg(1:cplex)*phiphj*one_over_rad2_(ir)**2
1699               end do
1700            end if
1701          end if
1702        end do
1703 
1704      end if
1705 
1706 !    -- End loop over ij channels
1707      jrhoij=jrhoij+pawrhoij%cplex_rhoij
1708    end do
1709 
1710 !  Compute tau1(r=0) and ttau1(r=0)
1711    if (cplex==2) then
1712      ABI_MALLOC(aa,(5))
1713      ABI_MALLOC(bb,(5))
1714    end if
1715    if (opt_dens==0.or.opt_dens==1) then
1716      do ll=0,pawtab%lcut_size-1
1717        do ilm=ll**2+1,(ll+1)**2
1718          if (lmselectin(ilm)) then
1719            if (cplex==1) then
1720              call pawrad_deducer0(tau1 (:,ilm,ispden),mesh_size,pawrad)
1721              call pawrad_deducer0(ttau1(:,ilm,ispden),mesh_size,pawrad)
1722            else
1723              do ii=0,1
1724                do ir=2,5
1725                  aa(ir)=tau1 (2*ir-ii,ilm,ispden)
1726                  bb(ir)=ttau1(2*ir-ii,ilm,ispden)
1727                end do
1728                call pawrad_deducer0(aa,5,pawrad)
1729                call pawrad_deducer0(bb,5,pawrad)
1730                tau1 (2-ii,ilm,ispden)=aa(1)
1731                ttau1(2-ii,ilm,ispden)=bb(1)
1732              end do
1733            end if
1734          end if
1735        end do
1736      end do
1737    else if (opt_dens==2) then
1738      do ll=0,pawtab%lcut_size-1
1739        do ilm=ll**2+1,(ll+1)**2
1740          if (lmselectin(ilm)) then
1741            if (cplex==1) then
1742              call pawrad_deducer0(tau1(:,ilm,ispden),mesh_size,pawrad)
1743            else
1744              do ii=0,1
1745                do ir=2,5
1746                  aa(ir)=tau1(2*ir-ii,ilm,ispden)
1747                end do
1748                call pawrad_deducer0(aa,5,pawrad)
1749                tau1(2-ii,ilm,ispden)=aa(1)
1750              end do
1751            end if
1752          end if
1753        end do
1754      end do
1755    end if
1756    if (cplex==2)  then
1757      ABI_FREE(aa)
1758      ABI_FREE(bb)
1759    end if
1760 
1761 !  ----- End loop over spin components
1762  end do
1763 
1764  if (.not.present(one_over_rad2))  then
1765    ABI_FREE(one_over_rad2_)
1766  end if
1767 
1768  DBG_EXIT("COLL")
1769 
1770 end subroutine pawkindensities