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

17 #if defined HAVE_CONFIG_H
18 #include "config.h"
19 #endif
20 
21 #include "abi_common.h"
22 
23 MODULE m_paw_denpot
24 
25  use defs_basis
26  use m_abicore
27  use m_errors
28  use m_xmpi
29  use m_time, only : timab
30 
31  use m_pawang,           only : pawang_type
32  use m_pawrad,           only : pawrad_type,pawrad_deducer0,poisson,simp_gen
33  use m_pawtab,           only : pawtab_type
34  use m_paw_an,           only : paw_an_type
35  use m_paw_ij,           only : paw_ij_type
36  use m_pawfgrtab,        only : pawfgrtab_type
37  use m_pawrhoij,         only : pawrhoij_type
38  use m_pawdij,           only : pawdijhartree,pawdiju_euijkl,pawdijnd,pawdijso,pawxpot,pawdijfock,symdij,symdij_all
39  use m_pawxc,            only : pawxc,pawxc_dfpt,pawxcm,pawxcm_dfpt,pawxcpositron,pawxcmpositron
40  use m_paw_finegrid,     only : pawgylm
41  use m_paral_atom,       only : get_my_atmtab,free_my_atmtab
42  use m_paw_correlations, only : pawuenergy,pawxenergy,setnoccmmp
43  use m_paral_atom,       only : get_my_atmtab,free_my_atmtab
44 
45  use m_crystal,          only : crystal_t
46  use m_electronpositron, only : electronpositron_type,electronpositron_calctype
47 
48  implicit none
49 
50  private
51 
52 !public procedures.
53  public :: pawdenpot    ! Compute different (PAW) energies, densities and potentials inside PAW spheres
54  public :: pawdensities ! Compute PAW on-site densities (all-electron ,pseudo and compensation)
55  public :: paw_mknewh0  ! Compute bare PAW on-site Hamiltonian (-> GW calculations)
56 
57 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

PARENTS

      calc_vhxc_me

CHILDREN

      free_my_atmtab,get_my_atmtab,pawgylm,symdij,symdij_all,wrtout

SOURCE

