TABLE OF CONTENTS


ABINIT/m_paw_nhat [ Modules ]

[ Top ] [ Modules ]

NAME

  m_paw_nhat

FUNCTION

  This module contains several routines related to the PAW compensation
    charge density (i.e. n^hat(r)).

COPYRIGHT

 Copyright (C) 2018-2018 ABINIT group (FJ, MT, MG, TRangel)
 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

18 #if defined HAVE_CONFIG_H
19 #include "config.h"
20 #endif
21 
22 #include "abi_common.h"
23 
24 MODULE m_paw_nhat
25 
26  use defs_basis
27  use m_abicore
28  use defs_abitypes
29  use m_errors
30  use m_xmpi
31 
32  use m_time,         only : timab
33  use m_pawang,       only : pawang_type
34  use m_pawtab,       only : pawtab_type
35  use m_pawfgrtab,    only : pawfgrtab_type
36  use m_pawrhoij,     only : pawrhoij_type
37  use m_pawcprj,      only : pawcprj_type
38  use m_paw_finegrid, only : pawgylm,pawrfgd_fft,pawrfgd_wvl,pawexpiqr
39  use m_paral_atom,   only : get_my_atmtab, free_my_atmtab
40  use m_distribfft,   only : distribfft_type,init_distribfft_seq,destroy_distribfft
41  use m_geometry,     only : xred2xcart
42  use m_cgtools,      only : mean_fftr
43  use m_mpinfo,       only : set_mpi_enreg_fft,unset_mpi_enreg_fft,initmpi_seq
44  use m_fft,          only : zerosym, fourwf, fourdp
45  use m_paw_lmn,      only : klmn2ijlmn
46 
47  implicit none
48 
49  private
50 
51 !public procedures.
52  public :: pawmknhat        ! Compute compensation charge density on the real space (fine) grid
53  public :: pawmknhat_psipsi ! Compute compensation charge density associated to the product of two WF
54  public :: pawnhatfr        ! Compute frozen part of 1st-order compensation charge density nhat^(1) (DFPT)
55  public :: pawsushat        ! Compute contrib. to the product of two WF from compensation charge density
56  public :: nhatgrid         ! Determine points of the (fine) grid that are located around atoms - PW version
57  public :: wvl_nhatgrid     ! Determine points of the (fine) grid that are located around atoms - WVL version
58 
59 CONTAINS  !========================================================================================

m_paw_nhat/nhatgrid [ Functions ]

[ Top ] [ m_paw_nhat ] [ Functions ]

NAME

 nhatgrid

FUNCTION

 Determine parts of the rectangular (fine) grid that are contained
 inside spheres around atoms (used to compute n_hat density).
 If corresponding option is selected, compute also g_l(r)*Y_lm(r)
 (and derivatives) on this grid (g_l=radial shape function).

INPUTS

  atindx1(natom)=index table for atoms, inverse of atindx (see gstate.f)
  distribfft<type(distribfft_type)>=--optional-- contains all the informations related
                                    to the FFT parallelism and plane sharing
  gmet(3,3)=reciprocal space metric tensor in bohr**-2
  mpi_atmtab(:)=--optional-- indexes of the atoms treated by current proc
  comm_atom=--optional-- MPI communicator over atoms
  comm_fft=--optional-- MPI communicator over FFT components
  my_natom=number of atoms treated by current processor
  natom=total number of atoms in cell
  nattyp(ntypat)= # atoms of each type.
  ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
  ntypat=number of types of atoms in unit cell
  optcut= option for the cut-off radius of spheres:
          if optcut=0, cut-off radius=pawtab%rshp=cut-off radius of compensation charge
          if optcut=1, cut-off radius=pawtab%rpaw=radius of PAW augmentation regions
  optgr0= 1 if g_l(r)*Y_lm(r) are computed
  optgr1= 1 if first derivatives of g_l(r)*Y_lm(r) are computed
  optgr2= 1 if second derivatives of g_l(r)*Y_lm(r) are computed
  optrad= 1 if vectors (r-r_atom) on the fine grid around atoms have to be stored
  pawfgrtab(natom) <type(pawfgrtab_type)>=atomic data given on fine rectangular grid
  pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data
  rprimd(3,3)=dimensional primitive translations in real space (bohr)
  typat(natom)=type (integer) for each atom
  typord=1 if the output is ordered by type of atoms, 0 otherwise
  ucvol=unit cell volume in bohr**3
  xred(3,natom)=reduced dimensionless atomic coordinates

OUTPUT

  pawfgrtab(natom)%ifftsph(nfgd)=FFT index (fine grid) of a points in paw spheres around each atom
  pawfgrtab(natom)%nfgd= number of (fine grid) FFT points in paw spheres around atoms
  if (optgr0==1)
    pawfgrtab(natom)%gylm(nfgd,l_size**2)= g_l(r)*Y_lm(r) around each atom
  if (optgr1==1)
    pawfgrtab(natom)%gylmgr(3,nfgd,l_size**2)= derivatives of g_l(r)*Y_lm(r) wrt cart. coordinates
  if (optgr2==1)
    pawfgrtab(natom)%gylmgr2(6,nfgd,l_size**2)= second derivatives of g_l(r)*Y_lm(r) wrt cart. coordinates
  if (optrad==1)
    pawfgrtab(natom)%rfgd(3,nfgd)= coordinates of r-r_atom around each atom

PARENTS

      afterscfloop,bethe_salpeter,classify_bands,denfgr,exc_plot,m_wfd
      pawmkaewf,respfn,scfcv,screening,sigma,wfk_analyze

CHILDREN

      free_my_atmtab,get_my_atmtab,pawgylm,pawrfgd_fft,timab

SOURCE

1740 subroutine nhatgrid(atindx1,gmet,my_natom,natom,nattyp,ngfft,ntypat,&
1741 & optcut,optgr0,optgr1,optgr2,optrad,pawfgrtab,pawtab,rprimd,typat,ucvol,xred, &
1742 & mpi_atmtab,comm_atom,comm_fft,distribfft,typord) ! optional arguments (parallelism)
1743 
1744 
1745 !This section has been created automatically by the script Abilint (TD).
1746 !Do not modify the following lines by hand.
1747 #undef ABI_FUNC
1748 #define ABI_FUNC 'nhatgrid'
1749 !End of the abilint section
1750 
1751  implicit none
1752 
1753 !Arguments ---------------------------------------------
1754 !scalars
1755  integer,intent(in) :: my_natom,natom,ntypat,optcut,optgr0,optgr1,optgr2,optrad
1756  integer,optional,intent(in) :: comm_atom,comm_fft,typord
1757  real(dp),intent(in) :: ucvol
1758  type(distribfft_type),optional,target,intent(in)  :: distribfft
1759 !arrays
1760  integer,intent(in) :: ngfft(18),typat(natom)
1761  integer,intent(in),target :: atindx1(natom),nattyp(ntypat)
1762  integer,optional,target,intent(in) :: mpi_atmtab(:)
1763  real(dp),intent(in) :: gmet(3,3),rprimd(3,3),xred(3,natom)
1764  type(pawfgrtab_type),intent(inout) :: pawfgrtab(my_natom)
1765  type(pawtab_type),intent(in) :: pawtab(ntypat)
1766 
1767 !Local variables ------------------------------
1768 !scalars
1769  integer :: i3,iat,iatm,iatom,iatom_,iatom_tot,itypat,lm_size,me_fft,my_comm_atom,n1,n2,n3,nfgd
1770  logical :: grid_found,my_atmtab_allocated,paral_atom
1771  real(dp) :: rcut
1772  character(len=500) :: msg
1773 !arrays
1774  integer,allocatable :: ifftsph_tmp(:)
1775  integer,pointer :: my_atindx1(:),my_atmtab(:),my_nattyp(:)
1776  integer, ABI_CONTIGUOUS pointer :: fftn3_distrib(:),ffti3_local(:)
1777  real(dp) :: tsec(2)
1778  real(dp),allocatable :: rfgd_tmp(:,:)
1779 
1780 ! *************************************************************************
1781 
1782  DBG_ENTER("COLL")
1783 
1784  call timab(559,1,tsec)
1785  if (my_natom==0) return
1786 
1787 !Set up parallelism over FFT
1788  me_fft=0
1789  if (present(comm_fft)) then
1790    me_fft=xmpi_comm_rank(comm_fft)
1791  end if
1792 
1793 !Set up parallelism over atoms
1794  paral_atom=(present(comm_atom).and.(my_natom/=natom))
1795  nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab
1796  my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom
1797  call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,&
1798 & my_natom_ref=my_natom)
1799  if (paral_atom) then
1800    ABI_ALLOCATE(my_atindx1,(natom))
1801    ABI_ALLOCATE(my_nattyp,(ntypat))
1802    my_atindx1(:)=0;my_nattyp(:)=0
1803    iat=1
1804    do itypat=1,ntypat
1805      if (my_natom>0) then
1806        do iatom=1,my_natom
1807          if(typat(my_atmtab(iatom))==itypat)then
1808            my_nattyp(itypat)=my_nattyp(itypat)+1
1809            my_atindx1(iat)=iatom
1810            iat=iat+1
1811          end if
1812        end do
1813      end if
1814    end do
1815  else
1816    my_atindx1 => atindx1
1817    my_nattyp => nattyp
1818  end if
1819 
1820 !Get the distrib associated with this fft_grid
1821  n1=ngfft(1);n2=ngfft(2);n3=ngfft(3)
1822  if (present(distribfft)) then
1823    grid_found=.false.
1824    if (n2 == distribfft%n2_coarse) then
1825      if (n3== size(distribfft%tab_fftdp3_distrib)) then
1826        fftn3_distrib => distribfft%tab_fftdp3_distrib
1827        ffti3_local => distribfft%tab_fftdp3_local
1828        grid_found=.true.
1829      end if
1830    end if
1831    if (n2 == distribfft%n2_fine) then
1832      if (n3 == size(distribfft%tab_fftdp3dg_distrib)) then
1833        fftn3_distrib => distribfft%tab_fftdp3dg_distrib
1834        ffti3_local => distribfft%tab_fftdp3dg_local
1835        grid_found = .true.
1836      end if
1837    end if
1838    if (.not.(grid_found)) then
1839      msg='Unable to find an allocated distrib for this fft grid!'
1840      MSG_BUG(msg)
1841    end if
1842  else
1843    ABI_ALLOCATE(fftn3_distrib,(n3))
1844    ABI_ALLOCATE(ffti3_local,(n3))
1845    fftn3_distrib=0;ffti3_local=(/(i3,i3=1,n3)/)
1846  end if
1847 
1848 !Loop over types of atom
1849 !-------------------------------------------
1850  iatm=0
1851  do itypat=1,ntypat
1852 
1853    if (optcut==1) then
1854      rcut=pawtab(itypat)%rpaw
1855    else
1856      rcut=pawtab(itypat)%rshp
1857    end if
1858 
1859 !  Loop over atoms
1860 !  -------------------------------------------
1861    do iat=1,my_nattyp(itypat)
1862      iatm=iatm+1;iatom=my_atindx1(iatm)
1863      iatom_tot=iatom;if (paral_atom) iatom_tot=my_atmtab(iatom)
1864      iatom_=iatom;if(present(typord)) iatom_=merge(iatm,iatom,typord==1)
1865      lm_size=pawfgrtab(iatom_)%l_size**2
1866 
1867 !    ------------------------------------------------------------------
1868 !    A-Determine FFT points and r-R vectors around the atom
1869 !    ------------------------------------------------------------------
1870 
1871      call pawrfgd_fft(ifftsph_tmp,gmet,n1,n2,n3,nfgd,rcut,rfgd_tmp,rprimd,ucvol,&
1872 &     xred(:,iatom_tot),fft_distrib=fftn3_distrib,fft_index=ffti3_local,me_fft=me_fft)
1873 
1874 !    Allocate arrays defining sphere (and related data) around current atom
1875      if (allocated(pawfgrtab(iatom_)%ifftsph)) then
1876        ABI_DEALLOCATE(pawfgrtab(iatom_)%ifftsph)
1877      end if
1878      ABI_ALLOCATE(pawfgrtab(iatom_)%ifftsph,(nfgd))
1879      pawfgrtab(iatom_)%nfgd=nfgd
1880      pawfgrtab(iatom_)%ifftsph(1:nfgd)=ifftsph_tmp(1:nfgd)
1881 
1882      if (optrad==1) then
1883        if (allocated(pawfgrtab(iatom_)%rfgd))  then
1884          ABI_DEALLOCATE(pawfgrtab(iatom_)%rfgd)
1885        end if
1886        ABI_ALLOCATE(pawfgrtab(iatom_)%rfgd,(3,nfgd))
1887        pawfgrtab(iatom_)%rfgd_allocated=1
1888        pawfgrtab(iatom_)%rfgd(1:3,1:nfgd)=rfgd_tmp(1:3,1:nfgd)
1889      end if
1890 
1891      if (optgr0==1) then
1892        if (allocated(pawfgrtab(iatom_)%gylm))  then
1893          ABI_DEALLOCATE(pawfgrtab(iatom_)%gylm)
1894        end if
1895        ABI_ALLOCATE(pawfgrtab(iatom_)%gylm,(nfgd,lm_size))
1896        pawfgrtab(iatom_)%gylm_allocated=1
1897      end if
1898 
1899      if (optgr1==1) then
1900        if (allocated(pawfgrtab(iatom_)%gylmgr))  then
1901          ABI_DEALLOCATE(pawfgrtab(iatom_)%gylmgr)
1902        end if
1903        ABI_ALLOCATE(pawfgrtab(iatom_)%gylmgr,(3,nfgd,lm_size))
1904        pawfgrtab(iatom_)%gylmgr_allocated=1
1905      end if
1906 
1907      if (optgr2==1) then
1908        if (allocated(pawfgrtab(iatom_)%gylmgr2))  then
1909          ABI_DEALLOCATE(pawfgrtab(iatom_)%gylmgr2)
1910        end if
1911        ABI_ALLOCATE(pawfgrtab(iatom_)%gylmgr2,(6,nfgd,lm_size))
1912        pawfgrtab(iatom_)%gylmgr2_allocated=1
1913      end if
1914 
1915 !    ------------------------------------------------------------------
1916 !    B-Calculate g_l(r-R)*Y_lm(r-R) for each r around the atom R
1917 !    ------------------------------------------------------------------
1918      if (optgr0+optgr1+optgr2>0) then
1919        call pawgylm(pawfgrtab(iatom_)%gylm,pawfgrtab(iatom_)%gylmgr,pawfgrtab(iatom_)%gylmgr2,&
1920 &       lm_size,nfgd,optgr0,optgr1,optgr2,pawtab(itypat),rfgd_tmp(:,1:nfgd))
1921      end if
1922 
1923 !    End loops over types/atoms
1924 !    -------------------------------------------
1925      ABI_DEALLOCATE(ifftsph_tmp)
1926      ABI_DEALLOCATE(rfgd_tmp)
1927    end do
1928  end do
1929 
1930 !Destroy atom tables used for parallelism
1931  call free_my_atmtab(my_atmtab,my_atmtab_allocated)
1932  if (paral_atom) then
1933    ABI_DEALLOCATE(my_atindx1)
1934    ABI_DEALLOCATE(my_nattyp)
1935  end if
1936 
1937  if (.not.present(distribfft)) then
1938    ABI_DEALLOCATE(fftn3_distrib)
1939    ABI_DEALLOCATE(ffti3_local)
1940  end if
1941 
1942  call timab(559,2,tsec)
1943 
1944  DBG_EXIT("COLL")
1945 
1946 end subroutine nhatgrid

m_paw_nhat/pawmknhat [ Functions ]

[ Top ] [ m_paw_nhat ] [ Functions ]

NAME

 pawmknhat

FUNCTION

 PAW only:
 Compute compensation charge density (and derivatives) on the fine FFT grid
 Can also compute first-order compensation charge density (RF calculations)

INPUTS

  cplex: if 1, real space 1-order functions on FFT grid are REAL, if 2, COMPLEX
  distribfft<type(distribfft_type)>=--optional-- contains infos related to FFT parallelism
  gprimd(3,3)=dimensional primitive translations for reciprocal space
  ider= 0: nhat(r) is computed
        1: cartesian derivatives of nhat(r) are computed
        2: nhat(r) and derivatives are computed
  idir=direction of atomic displacement (in case of atomic displ. perturb.)
  ipert=index of perturbation; must be 0 for ground-state calculations
  izero=if 1, unbalanced components of nhat(g) have to be set to zero
  me_g0=--optional-- 1 if the current process treat the g=0 plane-wave (only needed when comm_fft is present)
  mpi_atmtab(:)=--optional-- indexes of the atoms treated by current proc
  comm_atom=--optional-- MPI communicator over atoms
  comm_fft=--optional-- MPI communicator over FFT components
  my_natom=number of atoms treated by current processor
  natom=total number of atoms in cell
  nfft=number of point on the rectangular fft grid
  ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
  nhatgrdim= -PAW only- 0 if pawgrnhat array is not used ; 1 otherwise
  ntypat=number of types of atoms in unit cell.
  paral_kgb=--optional-- 1 if "band-FFT" parallelism is activated (only needed when comm_fft is present)
  pawang <type(pawang_type)>=paw angular mesh and related data
  pawfgrtab(my_natom) <type(pawfgrtab_type)>=atomic data given on fine rectangular grid
  pawrhoij(my_natom) <type(pawrhoij_type)>= paw rhoij occupancies and related data
                                         (1st-order occupancies if ipert>0)
  pawrhoij0(my_natom) <type(pawrhoij_type)>= GS paw rhoij occupancies and related data (used only if ipert>0)
                                          set equat to pawrhoij for GS calculations
  pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data
  qphon(3)=wavevector of the phonon (RF only)
  rprimd(3,3)=dimensional primitive translations for real space
  ucvol=volume of the unit cell
  xred(3,natom)= reduced atomic coordinates

OUTPUT

  === if ider=0 or 2
    compch_fft=compensation charge inside spheres computed over fine fft grid
    pawnhat(nfft,ispden)=nhat on fine rectangular grid
  === if ider=1 or 2
    pawgrnhat(nfft,ispden,3)=derivatives of nhat on fine rectangular grid (and derivatives)

PARENTS

      bethe_salpeter,dfpt_scfcv,energy,nres2vres,odamix,paw_qpscgw,pawmkrho
      respfn,scfcv,screening,setup_positron,sigma

CHILDREN

      destroy_distribfft,fourdp,free_my_atmtab,get_my_atmtab
      init_distribfft_seq,initmpi_seq,mean_fftr,pawexpiqr,pawgylm,pawnhatfr
      set_mpi_enreg_fft,timab,unset_mpi_enreg_fft,xmpi_sum,zerosym

SOURCE