1575 subroutine paw_mknewh0(my_natom,nsppol,nspden,nfftf,pawspnorb,pawprtvol,Cryst,&
1576 &          Pawtab,Paw_an,Paw_ij,Pawang,Pawfgrtab,vxc,vxc_val,vtrial,&
1577 &          mpi_atmtab,comm_atom) ! optional arguments (parallelism)
1578 
1579 
1580 !This section has been created automatically by the script Abilint (TD).
1581 !Do not modify the following lines by hand.
1582 #undef ABI_FUNC
1583 #define ABI_FUNC 'paw_mknewh0'
1584 !End of the abilint section
1585 
1586  implicit none
1587 
1588 !Arguments ------------------------------------
1589 !scalars
1590  integer,intent(in) :: my_natom,nsppol,nspden,nfftf,pawprtvol,pawspnorb
1591  integer,optional,intent(in) :: comm_atom
1592 !arrays
1593  integer,optional,target,intent(in) :: mpi_atmtab(:)
1594  real(dp),intent(in) :: vxc(nfftf,nspden),vxc_val(nfftf,nspden),vtrial(nfftf,nspden)
1595  type(crystal_t),intent(in) :: Cryst
1596  type(Pawang_type),intent(in) :: Pawang
1597  type(Pawtab_type),target,intent(in) :: Pawtab(Cryst%ntypat)
1598  type(Paw_an_type),intent(in) :: Paw_an(my_natom)
1599  type(Paw_ij_type),intent(inout) :: Paw_ij(my_natom)
1600  type(Pawfgrtab_type),intent(inout) :: Pawfgrtab(my_natom)
1601 
1602 !Local variables-------------------------------
1603 !scalars
1604  integer,parameter :: ipert0=0
1605  integer :: iat,iat_tot,idij,cplex_rf,ndij,option_dij
1606  integer :: itypat,lmn_size,j0lmn,jlmn,ilmn,klmn,klmn1,klm
1607  integer :: lmin,lmax,mm,isel,lm_size,lmn2_size,my_comm_atom,cplex_dij
1608  integer :: ils,ilslm,ic,lm0
1609  integer :: nsploop,is2fft
1610  real(dp) :: gylm,qijl
1611  logical :: ltest,my_atmtab_allocated,paral_atom
1612  character(len=500) :: msg
1613 !arrays
1614  integer, pointer :: indklmn_(:,:)
1615  integer,pointer :: my_atmtab(:)
1616  real(dp) :: rdum(1),rdum2(1)
1617  real(dp),allocatable :: prod_hloc(:,:),prodhxc_core(:,:)
1618  real(dp),allocatable :: dijhl_hat(:,:),dijhmxc_val(:,:)
1619 
1620 ! *************************************************************************
1621 
1622  DBG_ENTER("COLL")
1623 
1624  call wrtout(std_out,'Assembling PAW strengths for the bare Hamiltonian','COLL')
1625 
1626 !== Set up parallelism over atoms ===
1627  paral_atom=(present(comm_atom).and.(my_natom/=Cryst%natom))
1628  nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab
1629  my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom
1630  call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,Cryst%natom,my_natom_ref=my_natom)
1631 
1632  if (my_natom>0) then
1633 
1634 !  === Test if required pointers in paw_ij are allocated ===
1635    ltest = (allocated(Paw_ij(1)%dijxc).and.allocated(Paw_ij(1)%dijxc_val) )
1636    ABI_CHECK(ltest,'dijxc or dijxc_val not calculated')
1637 
1638    ltest=(allocated(Paw_ij(1)%dijhat)) !.and.Paw_ij(1)%has_dijhat==2)
1639    ABI_CHECK(ltest,'dijhat not calculated')
1640 
1641    ltest=(allocated(Paw_ij(1)%dijhartree)) !.and.Paw_ij(1)%has_dijhartree==2)
1642    ABI_CHECK(ltest,'dijhartree not calculated')
1643 
1644    if (ANY(Pawtab(:)%usepawu>0)) then
1645      do iat=1,my_natom
1646        iat_tot=iat;if (paral_atom) iat_tot=my_atmtab(iat)
1647        itypat=Cryst%typat(iat_tot)
1648        if (Pawtab(itypat)%usepawu>0) then
1649          ltest=(allocated(Paw_ij(iat)%dijU) ) !.and.Paw_ij(iat)%has_dijU==2)
1650          write(msg,'(a,i3,a)')" For atom no. ",iat," %dijU(iat) has not been calculated."
1651          ABI_CHECK(ltest,msg)
1652        end if
1653      end do
1654    end if
1655 
1656    if (pawspnorb>0) then
1657      do iat=1,my_natom
1658        ltest=(allocated(Paw_ij(iat)%dijso) ) !.and.Paw_ij(iat)%has_dijso==2)
1659        write(msg,'(a,i3,a)')" For atom no. ",iat," %dijso(iat) has not been calculated."
1660        ABI_CHECK(ltest,msg)
1661      end do
1662    end if
1663  end if ! my_natom>0
1664 
1665 !== Construct the new PAW H0 Hamiltonian ===
1666  do iat=1,my_natom
1667    iat_tot=iat;if (paral_atom) iat_tot=my_atmtab(iat)
1668 
1669    itypat    = Cryst%typat(iat_tot)
1670    lmn_size  = Pawtab(itypat)%lmn_size
1671    lmn2_size = Pawtab(itypat)%lmn2_size
1672    lm_size   = Paw_an(iat)%lm_size
1673    cplex_rf  = Paw_ij(iat)%cplex_rf
1674    cplex_dij = Paw_ij(iat)%cplex_dij
1675    ndij      = Paw_ij(iat)%ndij
1676 
1677    ABI_CHECK(cplex_rf==1,'cplex_rf/=1 not implemented')
1678    ABI_CHECK(cplex_dij==1,'cplex_dij/=1 not implemented')
1679 !
1680 !  Eventually compute g_l(r).Y_lm(r) factors for the current atom (if not already done)
1681    if (Pawfgrtab(iat)%gylm_allocated==0) then
1682      if (allocated(Pawfgrtab(iat)%gylm))  then
1683        ABI_DEALLOCATE(Pawfgrtab(iat)%gylm)
1684      end if
1685      ABI_ALLOCATE(Pawfgrtab(iat)%gylm,(Pawfgrtab(iat)%nfgd,lm_size))
1686      Pawfgrtab(iat)%gylm_allocated=2
1687 
1688      call pawgylm(Pawfgrtab(iat)%gylm,rdum,rdum2,lm_size,&
1689 &     Pawfgrtab(iat)%nfgd,1,0,0,Pawtab(itypat),Pawfgrtab(iat)%rfgd)
1690    end if
1691 
1692 !  === Calculate LM contribution to dijhmxc_val for this atom ===
1693 !  * Dijxc contains also the Hat term on the FFT mesh while Dijxc_val does not
1694 !  contain neither the hat term nor the LM sum of onsite terms (they should cancel each other)
1695 !  FIXME change paw_dij,  otherwise I miss tnc in vxc
1696 !  * prodhxc_core is used to assemble $\int g_l Ylm (vtrial - vxc_val[tn+nhat] dr$ on the FFT mesh ===
1697 !  * The following quantities do not depend on ij
1698    ABI_ALLOCATE(prod_hloc   ,(lm_size,ndij))
1699    ABI_ALLOCATE(prodhxc_core,(lm_size,ndij))
1700    prod_hloc   =zero
1701    prodhxc_core=zero
1702    do idij=1,ndij
1703      do ilslm=1,lm_size
1704        do ic=1,Pawfgrtab(iat)%nfgd
1705          is2fft=Pawfgrtab(iat)%ifftsph(ic)
1706          gylm=Pawfgrtab(iat)%gylm(ic,ilslm)
1707          prod_hloc (ilslm,idij)=prod_hloc (ilslm,idij) + (vtrial(is2fft,idij)-vxc(is2fft,idij))*gylm
1708 !        prodhxc_core(ilslm,idij)=prodhxc_core(ilslm,idij) + (vxc_val(is2fft,idij))*gylm
1709          prodhxc_core(ilslm,idij)=prodhxc_core(ilslm,idij) + (vtrial(is2fft,idij)-vxc_val(is2fft,idij))*gylm
1710        end do
1711      end do
1712    end do !idij
1713 
1714 !  === Assembly the "Hat" contribution for this atom ====
1715    ABI_ALLOCATE(dijhl_hat  ,(cplex_dij*lmn2_size,ndij))
1716    ABI_ALLOCATE(dijhmxc_val,(cplex_dij*lmn2_size,ndij))
1717    dijhl_hat  =zero
1718    dijhmxc_val=zero
1719    indklmn_ => Pawtab(itypat)%indklmn(1:6,1:lmn2_size)
1720 
1721    do idij=1,ndij
1722      do klmn=1,lmn2_size
1723        klm =indklmn_(1,klmn)
1724        lmin=indklmn_(3,klmn)
1725        lmax=indklmn_(4,klmn)
1726 
1727 !      === $\sum_lm q_ij^l prod* for each idij$ ===
1728        do ils=lmin,lmax,2
1729          lm0=ils**2+ils+1
1730          do mm=-ils,ils
1731            ilslm=lm0+mm
1732            isel=Pawang%gntselect(lm0+mm,klm)
1733            if (isel>0) then
1734              qijl=Pawtab(itypat)%qijl(ilslm,klmn)
1735              dijhl_hat  (klmn,idij)=dijhl_hat  (klmn,idij) +  prod_hloc (ilslm,idij)*qijl
1736              dijhmxc_val(klmn,idij)=dijhmxc_val(klmn,idij) +prodhxc_core(ilslm,idij)*qijl
1737            end if
1738          end do
1739        end do
1740      end do
1741    end do
1742 
1743    ABI_DEALLOCATE(prod_hloc)
1744    ABI_DEALLOCATE(prodhxc_core)
1745 
1746 !  * Normalization factor due to integration on the FFT mesh
1747    dijhl_hat  = dijhl_hat  *Cryst%ucvol/DBLE(nfftf)
1748    dijhmxc_val= dijhmxc_val*Cryst%ucvol/DBLE(nfftf)
1749 
1750 !  === Now assembly the bare Hamiltonian ===
1751 !  * Loop over density components overwriting %dij
1752    nsploop=nsppol; if (Paw_ij(iat)%ndij==4) nsploop=4
1753 
1754    do idij=1,nsploop
1755      klmn1=1
1756 
1757      do jlmn=1,lmn_size
1758        j0lmn=jlmn*(jlmn-1)/2
1759        do ilmn=1,jlmn
1760          klmn=j0lmn+ilmn
1761 
1762 !        The following gives back the input dij.
1763 !        since dijxc contains the hat term done on the FFT mesh
1764          if (.FALSE.) then
1765            Paw_ij(iat)%dij(klmn,idij) =        &
1766 &           Pawtab(itypat)%dij0    (klmn)      &
1767 &           +Paw_ij(iat)%dijhartree(klmn)      &
1768 &           +Paw_ij(iat)%dijxc     (klmn,idij) &
1769 &           +dijhl_hat   (klmn,idij)
1770 
1771          else
1772 !          === Make nonlocal part of h0 removing the valence contribution ===
1773 !          Remeber that XC contains already the Hat contribution
1774            Paw_ij(iat)%dij(klmn,idij) =        &
1775 &           Pawtab(itypat)%dij0      (klmn)    &
1776 &           +Paw_ij(iat)%dijhartree(klmn)      &
1777 &           +Paw_ij(iat)%dijxc     (klmn,idij) &  ! 2 lines to get the d1-dt1 XC core contribution + XC hat (core+val)
1778 &          -Paw_ij(iat)%dijxc_val (klmn,idij) &  ! I suppose that the "hat" term on the FFT mesh in included in both.
1779 &          +dijhmxc_val(klmn,idij)               ! Local + Hartree - XC val contribution to the "hat" term.
1780 
1781 !          Add the U contribution to the
1782 !          if (.FALSE. .and. Pawtab(itypat)%usepawu>0) then
1783            if (.TRUE. .and. Pawtab(itypat)%usepawu>0) then
1784              Paw_ij(iat)%dij(klmn,idij) = Paw_ij(iat)%dij(klmn,idij) + Paw_ij(iat)%dijU(klmn,idij)
1785            end if
1786          end if
1787 !        TODO dijso, dijU, vpawx?
1788 !        Just to be consistent, update some values.
1789 !$Paw_ij(iat)%dijhat(klmn,idij)=Paw_ij(iat)%dijhat(klmn,idij)-dijhmxc_val(klmn,idij)
1790 
1791        end do !ilmn
1792      end do !jlmn
1793    end do !idij
1794 
1795 !  this is to be consistent?
1796 !  deallocate(Paw_ij(iat)%dijvxc_val)
1797    ABI_DEALLOCATE(dijhl_hat)
1798    ABI_DEALLOCATE(dijhmxc_val)
1799  end do !iat
1800 
1801 !=== Symmetrize total Dij ===
1802  option_dij=0 ! For total Dij.
1803 #if 0
1804  if (paral_atom) then
1805    call symdij(Cryst%gprimd,Cryst%indsym,ipert0,my_natom,Cryst%natom,Cryst%nsym,Cryst%ntypat,option_dij,&
1806 &   Paw_ij,Pawang,pawprtvol,Pawtab,Cryst%rprimd,Cryst%symafm,Cryst%symrec,&
1807 &   comm_atom=my_comm_atom,mpi_atmtab=my_atmtab)
1808  else
1809    call symdij(Cryst%gprimd,,Cryst%indsym,ipert0,my_natom,Cryst%natom,Cryst%nsym,Cryst%ntypat,option_dij,&
1810 &   Paw_ij,Pawang,pawprtvol,Pawtab,Cryst%rprimd,Cryst%symafm,Cryst%symrec)
1811  end if
1812 #else
1813  if (paral_atom) then
1814    call symdij_all(Cryst%gprimd,Cryst%indsym,ipert0,my_natom,Cryst%natom,Cryst%nsym,Cryst%ntypat,&
1815 &   Paw_ij,Pawang,pawprtvol,Pawtab,Cryst%rprimd,Cryst%symafm,Cryst%symrec,&
1816 &   comm_atom=my_comm_atom,mpi_atmtab=my_atmtab)
1817  else
1818    call symdij_all(Cryst%gprimd,Cryst%indsym,ipert0,my_natom,Cryst%natom,Cryst%nsym,Cryst%ntypat,&
1819 &   Paw_ij,Pawang,pawprtvol,Pawtab,Cryst%rprimd,Cryst%symafm,Cryst%symrec)
1820  end if
1821 #endif
1822 
1823 !Destroy atom table used for parallelism
1824  call free_my_atmtab(my_atmtab,my_atmtab_allocated)
1825 
1826  DBG_EXIT("COLL")
1827 
1828 end subroutine paw_mknewh0

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,my_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 allowed 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(cplex*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[m](cplex*mesh_size,:,nspden)=XC potential calculated from "on-site" density
    paw_an(my_natom)%vxct1[m](cplex*mesh_size,:,nspden)=XC potential calculated from "on-site" pseudo density
    ==== if paw_an(iatom_tot)%has_vxcval==1 compute also XC potentials neglecting core charge
      paw_an(my_natom)%vxc1_val[m](cplex*mesh_size,:nspden)=XC potential calculated from spherical valence density
      paw_an(my_natom)%vxct1_val[m](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(cplex*cplex_dij*lmn2_size,nspden)=spin-orbit contribution to dij

NOTES

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

PARENTS

      bethe_salpeter,dfpt_scfcv,odamix,paw_qpscgw,respfn,scfcv,screening
      sigma

CHILDREN

      free_my_atmtab,get_my_atmtab,pawdensities,pawdijfock,pawdijhartree
      pawdijnd,pawdijso,pawrad_deducer0,pawuenergy,pawxc,pawxc_dfpt,pawxcm
      pawxcm_dfpt,pawxcmpositron,pawxcpositron,pawxenergy,pawxpot,poisson
      setnoccmmp,timab,wrtout,xmpi_sum

SOURCE

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

PARENTS

      make_efg_onsite,pawdenpot,pawdfptenergy,poslifetime,posratecore

CHILDREN

      pawrad_deducer0,simp_gen,wrtout

SOURCE

1178 subroutine pawdensities(compch_sph,cplex,iatom,lmselectin,lmselectout,lm_size,nhat1,nspden,nzlmopt,&
1179 &          opt_compch,opt_dens,opt_l,opt_print,pawang,pawprtvol,pawrad,pawrhoij,pawtab,rho1,trho1,&
1180 &          one_over_rad2) ! optional
1181 
1182 
1183 !This section has been created automatically by the script Abilint (TD).
1184 !Do not modify the following lines by hand.
1185 #undef ABI_FUNC
1186 #define ABI_FUNC 'pawdensities'
1187 !End of the abilint section
1188 
1189  implicit none
1190 
1191 !Arguments ---------------------------------------------
1192 !scalars
1193  integer,intent(in) :: cplex,iatom,lm_size,nspden,nzlmopt,opt_compch,opt_dens,opt_l,opt_print,pawprtvol
1194 ! jmb  real(dp),intent(out) :: compch_sph
1195  real(dp),intent(inout) :: compch_sph
1196  type(pawang_type),intent(in) :: pawang
1197  type(pawrad_type),intent(in) :: pawrad
1198  type(pawrhoij_type),intent(in) :: pawrhoij
1199  type(pawtab_type),intent(in) :: pawtab
1200 !arrays
1201  logical,intent(in) :: lmselectin(lm_size)
1202  logical,intent(inout) :: lmselectout(lm_size)
1203  real(dp),intent(in),target,optional :: one_over_rad2(pawtab%mesh_size)
1204  real(dp),intent(out) :: nhat1(cplex*pawtab%mesh_size,lm_size,nspden*(1-((opt_dens+1)/2)))
1205  real(dp),intent(out) ::  rho1(cplex*pawtab%mesh_size,lm_size,nspden)
1206  real(dp),intent(out) :: trho1(cplex*pawtab%mesh_size,lm_size,nspden*(1-(opt_dens/2)))
1207 
1208 !Local variables ---------------------------------------
1209 !scalars
1210  integer :: dplex,ii,ilm,iplex,ir,irhoij,isel,ispden,jrhoij
1211  integer :: klm,klmn,kln,ll,lmax,lmin,mesh_size
1212  real(dp) :: m1,mt1,rdum
1213  character(len=500) :: msg
1214 !arrays
1215  real(dp) :: compchspha(cplex),compchsphb(cplex),ro(cplex),ro_ql(cplex),ro_rg(cplex)
1216  real(dp),allocatable :: aa(:),bb(:)
1217  real(dp),pointer :: one_over_rad2_(:)
1218 
1219 ! *************************************************************************
1220 
1221  DBG_ENTER("COLL")
1222 
1223 !Compatibility tests
1224  if (opt_dens/=2.and.opt_l>=0) then
1225    msg='  opt_dens/=2 incompatible with opt_l>=0 !'
1226    MSG_BUG(msg)
1227  end if
1228  if(nzlmopt/=0.and.nzlmopt/=1.and.nzlmopt/=-1) then
1229    msg='  invalid value for variable "nzlmopt".'
1230    MSG_BUG(msg)
1231  end if
1232  if(nspden>pawrhoij%nspden) then
1233    msg='  nspden must be <= pawrhoij%nspden !'
1234    MSG_BUG(msg)
1235  end if
1236  if (cplex>pawrhoij%cplex) then
1237    msg='  cplex must be <= pawrhoij%cplex !'
1238    MSG_BUG(msg)
1239  end if
1240  if (nzlmopt/=1) then
1241    if (any(.not.lmselectin(1:lm_size))) then
1242      msg='  With nzlmopt/=1, lmselectin must be true !'
1243      MSG_BUG(msg)
1244    end if
1245  end if
1246  if (pawang%gnt_option==0) then
1247    msg='  pawang%gnt_option=0 !'
1248    MSG_BUG(msg)
1249  end if
1250 
1251 !Various inits
1252  rho1=zero
1253  if (opt_dens <2) trho1=zero
1254  if (opt_dens==0) nhat1=zero
1255  mesh_size=pawtab%mesh_size;dplex=cplex-1
1256  if (nzlmopt<1) lmselectout(1:lm_size)=.true.
1257  if (present(one_over_rad2)) then
1258    one_over_rad2_ => one_over_rad2
1259  else
1260    ABI_ALLOCATE(one_over_rad2_,(mesh_size))
1261    one_over_rad2_(1)=zero
1262    one_over_rad2_(2:mesh_size)=one/pawrad%rad(2:mesh_size)**2
1263  end if
1264 
1265 !===== Compute "on-site" densities (n1, ntild1, nhat1) =====
1266 !==========================================================
1267 
1268  do ispden=1,nspden
1269 
1270 !  -- Loop over ij channels (basis components)
1271    jrhoij=1
1272    do irhoij=1,pawrhoij%nrhoijsel
1273      klmn=pawrhoij%rhoijselect(irhoij)
1274      klm =pawtab%indklmn(1,klmn)
1275      kln =pawtab%indklmn(2,klmn)
1276      lmin=pawtab%indklmn(3,klmn)
1277      lmax=pawtab%indklmn(4,klmn)
1278 
1279 
1280 !    Retrieve rhoij
1281      if (pawrhoij%nspden/=2) then
1282        ro(1:cplex)=pawrhoij%rhoijp(jrhoij:jrhoij+dplex,ispden)
1283      else
1284        if (ispden==1) then
1285          ro(1:cplex)=pawrhoij%rhoijp(jrhoij:jrhoij+dplex,1)&
1286 &                   +pawrhoij%rhoijp(jrhoij:jrhoij+dplex,2)
1287        else if (ispden==2) then
1288          ro(1:cplex)=pawrhoij%rhoijp(jrhoij:jrhoij+dplex,1)
1289        end if
1290      end if
1291      ro(1:cplex)=pawtab%dltij(klmn)*ro(1:cplex)
1292 
1293 !    First option: all on-site densities are computed (opt_dens==0)
1294 !    --------------------------------------------------------------
1295      if (opt_dens==0) then
1296        do ll=lmin,lmax,2
1297          do ilm=ll**2+1,(ll+1)**2
1298            if (lmselectin(ilm)) then
1299              isel=pawang%gntselect(ilm,klm)
1300              if (isel>0) then
1301                ro_ql(1:cplex)=ro(1:cplex)*pawtab%qijl(ilm,klmn)
1302                ro_rg(1:cplex)=ro(1:cplex)*pawang%realgnt(isel)
1303 !              == nhat1(r=0)
1304                nhat1(1:cplex,ilm,ispden)=nhat1(1:cplex,ilm,ispden) &
1305 &               +ro_ql(1:cplex)*pawtab%shapefunc(1,ll+1)
1306 !              == rho1(r>0), trho1(r>0), nhat1(r>0)
1307                do ir=2,mesh_size
1308                  rho1(cplex*ir-dplex:ir*cplex,ilm,ispden) =rho1(cplex*ir-dplex:ir*cplex,ilm,ispden)&
1309 &                 +ro_rg(1:cplex)*pawtab%phiphj  (ir,kln)*one_over_rad2_(ir)
1310                  trho1(cplex*ir-dplex:ir*cplex,ilm,ispden)=trho1(cplex*ir-dplex:ir*cplex,ilm,ispden)&
1311 &                 +ro_rg(1:cplex)*pawtab%tphitphj(ir,kln)*one_over_rad2_(ir)
1312                  nhat1(cplex*ir-dplex:ir*cplex,ilm,ispden)=nhat1(cplex*ir-dplex:ir*cplex,ilm,ispden)&
1313 &                 +ro_ql(1:cplex)*pawtab%shapefunc(ir,ll+1)
1314                end do
1315              end if
1316            end if
1317          end do  ! End loops over ll,lm
1318        end do
1319 
1320 !      2nd option: AE and pseudo densities are computed (opt_dens==1)
1321 !      --------------------------------------------------------------
1322      else if (opt_dens==1) then
1323        do ll=lmin,lmax,2
1324          do ilm=ll**2+1,(ll+1)**2
1325            if (lmselectin(ilm)) then
1326              isel=pawang%gntselect(ilm,klm)
1327              if (isel>0) then
1328                ro_rg(1:cplex)=ro(1:cplex)*pawang%realgnt(isel)
1329 !              == rho1(r>0), trho1(r>0)
1330                do ir=2,mesh_size
1331                  rho1(cplex*ir-dplex:ir*cplex,ilm,ispden) =rho1(cplex*ir-dplex:ir*cplex,ilm,ispden)&
1332 &                 +ro_rg(1:cplex)*pawtab%phiphj  (ir,kln)*one_over_rad2_(ir)
1333                  trho1(cplex*ir-dplex:ir*cplex,ilm,ispden)=trho1(cplex*ir-dplex:ir*cplex,ilm,ispden)&
1334 &                 +ro_rg(1:cplex)*pawtab%tphitphj(ir,kln)*one_over_rad2_(ir)
1335                end do
1336              end if
1337            end if
1338          end do  ! End loops over ll,lm
1339        end do
1340 
1341 !      3rd option: only all-electron on-site density is computed (opt_dens==2)
1342 !      -----------------------------------------------------------------------
1343      else if (opt_dens==2) then
1344        if (opt_l<0.or.(pawtab%indklmn(3,klmn)==0.and.pawtab%indklmn(4,klmn)==2*opt_l)) then
1345          do ll=lmin,lmax,2
1346            do ilm=ll**2+1,(ll+1)**2
1347              if (lmselectin(ilm)) then
1348                isel=pawang%gntselect(ilm,klm)
1349                if (isel>0) then
1350                  ro_rg(1:cplex)=ro(1:cplex)*pawang%realgnt(isel)
1351 !                == rho1(r>0)
1352                  do ir=2,mesh_size
1353                    rho1(cplex*ir-dplex:ir*cplex,ilm,ispden) =rho1(cplex*ir-dplex:ir*cplex,ilm,ispden)&
1354 &                   +ro_rg(1:cplex)*pawtab%phiphj(ir,kln)*one_over_rad2_(ir)
1355                  end do
1356                end if
1357              end if
1358            end do  ! End loops over ll, lm
1359          end do
1360        end if
1361      end if
1362 
1363 !    -- End loop over ij channels
1364      jrhoij=jrhoij+pawrhoij%cplex
1365    end do
1366 
1367 !  Scale densities with 1/r**2 and compute rho1(r=0) and trho1(r=0)
1368    if (cplex==2)  then
1369      ABI_ALLOCATE(aa,(5))
1370      ABI_ALLOCATE(bb,(5))
1371    end if
1372    if (opt_dens==0.or.opt_dens==1) then
1373      do ll=0,pawtab%lcut_size-1
1374        do ilm=ll**2+1,(ll+1)**2
1375          if (lmselectin(ilm)) then
1376            if (cplex==1) then
1377              call pawrad_deducer0(rho1 (:,ilm,ispden),mesh_size,pawrad)
1378              call pawrad_deducer0(trho1(:,ilm,ispden),mesh_size,pawrad)
1379            else
1380              do ii=0,1
1381                do ir=2,5
1382                  aa(ir)=rho1 (2*ir-ii,ilm,ispden)
1383                  bb(ir)=trho1(2*ir-ii,ilm,ispden)
1384                end do
1385                call pawrad_deducer0(aa,5,pawrad)
1386                call pawrad_deducer0(bb,5,pawrad)
1387                rho1 (2-ii,ilm,ispden)=aa(1)
1388                trho1(2-ii,ilm,ispden)=bb(1)
1389              end do
1390            end if
1391          end if
1392        end do
1393      end do
1394    else
1395      do ll=0,pawtab%lcut_size-1
1396        do ilm=ll**2+1,(ll+1)**2
1397          if (lmselectin(ilm)) then
1398            if (cplex==1) then
1399              call pawrad_deducer0(rho1(:,ilm,ispden),mesh_size,pawrad)
1400            else
1401              do ii=0,1
1402                do ir=2,5
1403                  aa(ir)=rho1 (2*ir-ii,ilm,ispden)
1404                end do
1405                call pawrad_deducer0(aa,5,pawrad)
1406                rho1(2-ii,ilm,ispden)=aa(1)
1407              end do
1408            end if
1409          end if
1410        end do
1411      end do
1412    end if
1413    if (cplex==2)  then
1414      ABI_DEALLOCATE(aa)
1415      ABI_DEALLOCATE(bb)
1416    end if
1417 
1418 !  -- Test moments of densities and store non-zero ones
1419    if (nzlmopt==-1) then
1420      do ll=0,pawtab%lcut_size-1
1421        do ilm=ll**2+1,(ll+1)**2
1422          m1=zero;mt1=zero
1423          if (cplex==1) then
1424            m1=maxval(abs(rho1 (1:mesh_size,ilm,ispden)))
1425            if (opt_dens<2) mt1=maxval(abs(trho1(1:mesh_size,ilm,ispden)))
1426          else
1427            do ir=1,mesh_size
1428              rdum=sqrt(rho1(2*ir-1,ilm,ispden)**2+rho1(2*ir,ilm,ispden)**2)
1429              m1=max(m1,rdum)
1430            end do
1431            if (opt_dens<2) then
1432              do ir=1,mesh_size
1433                rdum=sqrt(trho1(2*ir-1,ilm,ispden)**2+trho1(2*ir,ilm,ispden)**2)
1434                mt1=max(mt1,rdum)
1435              end do
1436            end if
1437          end if
1438          if (ispden==1) then
1439            if ((ilm>1).and.(m1<tol16).and.(mt1<tol16)) then
1440              lmselectout(ilm)=.false.
1441            end if
1442          else if (.not.(lmselectout(ilm))) then
1443            lmselectout(ilm)=((m1>=tol16).or.(mt1>=tol16))
1444          end if
1445        end do
1446      end do
1447    end if
1448 
1449 !  -- Compute integral of (n1-tn1) inside spheres
1450    if (opt_compch==1.and.ispden==1.and.opt_dens<2) then
1451      ABI_ALLOCATE(aa,(mesh_size))
1452      aa(1:mesh_size)=(rho1(1:mesh_size,1,1)-trho1(1:mesh_size,1,1)) &
1453 &     *pawrad%rad(1:mesh_size)**2
1454 ! jmb
1455 !    compchspha(1)=zero
1456      call simp_gen(compchspha(1),aa,pawrad)
1457      compch_sph=compch_sph+compchspha(1)*sqrt(four_pi)
1458      ABI_DEALLOCATE(aa)
1459    end if
1460 
1461 !  -- Print out moments of densities (if requested)
1462    if (abs(pawprtvol)>=2.and.opt_print==1.and.opt_dens<2) then
1463      ABI_ALLOCATE(aa,(cplex*mesh_size))
1464      ABI_ALLOCATE(bb,(cplex*mesh_size))
1465      if (opt_dens==0) then
1466        write(msg,'(2a,i3,a,i1,3a)') ch10, &
1467 &       ' Atom ',iatom,' (ispden=',ispden,'):',ch10,&
1468 &       '  ******* Moment of (n1-tn1)  ** Moment of (n1-tn1-nhat1)'
1469      else
1470        write(msg,'(2a,i3,a,i1,3a)') ch10, &
1471 &       ' Atom ',iatom,' (ispden=',ispden,'):',ch10,&
1472 &       '  ******* Moment of (n1-tn1)'
1473      end if
1474      call wrtout(std_out,msg,'PERS')
1475      do ll=0,pawtab%lcut_size-1
1476        do ilm=ll**2+1,(ll+1)**2
1477          if (lmselectin(ilm)) then
1478            do iplex=1,cplex
1479              if (opt_dens==0) then
1480                do ir=1,mesh_size
1481                  ii=cplex*(ir-1)+iplex
1482                  ro(1)=pawrad%rad(ir)**(2+ll)
1483                  aa(ir)=ro(1)*(rho1(ii,ilm,ispden)-trho1(ii,ilm,ispden))
1484                  bb(ir)=ro(1)*nhat1(ii,ilm,ispden)
1485                end do
1486                call simp_gen(compchspha(iplex),aa,pawrad)
1487                call simp_gen(compchsphb(iplex),bb,pawrad)
1488              else
1489                do ir=1,mesh_size
1490                  ii=cplex*(ir-1)+iplex
1491                  ro(1)=pawrad%rad(ir)**(2+ll)
1492                  aa(ir)=ro(1)*(rho1(ii,ilm,ispden)-trho1(ii,ilm,ispden))
1493                end do
1494                call simp_gen(compchspha(iplex),aa,pawrad)
1495              end if
1496            end do
1497            if (opt_dens==0) then
1498              if (cplex==1) then
1499                write(msg,'(3x,a,2i2,2(a,es14.7))') &
1500 &               'l,m=',ll,ilm-(ll**2+ll+1),': M=',compchspha(1),&
1501 &               ' **    M=',compchspha(1)-compchsphb(1)
1502              else
1503                write(msg,'(3x,a,2i2,2(a,2es14.7))') &
1504 &               'l,m=',ll,ilm-(ll**2+ll+1),': M=',compchspha(1:2),&
1505 &               ' **    M=',compchspha(1:2)-compchsphb(1:2)
1506              end if
1507            else
1508              if (cplex==1) then
1509                write(msg,'(3x,a,2i2,a,es14.7)') &
1510 &               'l,m=',ll,ilm-(ll**2+ll+1),': M=',compchspha(1)
1511              else
1512                write(msg,'(3x,a,2i2,a,2es14.7)') &
1513 &               'l,m=',ll,ilm-(ll**2+ll+1),': M=',compchspha(1:2)
1514              end if
1515            end if
1516            call wrtout(std_out,msg,'PERS')
1517          end if
1518        end do
1519      end do
1520      ABI_DEALLOCATE(aa)
1521      ABI_DEALLOCATE(bb)
1522    end if
1523 
1524 !  ----- End loop over spin components
1525  end do
1526 
1527  if (.not.present(one_over_rad2))  then
1528    ABI_DEALLOCATE(one_over_rad2_)
1529  end if
1530 
1531  DBG_EXIT("COLL")
1532 
1533 end subroutine pawdensities