124 subroutine pawmknhat(compch_fft,cplex,ider,idir,ipert,izero,gprimd,&
125 &          my_natom,natom,nfft,ngfft,nhatgrdim,nspden,ntypat,pawang,pawfgrtab,&
126 &          pawgrnhat,pawnhat,pawrhoij,pawrhoij0,pawtab,qphon,rprimd,ucvol,usewvl,xred,&
127 &          mpi_atmtab,comm_atom,comm_fft,mpi_comm_wvl,me_g0,paral_kgb,distribfft) ! optional arguments
128 
129 
130 !This section has been created automatically by the script Abilint (TD).
131 !Do not modify the following lines by hand.
132 #undef ABI_FUNC
133 #define ABI_FUNC 'pawmknhat'
134 !End of the abilint section
135 
136  implicit none
137 
138 !Arguments ---------------------------------------------
139 !scalars
140  integer,intent(in) :: cplex,ider,idir,ipert,izero,my_natom,natom,nfft
141  integer,intent(in)  :: usewvl
142  integer,intent(in) :: nhatgrdim,nspden,ntypat
143  integer,optional,intent(in) :: me_g0,comm_atom,comm_fft,mpi_comm_wvl,paral_kgb
144  real(dp),intent(in) :: ucvol
145  real(dp),intent(out) :: compch_fft
146  type(distribfft_type),optional,intent(in),target :: distribfft
147  type(pawang_type),intent(in) :: pawang
148 !arrays
149  integer,intent(in) :: ngfft(18)
150  integer,optional,target,intent(in) :: mpi_atmtab(:)
151  real(dp),intent(in) :: gprimd(3,3),qphon(3),rprimd(3,3),xred(3,natom)
152  real(dp),intent(out) :: pawgrnhat(cplex*nfft,nspden,3*nhatgrdim)
153  real(dp),intent(inout) :: pawnhat(cplex*nfft,nspden) !vz_i
154  type(pawfgrtab_type),intent(inout) :: pawfgrtab(my_natom)
155  type(pawrhoij_type),intent(in) :: pawrhoij(my_natom),pawrhoij0(my_natom)
156  type(pawtab_type),intent(in) :: pawtab(ntypat)
157 
158 !Local variables ---------------------------------------
159 !scalars
160  integer :: dplex,iatom,iatom_tot,ic,ierr,ii,ils,ilslm,irhoij,ispden,itypat
161  integer :: jc,jrhoij,kc,klm,klmn,lmax,lmin,lm_size,mfgd,mm,mpi_comm_sphgrid
162  integer :: my_comm_atom,my_comm_fft,nfgd,nfftot,option,optgr0,optgr1,optgr2,paral_kgb_fft
163  logical :: compute_grad,compute_nhat,my_atmtab_allocated,need_frozen,paral_atom,qeq0
164  logical :: compute_phonons,has_phase
165  type(distribfft_type),pointer :: my_distribfft
166  type(mpi_type) :: mpi_enreg_fft
167 !arrays
168  integer,pointer :: my_atmtab(:)
169  real(dp) :: ro(cplex),ro_ql(cplex),tmp_compch_fft(nspden),tsec(2)
170  real(dp),allocatable :: pawgrnhat_atm(:,:),pawnhat_atm(:),work(:,:)
171 
172 ! *************************************************************************
173 
174  DBG_ENTER("COLL")
175 
176  compute_nhat=(ider==0.or.ider==2)
177  compute_grad=(ider==1.or.ider==2)
178  compute_phonons=(ipert>0.and.ipert<=natom)
179 
180 !Compatibility tests
181  qeq0=(qphon(1)**2+qphon(2)**2+qphon(3)**2<1.d-15)
182  if (present(comm_fft)) then
183    if ((.not.present(paral_kgb)).or.(.not.present(me_g0))) then
184      MSG_BUG('Need paral_kgb and me_g0 with comm_fft !')
185    end if
186  end if
187  if(ider>0.and.nhatgrdim==0) then
188    MSG_BUG(' Gradients of nhat required but not allocated !')
189  end if
190  if (my_natom>0) then
191    if(nspden>1.and.nspden/=pawrhoij(1)%nspden) then
192      MSG_BUG(' Wrong values for nspden and pawrhoij%nspden !')
193    end if
194    if(nspden>1.and.nspden/=pawfgrtab(1)%nspden) then
195      MSG_BUG(' Wrong values for nspden and pawfgrtab%nspden !')
196    end if
197    if(pawrhoij(1)%cplex<cplex) then
198      MSG_BUG('  Must have pawrhoij()%cplex >= cplex !')
199    end if
200    if (compute_phonons.and.(.not.qeq0)) then
201      if (pawfgrtab(1)%rfgd_allocated==0) then
202        MSG_BUG(' pawfgrtab()%rfgd array must be allocated  !')
203      end if
204      if (compute_grad.and.(.not.compute_nhat)) then
205        MSG_BUG(' When q<>0, nhat gradients need nhat !')
206      end if
207    end if
208  end if
209 
210 !nhat1 does not have to be computed for ddk or d2dk
211  if (ipert==natom+1.or.ipert==natom+10) return
212 
213 !Set up parallelism over atoms
214  paral_atom=(present(comm_atom).and.(my_natom/=natom))
215  nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab
216  my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom
217  call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,&
218 & my_natom_ref=my_natom)
219 
220 !Initialisations
221  if ((.not.compute_nhat).and.(.not.compute_grad)) return
222  mfgd=zero;if (my_natom>0) mfgd=maxval(pawfgrtab(1:my_natom)%nfgd)
223  dplex=cplex-1
224  if (compute_nhat) then
225    ABI_ALLOCATE(pawnhat_atm,(cplex*mfgd))
226    pawnhat=zero
227  end if
228  if (compute_grad) then
229    ABI_ALLOCATE(pawgrnhat_atm,(cplex*mfgd,3))
230    pawgrnhat=zero
231  end if
232 
233 !mpi communicators for spherical grid:
234  mpi_comm_sphgrid=xmpi_comm_self !no communicators passed
235  if(present(comm_fft) .and. usewvl==0) mpi_comm_sphgrid=comm_fft
236  if(present(mpi_comm_wvl) .and. usewvl==1) mpi_comm_sphgrid=mpi_comm_wvl
237 
238 !------------------------------------------------------------------------
239 !----- Loop over atoms
240 !------------------------------------------------------------------------
241 
242  do iatom=1,my_natom
243    iatom_tot=iatom;if (paral_atom) iatom_tot=my_atmtab(iatom)
244 
245    itypat=pawrhoij(iatom)%itypat
246    lm_size=pawfgrtab(iatom)%l_size**2
247    need_frozen=((compute_nhat).and.(ipert==iatom_tot.or.ipert==natom+3.or.ipert==natom+4))
248    nfgd=pawfgrtab(iatom)%nfgd
249 !  Eventually compute g_l(r).Y_lm(r) factors for the current atom (if not already done)
250    if (((compute_nhat).and.(pawfgrtab(iatom)%gylm_allocated==0)).or.&
251 &   ((compute_grad).and.(pawfgrtab(iatom)%gylmgr_allocated==0)).or.&
252 &   ((compute_grad.and.need_frozen).and.(pawfgrtab(iatom)%gylmgr2_allocated==0))) then
253      optgr0=0;optgr1=0;optgr2=0
254      if ((compute_nhat).and.(pawfgrtab(iatom)%gylm_allocated==0)) then
255        if (allocated(pawfgrtab(iatom)%gylm))  then
256          ABI_DEALLOCATE(pawfgrtab(iatom)%gylm)
257        end if
258        ABI_ALLOCATE(pawfgrtab(iatom)%gylm,(nfgd,pawfgrtab(iatom)%l_size**2))
259        pawfgrtab(iatom)%gylm_allocated=2;optgr0=1
260      end if
261      if ((compute_grad).and.(pawfgrtab(iatom)%gylmgr_allocated==0)) then
262        if (allocated(pawfgrtab(iatom)%gylmgr))  then
263          ABI_DEALLOCATE(pawfgrtab(iatom)%gylmgr)
264        end if
265        ABI_ALLOCATE(pawfgrtab(iatom)%gylmgr,(3,nfgd,pawfgrtab(iatom)%l_size**2))
266        pawfgrtab(iatom)%gylmgr_allocated=2;optgr1=1
267      end if
268      if ((compute_grad.and.need_frozen).and.(pawfgrtab(iatom)%gylmgr2_allocated==0)) then
269        if (allocated(pawfgrtab(iatom)%gylmgr2))  then
270          ABI_DEALLOCATE(pawfgrtab(iatom)%gylmgr2)
271        end if
272        ABI_ALLOCATE(pawfgrtab(iatom)%gylmgr2,(6,nfgd,pawfgrtab(iatom)%l_size**2))
273        pawfgrtab(iatom)%gylmgr2_allocated=2;optgr2=1
274      end if
275      if (optgr0+optgr1+optgr2>0) then
276        call pawgylm(pawfgrtab(iatom)%gylm,pawfgrtab(iatom)%gylmgr,pawfgrtab(iatom)%gylmgr2,&
277 &       lm_size,nfgd,optgr0,optgr1,optgr2,pawtab(itypat),pawfgrtab(iatom)%rfgd)
278      end if
279    end if
280 
281 
282 !  Eventually compute exp(-i.q.r) factors for the current atom (if not already done)
283    if (compute_phonons.and.(.not.qeq0).and.pawfgrtab(iatom)%expiqr_allocated==0) then
284      if (allocated(pawfgrtab(iatom)%expiqr))  then
285        ABI_DEALLOCATE(pawfgrtab(iatom)%expiqr)
286      end if
287      ABI_ALLOCATE(pawfgrtab(iatom)%expiqr,(2,nfgd))
288      call pawexpiqr(pawfgrtab(iatom)%expiqr,gprimd,nfgd,qphon,&
289 &     pawfgrtab(iatom)%rfgd,xred(:,iatom_tot))
290      pawfgrtab(iatom)%expiqr_allocated=2
291    end if
292    has_phase=(compute_phonons.and.pawfgrtab(iatom)%expiqr_allocated/=0)
293 
294 !  Eventually compute frozen part of nhat for the current atom (if not already done)
295    if ((need_frozen).and.((pawfgrtab(iatom)%nhatfr_allocated==0).or.&
296 &   (compute_grad.and.pawfgrtab(iatom)%nhatfrgr_allocated==0))) then
297      if (allocated(pawfgrtab(iatom)%nhatfr))  then
298        ABI_DEALLOCATE(pawfgrtab(iatom)%nhatfr)
299      end if
300      ABI_ALLOCATE(pawfgrtab(iatom)%nhatfr,(nfgd,pawfgrtab(iatom)%nspden))
301      option=0;pawfgrtab(iatom)%nhatfr_allocated=2
302      if (compute_grad) then
303        option=1
304        if (allocated(pawfgrtab(iatom)%nhatfrgr))  then
305          ABI_DEALLOCATE(pawfgrtab(iatom)%nhatfrgr)
306        end if
307        ABI_ALLOCATE(pawfgrtab(iatom)%nhatfrgr,(3,nfgd,pawfgrtab(iatom)%nspden))
308        pawfgrtab(iatom)%nhatfrgr_allocated=2
309      end if
310      call pawnhatfr(option,idir,ipert,1,natom,nspden,ntypat,pawang,pawfgrtab(iatom),&
311 &     pawrhoij0(iatom),pawtab,rprimd)
312    end if
313 
314 !  ------------------------------------------------------------------------
315 !  ----- Loop over density components
316 !  ------------------------------------------------------------------------
317 
318    do ispden=1,nspden
319 
320      if (compute_nhat) pawnhat_atm(1:cplex*nfgd)=zero
321      if (compute_grad) pawgrnhat_atm(1:cplex*nfgd,1:3)=zero
322 
323 !    ------------------------------------------------------------------------
324 !    ----- Loop over ij channels (basis components)
325 !    ------------------------------------------------------------------------
326      jrhoij=1
327      do irhoij=1,pawrhoij(iatom)%nrhoijsel
328        klmn=pawrhoij(iatom)%rhoijselect(irhoij)
329        klm =pawtab(itypat)%indklmn(1,klmn)
330        lmin=pawtab(itypat)%indklmn(3,klmn)
331        lmax=pawtab(itypat)%indklmn(4,klmn)
332 
333 !      Retrieve rhoij
334        if (nspden/=2) then
335          ro(1:cplex)=pawrhoij(iatom)%rhoijp(jrhoij:jrhoij+dplex,ispden)
336        else
337          if (ispden==1) then
338            ro(1:cplex)=pawrhoij(iatom)%rhoijp(jrhoij:jrhoij+dplex,1)&
339 &           +pawrhoij(iatom)%rhoijp(jrhoij:jrhoij+dplex,2)
340          else if (ispden==2) then
341            ro(1:cplex)=pawrhoij(iatom)%rhoijp(jrhoij:jrhoij+dplex,1)
342          end if
343        end if
344        ro(1:cplex)=pawtab(itypat)%dltij(klmn)*ro(1:cplex)
345        if (compute_nhat) then
346          if (cplex==1) then
347            do ils=lmin,lmax,2
348              do mm=-ils,ils
349                ilslm=ils*ils+ils+mm+1
350                if (pawang%gntselect(ilslm,klm)>0) then
351                  ro_ql(1)=ro(1)*pawtab(itypat)%qijl(ilslm,klmn)
352                  do ic=1,nfgd
353                    pawnhat_atm(ic)=pawnhat_atm(ic)+ro_ql(1)*pawfgrtab(iatom)%gylm(ic,ilslm)
354                  end do
355                end if
356              end do
357            end do
358          else
359            do ils=lmin,lmax,2
360              do mm=-ils,ils
361                ilslm=ils*ils+ils+mm+1
362                if (pawang%gntselect(ilslm,klm)>0) then
363                  ro_ql(1:2)=ro(1:2)*pawtab(itypat)%qijl(ilslm,klmn)
364                  do ic=1,nfgd
365                    jc=2*ic-1
366                    pawnhat_atm(jc:jc+1)=pawnhat_atm(jc:jc+1)+ro_ql(1:2)*pawfgrtab(iatom)%gylm(ic,ilslm)
367                  end do
368                end if
369              end do
370            end do
371          end if
372        end if
373 
374        if (compute_grad) then
375          if (cplex==1) then
376            do ils=lmin,lmax,2
377              do mm=-ils,ils
378                ilslm=ils*ils+ils+mm+1
379                if (pawang%gntselect(ilslm,klm)>0) then
380                  ro_ql(1)=ro(1)*pawtab(itypat)%qijl(ilslm,klmn)
381                  do ic=1,nfgd
382                    pawgrnhat_atm(ic,1)=pawgrnhat_atm(ic,1)+ro_ql(1)*pawfgrtab(iatom)%gylmgr(1,ic,ilslm)
383                    pawgrnhat_atm(ic,2)=pawgrnhat_atm(ic,2)+ro_ql(1)*pawfgrtab(iatom)%gylmgr(2,ic,ilslm)
384                    pawgrnhat_atm(ic,3)=pawgrnhat_atm(ic,3)+ro_ql(1)*pawfgrtab(iatom)%gylmgr(3,ic,ilslm)
385                  end do
386                end if
387              end do
388            end do
389          else
390            do ils=lmin,lmax,2
391              do mm=-ils,ils
392                ilslm=ils*ils+ils+mm+1
393                if (pawang%gntselect(ilslm,klm)>0) then
394                  ro_ql(1:2)=ro(1:2)*pawtab(itypat)%qijl(ilslm,klmn)
395                  do ic=1,nfgd
396                    jc=2*ic-1
397                    pawgrnhat_atm(jc:jc+1,1)=pawgrnhat_atm(jc:jc+1,1) &
398 &                   +ro_ql(1:2)*pawfgrtab(iatom)%gylmgr(1,ic,ilslm)
399                    pawgrnhat_atm(jc:jc+1,2)=pawgrnhat_atm(jc:jc+1,2) &
400 &                   +ro_ql(1:2)*pawfgrtab(iatom)%gylmgr(2,ic,ilslm)
401                    pawgrnhat_atm(jc:jc+1,3)=pawgrnhat_atm(jc:jc+1,3) &
402 &                   +ro_ql(1:2)*pawfgrtab(iatom)%gylmgr(3,ic,ilslm)
403                  end do
404                end if
405              end do
406            end do
407          end if
408        end if
409 
410 !      ------------------------------------------------------------------------
411 !      ----- End loop over ij channels
412 !      ------------------------------------------------------------------------
413        jrhoij=jrhoij+pawrhoij(iatom)%cplex
414      end do
415 
416 !    If RF calculation, add frozen part of 1st-order compensation density
417      if (need_frozen) then
418        if (cplex==1) then
419          do ic=1,nfgd
420            pawnhat_atm(ic)=pawnhat_atm(ic)+pawfgrtab(iatom)%nhatfr(ic,ispden)
421          end do
422        else
423          do ic=1,nfgd
424            jc=2*ic-1
425            pawnhat_atm(jc)=pawnhat_atm(jc)+pawfgrtab(iatom)%nhatfr(ic,ispden)
426          end do
427        end if
428        if (compute_grad) then
429          if (cplex==1) then
430            do ic=1,nfgd
431              pawgrnhat_atm(ic,1)=pawgrnhat_atm(ic,1)+pawfgrtab(iatom)%nhatfrgr(1,ic,ispden)
432              pawgrnhat_atm(ic,2)=pawgrnhat_atm(ic,2)+pawfgrtab(iatom)%nhatfrgr(2,ic,ispden)
433              pawgrnhat_atm(ic,3)=pawgrnhat_atm(ic,3)+pawfgrtab(iatom)%nhatfrgr(3,ic,ispden)
434            end do
435          else
436            do ic=1,nfgd
437              jc=2*ic-1
438              pawgrnhat_atm(jc,1)=pawgrnhat_atm(jc,1)+pawfgrtab(iatom)%nhatfrgr(1,ic,ispden)
439              pawgrnhat_atm(jc,2)=pawgrnhat_atm(jc,2)+pawfgrtab(iatom)%nhatfrgr(2,ic,ispden)
440              pawgrnhat_atm(jc,3)=pawgrnhat_atm(jc,3)+pawfgrtab(iatom)%nhatfrgr(3,ic,ispden)
441            end do
442          end if
443        end if
444      end if
445 
446 !    If needed, multiply eventually by exp(-i.q.r) phase
447      if (has_phase) then
448        if (cplex==1) then
449          do ic=1,nfgd
450            pawnhat_atm(ic)=pawnhat_atm(ic)*pawfgrtab(iatom)%expiqr(1,ic)
451          end do
452        else
453          do ic=1,nfgd
454            jc=2*ic-1
455            ro_ql(1)= pawfgrtab(iatom)%expiqr(1,ic)
456            ro_ql(2)=-pawfgrtab(iatom)%expiqr(2,ic)
457            ro(1:2)=pawnhat_atm(jc:jc+1)
458            pawnhat_atm(jc  )=ro(1)*ro_ql(1)-ro(2)*ro_ql(2)
459            pawnhat_atm(jc+1)=ro(2)*ro_ql(1)+ro(1)*ro_ql(2)
460          end do
461        end if
462        if (compute_grad) then
463          if (cplex==1) then
464            do ic=1,nfgd
465              pawgrnhat_atm(ic,1:3)=pawgrnhat_atm(ic,1:3)*pawfgrtab(iatom)%expiqr(1,ic)
466            end do
467          else
468            do ic=1,nfgd
469              jc=2*ic-1
470 !            dn^hat(r)/dr_i * exp(-i.q.r)
471              ro_ql(1)= pawfgrtab(iatom)%expiqr(1,ic)
472              ro_ql(2)=-pawfgrtab(iatom)%expiqr(2,ic)
473              do ii=1,3
474                ro(1:2)=pawgrnhat_atm(jc:jc+1,ii)
475                pawgrnhat_atm(jc  ,ii)=ro(1)*ro_ql(1)-ro(2)*ro_ql(2)
476                pawgrnhat_atm(jc+1,ii)=ro(2)*ro_ql(1)+ro(1)*ro_ql(2)
477              end do
478 !            -i.q_i * [n^hat(r).exp(-i.q.r)]
479              ro(1:2)=pawnhat_atm(jc:jc+1)
480              do ii=1,3
481                pawgrnhat_atm(jc  ,ii)=pawgrnhat_atm(kc  ,ii)+qphon(ii)*ro(2)
482                pawgrnhat_atm(jc+1,ii)=pawgrnhat_atm(kc+1,ii)-qphon(ii)*ro(1)
483              end do
484            end do
485          end if
486        end if
487      end if
488 
489 !    Add the contribution of the atom to the compensation charge
490      if (compute_nhat) then
491        if (cplex==1) then
492          do ic=1,nfgd
493            kc=pawfgrtab(iatom)%ifftsph(ic)
494            pawnhat(kc,ispden)=pawnhat(kc,ispden)+pawnhat_atm(ic)
495          end do
496        else
497          do ic=1,nfgd
498            jc=2*ic-1;kc=2*pawfgrtab(iatom)%ifftsph(ic)-1
499            pawnhat(kc:kc+1,ispden)=pawnhat(kc:kc+1,ispden)+pawnhat_atm(jc:jc+1)
500          end do
501        end if
502      end if
503      if (compute_grad) then
504        if (cplex==1) then
505          do ic=1,nfgd
506            kc=pawfgrtab(iatom)%ifftsph(ic)
507            pawgrnhat(kc,ispden,1:3)=pawgrnhat(kc,ispden,1:3)+pawgrnhat_atm(ic,1:3)
508          end do
509        else
510          do ic=1,nfgd
511            jc=2*ic-1;kc=2*pawfgrtab(iatom)%ifftsph(ic)-1
512            do ii=1,3
513              pawgrnhat(kc:kc+1,ispden,ii)=pawgrnhat(kc:kc+1,ispden,ii)+pawgrnhat_atm(jc:jc+1,ii)
514            end do
515          end do
516        end if
517      end if
518 !    ------------------------------------------------------------------------
519 !    ----- End loop over density components
520 !    ------------------------------------------------------------------------
521    end do
522 
523    if (pawfgrtab(iatom)%gylm_allocated==2) then
524      ABI_DEALLOCATE(pawfgrtab(iatom)%gylm)
525      ABI_ALLOCATE(pawfgrtab(iatom)%gylm,(0,0))
526      pawfgrtab(iatom)%gylm_allocated=0
527    end if
528    if (pawfgrtab(iatom)%gylmgr_allocated==2) then
529      ABI_DEALLOCATE(pawfgrtab(iatom)%gylmgr)
530      ABI_ALLOCATE(pawfgrtab(iatom)%gylmgr,(0,0,0))
531      pawfgrtab(iatom)%gylmgr_allocated=0
532    end if
533    if (pawfgrtab(iatom)%gylmgr2_allocated==2) then
534      ABI_DEALLOCATE(pawfgrtab(iatom)%gylmgr2)
535      ABI_ALLOCATE(pawfgrtab(iatom)%gylmgr2,(0,0,0))
536      pawfgrtab(iatom)%gylmgr2_allocated=0
537    end if
538    if (pawfgrtab(iatom)%nhatfr_allocated==2) then
539      ABI_DEALLOCATE(pawfgrtab(iatom)%nhatfr)
540      ABI_ALLOCATE(pawfgrtab(iatom)%nhatfr,(0,0))
541      pawfgrtab(iatom)%nhatfr_allocated=0
542    end if
543    if (pawfgrtab(iatom)%nhatfrgr_allocated==2) then
544      ABI_DEALLOCATE(pawfgrtab(iatom)%nhatfrgr)
545      ABI_ALLOCATE(pawfgrtab(iatom)%nhatfrgr,(0,0,0))
546      pawfgrtab(iatom)%nhatfrgr_allocated=0
547    end if
548    if (pawfgrtab(iatom)%expiqr_allocated==2) then
549      ABI_DEALLOCATE(pawfgrtab(iatom)%expiqr)
550      ABI_ALLOCATE(pawfgrtab(iatom)%expiqr,(0,0))
551      pawfgrtab(iatom)%expiqr_allocated=0
552    end if
553 
554 !  ------------------------------------------------------------------------
555 !  ----- End loop over atoms
556 !  ------------------------------------------------------------------------
557  end do
558 
559 !----- Free some memory
560  if (compute_nhat) then
561    ABI_DEALLOCATE(pawnhat_atm)
562  end if
563  if (compute_grad) then
564    ABI_DEALLOCATE(pawgrnhat_atm)
565  end if
566 
567 !----- Reduction in case of parallelism
568  if (paral_atom) then
569    call timab(48,1,tsec)
570    if (compute_nhat) then
571      call xmpi_sum(pawnhat,my_comm_atom,ierr)
572    end if
573    if (compute_grad) then
574      call xmpi_sum(pawgrnhat,my_comm_atom,ierr)
575    end if
576    call timab(48,2,tsec)
577  end if
578 
579 !----- Avoid unbalanced g-components numerical errors
580  if (izero==1.and.compute_nhat.and.usewvl==0) then
581 !  Create fake mpi_enreg to wrap fourdp
582    if (present(distribfft)) then
583      my_distribfft => distribfft
584    else
585      ABI_DATATYPE_ALLOCATE(my_distribfft,)
586      call init_distribfft_seq(my_distribfft,'f',ngfft(2),ngfft(3),'fourdp')
587    end if
588    call initmpi_seq(mpi_enreg_fft)
589    ABI_DATATYPE_DEALLOCATE(mpi_enreg_fft%distribfft)
590    if (present(comm_fft)) then
591      call set_mpi_enreg_fft(mpi_enreg_fft,comm_fft,my_distribfft,me_g0,paral_kgb)
592      my_comm_fft=comm_fft;paral_kgb_fft=paral_kgb
593    else
594      my_comm_fft=xmpi_comm_self;paral_kgb_fft=0;
595      mpi_enreg_fft%distribfft => my_distribfft
596    end if
597 !  do FFT
598    ABI_ALLOCATE(work,(2,nfft))
599    do ispden=1,min(2,nspden)
600      call fourdp(cplex,work,pawnhat(:,ispden),-1,mpi_enreg_fft,nfft,ngfft,paral_kgb_fft,0)
601      call zerosym(work,2,ngfft(1),ngfft(2),ngfft(3),comm_fft=my_comm_fft,distribfft=my_distribfft)
602      call fourdp(cplex,work,pawnhat(:,ispden),+1,mpi_enreg_fft,nfft,ngfft,paral_kgb_fft,0)
603    end do
604    ABI_DEALLOCATE(work)
605 !  Destroy fake mpi_enreg
606    call unset_mpi_enreg_fft(mpi_enreg_fft)
607    if (.not.present(distribfft)) then
608      call destroy_distribfft(my_distribfft)
609      ABI_DATATYPE_DEALLOCATE(my_distribfft)
610    end if
611  end if
612 
613 !----- Computation of compensation charge over real space grid
614  if (compute_nhat.and.ipert==0) then
615    nfftot=PRODUCT(ngfft(1:3))
616    call mean_fftr(pawnhat,tmp_compch_fft,nfft,nfftot,1,mpi_comm_sphgrid)
617    compch_fft = tmp_compch_fft(1)
618    compch_fft=compch_fft*ucvol
619  end if
620 
621 !Destroy atom table used for parallelism
622  call free_my_atmtab(my_atmtab,my_atmtab_allocated)
623 
624  DBG_EXIT("COLL")
625 
626 end subroutine pawmknhat

m_paw_nhat/pawmknhat_psipsi [ Functions ]

[ Top ] [ m_paw_nhat ] [ Functions ]

NAME

 pawmknhat_psipsi

FUNCTION

 PAW only:
 Compute on the fine FFT grid the compensation charge density (and derivatives) associated
 to the product of two wavefunctions n_{12}(r) = \Psi_1* \Psi_2. Based on pawmknhat.

INPUTS

  cprj1(natom,nspinor), cprj2(natom,nspinor) <type(pawcprj_type)>=
   projected input wave functions <Proj_i|Cnk> with all NL projectors corresponding to
   the \Psi_1 and \Psi_2, respectively.
  distribfft<type(distribfft_type)>=--optional-- contains infos related to FFT parallelism
  ider= 0: nhat(r) is computed
        1: cartesian derivatives of nhat(r) are computed
        2: nhat(r) and derivatives are computed
        3: nhat(r) and gradients of nhat  wrt atomic coordinates are computed
        Note: ider>0 not compatible with ipert>0
  izero=if 1, unbalanced components of nhat(g) have to be set to zero
  me_g0=--optional-- 1 if the current process treat the g=0 plane-wave (only needed when comm_fft is present)
  mpi_atmtab(:)=--optional-- indexes of the atoms treated by current proc
  comm_atom=--optional-- MPI communicator over atoms
  comm_fft=--optional-- MPI communicator over FFT components
  my_natom=number of atoms treated by current processor
  natom=total number of atoms in cell
  nfft=number of point on the rectangular fft grid
  ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
  nhat12_grdim= 0 if grnhat12 array is not used ; 1 otherwise
  ntypat=number of types of atoms in unit cell.
  paral_kgb=--optional-- 1 if "band-FFT" parallelism is activated (only needed when comm_fft is present)
  pawang <type(pawang_type)>=paw angular mesh and related data
  pawfgrtab(my_natom) <type(pawfgrtab_type)>=atomic data given on fine rectangular grid
  pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data

OUTPUT

  === if ider=0 or 2
    nhat12(2,nfft,nspinor**2)=nhat on fine rectangular grid*exp(iqr)
  === if ider=1 or 2
    grnhat12(nfft,nspinor**2,3)=gradient of (nhat*exp(iqr)) on fine rectangular grid (derivative versus r)
  === if ider=3
    grnhat_12(nfft,nspinor**2,3,natom*(ider/3))=derivatives of nhat on fine rectangular grid versus R*exp(iqr)

PARENTS

      calc_sigx_me,fock_getghc,prep_calc_ucrpa

CHILDREN

      destroy_distribfft,fourdp,free_my_atmtab,get_my_atmtab
      init_distribfft_seq,initmpi_seq,pawexpiqr,pawgylm,set_mpi_enreg_fft
      timab,unset_mpi_enreg_fft,xmpi_sum,zerosym

SOURCE

 684 subroutine pawmknhat_psipsi(cprj1,cprj2,ider,izero,my_natom,natom,nfft,ngfft,nhat12_grdim,&
 685 &          nspinor,ntypat,pawang,pawfgrtab,grnhat12,nhat12,pawtab, &
 686 &          gprimd,grnhat_12,qphon,xred,atindx,mpi_atmtab,comm_atom,comm_fft,me_g0,paral_kgb,distribfft) ! optional arguments
 687 
 688 
 689 !This section has been created automatically by the script Abilint (TD).
 690 !Do not modify the following lines by hand.
 691 #undef ABI_FUNC
 692 #define ABI_FUNC 'pawmknhat_psipsi'
 693 !End of the abilint section
 694 
 695  implicit none
 696 
 697 !Arguments ---------------------------------------------
 698 !scalars
 699  integer,intent(in) :: ider,izero,my_natom,natom,nfft,nhat12_grdim,ntypat,nspinor
 700  integer,optional,intent(in) :: me_g0,comm_fft,paral_kgb
 701  integer,optional,intent(in) :: comm_atom
 702  type(distribfft_type),optional,intent(in),target :: distribfft
 703  type(pawang_type),intent(in) :: pawang
 704 !arrays
 705  integer,intent(in) :: ngfft(18)
 706  integer,optional,intent(in) ::atindx(natom)
 707  integer,optional,target,intent(in) :: mpi_atmtab(:)
 708  real(dp),optional, intent(in) ::gprimd(3,3),qphon(3),xred(3,natom)
 709  real(dp),intent(out) :: grnhat12(2,nfft,nspinor**2,3*nhat12_grdim)
 710  real(dp),optional,intent(out) :: grnhat_12(2,nfft,nspinor**2,3,natom*(ider/3))
 711  real(dp),intent(out) :: nhat12(2,nfft,nspinor**2)
 712  type(pawfgrtab_type),intent(inout) :: pawfgrtab(my_natom)
 713  type(pawtab_type),intent(in) :: pawtab(ntypat)
 714  type(pawcprj_type),intent(in) :: cprj1(natom,nspinor),cprj2(natom,nspinor)
 715 
 716 !Local variables ---------------------------------------
 717 !scalars
 718  integer :: iatm,iatom,iatom_tot,ic,ierr,ils,ilslm,isp1,isp2,isploop,itypat,jc,klm,klmn
 719  integer :: lmax,lmin,lm_size,mm,my_comm_atom,my_comm_fft,optgr0,optgr1,paral_kgb_fft
 720  integer :: cplex,ilmn,jlmn,lmn_size,lmn2_size
 721  real(dp) :: re_p,im_p
 722  logical :: compute_grad,compute_grad1,compute_nhat,my_atmtab_allocated,paral_atom,qeq0,compute_phonon,order
 723  type(distribfft_type),pointer :: my_distribfft
 724  type(mpi_type) :: mpi_enreg_fft
 725 !arrays
 726  integer,parameter :: spinor_idxs(2,4)=RESHAPE((/1,1,2,2,1,2,2,1/),(/2,4/))
 727  integer,pointer :: my_atmtab(:)
 728  real(dp) :: rdum(1),cpf(2),cpf_ql(2),tsec(2),ro(2),ro_ql(2),nhat12_atm(2,nfft,nspinor**2)
 729  real(dp),allocatable :: work(:,:), qijl(:,:)
 730 
 731 ! *************************************************************************
 732 
 733  DBG_ENTER("COLL")
 734 
 735 !Compatibility tests
 736  if (present(comm_fft)) then
 737    if ((.not.present(paral_kgb)).or.(.not.present(me_g0))) then
 738      MSG_BUG('Need paral_kgb and me_g0 with comm_fft !')
 739    end if
 740    if (present(paral_kgb)) then
 741      if (paral_kgb/=0) then
 742        MSG_BUG('paral_kgb/=0 not coded!')
 743      end if
 744    end if
 745  end if
 746  if (ider>0.and.nhat12_grdim==0) then
 747 !   MSG_BUG('Gradients of nhat required but not allocated !')
 748  end if
 749  if (nspinor==2) then
 750    MSG_BUG('nspinor==2 not coded!')
 751  end if
 752 
 753  compute_phonon=.false.;qeq0=.false.
 754  if (present(gprimd).and.present(qphon).and.present(xred)) compute_phonon=.true.
 755  if (compute_phonon) qeq0=(qphon(1)**2+qphon(2)**2+qphon(3)**2<1.d-15)
 756  if (present(atindx)) order=.true.
 757 !Set up parallelism over atoms
 758  paral_atom=(present(comm_atom).and.(my_natom/=natom))
 759  nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab
 760  my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom
 761  call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,my_natom_ref=my_natom)
 762 
 763 !Initialisations
 764  compute_nhat=(ider==0.or.ider==2.or.ider==3)
 765  compute_grad=(ider==1.or.ider==2)
 766  compute_grad1=(ider==3)
 767  if ((.not.compute_nhat).and.(.not.compute_grad)) return
 768 
 769  if (compute_nhat) nhat12=zero
 770  if (compute_grad) grnhat12=zero
 771  if (compute_grad1) grnhat_12=zero
 772 
 773  if (compute_grad) then
 774 !   MSG_BUG('compute_grad not tested!')
 775  end if
 776 
 777 !------------------------------------------------------------------------
 778 !----- Loop over atoms
 779 !------------------------------------------------------------------------
 780  do iatom=1,my_natom
 781    iatom_tot=iatom;if (paral_atom) iatom_tot=my_atmtab(iatom)
 782    iatm=iatom_tot
 783    if (order) iatm=atindx(iatom_tot)
 784    itypat    = pawfgrtab(iatom)%itypat
 785    lm_size   = pawfgrtab(iatom)%l_size**2
 786    lmn_size  = pawtab(itypat)%lmn_size
 787    lmn2_size = pawtab(itypat)%lmn2_size
 788    ABI_ALLOCATE(qijl,(lm_size,lmn2_size))
 789    qijl=zero
 790    qijl=pawtab(itypat)%qijl
 791    if (compute_nhat) nhat12_atm=zero
 792 !  Eventually compute g_l(r).Y_lm(r) factors for the current atom (if not already done)
 793    if (((compute_nhat).and.(pawfgrtab(iatom)%gylm_allocated==0)).or.&
 794 &   (((compute_grad).or.(compute_grad1)).and.(pawfgrtab(iatom)%gylmgr_allocated==0))) then
 795      optgr0=0; optgr1=0
 796      if ((compute_nhat).and.(pawfgrtab(iatom)%gylm_allocated==0)) then
 797        if (allocated(pawfgrtab(iatom)%gylm))  then
 798          ABI_DEALLOCATE(pawfgrtab(iatom)%gylm)
 799        end if
 800        ABI_ALLOCATE(pawfgrtab(iatom)%gylm,(pawfgrtab(iatom)%nfgd,pawfgrtab(iatom)%l_size**2))
 801        pawfgrtab(iatom)%gylm_allocated=2;optgr0=1
 802      end if
 803 
 804      if (((compute_grad).or.(compute_grad1)).and.(pawfgrtab(iatom)%gylmgr_allocated==0)) then
 805        if (allocated(pawfgrtab(iatom)%gylmgr))  then
 806          ABI_DEALLOCATE(pawfgrtab(iatom)%gylmgr)
 807        end if
 808        ABI_ALLOCATE(pawfgrtab(iatom)%gylmgr,(3,pawfgrtab(iatom)%nfgd,pawfgrtab(iatom)%l_size**2))
 809        pawfgrtab(iatom)%gylmgr_allocated=2;optgr1=1
 810      end if
 811 
 812      if (optgr0+optgr1>0) then
 813        call pawgylm(pawfgrtab(iatom)%gylm,pawfgrtab(iatom)%gylmgr,rdum,&
 814 &       lm_size,pawfgrtab(iatom)%nfgd,optgr0,optgr1,0,pawtab(itypat),&
 815 &       pawfgrtab(iatom)%rfgd)
 816      end if
 817 
 818    end if
 819    if (compute_phonon.and.(.not.qeq0).and.(pawfgrtab(iatom)%expiqr_allocated==0)) then
 820      if (allocated(pawfgrtab(iatom)%expiqr))  then
 821        ABI_DEALLOCATE(pawfgrtab(iatom)%expiqr)
 822      end if
 823      ABI_ALLOCATE(pawfgrtab(iatom)%expiqr,(2,pawfgrtab(iatom)%nfgd))
 824      call pawexpiqr(pawfgrtab(iatom)%expiqr,gprimd,pawfgrtab(iatom)%nfgd,qphon,&
 825 &     pawfgrtab(iatom)%rfgd,xred(:,iatom_tot))
 826      pawfgrtab(iatom)%expiqr_allocated=2
 827    end if
 828 
 829    do isploop=1,nspinor**2    ! Loop over density components of the compensation charge.
 830 !    TODO Here we might take advantage of symmetry relations between the four components if nspinor==2
 831      isp1=spinor_idxs(1,isploop)
 832      isp2=spinor_idxs(2,isploop)
 833 
 834      do klmn=1,lmn2_size  ! Loop over ij channels of this atom type.
 835        klm =pawtab(itypat)%indklmn(1,klmn)
 836        lmin=pawtab(itypat)%indklmn(3,klmn)  ! abs(il-jl)
 837        lmax=pawtab(itypat)%indklmn(4,klmn)  ! il+jl
 838        ilmn=pawtab(itypat)%indklmn(7,klmn)
 839        jlmn=pawtab(itypat)%indklmn(8,klmn)
 840 !       call klmn2ijlmn(klmn,lmn_size,ilmn,jlmn)  ! This mapping should be stored in pawtab_type
 841 
 842 !      Retrieve the factor due to the PAW projections.
 843        re_p =  cprj1(iatm,isp1)%cp(1,ilmn) * cprj2(iatm,isp2)%cp(1,jlmn) &
 844 &       +cprj1(iatm,isp1)%cp(2,ilmn) * cprj2(iatm,isp2)%cp(2,jlmn) &
 845 &       +cprj1(iatm,isp1)%cp(1,jlmn) * cprj2(iatm,isp2)%cp(1,ilmn) &
 846 &       +cprj1(iatm,isp1)%cp(2,jlmn) * cprj2(iatm,isp2)%cp(2,ilmn)
 847 
 848        im_p =  cprj1(iatm,isp1)%cp(1,ilmn) * cprj2(iatm,isp2)%cp(2,jlmn) &
 849 &       -cprj1(iatm,isp1)%cp(2,ilmn) * cprj2(iatm,isp2)%cp(1,jlmn) &
 850 &       +cprj1(iatm,isp1)%cp(1,jlmn) * cprj2(iatm,isp2)%cp(2,ilmn) &
 851 &       -cprj1(iatm,isp1)%cp(2,jlmn) * cprj2(iatm,isp2)%cp(1,ilmn)
 852 
 853        cpf(1)=re_p*pawtab(itypat)%dltij(klmn)*half
 854        cpf(2)=im_p*pawtab(itypat)%dltij(klmn)*half
 855 
 856        if (compute_nhat) then
 857          do ils=lmin,lmax,2   ! Sum over (L,M)
 858            do mm=-ils,ils
 859              ilslm=ils*ils+ils+mm+1
 860              if (pawang%gntselect(ilslm,klm)>0) then
 861                cpf_ql(1)=cpf(1)*qijl(ilslm,klmn)
 862                cpf_ql(2)=cpf(2)*qijl(ilslm,klmn)
 863                do ic=1,pawfgrtab(iatom)%nfgd
 864                  jc=pawfgrtab(iatom)%ifftsph(ic)
 865                  nhat12_atm(1,jc,isploop)=nhat12_atm(1,jc,isploop)+cpf_ql(1)*pawfgrtab(iatom)%gylm(ic,ilslm)
 866                  nhat12_atm(2,jc,isploop)=nhat12_atm(2,jc,isploop)+cpf_ql(2)*pawfgrtab(iatom)%gylm(ic,ilslm)
 867                end do
 868              end if
 869            end do
 870          end do
 871        end if ! compute_nhat
 872 
 873        if (compute_grad) then
 874          do ils=lmin,lmax,2  ! Sum over (L,M)
 875            do mm=-ils,ils
 876              ilslm=ils*ils+ils+mm+1
 877              if (pawang%gntselect(ilslm,klm)>0) then
 878                cpf_ql(1)=cpf(1)*qijl(ilslm,klmn)
 879                cpf_ql(2)=cpf(2)*qijl(ilslm,klmn)
 880                do ic=1,pawfgrtab(iatom)%nfgd
 881                  jc=pawfgrtab(iatom)%ifftsph(ic)
 882                  grnhat12(1,jc,isploop,1)=grnhat12(1,jc,isploop,1)+cpf_ql(1)*pawfgrtab(iatom)%gylmgr(1,ic,ilslm)
 883                  grnhat12(1,jc,isploop,2)=grnhat12(1,jc,isploop,2)+cpf_ql(1)*pawfgrtab(iatom)%gylmgr(2,ic,ilslm)
 884                  grnhat12(1,jc,isploop,3)=grnhat12(1,jc,isploop,3)+cpf_ql(1)*pawfgrtab(iatom)%gylmgr(3,ic,ilslm)
 885 
 886                  grnhat12(2,jc,isploop,1)=grnhat12(2,jc,isploop,1)+cpf_ql(2)*pawfgrtab(iatom)%gylmgr(1,ic,ilslm)
 887                  grnhat12(2,jc,isploop,2)=grnhat12(2,jc,isploop,2)+cpf_ql(2)*pawfgrtab(iatom)%gylmgr(2,ic,ilslm)
 888                  grnhat12(2,jc,isploop,3)=grnhat12(2,jc,isploop,3)+cpf_ql(2)*pawfgrtab(iatom)%gylmgr(3,ic,ilslm)
 889                end do
 890              end if
 891            end do
 892          end do
 893        end if ! compute_grad
 894        if (compute_grad1) then
 895          do ils=lmin,lmax,2  ! Sum over (L,M)
 896            do mm=-ils,ils
 897              ilslm=ils*ils+ils+mm+1
 898              if (pawang%gntselect(ilslm,klm)>0) then
 899                cpf_ql(1)=cpf(1)*qijl(ilslm,klmn)
 900                cpf_ql(2)=cpf(2)*qijl(ilslm,klmn)
 901                do ic=1,pawfgrtab(iatom)%nfgd
 902                  jc=pawfgrtab(iatom)%ifftsph(ic)
 903                  grnhat_12(1,jc,isploop,1,iatom)=grnhat_12(1,jc,isploop,1,iatom)+cpf_ql(1)*pawfgrtab(iatom)%gylmgr(1,ic,ilslm)
 904                  grnhat_12(1,jc,isploop,2,iatom)=grnhat_12(1,jc,isploop,2,iatom)+cpf_ql(1)*pawfgrtab(iatom)%gylmgr(2,ic,ilslm)
 905                  grnhat_12(1,jc,isploop,3,iatom)=grnhat_12(1,jc,isploop,3,iatom)+cpf_ql(1)*pawfgrtab(iatom)%gylmgr(3,ic,ilslm)
 906 
 907                  grnhat_12(2,jc,isploop,1,iatom)=grnhat_12(2,jc,isploop,1,iatom)+cpf_ql(2)*pawfgrtab(iatom)%gylmgr(1,ic,ilslm)
 908                  grnhat_12(2,jc,isploop,2,iatom)=grnhat_12(2,jc,isploop,2,iatom)+cpf_ql(2)*pawfgrtab(iatom)%gylmgr(2,ic,ilslm)
 909                  grnhat_12(2,jc,isploop,3,iatom)=grnhat_12(2,jc,isploop,3,iatom)+cpf_ql(2)*pawfgrtab(iatom)%gylmgr(3,ic,ilslm)
 910                end do
 911              end if
 912            end do
 913          end do
 914        end if ! compute_grad1
 915      end do  ! klmn (ij channels)
 916 !    If needed, multiply eventually by exp(-i.q.r) phase
 917      if (compute_nhat) then
 918        if(compute_phonon.and.(.not.qeq0).and.pawfgrtab(iatom)%expiqr_allocated/=0) then
 919          do ic=1,pawfgrtab(iatom)%nfgd
 920            jc=pawfgrtab(iatom)%ifftsph(ic)
 921            ro_ql(1)= pawfgrtab(iatom)%expiqr(1,ic)
 922            ro_ql(2)= pawfgrtab(iatom)%expiqr(2,ic)
 923            ro(1:2)=nhat12_atm(1:2,jc,isploop)
 924            nhat12_atm(1,jc,isploop)=ro(1)*ro_ql(1)-ro(2)*ro_ql(2)
 925            nhat12_atm(2,jc,isploop)=ro(2)*ro_ql(1)+ro(1)*ro_ql(2)
 926          end do
 927        end if
 928      end if
 929 
 930      if (compute_grad) then
 931        if(compute_phonon.and.(.not.qeq0).and.pawfgrtab(iatom)%expiqr_allocated/=0) then
 932          do ic=1,pawfgrtab(iatom)%nfgd
 933            jc=pawfgrtab(iatom)%ifftsph(ic)
 934            ro_ql(1)= pawfgrtab(iatom)%expiqr(1,ic)
 935            ro_ql(2)= pawfgrtab(iatom)%expiqr(2,ic)
 936            ro(1)=grnhat12(1,jc,isploop,1)-qphon(1)*nhat12_atm(2,jc,isploop)
 937            ro(2)=grnhat12(2,jc,isploop,1)+qphon(1)*nhat12_atm(1,jc,isploop)
 938            grnhat12(1,jc,isploop,1)=ro(1)*ro_ql(1)-ro(2)*ro_ql(2)
 939            grnhat12(2,jc,isploop,1)=ro(2)*ro_ql(1)+ro(1)*ro_ql(2)
 940            ro(1)=grnhat12(1,jc,isploop,2)-qphon(2)*nhat12_atm(2,jc,isploop)
 941            ro(2)=grnhat12(2,jc,isploop,2)+qphon(2)*nhat12_atm(1,jc,isploop)
 942            grnhat12(1,jc,isploop,2)=ro(1)*ro_ql(1)-ro(2)*ro_ql(2)
 943            grnhat12(2,jc,isploop,2)=ro(2)*ro_ql(1)+ro(1)*ro_ql(2)
 944            ro(1)=grnhat12(1,jc,isploop,3)-qphon(3)*nhat12_atm(2,jc,isploop)
 945            ro(2)=grnhat12(2,jc,isploop,3)+qphon(3)*nhat12_atm(1,jc,isploop)
 946            grnhat12(1,jc,isploop,3)=ro(1)*ro_ql(1)-ro(2)*ro_ql(2)
 947            grnhat12(2,jc,isploop,3)=ro(2)*ro_ql(1)+ro(1)*ro_ql(2)
 948          end do
 949        end if
 950      end if
 951      if (compute_grad1) then
 952        if(compute_phonon.and.(.not.qeq0).and.pawfgrtab(iatom)%expiqr_allocated/=0) then
 953          do ic=1,pawfgrtab(iatom)%nfgd
 954            jc=pawfgrtab(iatom)%ifftsph(ic)
 955            ro_ql(1)= pawfgrtab(iatom)%expiqr(1,ic)
 956            ro_ql(2)= pawfgrtab(iatom)%expiqr(2,ic)
 957            ro(1)=grnhat_12(1,jc,isploop,1,iatom)
 958            ro(2)=grnhat_12(2,jc,isploop,1,iatom)
 959            grnhat_12(1,jc,isploop,1,iatom)=ro(1)*ro_ql(1)-ro(2)*ro_ql(2)
 960            grnhat_12(2,jc,isploop,1,iatom)=ro(2)*ro_ql(1)+ro(1)*ro_ql(2)
 961            ro(1)=grnhat_12(1,jc,isploop,2,iatom)
 962            ro(2)=grnhat_12(2,jc,isploop,2,iatom)
 963            grnhat_12(1,jc,isploop,2,iatom)=ro(1)*ro_ql(1)-ro(2)*ro_ql(2)
 964            grnhat_12(2,jc,isploop,2,iatom)=ro(2)*ro_ql(1)+ro(1)*ro_ql(2)
 965            ro(1)=grnhat_12(1,jc,isploop,3,iatom)
 966            ro(2)=grnhat_12(2,jc,isploop,3,iatom)
 967            grnhat_12(1,jc,isploop,3,iatom)=ro(1)*ro_ql(1)-ro(2)*ro_ql(2)
 968            grnhat_12(2,jc,isploop,3,iatom)=ro(2)*ro_ql(1)+ro(1)*ro_ql(2)
 969          end do
 970        end if
 971      end if
 972    end do ! isploop (density components of the compensation charge)
 973 ! accumlate nhat12 for all the atoms
 974    if (compute_nhat) nhat12=nhat12+nhat12_atm
 975 
 976    if (pawfgrtab(iatom)%gylm_allocated==2) then
 977      ABI_DEALLOCATE(pawfgrtab(iatom)%gylm)
 978      ABI_ALLOCATE(pawfgrtab(iatom)%gylm,(0,0))
 979      pawfgrtab(iatom)%gylm_allocated=0
 980    end if
 981    if (pawfgrtab(iatom)%gylmgr_allocated==2) then
 982      ABI_DEALLOCATE(pawfgrtab(iatom)%gylmgr)
 983      ABI_ALLOCATE(pawfgrtab(iatom)%gylmgr,(0,0,0))
 984      pawfgrtab(iatom)%gylmgr_allocated=0
 985    end if
 986    ABI_DEALLOCATE(qijl)
 987    if (pawfgrtab(iatom)%expiqr_allocated==2) then
 988      ABI_DEALLOCATE(pawfgrtab(iatom)%expiqr)
 989      ABI_ALLOCATE(pawfgrtab(iatom)%expiqr,(0,0))
 990      pawfgrtab(iatom)%expiqr_allocated=0
 991    end if
 992 
 993  end do ! iatom
 994 
 995  if (compute_grad1) grnhat_12=-grnhat_12
 996 
 997 !----- Reduction in case of parallelism -----!
 998  if (paral_atom)then
 999    call timab(48,1,tsec)
1000    if (compute_nhat) then
1001      call xmpi_sum(nhat12,my_comm_atom,ierr)
1002    end if
1003    if (compute_grad) then
1004      call xmpi_sum(grnhat12,my_comm_atom,ierr)
1005    end if
1006    if (compute_grad1) then
1007      call xmpi_sum(grnhat_12,my_comm_atom,ierr)
1008    end if
1009    call timab(48,2,tsec)
1010  end if
1011 
1012 !----- Avoid unbalanced g-components numerical errors -----!
1013 
1014  if (izero==1.and.compute_nhat) then
1015 !  Create fake mpi_enreg to wrap fourdp
1016    if (present(distribfft)) then
1017      my_distribfft => distribfft
1018    else
1019      ABI_DATATYPE_ALLOCATE(my_distribfft,)
1020      call init_distribfft_seq(my_distribfft,'f',ngfft(2),ngfft(3),'fourdp')
1021    end if
1022    call initmpi_seq(mpi_enreg_fft)
1023    ABI_DATATYPE_DEALLOCATE(mpi_enreg_fft%distribfft)
1024    if (present(comm_fft)) then
1025      call set_mpi_enreg_fft(mpi_enreg_fft,comm_fft,my_distribfft,me_g0,paral_kgb)
1026      my_comm_fft=comm_fft;paral_kgb_fft=paral_kgb
1027    else
1028      my_comm_fft=xmpi_comm_self;paral_kgb_fft=0;
1029      mpi_enreg_fft%distribfft => my_distribfft
1030    end if
1031 !  Do FFT
1032    ABI_ALLOCATE(work,(2,nfft))
1033    cplex=2
1034    do isp1=1,MIN(2,nspinor**2)
1035      call fourdp(cplex,work,nhat12(:,:,isp1),-1,mpi_enreg_fft,nfft,ngfft,paral_kgb_fft,0)
1036      call zerosym(work,cplex,ngfft(1),ngfft(2),ngfft(3),comm_fft=my_comm_fft,distribfft=my_distribfft)
1037      call fourdp(cplex,work,nhat12(:,:,isp1),+1,mpi_enreg_fft,nfft,ngfft,paral_kgb_fft,0)
1038    end do
1039    ABI_DEALLOCATE(work)
1040 !  Destroy fake mpi_enreg
1041    call unset_mpi_enreg_fft(mpi_enreg_fft)
1042    if (.not.present(distribfft)) then
1043      call destroy_distribfft(my_distribfft)
1044      ABI_DATATYPE_DEALLOCATE(my_distribfft)
1045    end if
1046  end if
1047 
1048 !Destroy atom table used for parallelism
1049  call free_my_atmtab(my_atmtab,my_atmtab_allocated)
1050 
1051  DBG_EXIT("COLL")
1052 
1053 end subroutine pawmknhat_psipsi

m_paw_nhat/pawnhatfr [ Functions ]

[ Top ] [ m_paw_nhat ] [ Functions ]

NAME

 pawnhatfr

FUNCTION

 PAW: Compute frozen part of 1st-order compensation charge density nhat^(1)
      nhatfr(r)=Sum_ij,lm[rhoij_ij.q_ij^l.(g_l(r).Y_lm(r))^(1)]
      Depends on q wave vector but not on first-order wave-function.

INPUTS

  ider=0: computes frozen part of compensation density
       1: computes frozen part of compensation density and cartesian gradients
  idir=direction of atomic displacement (in case of phonons perturb.)
  ipert=nindex of perturbation
  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
  pawang <type(pawang_type)>=paw angular mesh and related data
  pawfgrtab(my_natom) <type(pawfgrtab_type)>=atomic data given on fine rectangular grid
  pawrhoij(my_natom) <type(pawrhoij_type)>= Ground-State paw rhoij occupancies and related data
  pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data
  rprimd(3,3)=dimensional primitive translations for real space

OUTPUT

  pawfgrtab(iatom)%nhatfr(nfgd,nspden)
                  frozen part of charge compensation density (inside PAW spheres)
                  =Sum_ij,lm[rhoij_ij.q_ij^l.(g_l(r).Y_lm(r))^(1)]
  === If ider==1
  pawfgrtab(iatom)%nhatfrgr(3,nfgd,nspden)
                  gradients of frozen part of charge compensation density (inside PAW spheres)
                  =Sum_ij,lm[rhoij_ij.q_ij^l . d/dr((g_l(r).Y_lm(r))^(1))]

PARENTS

      dfpt_nstpaw,dfpt_scfcv,pawmknhat

CHILDREN

      free_my_atmtab,get_my_atmtab,pawgylm

SOURCE

1102 subroutine pawnhatfr(ider,idir,ipert,my_natom,natom,nspden,ntypat,&
1103 &                    pawang,pawfgrtab,pawrhoij,pawtab,rprimd, &
1104 &                    mpi_atmtab,comm_atom) ! optional arguments (parallelism)
1105 
1106 
1107 !This section has been created automatically by the script Abilint (TD).
1108 !Do not modify the following lines by hand.
1109 #undef ABI_FUNC
1110 #define ABI_FUNC 'pawnhatfr'
1111 !End of the abilint section
1112 
1113  implicit none
1114 
1115 !Arguments ------------------------------------
1116 !scalars
1117  integer,intent(in) :: ider,idir,ipert,my_natom,natom,nspden,ntypat
1118  integer,optional,intent(in) :: comm_atom
1119  type(pawang_type),intent(in) :: pawang
1120 !arrays
1121  integer,optional,target,intent(in) :: mpi_atmtab(:)
1122  real(dp),intent(in) :: rprimd(3,3)
1123  type(pawfgrtab_type),intent(inout) :: pawfgrtab(my_natom)
1124  type(pawrhoij_type),intent(in) :: pawrhoij(my_natom)
1125  type(pawtab_type),intent(in) :: pawtab(ntypat)
1126 
1127 !Local variables-------------------------------
1128 !scalars
1129  integer :: iatom,iatom_tot,ic,ils,ilslm,irhoij,isel,ispden,istr,itypat,jrhoij
1130  integer :: klm,klmn,lm_size,lmn2_size,lm0,lmax,lmin,mua,mub,mm,mu,my_comm_atom,nfgd,nu,optgr0,optgr1,optgr2
1131  logical :: my_atmtab_allocated,my_pert,paral_atom
1132  real(dp) :: contrib,ro
1133 !arrays
1134  integer,parameter :: voigt(3,3)=reshape((/1,6,5,6,2,4,5,4,3/),(/3,3/))
1135  integer,parameter :: alpha(9)=(/1,2,3,3,3,2,2,1,1/),beta(9)=(/1,2,3,2,1,1,3,3,2/)
1136  integer,pointer :: my_atmtab(:)
1137  real(dp),allocatable :: nhatfr_tmp(:,:),nhatfrgr_tmp(:,:,:)
1138 
1139 ! *************************************************************************
1140 
1141  DBG_ENTER("COLL")
1142 
1143 !Only relevant for atomic displacement and strain perturbation
1144  if (ipert>natom.and.ipert/=natom+3.and.ipert/=natom+4) return
1145 
1146 !Set up parallelism over atoms
1147  paral_atom=(present(comm_atom).and.(my_natom/=natom))
1148  nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab
1149  my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom
1150  call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,my_natom_ref=my_natom)
1151 
1152 !Compatibility tests
1153  if (my_natom>0) then
1154    if ((pawfgrtab(1)%gylm_allocated==0.or.pawfgrtab(1)%gylmgr_allocated==0).and. &
1155 &   pawfgrtab(1)%rfgd_allocated==0) then
1156      MSG_BUG('  pawfgrtab()%rfgd array must be allocated  !')
1157    end if
1158  end if
1159 
1160  my_pert = (ipert<=natom).or.ipert==natom+3.or.ipert==natom+4
1161 
1162 !Get correct index of strain pertubation
1163  if (ipert==natom+3) istr = idir
1164  if (ipert==natom+4) istr = idir + 3
1165 
1166 !Loops over  atoms
1167  do iatom=1,my_natom
1168    iatom_tot=iatom;if (paral_atom) iatom_tot=my_atmtab(iatom)
1169 
1170 !  Eventually allocate frozen nhat points
1171    if (my_pert) then
1172      if (pawfgrtab(iatom)%nhatfr_allocated==0) then
1173        if (allocated(pawfgrtab(iatom)%nhatfr))  then
1174          ABI_DEALLOCATE(pawfgrtab(iatom)%nhatfr)
1175        end if
1176        ABI_ALLOCATE(pawfgrtab(iatom)%nhatfr,(pawfgrtab(iatom)%nfgd,nspden))
1177        pawfgrtab(iatom)%nhatfr_allocated=1
1178      end if
1179      if (ider==1.and.pawfgrtab(iatom)%nhatfrgr_allocated==0) then
1180        if (allocated(pawfgrtab(iatom)%nhatfrgr))  then
1181          ABI_DEALLOCATE(pawfgrtab(iatom)%nhatfrgr)
1182        end if
1183        ABI_ALLOCATE(pawfgrtab(iatom)%nhatfrgr,(3,pawfgrtab(iatom)%nfgd,nspden))
1184        pawfgrtab(iatom)%nhatfrgr_allocated=1
1185      end if
1186    end if
1187 
1188 !  Select if frozen part of nhat exists for the current perturbation
1189    if ((.not.my_pert).or.(pawfgrtab(iatom)%nhatfr_allocated==0)) cycle
1190 
1191 !  Some atom-dependent quantities
1192    itypat=pawfgrtab(iatom)%itypat
1193    lm_size=pawfgrtab(iatom)%l_size**2
1194    lmn2_size=pawtab(itypat)%lmn2_size
1195 
1196 !  Eventually compute g_l(r).Y_lm(r) factors for the current atom (if not already done)
1197    nfgd=pawfgrtab(iatom)%nfgd
1198    if ((pawfgrtab(iatom)%gylmgr_allocated==0).or. &
1199 &   (pawfgrtab(iatom)%gylmgr2_allocated==0.and.ider==1)) then
1200      optgr0=0;optgr1=0;optgr2=0
1201      if(ipert==natom+3.or.ipert==natom+4)then
1202        if (pawfgrtab(iatom)%gylm_allocated==0) then
1203          if (allocated(pawfgrtab(iatom)%gylm))  then
1204            ABI_DEALLOCATE(pawfgrtab(iatom)%gylm)
1205          end if
1206          ABI_ALLOCATE(pawfgrtab(iatom)%gylm,(nfgd,lm_size))
1207          pawfgrtab(iatom)%gylm_allocated=2;optgr0=1
1208        end if
1209      end if
1210      if (pawfgrtab(iatom)%gylmgr_allocated==0) then
1211        if (allocated(pawfgrtab(iatom)%gylmgr))  then
1212          ABI_DEALLOCATE(pawfgrtab(iatom)%gylmgr)
1213        end if
1214        ABI_ALLOCATE(pawfgrtab(iatom)%gylmgr,(3,nfgd,lm_size))
1215        pawfgrtab(iatom)%gylmgr_allocated=2;optgr1=1
1216      end if
1217      if (ider==1.and.pawfgrtab(iatom)%gylmgr2_allocated==0) then
1218        if (allocated(pawfgrtab(iatom)%gylmgr2))  then
1219          ABI_DEALLOCATE(pawfgrtab(iatom)%gylmgr2)
1220        end if
1221        ABI_ALLOCATE(pawfgrtab(iatom)%gylmgr2,(6,nfgd,lm_size))
1222        pawfgrtab(iatom)%gylmgr2_allocated=2;optgr2=1
1223      end if
1224      if (optgr0+optgr1+optgr2>0) then
1225        call pawgylm(pawfgrtab(iatom)%gylm,pawfgrtab(iatom)%gylmgr,pawfgrtab(iatom)%gylmgr2,&
1226 &       lm_size,nfgd,optgr0,optgr1,optgr2,pawtab(itypat),pawfgrtab(iatom)%rfgd)
1227      end if
1228    end if
1229 
1230 
1231 !  ============ Phonons ====================================
1232    if (ipert<=natom) then
1233 
1234 !    Loop over spin components
1235      do ispden=1,nspden
1236 
1237        ABI_ALLOCATE(nhatfr_tmp,(3,nfgd))
1238        nhatfr_tmp=zero
1239        if (ider==1) then
1240          ABI_ALLOCATE(nhatfrgr_tmp,(3,nfgd,3))
1241          nhatfrgr_tmp=zero
1242        end if
1243 
1244        jrhoij=1
1245        do irhoij=1,pawrhoij(iatom)%nrhoijsel
1246          klmn=pawrhoij(iatom)%rhoijselect(irhoij)
1247          klm =pawtab(itypat)%indklmn(1,klmn)
1248          lmin=pawtab(itypat)%indklmn(3,klmn)
1249          lmax=pawtab(itypat)%indklmn(4,klmn)
1250          if (nspden/=2) then
1251            ro=pawrhoij(iatom)%rhoijp(jrhoij,ispden)
1252          else
1253            if (ispden==1) then
1254              ro=pawrhoij(iatom)%rhoijp(jrhoij,1)+pawrhoij(iatom)%rhoijp(jrhoij,2)
1255            else if (ispden==2) then
1256              ro=pawrhoij(iatom)%rhoijp(jrhoij,1)
1257            end if
1258          end if
1259          ro=pawtab(itypat)%dltij(klmn)*ro
1260          do ils=lmin,lmax,2
1261            lm0=ils**2+ils+1
1262            do mm=-ils,ils
1263              ilslm=lm0+mm;isel=pawang%gntselect(lm0+mm,klm)
1264              if (isel>0) then
1265                do ic=1,nfgd
1266                  do mu=1,3
1267                    contrib=-ro*pawtab(itypat)%qijl(ilslm,klmn)&
1268 &                   *pawfgrtab(iatom)%gylmgr(mu,ic,ilslm)
1269                    nhatfr_tmp(mu,ic)=nhatfr_tmp(mu,ic)+contrib
1270                  end do
1271                end do
1272                if (ider==1) then
1273                  do ic=1,nfgd
1274                    do nu=1,3
1275                      do mu=1,3
1276                        contrib=-ro*pawtab(itypat)%qijl(ilslm,klmn) &
1277 &                       *pawfgrtab(iatom)%gylmgr2(voigt(mu,nu),ic,ilslm)
1278                        nhatfrgr_tmp(mu,ic,nu)=nhatfrgr_tmp(mu,ic,nu)+contrib
1279                      end do
1280                    end do
1281                  end do
1282                end if
1283              end if
1284            end do
1285          end do
1286          jrhoij=jrhoij+pawrhoij(iatom)%cplex
1287        end do
1288 
1289 !      Convert from cartesian to reduced coordinates
1290        do ic=1,nfgd
1291          pawfgrtab(iatom)%nhatfr(ic,ispden)= &
1292 &         rprimd(1,idir)*nhatfr_tmp(1,ic) &
1293 &         +rprimd(2,idir)*nhatfr_tmp(2,ic) &
1294 &         +rprimd(3,idir)*nhatfr_tmp(3,ic)
1295        end do
1296        if (ider==1) then
1297          do nu=1,3
1298            do ic=1,nfgd
1299              pawfgrtab(iatom)%nhatfrgr(nu,ic,ispden)= &
1300 &             rprimd(1,idir)*nhatfrgr_tmp(1,ic,nu) &
1301 &             +rprimd(2,idir)*nhatfrgr_tmp(2,ic,nu) &
1302 &             +rprimd(3,idir)*nhatfrgr_tmp(3,ic,nu)
1303            end do
1304          end do
1305        end if
1306        ABI_DEALLOCATE(nhatfr_tmp)
1307        if (ider==1) then
1308          ABI_DEALLOCATE(nhatfrgr_tmp)
1309        end if
1310 !      End loop over spin components
1311      end do ! ispden
1312 
1313 
1314 !  ============ Elastic tensor ===============================
1315    else if (ipert==natom+3.or.ipert==natom+4) then
1316 !    Loop over spin components
1317      pawfgrtab(iatom)%nhatfr(:,:) = zero
1318      do ispden=1,nspden
1319        jrhoij=1
1320        do irhoij=1,pawrhoij(iatom)%nrhoijsel
1321          klmn=pawrhoij(iatom)%rhoijselect(irhoij)
1322          klm =pawtab(itypat)%indklmn(1,klmn)
1323          lmin=pawtab(itypat)%indklmn(3,klmn)
1324          lmax=pawtab(itypat)%indklmn(4,klmn)
1325          if (nspden/=2) then
1326            ro=pawrhoij(iatom)%rhoijp(jrhoij,ispden)
1327          else
1328            if (ispden==1) then
1329              ro=pawrhoij(iatom)%rhoijp(jrhoij,1)+pawrhoij(iatom)%rhoijp(jrhoij,2)
1330            else if (ispden==2) then
1331              ro=pawrhoij(iatom)%rhoijp(jrhoij,1)
1332            end if
1333          end if
1334          ro=pawtab(itypat)%dltij(klmn)*ro
1335          do ils=lmin,lmax,2
1336            lm0=ils**2+ils+1
1337            do mm=-ils,ils
1338              ilslm=lm0+mm;isel=pawang%gntselect(lm0+mm,klm)
1339              if (isel>0) then
1340 !              Sum{[Q_ij_q^LM^(1)]}
1341                do ic=1,nfgd
1342                  mua=alpha(istr);mub=beta(istr)
1343                  pawfgrtab(iatom)%nhatfr(ic,ispden) = pawfgrtab(iatom)%nhatfr(ic,ispden)+&
1344 &                 ro*pawtab(itypat)%qijl(ilslm,klmn)*half*(&
1345 &                 pawfgrtab(iatom)%gylmgr(mua,ic,ilslm)*pawfgrtab(iatom)%rfgd(mub,ic)&
1346 &                 +pawfgrtab(iatom)%gylmgr(mub,ic,ilslm)*pawfgrtab(iatom)%rfgd(mua,ic))
1347                end do
1348 !              Add volume contribution
1349                if(istr<=3)then
1350                  do ic=1,nfgd
1351                    pawfgrtab(iatom)%nhatfr(ic,ispden) = pawfgrtab(iatom)%nhatfr(ic,ispden)+&
1352 &                   ro*pawtab(itypat)%qijl(ilslm,klmn)*pawfgrtab(iatom)%gylm(ic,ilslm)
1353                  end do
1354                end if
1355                if (ider==1) then
1356                  MSG_ERROR("nhatgr not implemented for strain perturbationxs")
1357 !                 do ic=1,nfgd
1358 !                   do nu=1,6
1359 !                     do mu=1,6
1360 !                       contrib=-ro*pawtab(itypat)%qijl(ilslm,klmn) &
1361 !&                       *pawfgrtab(iatom)%gylmgr2(voigt(mu,nu),ic,ilslm)
1362 !                       nhatfrgr_tmp(mu,ic,nu)=nhatfrgr_tmp(mu,ic,nu)+contrib
1363 !                     end do
1364 !                   end do
1365 !                 end do
1366                end if
1367              end if
1368            end do
1369          end do
1370          jrhoij=jrhoij+pawrhoij(iatom)%cplex
1371        end do
1372      end do ! ispden
1373    end if
1374 
1375 !  Eventually free temporary space for g_l(r).Y_lm(r) gradients and exp(-i.q.r)
1376    if (pawfgrtab(iatom)%gylmgr_allocated==2) then
1377      ABI_DEALLOCATE(pawfgrtab(iatom)%gylmgr)
1378      ABI_ALLOCATE(pawfgrtab(iatom)%gylmgr,(0,0,0))
1379      pawfgrtab(iatom)%gylmgr_allocated=0
1380    end if
1381    if (pawfgrtab(iatom)%gylmgr2_allocated==2) then
1382      ABI_DEALLOCATE(pawfgrtab(iatom)%gylmgr2)
1383      ABI_ALLOCATE(pawfgrtab(iatom)%gylmgr2,(0,0,0))
1384      pawfgrtab(iatom)%gylmgr2_allocated=0
1385    end if
1386 
1387 !  End loop on atoms
1388  end do
1389 
1390 
1391 !Destroy atom table used for parallelism
1392  call free_my_atmtab(my_atmtab,my_atmtab_allocated)
1393 
1394  DBG_EXIT("COLL")
1395 
1396 end subroutine pawnhatfr

m_paw_nhat/pawsushat [ Functions ]

[ Top ] [ m_paw_nhat ] [ Functions ]

NAME

 pawsushat

FUNCTION

 PAW only, for susceptibility matrix:
 Compute contribution to the product of two wavefunctions (exchange charge density)
 from hat (compensation charge) density (in reciprocal space and eventually in real space):
    sushat_{ij,R}(g)=Sum_{L}[Q^L_ijR(g)]

INPUTS

  atindx(natom)=index table for atoms, inverse of atindx
  cprj_k(natom,nspinor*nband_k)= wave functions projected with non-local projectors:
                                 cprj_k=<p_i|Cnk> where p_i is a non-local projector.
                                 WARNING: cprj(iatom,:) ARE SORTED BY ATOM TYPE !!!
  distribfft<type(distribfft_type)>=--optional-- contains infos related to FFT parallelism
  gbound_diel(2*mgfftdiel+8,2)=G sphere boundary for small FFT sphere.
  gylmg_diel(npwdiel,lmax_diel**2,ntypat)= -PAW only- Fourier transform of g_l(r).Y_ml(r) shape functions
  iband1,iband2= indices of the bands concerned with
  ispinor1,ispinor2= indices of spinorial components concerned with
  istwf_k=input option parameter that describes the storage of wfs
  kg_diel(3,npwdiel)=reduced planewave coordinates for the dielectric matrix.
  lmax_diel=1+max. value of l angular momentum used for dielectric matrix
  me_g0=--optional-- 1 if the current process treat the g=0 plane-wave (only needed when comm_fft is present)
  mgfftdiel=maximum size of 1D FFTs, for the computation of the dielectric matrix
  mpi_atmtab(:)=--optional-- indexes of the atoms treated by current proc
  comm_atom=--optional-- MPI communicator over atoms
  comm_fft=--optional-- MPI communicator over FT components
  natom=number of atoms in cell
  nband=number of bands at this k point for that spin polarization
  ndiel4,ndiel5,ndiel6= FFT dimensions, modified to avoid cache trashing
  nfftdiel=number of FFT grid points for the small (diel) grid
  ngfftdiel(18)=contain all needed information about 3D FFT, for dielectric matrix
  nspinor=number of spinorial components of the wavefunctions
  ntypat=number of types of atoms in unit cell.
  optreal=0 if WF product has to be output in reciprocal space
          1 if WF product has to be output in real space
  paral_kgb=--optional-- 1 if "band-FFT" parallelism is activated (only needed when comm_fft is present)
  pawang <type(pawang_type)>=paw angular mesh and related data
  pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data
  ph3d_diel(2,npwdiel,natom*usepaw)=3-dim structure factors, for each atom and plane wave, for dielectric matrix
  typat(natom)=type (integer) for each atom

SIDE EFFECTS

  === if optreal=0
  wfprod(2,npwdiel)=PAW contrib. to product of two wavefunctions (iband1,iband2):
                    is added (in reciprocal space)
  === if optreal=1
  wfraug(2,ndiel4,ndiel5,ndiel6)=PAW contrib. to product of two wavefunctions (iband1,iband2)
                                 is added (in real space)

PARENTS

      susk,suskmm

CHILDREN

      destroy_distribfft,fourwf,free_my_atmtab,get_my_atmtab
      init_distribfft_seq,initmpi_seq,set_mpi_enreg_fft,unset_mpi_enreg_fft
      xmpi_sum

SOURCE

1462 subroutine pawsushat(atindx,cprj_k,gbound_diel,gylmg_diel,iband1,iband2,ispinor1,ispinor2,istwf_k,kg_diel,&
1463 &                    lmax_diel,mgfftdiel,natom,nband,ndiel4,ndiel5,ndiel6,&
1464 &                    ngfftdiel,npwdiel,nspinor,ntypat,optreal,&
1465 &                    pawang,pawtab,ph3d_diel,typat,wfprod,wfraug, &
1466 &                    mpi_atmtab,comm_atom,comm_fft,me_g0,paral_kgb,distribfft) ! optional arguments (parallelism)
1467 
1468 
1469 !This section has been created automatically by the script Abilint (TD).
1470 !Do not modify the following lines by hand.
1471 #undef ABI_FUNC
1472 #define ABI_FUNC 'pawsushat'
1473 !End of the abilint section
1474 
1475  implicit none
1476 
1477 !Arguments ---------------------------------------------
1478 !scalars
1479  integer,intent(in) :: iband1,iband2,ispinor1,ispinor2,istwf_k,lmax_diel,mgfftdiel
1480  integer,intent(in) :: natom,nband,ndiel4,ndiel5,ndiel6,npwdiel,nspinor
1481  integer,intent(in) :: ntypat,optreal
1482  integer,optional,intent(in) :: me_g0,comm_atom,comm_fft,paral_kgb
1483  type(distribfft_type),optional,intent(in),target :: distribfft
1484  type(pawang_type),intent(in) :: pawang
1485 !arrays
1486  integer,intent(in) :: atindx(natom),gbound_diel(2*mgfftdiel+8,2)
1487  integer,intent(in) :: kg_diel(3,npwdiel),ngfftdiel(18),typat(natom)
1488  integer,optional,target,intent(in) :: mpi_atmtab(:)
1489  real(dp),intent(in) :: gylmg_diel(npwdiel,lmax_diel**2,ntypat)
1490  real(dp),intent(in) :: ph3d_diel(2,npwdiel,natom)
1491  real(dp),intent(inout) :: wfprod(2,npwdiel*(1-optreal))
1492  real(dp),intent(inout) :: wfraug(2,ndiel4,ndiel5,ndiel6*optreal)
1493  type(pawcprj_type),intent(in) :: cprj_k(natom,nspinor*nband)
1494  type(pawtab_type),intent(in) :: pawtab(ntypat)
1495 
1496 !Local variables ---------------------------------------
1497 !scalars
1498  integer :: cplex,iatm,iatom,iatom_tot,ibsp1,ibsp2,ierr,il,ilmn,ils,ilslm,ipw
1499  integer :: itypat,j0lmn,jlmn,klm,klmn,lmax,lmin,mm,my_comm_atom,my_comm_fft,my_natom,paral_kgb_fft,tim_fourwf
1500  real(dp) :: phil1,phil2,sgn,weight_dum,wf1,wf2
1501  logical :: my_atmtab_allocated,parity,paral_atom
1502  type(distribfft_type),pointer :: my_distribfft
1503  type(mpi_type) :: mpi_enreg_fft
1504 !arrays
1505  integer,pointer :: my_atmtab(:)
1506  real(dp) :: ro(2),ro_ql(2)
1507  real(dp),allocatable :: dummy(:,:),wfprod_paw(:,:),wfraug_paw(:,:,:,:)
1508 
1509 ! *************************************************************************
1510 
1511  DBG_ENTER("COLL")
1512 
1513  if (present(comm_fft)) then
1514    if ((.not.present(paral_kgb)).or.(.not.present(me_g0))) then
1515      MSG_BUG('Need paral_kgb and me_g0 with comm_fft !')
1516    end if
1517  end if
1518 
1519 !Set up parallelism over atoms
1520  paral_atom=(present(comm_atom).and.(my_natom/=natom))
1521  nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab
1522  my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom
1523  call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom)
1524  my_natom=natom;if (paral_atom) my_natom=size(my_atmtab)
1525 
1526  cplex=1;if (istwf_k>1) cplex=2
1527  ABI_ALLOCATE(wfprod_paw,(2,npwdiel))
1528  wfprod_paw(:,:)=zero
1529  ibsp1=(iband1-1)*nspinor+ispinor1
1530  ibsp2=(iband2-1)*nspinor+ispinor2
1531 
1532 !------------------------------------------------------------------------
1533 !----- Loop over atoms
1534 !------------------------------------------------------------------------
1535  do iatom=1,my_natom
1536    iatom_tot=iatom;if (paral_atom) iatom_tot=my_atmtab(iatom)
1537    iatm=atindx(iatom_tot)
1538    itypat=typat(iatom_tot)
1539 
1540 !  ------------------------------------------------------------------------
1541 !  ----- Loop over ij channels (basis components)
1542 !  ------------------------------------------------------------------------
1543    do jlmn=1,pawtab(itypat)%lmn_size
1544      j0lmn=jlmn*(jlmn-1)/2
1545      do ilmn=1,jlmn
1546        klmn=j0lmn+ilmn
1547        klm =pawtab(itypat)%indklmn(1,klmn)
1548        lmin=pawtab(itypat)%indklmn(3,klmn)
1549        lmax=pawtab(itypat)%indklmn(4,klmn)
1550 
1551        ro(1)=cprj_k(iatm,ibsp1)%cp(1,ilmn)*cprj_k(iatm,ibsp2)%cp(1,jlmn)
1552        if (cplex==2) then
1553          ro(1)=ro(1)+cprj_k(iatm,ibsp1)%cp(2,ilmn)*cprj_k(iatm,ibsp2)%cp(2,jlmn)
1554          ro(2)=cprj_k(iatm,ibsp1)%cp(2,ilmn)*cprj_k(iatm,ibsp2)%cp(1,jlmn) &
1555 &         -cprj_k(iatm,ibsp1)%cp(1,ilmn)*cprj_k(iatm,ibsp2)%cp(2,jlmn)
1556        end if
1557        ro(1:cplex)=ro(1:cplex)*pawtab(itypat)%dltij(klmn)
1558 
1559        do ils=lmin,lmax,2
1560          il=mod(ils,4);parity=(mod(il,2)==0)
1561          sgn=one;if (il>1) sgn=-one
1562 
1563          do mm=-ils,ils
1564            ilslm=ils*ils+ils+mm+1
1565            if (pawang%gntselect(ilslm,klm)>0) then
1566 
1567              ro_ql(1:cplex)=pawtab(itypat)%qijl(ilslm,klmn)*ro(1:cplex)
1568 
1569 !            Compute: Sum_{ijR} [ cpi* cpj qij^l (-i)^l g_l(g) S_lm(g) ]
1570 
1571              if (cplex==1) then
1572                if (parity) then
1573                  do ipw=1,npwdiel
1574                    phil1= sgn*ph3d_diel(1,ipw,iatm)     ! (i)^l.exp(i.g.R)
1575                    phil2= sgn*ph3d_diel(2,ipw,iatm)
1576                    wf1= phil1*ro_ql(1)                  ! cpi* cpj qij^l (-i)^l.exp(-i.g.R)
1577                    wf2=-phil2*ro_ql(1)
1578                    wfprod_paw(1,ipw)=wfprod_paw(1,ipw)+wf1*gylmg_diel(ipw,ilslm,itypat)
1579                    wfprod_paw(2,ipw)=wfprod_paw(2,ipw)+wf2*gylmg_diel(ipw,ilslm,itypat)
1580                  end do
1581                else
1582                  do ipw=1,npwdiel
1583                    phil1=-sgn*ph3d_diel(2,ipw,iatm)  ! (i)^l.exp(i.g.R)
1584                    phil2= sgn*ph3d_diel(1,ipw,iatm)
1585                    wf1= phil1*ro_ql(1)               ! cpi* cpj qij^l (-i)^l.exp(-i.g.R)
1586                    wf2=-phil2*ro_ql(1)
1587                    wfprod_paw(1,ipw)=wfprod_paw(1,ipw)+wf1*gylmg_diel(ipw,ilslm,itypat)
1588                    wfprod_paw(2,ipw)=wfprod_paw(2,ipw)+wf2*gylmg_diel(ipw,ilslm,itypat)
1589                  end do
1590                end if
1591 
1592              else
1593 
1594                if (parity) then
1595                  do ipw=1,npwdiel
1596                    phil1= sgn*ph3d_diel(1,ipw,iatm)     ! (i)^l.exp(i.g.R)
1597                    phil2= sgn*ph3d_diel(2,ipw,iatm)
1598                    wf1=phil1*ro_ql(1)+phil2*ro_ql(2)    ! cpi* cpj qij^l (-i)^l.exp(-i.g.R)
1599                    wf2=phil1*ro_ql(2)-phil2*ro_ql(1)
1600                    wfprod_paw(1,ipw)=wfprod_paw(1,ipw)+wf1*gylmg_diel(ipw,ilslm,itypat)
1601                    wfprod_paw(2,ipw)=wfprod_paw(2,ipw)+wf2*gylmg_diel(ipw,ilslm,itypat)
1602                  end do
1603                else
1604                  do ipw=1,npwdiel
1605                    phil1=-sgn*ph3d_diel(2,ipw,iatm)     ! (i)^l.exp(i.g.R)
1606                    phil2= sgn*ph3d_diel(1,ipw,iatm)
1607                    wf1=phil1*ro_ql(1)+phil2*ro_ql(2)    ! cpi* cpj qij^l (-i)^l.exp(-i.g.R)
1608                    wf2=phil1*ro_ql(2)-phil2*ro_ql(1)
1609                    wfprod_paw(1,ipw)=wfprod_paw(1,ipw)+wf1*gylmg_diel(ipw,ilslm,itypat)
1610                    wfprod_paw(2,ipw)=wfprod_paw(2,ipw)+wf2*gylmg_diel(ipw,ilslm,itypat)
1611                  end do
1612                end if
1613 
1614              end if
1615            end if
1616          end do
1617        end do
1618 
1619 !      ----- End loop over ij channels
1620      end do
1621    end do
1622 
1623 !  ----- End loop over atoms
1624  end do
1625 
1626 !Reduction in case of parallelism over atoms
1627  if (paral_atom) then
1628    call xmpi_sum(wfprod_paw,my_comm_atom,ierr)
1629  end if
1630 
1631  if (optreal==0) then
1632 
1633 !  === Output in reciprocal space
1634    wfprod(:,:)=wfprod(:,:)+wfprod_paw(:,:)
1635 
1636  else
1637 !  === Output in reciprocal space
1638    tim_fourwf=17;weight_dum=0
1639 !  Create fake mpi_enreg to wrap fourdp
1640    if (present(distribfft)) then
1641      my_distribfft => distribfft
1642    else
1643      ABI_DATATYPE_ALLOCATE(my_distribfft,)
1644      call init_distribfft_seq(my_distribfft,'c',ngfftdiel(2),ngfftdiel(3),'fourwf')
1645    end if
1646    call initmpi_seq(mpi_enreg_fft)
1647    ABI_DATATYPE_DEALLOCATE(mpi_enreg_fft%distribfft)
1648    if (present(comm_fft)) then
1649      call set_mpi_enreg_fft(mpi_enreg_fft,comm_fft,my_distribfft,me_g0,paral_kgb)
1650      my_comm_fft=comm_fft;paral_kgb_fft=paral_kgb
1651    else
1652      my_comm_fft=xmpi_comm_self;paral_kgb_fft=0;
1653      mpi_enreg_fft%distribfft => my_distribfft
1654    end if
1655 !  do FFT
1656    ABI_ALLOCATE(wfraug_paw,(2,ndiel4,ndiel5,ndiel6))
1657    call fourwf(1,dummy,wfprod_paw,dummy,wfraug_paw,gbound_diel,gbound_diel,&
1658 &   istwf_k,kg_diel,kg_diel,mgfftdiel,mpi_enreg_fft,1,ngfftdiel,1,npwdiel,&
1659 &   ndiel4,ndiel5,ndiel6,0,paral_kgb_fft,tim_fourwf,weight_dum,weight_dum)
1660    wfraug(:,:,:,:)=wfraug(:,:,:,:)+wfraug_paw(:,:,:,:)
1661    ABI_DEALLOCATE(wfraug_paw)
1662    call unset_mpi_enreg_fft(mpi_enreg_fft)
1663    if (.not.present(distribfft)) then
1664      call destroy_distribfft(my_distribfft)
1665      ABI_DATATYPE_DEALLOCATE(my_distribfft)
1666    end if
1667  end if
1668 
1669  ABI_DEALLOCATE(wfprod_paw)
1670 
1671 !Destroy atom table used for parallelism
1672  call free_my_atmtab(my_atmtab,my_atmtab_allocated)
1673 
1674  DBG_EXIT("COLL")
1675 
1676 end subroutine pawsushat

m_paw_nhat/wvl_nhatgrid [ Functions ]

[ Top ] [ m_paw_nhat ] [ Functions ]

NAME

 wvl_nhatgrid

FUNCTION

 Determine parts of the rectangular (fine) grid that are contained
 inside spheres around atoms (used to compute n_hat density).
 If corresponding option is selected, compute also g_l(r)*Y_lm(r)
 (and derivatives) on this grid (g_l=radial shape function).

INPUTS

OUTPUT

  pawfgrtab(natom)%ifftsph(nfgd)=FFT index (fine grid) of a points in paw spheres around each atom
  pawfgrtab(natom)%nfgd= number of (fine grid) FFT points in paw spheres around atoms
  if (optgr0==1)
    pawfgrtab(natom)%gylm(nfgd,l_size**2)= g_l(r)*Y_lm(r) around each atom
  if (optgr1==1)
    pawfgrtab(natom)%gylmgr(3,nfgd,l_size**2)= derivatives of g_l(r)*Y_lm(r) wrt cart. coordinates
  if (optgr2==1)
    pawfgrtab(natom)%gylmgr2(6,nfgd,l_size**2)= second derivatives of g_l(r)*Y_lm(r) wrt cart. coordinates
  if (optrad==1)
    pawfgrtab(natom)%rfgd(3,nfgd)= coordinates of r-r_atom around each atom

NOTES

   PENDING: ADD PARALELLISM OVER ATOMS:
   COPY NHATGRID

PARENTS

      afterscfloop,scfcv

CHILDREN

      pawgylm,pawrfgd_wvl,timab,xred2xcart

SOURCE

1987 subroutine wvl_nhatgrid(atindx1,geocode,h,i3s,natom,natom_tot,&
1988 & nattyp,ntypat,n1,n1i,n2,n2i,n3,n3pi,optcut,optgr0,optgr1,optgr2,optrad,&
1989 & pawfgrtab,pawtab,psppar,rprimd,shift,xred)
1990 
1991 
1992 !This section has been created automatically by the script Abilint (TD).
1993 !Do not modify the following lines by hand.
1994 #undef ABI_FUNC
1995 #define ABI_FUNC 'wvl_nhatgrid'
1996 !End of the abilint section
1997 
1998  implicit none
1999 
2000 !Arguments ---------------------------------------------
2001 !scalars
2002  integer,intent(in) :: i3s,natom,natom_tot,ntypat,optcut,optgr0,optgr1,optgr2,optrad
2003  integer,intent(in) :: n1,n2,n3,n1i,n2i,n3pi,shift
2004  real(dp),intent(in) :: h(3)
2005  character(1),intent(in) :: geocode
2006 !integer,intent(in),optional :: mpi_comm_wvl
2007 !arrays
2008  integer,intent(in) :: atindx1(natom),nattyp(ntypat)
2009  real(dp),intent(in) :: psppar(0:4,0:6,ntypat),rprimd(3,3)
2010  real(dp),intent(inout) :: xred(3,natom)
2011  type(pawfgrtab_type),intent(inout) :: pawfgrtab(natom)
2012  type(pawtab_type),intent(in) :: pawtab(ntypat)
2013 
2014 !Local variables ------------------------------
2015 !scalars
2016 !buffer to be added at the end of the last dimension of an array to control bounds_check
2017  integer :: iat,iatm,iatom,iatom_tot,itypat,lm_size,nfgd
2018  real(dp) :: rloc,rshp,xcart(3,natom)
2019 !arrays
2020  integer,allocatable :: ifftsph_tmp(:)
2021  real(dp) :: hh(3) !fine grid spacing for wavelets
2022  real(dp) :: tsec(2)
2023  real(dp),allocatable :: rfgd_tmp(:,:)
2024 
2025 ! *************************************************************************
2026 
2027  DBG_ENTER("COLL")
2028 
2029 #if !defined HAVE_BIGDFT
2030  BIGDFT_NOTENABLED_ERROR()
2031 #endif
2032 
2033  call timab(559,1,tsec)
2034 
2035 !Set up parallelism for wvl
2036 !for debug: use me_wvl=xmpi_comm_rank(MPI_COMM_WORLD)
2037 !if (present(mpi_comm_wvl)) then
2038 !me_wvl=xmpi_comm_rank(mpi_comm_wvl)
2039 !nproc_wvl=xmpi_comm_size(mpi_comm_wvl)
2040 !else
2041 !me_wvl=0;nproc_wvl=1
2042 !end if
2043 !Pending: parallelism over atoms: see nhatgrid
2044 
2045  if (natom_tot<natom) then   ! This test has to be remove when natom_tot is used
2046    MSG_BUG(' natom_tot<natom !')
2047  end if
2048 
2049 !Fine grid
2050  hh(:)=0.5d0*h(:)
2051 
2052 !Compute xcart from xred
2053  call xred2xcart(natom,rprimd,xcart,xred)
2054 
2055 !Loop over types of atom
2056  iatm=0
2057  do itypat=1,ntypat
2058 
2059    rloc=psppar(0,0,itypat)
2060    if (optcut==1) then
2061      rshp=pawtab(itypat)%rpaw
2062    else
2063      rshp=pawtab(itypat)%rshp
2064    end if
2065 
2066 !  Loop over atoms
2067    do iat=1,nattyp(itypat)
2068      iatm=iatm+1;iatom=atindx1(iatm)
2069      iatom_tot=iatom; !if (paral_atom) iatom_tot=my_atmtab(iatom)
2070      lm_size=pawfgrtab(iatom)%l_size**2
2071 
2072 !    Determine FFT points and r-R vectors around the atom
2073      call pawrfgd_wvl(geocode,hh,ifftsph_tmp,i3s,n1,n1i,n2,n2i,n3,n3pi,nfgd,rshp,rloc,&
2074 &     rfgd_tmp,shift,xcart(:,iatom_tot))
2075 
2076 !    Allocate arrays defining sphere (and related data) around current atom
2077      if (allocated(pawfgrtab(iatom)%ifftsph)) then
2078        ABI_DEALLOCATE(pawfgrtab(iatom)%ifftsph)
2079      end if
2080      ABI_ALLOCATE(pawfgrtab(iatom)%ifftsph,(nfgd))
2081      pawfgrtab(iatom)%nfgd=nfgd
2082      pawfgrtab(iatom)%ifftsph(1:nfgd)=ifftsph_tmp(1:nfgd)
2083 
2084      if (optrad==1) then
2085        if (allocated(pawfgrtab(iatom)%rfgd)) then
2086          ABI_DEALLOCATE(pawfgrtab(iatom)%rfgd)
2087        end if
2088        ABI_ALLOCATE(pawfgrtab(iatom)%rfgd,(3,nfgd))
2089        pawfgrtab(iatom)%rfgd_allocated=1
2090        pawfgrtab(iatom)%rfgd(1:3,1:nfgd)=rfgd_tmp(1:3,1:nfgd)
2091      end if
2092 
2093      if (optgr0==1) then
2094        if (allocated(pawfgrtab(iatom)%gylm)) then
2095          ABI_DEALLOCATE(pawfgrtab(iatom)%gylm)
2096        end if
2097        ABI_ALLOCATE(pawfgrtab(iatom)%gylm,(nfgd,lm_size))
2098        pawfgrtab(iatom)%gylm_allocated=1
2099      end if
2100 
2101      if (optgr1==1) then
2102        if (allocated(pawfgrtab(iatom)%gylmgr)) then
2103          ABI_DEALLOCATE(pawfgrtab(iatom)%gylmgr)
2104        end if
2105        ABI_ALLOCATE(pawfgrtab(iatom)%gylmgr,(3,nfgd,lm_size))
2106        pawfgrtab(iatom)%gylmgr_allocated=1
2107      end if
2108 
2109      if (optgr2==1) then
2110        if (allocated(pawfgrtab(iatom)%gylmgr2)) then
2111          ABI_DEALLOCATE(pawfgrtab(iatom)%gylmgr2)
2112        end if
2113        ABI_ALLOCATE(pawfgrtab(iatom)%gylmgr2,(6,nfgd,lm_size))
2114        pawfgrtab(iatom)%gylmgr2_allocated=1
2115      end if
2116 
2117 !    Calculate g_l(r-R)*Y_lm(r-R) for each r around the atom R
2118      if (optgr0+optgr1+optgr2>0) then
2119        call pawgylm(pawfgrtab(iatom)%gylm,pawfgrtab(iatom)%gylmgr,pawfgrtab(iatom)%gylmgr2,&
2120 &       lm_size,nfgd,optgr0,optgr1,optgr2,pawtab(itypat),rfgd_tmp(:,1:nfgd))
2121      end if
2122 
2123 !    End loops over types/atoms
2124      ABI_DEALLOCATE(ifftsph_tmp)
2125      ABI_DEALLOCATE(rfgd_tmp)
2126    end do
2127  end do
2128 
2129  call timab(559,2,tsec)
2130 
2131  DBG_EXIT("COLL")
2132 
2133 end subroutine wvl_nhatgrid