TABLE OF CONTENTS
- ABINIT/m_paw_nhat
- m_paw_nhat/nhatgrid
- m_paw_nhat/pawmknhat
- m_paw_nhat/pawmknhat_psipsi
- m_paw_nhat/pawnhatfr
- m_paw_nhat/pawsushat
- m_paw_nhat/wvl_nhatgrid
ABINIT/m_paw_nhat [ 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-2024 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
17 #if defined HAVE_CONFIG_H 18 #include "config.h" 19 #endif 20 21 #include "abi_common.h" 22 23 MODULE m_paw_nhat 24 25 use defs_basis 26 use m_abicore 27 use m_errors 28 use m_xmpi 29 30 use defs_abitypes, only : MPI_type 31 use m_time, only : timab 32 use m_pawang, only : pawang_type 33 use m_pawtab, only : pawtab_type 34 use m_pawfgrtab, only : pawfgrtab_type 35 use m_pawrhoij, only : pawrhoij_type 36 use m_pawcprj, only : pawcprj_type 37 use m_paw_finegrid, only : pawgylm,pawrfgd_fft,pawrfgd_wvl,pawexpiqr 38 use m_paral_atom, only : get_my_atmtab, free_my_atmtab 39 use m_distribfft, only : distribfft_type,init_distribfft_seq,destroy_distribfft 40 use m_geometry, only : xred2xcart 41 use m_cgtools, only : mean_fftr 42 use m_mpinfo, only : set_mpi_enreg_fft,unset_mpi_enreg_fft,initmpi_seq 43 use m_fft, only : zerosym, fourwf, fourdp 44 use m_paw_lmn, only : klmn2ijlmn 45 46 implicit none 47 48 private 49 50 !public procedures. 51 public :: pawmknhat ! Compute compensation charge density on the real space (fine) grid 52 public :: pawmknhat_psipsi ! Compute compensation charge density associated to the product of two WF 53 public :: pawnhatfr ! Compute frozen part of 1st-order compensation charge density nhat^(1) (DFPT) 54 public :: pawsushat ! Compute contrib. to the product of two WF from compensation charge density 55 public :: nhatgrid ! Determine points of the (fine) grid that are located around atoms - PW version 56 public :: wvl_nhatgrid ! Determine points of the (fine) grid that are located around atoms - WVL version 57 58 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 information 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
SOURCE
1683 subroutine nhatgrid(atindx1,gmet,my_natom,natom,nattyp,ngfft,ntypat,& 1684 & optcut,optgr0,optgr1,optgr2,optrad,pawfgrtab,pawtab,rprimd,typat,ucvol,xred, & 1685 & mpi_atmtab,comm_atom,comm_fft,distribfft,typord) ! optional arguments (parallelism) 1686 1687 implicit none 1688 1689 !Arguments --------------------------------------------- 1690 !scalars 1691 integer,intent(in) :: my_natom,natom,ntypat,optcut,optgr0,optgr1,optgr2,optrad 1692 integer,optional,intent(in) :: comm_atom,comm_fft,typord 1693 real(dp),intent(in) :: ucvol 1694 type(distribfft_type),optional,target,intent(in) :: distribfft 1695 !arrays 1696 integer,intent(in) :: ngfft(18),typat(natom) 1697 integer,intent(in),target :: atindx1(natom),nattyp(ntypat) 1698 integer,optional,target,intent(in) :: mpi_atmtab(:) 1699 real(dp),intent(in) :: gmet(3,3),rprimd(3,3),xred(3,natom) 1700 type(pawfgrtab_type),intent(inout) :: pawfgrtab(my_natom) 1701 type(pawtab_type),intent(in) :: pawtab(ntypat) 1702 1703 !Local variables ------------------------------ 1704 !scalars 1705 integer :: i3,iat,iatm,iatom,iatom_,iatom_tot,itypat,lm_size,me_fft,my_comm_atom,n1,n2,n3,nfgd 1706 logical :: grid_found,my_atmtab_allocated,paral_atom 1707 real(dp) :: rcut 1708 character(len=500) :: msg 1709 !arrays 1710 integer,allocatable :: ifftsph_tmp(:) 1711 integer,pointer :: my_atindx1(:),my_atmtab(:),my_nattyp(:) 1712 integer, ABI_CONTIGUOUS pointer :: fftn3_distrib(:),ffti3_local(:) 1713 real(dp) :: tsec(2) 1714 real(dp),allocatable :: rfgd_tmp(:,:) 1715 1716 ! ************************************************************************* 1717 1718 DBG_ENTER("COLL") 1719 1720 call timab(559,1,tsec) 1721 if (my_natom==0) return 1722 1723 !Set up parallelism over FFT 1724 me_fft=0 1725 if (present(comm_fft)) then 1726 me_fft=xmpi_comm_rank(comm_fft) 1727 end if 1728 1729 !Set up parallelism over atoms 1730 paral_atom=(present(comm_atom).and.(my_natom/=natom)) 1731 nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab 1732 my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom 1733 call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,& 1734 & my_natom_ref=my_natom) 1735 if (paral_atom) then 1736 ABI_MALLOC(my_atindx1,(natom)) 1737 ABI_MALLOC(my_nattyp,(ntypat)) 1738 my_atindx1(:)=0;my_nattyp(:)=0 1739 iat=1 1740 do itypat=1,ntypat 1741 if (my_natom>0) then 1742 do iatom=1,my_natom 1743 if(typat(my_atmtab(iatom))==itypat)then 1744 my_nattyp(itypat)=my_nattyp(itypat)+1 1745 my_atindx1(iat)=iatom 1746 iat=iat+1 1747 end if 1748 end do 1749 end if 1750 end do 1751 else 1752 my_atindx1 => atindx1 1753 my_nattyp => nattyp 1754 end if 1755 1756 !Get the distrib associated with this fft_grid 1757 n1=ngfft(1);n2=ngfft(2);n3=ngfft(3) 1758 if (present(distribfft)) then 1759 grid_found=.false. 1760 if (n2 == distribfft%n2_coarse) then 1761 if (n3== size(distribfft%tab_fftdp3_distrib)) then 1762 fftn3_distrib => distribfft%tab_fftdp3_distrib 1763 ffti3_local => distribfft%tab_fftdp3_local 1764 grid_found=.true. 1765 end if 1766 end if 1767 if (n2 == distribfft%n2_fine) then 1768 if (n3 == size(distribfft%tab_fftdp3dg_distrib)) then 1769 fftn3_distrib => distribfft%tab_fftdp3dg_distrib 1770 ffti3_local => distribfft%tab_fftdp3dg_local 1771 grid_found = .true. 1772 end if 1773 end if 1774 if (.not.(grid_found)) then 1775 msg='Unable to find an allocated distrib for this fft grid!' 1776 ABI_BUG(msg) 1777 end if 1778 else 1779 ABI_MALLOC(fftn3_distrib,(n3)) 1780 ABI_MALLOC(ffti3_local,(n3)) 1781 fftn3_distrib=0;ffti3_local=(/(i3,i3=1,n3)/) 1782 end if 1783 1784 !Loop over types of atom 1785 !------------------------------------------- 1786 iatm=0 1787 do itypat=1,ntypat 1788 1789 if (optcut==1) then 1790 rcut=pawtab(itypat)%rpaw 1791 else 1792 rcut=pawtab(itypat)%rshp 1793 end if 1794 1795 ! Loop over atoms 1796 ! ------------------------------------------- 1797 do iat=1,my_nattyp(itypat) 1798 iatm=iatm+1;iatom=my_atindx1(iatm) 1799 iatom_tot=iatom;if (paral_atom) iatom_tot=my_atmtab(iatom) 1800 iatom_=iatom;if(present(typord)) iatom_=merge(iatm,iatom,typord==1) 1801 lm_size=pawfgrtab(iatom_)%l_size**2 1802 1803 ! ------------------------------------------------------------------ 1804 ! A-Determine FFT points and r-R vectors around the atom 1805 ! ------------------------------------------------------------------ 1806 1807 call pawrfgd_fft(ifftsph_tmp,gmet,n1,n2,n3,nfgd,rcut,rfgd_tmp,rprimd,ucvol,& 1808 & xred(:,iatom_tot),fft_distrib=fftn3_distrib,fft_index=ffti3_local,me_fft=me_fft) 1809 1810 ! Allocate arrays defining sphere (and related data) around current atom 1811 if (allocated(pawfgrtab(iatom_)%ifftsph)) then 1812 ABI_FREE(pawfgrtab(iatom_)%ifftsph) 1813 end if 1814 ABI_MALLOC(pawfgrtab(iatom_)%ifftsph,(nfgd)) 1815 pawfgrtab(iatom_)%nfgd=nfgd 1816 pawfgrtab(iatom_)%ifftsph(1:nfgd)=ifftsph_tmp(1:nfgd) 1817 1818 if (optrad==1) then 1819 if (allocated(pawfgrtab(iatom_)%rfgd)) then 1820 ABI_FREE(pawfgrtab(iatom_)%rfgd) 1821 end if 1822 ABI_MALLOC(pawfgrtab(iatom_)%rfgd,(3,nfgd)) 1823 pawfgrtab(iatom_)%rfgd_allocated=1 1824 pawfgrtab(iatom_)%rfgd(1:3,1:nfgd)=rfgd_tmp(1:3,1:nfgd) 1825 end if 1826 1827 if (optgr0==1) then 1828 if (allocated(pawfgrtab(iatom_)%gylm)) then 1829 ABI_FREE(pawfgrtab(iatom_)%gylm) 1830 end if 1831 ABI_MALLOC(pawfgrtab(iatom_)%gylm,(nfgd,lm_size)) 1832 pawfgrtab(iatom_)%gylm_allocated=1 1833 end if 1834 1835 if (optgr1==1) then 1836 if (allocated(pawfgrtab(iatom_)%gylmgr)) then 1837 ABI_FREE(pawfgrtab(iatom_)%gylmgr) 1838 end if 1839 ABI_MALLOC(pawfgrtab(iatom_)%gylmgr,(3,nfgd,lm_size)) 1840 pawfgrtab(iatom_)%gylmgr_allocated=1 1841 end if 1842 1843 if (optgr2==1) then 1844 if (allocated(pawfgrtab(iatom_)%gylmgr2)) then 1845 ABI_FREE(pawfgrtab(iatom_)%gylmgr2) 1846 end if 1847 ABI_MALLOC(pawfgrtab(iatom_)%gylmgr2,(6,nfgd,lm_size)) 1848 pawfgrtab(iatom_)%gylmgr2_allocated=1 1849 end if 1850 1851 ! ------------------------------------------------------------------ 1852 ! B-Calculate g_l(r-R)*Y_lm(r-R) for each r around the atom R 1853 ! ------------------------------------------------------------------ 1854 if (optgr0+optgr1+optgr2>0) then 1855 call pawgylm(pawfgrtab(iatom_)%gylm,pawfgrtab(iatom_)%gylmgr,pawfgrtab(iatom_)%gylmgr2,& 1856 & lm_size,nfgd,optgr0,optgr1,optgr2,pawtab(itypat),rfgd_tmp(:,1:nfgd)) 1857 end if 1858 1859 ! End loops over types/atoms 1860 ! ------------------------------------------- 1861 ABI_FREE(ifftsph_tmp) 1862 ABI_FREE(rfgd_tmp) 1863 end do 1864 end do 1865 1866 !Destroy atom tables used for parallelism 1867 call free_my_atmtab(my_atmtab,my_atmtab_allocated) 1868 if (paral_atom) then 1869 ABI_FREE(my_atindx1) 1870 ABI_FREE(my_nattyp) 1871 end if 1872 1873 if (.not.present(distribfft)) then 1874 ABI_FREE(fftn3_distrib) 1875 ABI_FREE(ffti3_local) 1876 end if 1877 1878 call timab(559,2,tsec) 1879 1880 DBG_EXIT("COLL") 1881 1882 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)
SOURCE
114 subroutine pawmknhat(compch_fft,cplex,ider,idir,ipert,izero,gprimd,& 115 & my_natom,natom,nfft,ngfft,nhatgrdim,nspden,ntypat,pawang,pawfgrtab,& 116 & pawgrnhat,pawnhat,pawrhoij,pawrhoij0,pawtab,qphon,rprimd,ucvol,usewvl,xred,& 117 & mpi_atmtab,comm_atom,comm_fft,mpi_comm_wvl,me_g0,paral_kgb,distribfft) ! optional arguments 118 119 implicit none 120 121 !Arguments --------------------------------------------- 122 !scalars 123 integer,intent(in) :: cplex,ider,idir,ipert,izero,my_natom,natom,nfft 124 integer,intent(in) :: usewvl 125 integer,intent(in) :: nhatgrdim,nspden,ntypat 126 integer,optional,intent(in) :: me_g0,comm_atom,comm_fft,mpi_comm_wvl,paral_kgb 127 real(dp),intent(in) :: ucvol 128 real(dp),intent(inout) :: compch_fft 129 type(distribfft_type),optional,intent(in),target :: distribfft 130 type(pawang_type),intent(in) :: pawang 131 !arrays 132 integer,intent(in) :: ngfft(18) 133 integer,optional,target,intent(in) :: mpi_atmtab(:) 134 real(dp),intent(in) :: gprimd(3,3),qphon(3),rprimd(3,3),xred(3,natom) 135 real(dp),intent(out) :: pawgrnhat(cplex*nfft,nspden,3*nhatgrdim) 136 real(dp),intent(inout) :: pawnhat(cplex*nfft,nspden) !vz_i 137 type(pawfgrtab_type),intent(inout) :: pawfgrtab(my_natom) 138 type(pawrhoij_type),intent(in) :: pawrhoij(my_natom),pawrhoij0(my_natom) 139 type(pawtab_type),intent(in) :: pawtab(ntypat) 140 141 !Local variables --------------------------------------- 142 !scalars 143 integer :: cplex_rhoij,iatom,iatom_tot,ic,ierr,ii,ils,ilslm,iq0,irhoij,ispden,itypat 144 integer :: jc,jrhoij,kc,klm,klmn,lmax,lmin,lm_size,mfgd,mm,mpi_comm_sphgrid 145 integer :: my_comm_atom,my_comm_fft,nfgd,nfftot,option,optgr0,optgr1,optgr2,paral_kgb_fft 146 logical :: compute_grad,compute_nhat,my_atmtab_allocated,need_frozen,paral_atom,qeq0 147 logical :: compute_phonons,has_phase 148 type(distribfft_type),pointer :: my_distribfft 149 type(mpi_type) :: mpi_enreg_fft 150 !arrays 151 integer,pointer :: my_atmtab(:) 152 real(dp) :: ro(cplex),ro_ql(cplex),tmp_compch_fft(nspden),tsec(2) 153 real(dp),allocatable :: pawgrnhat_atm(:,:),pawnhat_atm(:),work(:,:) 154 155 ! ************************************************************************* 156 157 DBG_ENTER("COLL") 158 159 compute_nhat=(ider==0.or.ider==2) 160 compute_grad=(ider==1.or.ider==2) 161 compute_phonons=(ipert>0.and.ipert<=natom) 162 163 !Compatibility tests 164 qeq0=(qphon(1)**2+qphon(2)**2+qphon(3)**2<1.d-15) 165 if (present(comm_fft)) then 166 if ((.not.present(paral_kgb)).or.(.not.present(me_g0))) then 167 ABI_BUG('Need paral_kgb and me_g0 with comm_fft !') 168 end if 169 end if 170 if(ider>0.and.nhatgrdim==0) then 171 ABI_BUG('Gradients of nhat required but not allocated!') 172 end if 173 if (my_natom>0) then 174 if(nspden>1.and.nspden/=pawrhoij(1)%nspden) then 175 ABI_BUG('Wrong values for nspden and pawrhoij%nspden!') 176 end if 177 if(nspden>1.and.nspden/=pawfgrtab(1)%nspden) then 178 ABI_BUG('Wrong values for nspden and pawfgrtab%nspden!') 179 end if 180 if(pawrhoij(1)%qphase<cplex) then 181 ABI_BUG('Must have pawrhoij()%qphase >= cplex!') 182 end if 183 if (compute_phonons.and.(.not.qeq0)) then 184 if (pawfgrtab(1)%rfgd_allocated==0) then 185 ABI_BUG('pawfgrtab()%rfgd array must be allocated!') 186 end if 187 if (compute_grad.and.(.not.compute_nhat)) then 188 ABI_BUG('When q<>0, nhat gradients need nhat!') 189 end if 190 end if 191 end if 192 193 !nhat1 does not have to be computed for ddk or d2dk 194 if (ipert==natom+1.or.ipert==natom+10) return 195 196 !Set up parallelism over atoms 197 paral_atom=(present(comm_atom).and.(my_natom/=natom)) 198 nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab 199 my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom 200 call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,& 201 & my_natom_ref=my_natom) 202 203 !Initialisations 204 if ((.not.compute_nhat).and.(.not.compute_grad)) return 205 mfgd=zero;if (my_natom>0) mfgd=maxval(pawfgrtab(1:my_natom)%nfgd) 206 if (compute_nhat) then 207 ABI_MALLOC(pawnhat_atm,(cplex*mfgd)) 208 pawnhat=zero 209 end if 210 if (compute_grad) then 211 ABI_MALLOC(pawgrnhat_atm,(cplex*mfgd,3)) 212 pawgrnhat=zero 213 end if 214 215 !mpi communicators for spherical grid: 216 mpi_comm_sphgrid=xmpi_comm_self !no communicators passed 217 if(present(comm_fft) .and. usewvl==0) mpi_comm_sphgrid=comm_fft 218 if(present(mpi_comm_wvl) .and. usewvl==1) mpi_comm_sphgrid=mpi_comm_wvl 219 220 !------------------------------------------------------------------------ 221 !----- Loop over atoms 222 !------------------------------------------------------------------------ 223 224 do iatom=1,my_natom 225 iatom_tot=iatom;if (paral_atom) iatom_tot=my_atmtab(iatom) 226 227 itypat=pawrhoij(iatom)%itypat 228 lm_size=pawfgrtab(iatom)%l_size**2 229 need_frozen=((compute_nhat).and.(ipert==iatom_tot.or.ipert==natom+3.or.ipert==natom+4)) 230 nfgd=pawfgrtab(iatom)%nfgd 231 cplex_rhoij=pawrhoij(iatom)%cplex_rhoij 232 iq0=cplex_rhoij*pawrhoij(iatom)%lmn2_size 233 234 ! Eventually compute g_l(r).Y_lm(r) factors for the current atom (if not already done) 235 if (((compute_nhat).and.(pawfgrtab(iatom)%gylm_allocated==0)).or.& 236 & ((compute_grad).and.(pawfgrtab(iatom)%gylmgr_allocated==0)).or.& 237 & ((compute_grad.and.need_frozen).and.(pawfgrtab(iatom)%gylmgr2_allocated==0))) then 238 optgr0=0;optgr1=0;optgr2=0 239 if ((compute_nhat).and.(pawfgrtab(iatom)%gylm_allocated==0)) then 240 if (allocated(pawfgrtab(iatom)%gylm)) then 241 ABI_FREE(pawfgrtab(iatom)%gylm) 242 end if 243 ABI_MALLOC(pawfgrtab(iatom)%gylm,(nfgd,pawfgrtab(iatom)%l_size**2)) 244 pawfgrtab(iatom)%gylm_allocated=2;optgr0=1 245 end if 246 if ((compute_grad).and.(pawfgrtab(iatom)%gylmgr_allocated==0)) then 247 if (allocated(pawfgrtab(iatom)%gylmgr)) then 248 ABI_FREE(pawfgrtab(iatom)%gylmgr) 249 end if 250 ABI_MALLOC(pawfgrtab(iatom)%gylmgr,(3,nfgd,pawfgrtab(iatom)%l_size**2)) 251 pawfgrtab(iatom)%gylmgr_allocated=2;optgr1=1 252 end if 253 if ((compute_grad.and.need_frozen).and.(pawfgrtab(iatom)%gylmgr2_allocated==0)) then 254 if (allocated(pawfgrtab(iatom)%gylmgr2)) then 255 ABI_FREE(pawfgrtab(iatom)%gylmgr2) 256 end if 257 ABI_MALLOC(pawfgrtab(iatom)%gylmgr2,(6,nfgd,pawfgrtab(iatom)%l_size**2)) 258 pawfgrtab(iatom)%gylmgr2_allocated=2;optgr2=1 259 end if 260 if (optgr0+optgr1+optgr2>0) then 261 call pawgylm(pawfgrtab(iatom)%gylm,pawfgrtab(iatom)%gylmgr,pawfgrtab(iatom)%gylmgr2,& 262 & lm_size,nfgd,optgr0,optgr1,optgr2,pawtab(itypat),pawfgrtab(iatom)%rfgd) 263 end if 264 end if 265 266 267 ! Eventually compute exp(-i.q.r) factors for the current atom (if not already done) 268 if (compute_phonons.and.(.not.qeq0).and.pawfgrtab(iatom)%expiqr_allocated==0) then 269 if (allocated(pawfgrtab(iatom)%expiqr)) then 270 ABI_FREE(pawfgrtab(iatom)%expiqr) 271 end if 272 ABI_MALLOC(pawfgrtab(iatom)%expiqr,(2,nfgd)) 273 call pawexpiqr(pawfgrtab(iatom)%expiqr,gprimd,nfgd,qphon,& 274 & pawfgrtab(iatom)%rfgd,xred(:,iatom_tot)) 275 pawfgrtab(iatom)%expiqr_allocated=2 276 end if 277 has_phase=(compute_phonons.and.pawfgrtab(iatom)%expiqr_allocated/=0) 278 279 ! Eventually compute frozen part of nhat for the current atom (if not already done) 280 if ((need_frozen).and.((pawfgrtab(iatom)%nhatfr_allocated==0).or.& 281 & (compute_grad.and.pawfgrtab(iatom)%nhatfrgr_allocated==0))) then 282 if (allocated(pawfgrtab(iatom)%nhatfr)) then 283 ABI_FREE(pawfgrtab(iatom)%nhatfr) 284 end if 285 ABI_MALLOC(pawfgrtab(iatom)%nhatfr,(nfgd,pawfgrtab(iatom)%nspden)) 286 option=0;pawfgrtab(iatom)%nhatfr_allocated=2 287 if (compute_grad) then 288 option=1 289 if (allocated(pawfgrtab(iatom)%nhatfrgr)) then 290 ABI_FREE(pawfgrtab(iatom)%nhatfrgr) 291 end if 292 ABI_MALLOC(pawfgrtab(iatom)%nhatfrgr,(3,nfgd,pawfgrtab(iatom)%nspden)) 293 pawfgrtab(iatom)%nhatfrgr_allocated=2 294 end if 295 call pawnhatfr(option,idir,ipert,1,natom,nspden,ntypat,pawang,pawfgrtab(iatom),& 296 & pawrhoij0(iatom),pawtab,rprimd) 297 end if 298 299 ! ------------------------------------------------------------------------ 300 ! ----- Loop over density components 301 ! ------------------------------------------------------------------------ 302 303 do ispden=1,nspden 304 305 if (compute_nhat) pawnhat_atm(1:cplex*nfgd)=zero 306 if (compute_grad) pawgrnhat_atm(1:cplex*nfgd,1:3)=zero 307 308 ! ------------------------------------------------------------------------ 309 ! ----- Loop over ij channels (basis components) 310 ! ------------------------------------------------------------------------ 311 jrhoij=1 312 do irhoij=1,pawrhoij(iatom)%nrhoijsel 313 klmn=pawrhoij(iatom)%rhoijselect(irhoij) 314 klm =pawtab(itypat)%indklmn(1,klmn) 315 lmin=pawtab(itypat)%indklmn(3,klmn) 316 lmax=pawtab(itypat)%indklmn(4,klmn) 317 318 ! Retrieve rhoij 319 if (pawrhoij(iatom)%nspden/=2) then 320 ro(1)=pawrhoij(iatom)%rhoijp(jrhoij,ispden) 321 if (cplex==2) ro(2)=pawrhoij(iatom)%rhoijp(iq0+jrhoij,ispden) 322 else 323 if (ispden==1) then 324 ro(1)=pawrhoij(iatom)%rhoijp(jrhoij,1)+pawrhoij(iatom)%rhoijp(jrhoij,2) 325 if (cplex==2) ro(2)=pawrhoij(iatom)%rhoijp(iq0+jrhoij,1)+pawrhoij(iatom)%rhoijp(iq0+jrhoij,2) 326 else if (ispden==2) then 327 ro(1)=pawrhoij(iatom)%rhoijp(jrhoij,1) 328 if (cplex==2) ro(2)=pawrhoij(iatom)%rhoijp(iq0+jrhoij,1) 329 end if 330 end if 331 ro(1:cplex)=pawtab(itypat)%dltij(klmn)*ro(1:cplex) 332 333 if (compute_nhat) then 334 if (cplex==1) then 335 do ils=lmin,lmax,2 336 do mm=-ils,ils 337 ilslm=ils*ils+ils+mm+1 338 if (pawang%gntselect(ilslm,klm)>0) then 339 ro_ql(1)=ro(1)*pawtab(itypat)%qijl(ilslm,klmn) 340 do ic=1,nfgd 341 pawnhat_atm(ic)=pawnhat_atm(ic)+ro_ql(1)*pawfgrtab(iatom)%gylm(ic,ilslm) 342 end do 343 end if 344 end do 345 end do 346 else 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:2)=ro(1:2)*pawtab(itypat)%qijl(ilslm,klmn) 352 do ic=1,nfgd 353 jc=2*ic-1 354 pawnhat_atm(jc:jc+1)=pawnhat_atm(jc:jc+1)+ro_ql(1:2)*pawfgrtab(iatom)%gylm(ic,ilslm) 355 end do 356 end if 357 end do 358 end do 359 end if 360 end if 361 362 if (compute_grad) then 363 if (cplex==1) then 364 do ils=lmin,lmax,2 365 do mm=-ils,ils 366 ilslm=ils*ils+ils+mm+1 367 if (pawang%gntselect(ilslm,klm)>0) then 368 ro_ql(1)=ro(1)*pawtab(itypat)%qijl(ilslm,klmn) 369 do ic=1,nfgd 370 pawgrnhat_atm(ic,1)=pawgrnhat_atm(ic,1)+ro_ql(1)*pawfgrtab(iatom)%gylmgr(1,ic,ilslm) 371 pawgrnhat_atm(ic,2)=pawgrnhat_atm(ic,2)+ro_ql(1)*pawfgrtab(iatom)%gylmgr(2,ic,ilslm) 372 pawgrnhat_atm(ic,3)=pawgrnhat_atm(ic,3)+ro_ql(1)*pawfgrtab(iatom)%gylmgr(3,ic,ilslm) 373 end do 374 end if 375 end do 376 end do 377 else 378 do ils=lmin,lmax,2 379 do mm=-ils,ils 380 ilslm=ils*ils+ils+mm+1 381 if (pawang%gntselect(ilslm,klm)>0) then 382 ro_ql(1:2)=ro(1:2)*pawtab(itypat)%qijl(ilslm,klmn) 383 do ic=1,nfgd 384 jc=2*ic-1 385 pawgrnhat_atm(jc:jc+1,1)=pawgrnhat_atm(jc:jc+1,1) & 386 & +ro_ql(1:2)*pawfgrtab(iatom)%gylmgr(1,ic,ilslm) 387 pawgrnhat_atm(jc:jc+1,2)=pawgrnhat_atm(jc:jc+1,2) & 388 & +ro_ql(1:2)*pawfgrtab(iatom)%gylmgr(2,ic,ilslm) 389 pawgrnhat_atm(jc:jc+1,3)=pawgrnhat_atm(jc:jc+1,3) & 390 & +ro_ql(1:2)*pawfgrtab(iatom)%gylmgr(3,ic,ilslm) 391 end do 392 end if 393 end do 394 end do 395 end if 396 end if 397 398 ! ------------------------------------------------------------------------ 399 ! ----- End loop over ij channels 400 ! ------------------------------------------------------------------------ 401 jrhoij=jrhoij+cplex_rhoij 402 end do 403 404 ! If RF calculation, add frozen part of 1st-order compensation density 405 if (need_frozen) then 406 if (cplex==1) then 407 do ic=1,nfgd 408 pawnhat_atm(ic)=pawnhat_atm(ic)+pawfgrtab(iatom)%nhatfr(ic,ispden) 409 end do 410 else 411 do ic=1,nfgd 412 jc=2*ic-1 413 pawnhat_atm(jc)=pawnhat_atm(jc)+pawfgrtab(iatom)%nhatfr(ic,ispden) 414 end do 415 end if 416 if (compute_grad) then 417 if (cplex==1) then 418 do ic=1,nfgd 419 pawgrnhat_atm(ic,1)=pawgrnhat_atm(ic,1)+pawfgrtab(iatom)%nhatfrgr(1,ic,ispden) 420 pawgrnhat_atm(ic,2)=pawgrnhat_atm(ic,2)+pawfgrtab(iatom)%nhatfrgr(2,ic,ispden) 421 pawgrnhat_atm(ic,3)=pawgrnhat_atm(ic,3)+pawfgrtab(iatom)%nhatfrgr(3,ic,ispden) 422 end do 423 else 424 do ic=1,nfgd 425 jc=2*ic-1 426 pawgrnhat_atm(jc,1)=pawgrnhat_atm(jc,1)+pawfgrtab(iatom)%nhatfrgr(1,ic,ispden) 427 pawgrnhat_atm(jc,2)=pawgrnhat_atm(jc,2)+pawfgrtab(iatom)%nhatfrgr(2,ic,ispden) 428 pawgrnhat_atm(jc,3)=pawgrnhat_atm(jc,3)+pawfgrtab(iatom)%nhatfrgr(3,ic,ispden) 429 end do 430 end if 431 end if 432 end if 433 434 ! If needed, multiply eventually by exp(-i.q.r) phase 435 if (has_phase) then 436 if (cplex==1) then 437 do ic=1,nfgd 438 pawnhat_atm(ic)=pawnhat_atm(ic)*pawfgrtab(iatom)%expiqr(1,ic) 439 end do 440 else 441 do ic=1,nfgd 442 jc=2*ic-1 443 ro_ql(1)= pawfgrtab(iatom)%expiqr(1,ic) 444 ro_ql(2)=-pawfgrtab(iatom)%expiqr(2,ic) 445 ro(1:2)=pawnhat_atm(jc:jc+1) 446 pawnhat_atm(jc )=ro(1)*ro_ql(1)-ro(2)*ro_ql(2) 447 pawnhat_atm(jc+1)=ro(2)*ro_ql(1)+ro(1)*ro_ql(2) 448 end do 449 end if 450 if (compute_grad) then 451 if (cplex==1) then 452 do ic=1,nfgd 453 pawgrnhat_atm(ic,1:3)=pawgrnhat_atm(ic,1:3)*pawfgrtab(iatom)%expiqr(1,ic) 454 end do 455 else 456 do ic=1,nfgd 457 jc=2*ic-1 458 ! dn^hat(r)/dr_i * exp(-i.q.r) 459 ro_ql(1)= pawfgrtab(iatom)%expiqr(1,ic) 460 ro_ql(2)=-pawfgrtab(iatom)%expiqr(2,ic) 461 do ii=1,3 462 ro(1:2)=pawgrnhat_atm(jc:jc+1,ii) 463 pawgrnhat_atm(jc ,ii)=ro(1)*ro_ql(1)-ro(2)*ro_ql(2) 464 pawgrnhat_atm(jc+1,ii)=ro(2)*ro_ql(1)+ro(1)*ro_ql(2) 465 end do 466 ! -i.q_i * [n^hat(r).exp(-i.q.r)] 467 ro(1:2)=pawnhat_atm(jc:jc+1) 468 do ii=1,3 469 pawgrnhat_atm(jc ,ii)=pawgrnhat_atm(jc ,ii)+qphon(ii)*ro(2) 470 pawgrnhat_atm(jc+1,ii)=pawgrnhat_atm(jc+1,ii)-qphon(ii)*ro(1) 471 end do 472 end do 473 end if 474 end if 475 end if 476 477 ! Add the contribution of the atom to the compensation charge 478 if (compute_nhat) then 479 if (cplex==1) then 480 do ic=1,nfgd 481 kc=pawfgrtab(iatom)%ifftsph(ic) 482 pawnhat(kc,ispden)=pawnhat(kc,ispden)+pawnhat_atm(ic) 483 end do 484 else 485 do ic=1,nfgd 486 jc=2*ic-1;kc=2*pawfgrtab(iatom)%ifftsph(ic)-1 487 pawnhat(kc:kc+1,ispden)=pawnhat(kc:kc+1,ispden)+pawnhat_atm(jc:jc+1) 488 end do 489 end if 490 end if 491 if (compute_grad) then 492 if (cplex==1) then 493 do ic=1,nfgd 494 kc=pawfgrtab(iatom)%ifftsph(ic) 495 pawgrnhat(kc,ispden,1:3)=pawgrnhat(kc,ispden,1:3)+pawgrnhat_atm(ic,1:3) 496 end do 497 else 498 do ic=1,nfgd 499 jc=2*ic-1;kc=2*pawfgrtab(iatom)%ifftsph(ic)-1 500 do ii=1,3 501 pawgrnhat(kc:kc+1,ispden,ii)=pawgrnhat(kc:kc+1,ispden,ii)+pawgrnhat_atm(jc:jc+1,ii) 502 end do 503 end do 504 end if 505 end if 506 ! ------------------------------------------------------------------------ 507 ! ----- End loop over density components 508 ! ------------------------------------------------------------------------ 509 end do 510 511 if (pawfgrtab(iatom)%gylm_allocated==2) then 512 ABI_FREE(pawfgrtab(iatom)%gylm) 513 ABI_MALLOC(pawfgrtab(iatom)%gylm,(0,0)) 514 pawfgrtab(iatom)%gylm_allocated=0 515 end if 516 if (pawfgrtab(iatom)%gylmgr_allocated==2) then 517 ABI_FREE(pawfgrtab(iatom)%gylmgr) 518 ABI_MALLOC(pawfgrtab(iatom)%gylmgr,(0,0,0)) 519 pawfgrtab(iatom)%gylmgr_allocated=0 520 end if 521 if (pawfgrtab(iatom)%gylmgr2_allocated==2) then 522 ABI_FREE(pawfgrtab(iatom)%gylmgr2) 523 ABI_MALLOC(pawfgrtab(iatom)%gylmgr2,(0,0,0)) 524 pawfgrtab(iatom)%gylmgr2_allocated=0 525 end if 526 if (pawfgrtab(iatom)%nhatfr_allocated==2) then 527 ABI_FREE(pawfgrtab(iatom)%nhatfr) 528 ABI_MALLOC(pawfgrtab(iatom)%nhatfr,(0,0)) 529 pawfgrtab(iatom)%nhatfr_allocated=0 530 end if 531 if (pawfgrtab(iatom)%nhatfrgr_allocated==2) then 532 ABI_FREE(pawfgrtab(iatom)%nhatfrgr) 533 ABI_MALLOC(pawfgrtab(iatom)%nhatfrgr,(0,0,0)) 534 pawfgrtab(iatom)%nhatfrgr_allocated=0 535 end if 536 if (pawfgrtab(iatom)%expiqr_allocated==2) then 537 ABI_FREE(pawfgrtab(iatom)%expiqr) 538 ABI_MALLOC(pawfgrtab(iatom)%expiqr,(0,0)) 539 pawfgrtab(iatom)%expiqr_allocated=0 540 end if 541 542 ! ------------------------------------------------------------------------ 543 ! ----- End loop over atoms 544 ! ------------------------------------------------------------------------ 545 end do 546 547 !----- Free some memory 548 if (compute_nhat) then 549 ABI_FREE(pawnhat_atm) 550 end if 551 if (compute_grad) then 552 ABI_FREE(pawgrnhat_atm) 553 end if 554 555 !----- Reduction in case of parallelism 556 if (paral_atom) then 557 call timab(48,1,tsec) 558 if (compute_nhat) then 559 call xmpi_sum(pawnhat,my_comm_atom,ierr) 560 end if 561 if (compute_grad) then 562 call xmpi_sum(pawgrnhat,my_comm_atom,ierr) 563 end if 564 call timab(48,2,tsec) 565 end if 566 567 !----- Avoid unbalanced g-components numerical errors 568 if (izero==1.and.compute_nhat.and.usewvl==0) then 569 ! Create fake mpi_enreg to wrap fourdp 570 if (present(distribfft)) then 571 my_distribfft => distribfft 572 else 573 ABI_MALLOC(my_distribfft,) 574 call init_distribfft_seq(my_distribfft,'f',ngfft(2),ngfft(3),'fourdp') 575 end if 576 call initmpi_seq(mpi_enreg_fft) 577 ABI_FREE(mpi_enreg_fft%distribfft) 578 if (present(comm_fft)) then 579 call set_mpi_enreg_fft(mpi_enreg_fft,comm_fft,my_distribfft,me_g0,paral_kgb) 580 my_comm_fft=comm_fft;paral_kgb_fft=paral_kgb 581 else 582 my_comm_fft=xmpi_comm_self;paral_kgb_fft=0; 583 mpi_enreg_fft%distribfft => my_distribfft 584 end if 585 ! do FFT 586 ABI_MALLOC(work,(2,nfft)) 587 do ispden=1,min(2,nspden) 588 call fourdp(cplex,work,pawnhat(:,ispden),-1,mpi_enreg_fft,nfft,1,ngfft,0) 589 call zerosym(work,2,ngfft(1),ngfft(2),ngfft(3),comm_fft=my_comm_fft,distribfft=my_distribfft) 590 call fourdp(cplex,work,pawnhat(:,ispden),+1,mpi_enreg_fft,nfft,1,ngfft,0) 591 end do 592 ABI_FREE(work) 593 ! Destroy fake mpi_enreg 594 call unset_mpi_enreg_fft(mpi_enreg_fft) 595 if (.not.present(distribfft)) then 596 call destroy_distribfft(my_distribfft) 597 ABI_FREE(my_distribfft) 598 end if 599 end if 600 601 !----- Computation of compensation charge over real space grid 602 if (compute_nhat.and.ipert==0) then 603 nfftot=PRODUCT(ngfft(1:3)) 604 call mean_fftr(pawnhat,tmp_compch_fft,nfft,nfftot,1,mpi_comm_sphgrid) 605 compch_fft = tmp_compch_fft(1) 606 compch_fft=compch_fft*ucvol 607 end if 608 609 !Destroy atom table used for parallelism 610 call free_my_atmtab(my_atmtab,my_atmtab_allocated) 611 612 DBG_EXIT("COLL") 613 614 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)
SOURCE
664 subroutine pawmknhat_psipsi(cprj1,cprj2,ider,izero,my_natom,natom,nfft,ngfft,nhat12_grdim,& 665 & nspinor,ntypat,pawang,pawfgrtab,grnhat12,nhat12,pawtab, & 666 & gprimd,grnhat_12,qphon,xred,atindx,mpi_atmtab,comm_atom,comm_fft,me_g0,paral_kgb,distribfft) ! optional arguments 667 668 implicit none 669 670 !Arguments --------------------------------------------- 671 !scalars 672 integer,intent(in) :: ider,izero,my_natom,natom,nfft,nhat12_grdim,ntypat,nspinor 673 integer,optional,intent(in) :: me_g0,comm_fft,paral_kgb 674 integer,optional,intent(in) :: comm_atom 675 type(distribfft_type),optional,intent(in),target :: distribfft 676 type(pawang_type),intent(in) :: pawang 677 !arrays 678 integer,intent(in) :: ngfft(18) 679 integer,optional,intent(in) ::atindx(natom) 680 integer,optional,target,intent(in) :: mpi_atmtab(:) 681 real(dp),optional, intent(in) ::gprimd(3,3),qphon(3),xred(3,natom) 682 real(dp),intent(out) :: grnhat12(2,nfft,nspinor**2,3*nhat12_grdim) 683 real(dp),optional,intent(out) :: grnhat_12(2,nfft,nspinor**2,3,natom*(ider/3)) 684 real(dp),intent(out) :: nhat12(2,nfft,nspinor**2) 685 type(pawfgrtab_type),intent(inout) :: pawfgrtab(my_natom) 686 type(pawtab_type),intent(in) :: pawtab(ntypat) 687 type(pawcprj_type),intent(in) :: cprj1(natom,nspinor),cprj2(natom,nspinor) 688 689 !Local variables --------------------------------------- 690 !scalars 691 integer :: iatm,iatom,iatom_tot,ic,ierr,ils,ilslm,isp1,isp2,isploop,itypat,jc,klm,klmn 692 integer :: lmax,lmin,lm_size,mm,my_comm_atom,my_comm_fft,optgr0,optgr1,paral_kgb_fft 693 integer :: cplex,ilmn,jlmn,lmn_size,lmn2_size 694 real(dp) :: re_p,im_p 695 logical :: compute_grad,compute_grad1,compute_nhat,my_atmtab_allocated,paral_atom,qeq0,compute_phonon,order 696 type(distribfft_type),pointer :: my_distribfft 697 type(mpi_type) :: mpi_enreg_fft 698 !arrays 699 integer,parameter :: spinor_idxs(2,4)=RESHAPE((/1,1,2,2,1,2,2,1/),(/2,4/)) 700 integer,pointer :: my_atmtab(:) 701 real(dp) :: rdum(1),cpf(2),cpf_ql(2),tsec(2),ro(2),ro_ql(2),nhat12_atm(2,nfft,nspinor**2) 702 real(dp),allocatable :: work(:,:), qijl(:,:) 703 704 ! ************************************************************************* 705 706 DBG_ENTER("COLL") 707 708 !Compatibility tests 709 if (present(comm_fft)) then 710 if ((.not.present(paral_kgb)).or.(.not.present(me_g0))) then 711 ABI_BUG('Need paral_kgb and me_g0 with comm_fft!') 712 end if 713 if (present(paral_kgb)) then 714 if (paral_kgb/=0) then 715 ABI_BUG('paral_kgb/=0 not coded!') 716 end if 717 end if 718 end if 719 if (ider>0.and.nhat12_grdim==0) then 720 ! ABI_BUG('Gradients of nhat required but not allocated !') 721 end if 722 if (nspinor==2) then 723 ABI_BUG('nspinor==2 not coded!') 724 end if 725 726 compute_phonon=.false.;qeq0=.false. 727 if (present(gprimd).and.present(qphon).and.present(xred)) compute_phonon=.true. 728 if (compute_phonon) qeq0=(qphon(1)**2+qphon(2)**2+qphon(3)**2<1.d-15) 729 if (present(atindx)) order=.true. 730 !Set up parallelism over atoms 731 paral_atom=(present(comm_atom).and.(my_natom/=natom)) 732 nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab 733 my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom 734 call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,my_natom_ref=my_natom) 735 736 !Initialisations 737 compute_nhat=(ider==0.or.ider==2.or.ider==3) 738 compute_grad=(ider==1.or.ider==2) 739 compute_grad1=(ider==3) 740 if ((.not.compute_nhat).and.(.not.compute_grad)) return 741 742 if (compute_nhat) nhat12=zero 743 if (compute_grad) grnhat12=zero 744 if (compute_grad1) grnhat_12=zero 745 746 if (compute_grad) then 747 ! ABI_BUG('compute_grad not tested!') 748 end if 749 750 !------------------------------------------------------------------------ 751 !----- Loop over atoms 752 !------------------------------------------------------------------------ 753 do iatom=1,my_natom 754 iatom_tot=iatom;if (paral_atom) iatom_tot=my_atmtab(iatom) 755 iatm=iatom_tot 756 if (order) iatm=atindx(iatom_tot) 757 itypat = pawfgrtab(iatom)%itypat 758 lm_size = pawfgrtab(iatom)%l_size**2 759 lmn_size = pawtab(itypat)%lmn_size 760 lmn2_size = pawtab(itypat)%lmn2_size 761 ABI_MALLOC(qijl,(lm_size,lmn2_size)) 762 qijl=zero 763 qijl=pawtab(itypat)%qijl 764 if (compute_nhat) nhat12_atm=zero 765 766 ! Eventually compute g_l(r).Y_lm(r) factors for the current atom (if not already done) 767 if (((compute_nhat).and.(pawfgrtab(iatom)%gylm_allocated==0)).or.& 768 & (((compute_grad).or.(compute_grad1)).and.(pawfgrtab(iatom)%gylmgr_allocated==0))) then 769 optgr0=0; optgr1=0 770 if ((compute_nhat).and.(pawfgrtab(iatom)%gylm_allocated==0)) then 771 if (allocated(pawfgrtab(iatom)%gylm)) then 772 ABI_FREE(pawfgrtab(iatom)%gylm) 773 end if 774 ABI_MALLOC(pawfgrtab(iatom)%gylm,(pawfgrtab(iatom)%nfgd,pawfgrtab(iatom)%l_size**2)) 775 pawfgrtab(iatom)%gylm_allocated=2;optgr0=1 776 end if 777 if (((compute_grad).or.(compute_grad1)).and.(pawfgrtab(iatom)%gylmgr_allocated==0)) then 778 if (allocated(pawfgrtab(iatom)%gylmgr)) then 779 ABI_FREE(pawfgrtab(iatom)%gylmgr) 780 end if 781 ABI_MALLOC(pawfgrtab(iatom)%gylmgr,(3,pawfgrtab(iatom)%nfgd,pawfgrtab(iatom)%l_size**2)) 782 pawfgrtab(iatom)%gylmgr_allocated=2;optgr1=1 783 end if 784 if (optgr0+optgr1>0) then 785 call pawgylm(pawfgrtab(iatom)%gylm,pawfgrtab(iatom)%gylmgr,rdum,& 786 & lm_size,pawfgrtab(iatom)%nfgd,optgr0,optgr1,0,pawtab(itypat),& 787 & pawfgrtab(iatom)%rfgd) 788 end if 789 790 end if 791 if (compute_phonon.and.(.not.qeq0).and.(pawfgrtab(iatom)%expiqr_allocated==0)) then 792 if (allocated(pawfgrtab(iatom)%expiqr)) then 793 ABI_FREE(pawfgrtab(iatom)%expiqr) 794 end if 795 ABI_MALLOC(pawfgrtab(iatom)%expiqr,(2,pawfgrtab(iatom)%nfgd)) 796 call pawexpiqr(pawfgrtab(iatom)%expiqr,gprimd,pawfgrtab(iatom)%nfgd,qphon,& 797 & pawfgrtab(iatom)%rfgd,xred(:,iatom_tot)) 798 pawfgrtab(iatom)%expiqr_allocated=2 799 end if 800 801 do isploop=1,nspinor**2 ! Loop over density components of the compensation charge. 802 ! TODO Here we might take advantage of symmetry relations between the four components if nspinor==2 803 isp1=spinor_idxs(1,isploop) 804 isp2=spinor_idxs(2,isploop) 805 806 do klmn=1,lmn2_size ! Loop over ij channels of this atom type. 807 klm =pawtab(itypat)%indklmn(1,klmn) 808 lmin=pawtab(itypat)%indklmn(3,klmn) ! abs(il-jl) 809 lmax=pawtab(itypat)%indklmn(4,klmn) ! il+jl 810 ilmn=pawtab(itypat)%indklmn(7,klmn) 811 jlmn=pawtab(itypat)%indklmn(8,klmn) 812 ! call klmn2ijlmn(klmn,lmn_size,ilmn,jlmn) ! This mapping should be stored in pawtab_type 813 814 ! Retrieve the factor due to the PAW projections. 815 re_p = cprj1(iatm,isp1)%cp(1,ilmn) * cprj2(iatm,isp2)%cp(1,jlmn) & 816 & +cprj1(iatm,isp1)%cp(2,ilmn) * cprj2(iatm,isp2)%cp(2,jlmn) & 817 & +cprj1(iatm,isp1)%cp(1,jlmn) * cprj2(iatm,isp2)%cp(1,ilmn) & 818 & +cprj1(iatm,isp1)%cp(2,jlmn) * cprj2(iatm,isp2)%cp(2,ilmn) 819 820 im_p = cprj1(iatm,isp1)%cp(1,ilmn) * cprj2(iatm,isp2)%cp(2,jlmn) & 821 & -cprj1(iatm,isp1)%cp(2,ilmn) * cprj2(iatm,isp2)%cp(1,jlmn) & 822 & +cprj1(iatm,isp1)%cp(1,jlmn) * cprj2(iatm,isp2)%cp(2,ilmn) & 823 & -cprj1(iatm,isp1)%cp(2,jlmn) * cprj2(iatm,isp2)%cp(1,ilmn) 824 825 cpf(1)=re_p*pawtab(itypat)%dltij(klmn)*half 826 cpf(2)=im_p*pawtab(itypat)%dltij(klmn)*half 827 828 if (compute_nhat) then 829 do ils=lmin,lmax,2 ! Sum over (L,M) 830 do mm=-ils,ils 831 ilslm=ils*ils+ils+mm+1 832 if (pawang%gntselect(ilslm,klm)>0) then 833 cpf_ql(1)=cpf(1)*qijl(ilslm,klmn) 834 cpf_ql(2)=cpf(2)*qijl(ilslm,klmn) 835 do ic=1,pawfgrtab(iatom)%nfgd 836 jc=pawfgrtab(iatom)%ifftsph(ic) 837 nhat12_atm(1,jc,isploop)=nhat12_atm(1,jc,isploop)+cpf_ql(1)*pawfgrtab(iatom)%gylm(ic,ilslm) 838 nhat12_atm(2,jc,isploop)=nhat12_atm(2,jc,isploop)+cpf_ql(2)*pawfgrtab(iatom)%gylm(ic,ilslm) 839 end do 840 end if 841 end do 842 end do 843 end if ! compute_nhat 844 845 if (compute_grad) then 846 do ils=lmin,lmax,2 ! Sum over (L,M) 847 do mm=-ils,ils 848 ilslm=ils*ils+ils+mm+1 849 if (pawang%gntselect(ilslm,klm)>0) then 850 cpf_ql(1)=cpf(1)*qijl(ilslm,klmn) 851 cpf_ql(2)=cpf(2)*qijl(ilslm,klmn) 852 do ic=1,pawfgrtab(iatom)%nfgd 853 jc=pawfgrtab(iatom)%ifftsph(ic) 854 grnhat12(1,jc,isploop,1)=grnhat12(1,jc,isploop,1)+cpf_ql(1)*pawfgrtab(iatom)%gylmgr(1,ic,ilslm) 855 grnhat12(1,jc,isploop,2)=grnhat12(1,jc,isploop,2)+cpf_ql(1)*pawfgrtab(iatom)%gylmgr(2,ic,ilslm) 856 grnhat12(1,jc,isploop,3)=grnhat12(1,jc,isploop,3)+cpf_ql(1)*pawfgrtab(iatom)%gylmgr(3,ic,ilslm) 857 858 grnhat12(2,jc,isploop,1)=grnhat12(2,jc,isploop,1)+cpf_ql(2)*pawfgrtab(iatom)%gylmgr(1,ic,ilslm) 859 grnhat12(2,jc,isploop,2)=grnhat12(2,jc,isploop,2)+cpf_ql(2)*pawfgrtab(iatom)%gylmgr(2,ic,ilslm) 860 grnhat12(2,jc,isploop,3)=grnhat12(2,jc,isploop,3)+cpf_ql(2)*pawfgrtab(iatom)%gylmgr(3,ic,ilslm) 861 end do 862 end if 863 end do 864 end do 865 end if ! compute_grad 866 if (compute_grad1) then 867 do ils=lmin,lmax,2 ! Sum over (L,M) 868 do mm=-ils,ils 869 ilslm=ils*ils+ils+mm+1 870 if (pawang%gntselect(ilslm,klm)>0) then 871 cpf_ql(1)=cpf(1)*qijl(ilslm,klmn) 872 cpf_ql(2)=cpf(2)*qijl(ilslm,klmn) 873 do ic=1,pawfgrtab(iatom)%nfgd 874 jc=pawfgrtab(iatom)%ifftsph(ic) 875 grnhat_12(1,jc,isploop,1,iatom)=grnhat_12(1,jc,isploop,1,iatom)+cpf_ql(1)*pawfgrtab(iatom)%gylmgr(1,ic,ilslm) 876 grnhat_12(1,jc,isploop,2,iatom)=grnhat_12(1,jc,isploop,2,iatom)+cpf_ql(1)*pawfgrtab(iatom)%gylmgr(2,ic,ilslm) 877 grnhat_12(1,jc,isploop,3,iatom)=grnhat_12(1,jc,isploop,3,iatom)+cpf_ql(1)*pawfgrtab(iatom)%gylmgr(3,ic,ilslm) 878 879 grnhat_12(2,jc,isploop,1,iatom)=grnhat_12(2,jc,isploop,1,iatom)+cpf_ql(2)*pawfgrtab(iatom)%gylmgr(1,ic,ilslm) 880 grnhat_12(2,jc,isploop,2,iatom)=grnhat_12(2,jc,isploop,2,iatom)+cpf_ql(2)*pawfgrtab(iatom)%gylmgr(2,ic,ilslm) 881 grnhat_12(2,jc,isploop,3,iatom)=grnhat_12(2,jc,isploop,3,iatom)+cpf_ql(2)*pawfgrtab(iatom)%gylmgr(3,ic,ilslm) 882 end do 883 end if 884 end do 885 end do 886 end if ! compute_grad1 887 end do ! klmn (ij channels) 888 ! If needed, multiply eventually by exp(-i.q.r) phase 889 if (compute_nhat) then 890 if(compute_phonon.and.(.not.qeq0).and.pawfgrtab(iatom)%expiqr_allocated/=0) then 891 do ic=1,pawfgrtab(iatom)%nfgd 892 jc=pawfgrtab(iatom)%ifftsph(ic) 893 ro_ql(1)= pawfgrtab(iatom)%expiqr(1,ic) 894 ro_ql(2)= pawfgrtab(iatom)%expiqr(2,ic) 895 ro(1:2)=nhat12_atm(1:2,jc,isploop) 896 nhat12_atm(1,jc,isploop)=ro(1)*ro_ql(1)-ro(2)*ro_ql(2) 897 nhat12_atm(2,jc,isploop)=ro(2)*ro_ql(1)+ro(1)*ro_ql(2) 898 end do 899 end if 900 end if 901 902 if (compute_grad) then 903 if(compute_phonon.and.(.not.qeq0).and.pawfgrtab(iatom)%expiqr_allocated/=0) then 904 do ic=1,pawfgrtab(iatom)%nfgd 905 jc=pawfgrtab(iatom)%ifftsph(ic) 906 ro_ql(1)= pawfgrtab(iatom)%expiqr(1,ic) 907 ro_ql(2)= pawfgrtab(iatom)%expiqr(2,ic) 908 ro(1)=grnhat12(1,jc,isploop,1)-qphon(1)*nhat12_atm(2,jc,isploop) 909 ro(2)=grnhat12(2,jc,isploop,1)+qphon(1)*nhat12_atm(1,jc,isploop) 910 grnhat12(1,jc,isploop,1)=ro(1)*ro_ql(1)-ro(2)*ro_ql(2) 911 grnhat12(2,jc,isploop,1)=ro(2)*ro_ql(1)+ro(1)*ro_ql(2) 912 ro(1)=grnhat12(1,jc,isploop,2)-qphon(2)*nhat12_atm(2,jc,isploop) 913 ro(2)=grnhat12(2,jc,isploop,2)+qphon(2)*nhat12_atm(1,jc,isploop) 914 grnhat12(1,jc,isploop,2)=ro(1)*ro_ql(1)-ro(2)*ro_ql(2) 915 grnhat12(2,jc,isploop,2)=ro(2)*ro_ql(1)+ro(1)*ro_ql(2) 916 ro(1)=grnhat12(1,jc,isploop,3)-qphon(3)*nhat12_atm(2,jc,isploop) 917 ro(2)=grnhat12(2,jc,isploop,3)+qphon(3)*nhat12_atm(1,jc,isploop) 918 grnhat12(1,jc,isploop,3)=ro(1)*ro_ql(1)-ro(2)*ro_ql(2) 919 grnhat12(2,jc,isploop,3)=ro(2)*ro_ql(1)+ro(1)*ro_ql(2) 920 end do 921 end if 922 end if 923 if (compute_grad1) then 924 if(compute_phonon.and.(.not.qeq0).and.pawfgrtab(iatom)%expiqr_allocated/=0) then 925 do ic=1,pawfgrtab(iatom)%nfgd 926 jc=pawfgrtab(iatom)%ifftsph(ic) 927 ro_ql(1)= pawfgrtab(iatom)%expiqr(1,ic) 928 ro_ql(2)= pawfgrtab(iatom)%expiqr(2,ic) 929 ro(1)=grnhat_12(1,jc,isploop,1,iatom) 930 ro(2)=grnhat_12(2,jc,isploop,1,iatom) 931 grnhat_12(1,jc,isploop,1,iatom)=ro(1)*ro_ql(1)-ro(2)*ro_ql(2) 932 grnhat_12(2,jc,isploop,1,iatom)=ro(2)*ro_ql(1)+ro(1)*ro_ql(2) 933 ro(1)=grnhat_12(1,jc,isploop,2,iatom) 934 ro(2)=grnhat_12(2,jc,isploop,2,iatom) 935 grnhat_12(1,jc,isploop,2,iatom)=ro(1)*ro_ql(1)-ro(2)*ro_ql(2) 936 grnhat_12(2,jc,isploop,2,iatom)=ro(2)*ro_ql(1)+ro(1)*ro_ql(2) 937 ro(1)=grnhat_12(1,jc,isploop,3,iatom) 938 ro(2)=grnhat_12(2,jc,isploop,3,iatom) 939 grnhat_12(1,jc,isploop,3,iatom)=ro(1)*ro_ql(1)-ro(2)*ro_ql(2) 940 grnhat_12(2,jc,isploop,3,iatom)=ro(2)*ro_ql(1)+ro(1)*ro_ql(2) 941 end do 942 end if 943 end if 944 end do ! isploop (density components of the compensation charge) 945 ! accumlate nhat12 for all the atoms 946 if (compute_nhat) nhat12=nhat12+nhat12_atm 947 948 if (pawfgrtab(iatom)%gylm_allocated==2) then 949 ABI_FREE(pawfgrtab(iatom)%gylm) 950 ABI_MALLOC(pawfgrtab(iatom)%gylm,(0,0)) 951 pawfgrtab(iatom)%gylm_allocated=0 952 end if 953 if (pawfgrtab(iatom)%gylmgr_allocated==2) then 954 ABI_FREE(pawfgrtab(iatom)%gylmgr) 955 ABI_MALLOC(pawfgrtab(iatom)%gylmgr,(0,0,0)) 956 pawfgrtab(iatom)%gylmgr_allocated=0 957 end if 958 ABI_FREE(qijl) 959 if (pawfgrtab(iatom)%expiqr_allocated==2) then 960 ABI_FREE(pawfgrtab(iatom)%expiqr) 961 ABI_MALLOC(pawfgrtab(iatom)%expiqr,(0,0)) 962 pawfgrtab(iatom)%expiqr_allocated=0 963 end if 964 965 end do ! iatom 966 967 if (compute_grad1) grnhat_12=-grnhat_12 968 969 !----- Reduction in case of parallelism -----! 970 if (paral_atom)then 971 call timab(48,1,tsec) 972 if (compute_nhat) then 973 call xmpi_sum(nhat12,my_comm_atom,ierr) 974 end if 975 if (compute_grad) then 976 call xmpi_sum(grnhat12,my_comm_atom,ierr) 977 end if 978 if (compute_grad1) then 979 call xmpi_sum(grnhat_12,my_comm_atom,ierr) 980 end if 981 call timab(48,2,tsec) 982 end if 983 984 !----- Avoid unbalanced g-components numerical errors -----! 985 986 if (izero==1.and.compute_nhat) then 987 ! Create fake mpi_enreg to wrap fourdp 988 if (present(distribfft)) then 989 my_distribfft => distribfft 990 else 991 ABI_MALLOC(my_distribfft,) 992 call init_distribfft_seq(my_distribfft,'f',ngfft(2),ngfft(3),'fourdp') 993 end if 994 call initmpi_seq(mpi_enreg_fft) 995 ABI_FREE(mpi_enreg_fft%distribfft) 996 if (present(comm_fft)) then 997 call set_mpi_enreg_fft(mpi_enreg_fft,comm_fft,my_distribfft,me_g0,paral_kgb) 998 my_comm_fft=comm_fft;paral_kgb_fft=paral_kgb 999 else 1000 my_comm_fft=xmpi_comm_self;paral_kgb_fft=0; 1001 mpi_enreg_fft%distribfft => my_distribfft 1002 end if 1003 ! Do FFT 1004 ABI_MALLOC(work,(2,nfft)) 1005 cplex=2 1006 do isp1=1,MIN(2,nspinor**2) 1007 call fourdp(cplex,work,nhat12(:,:,isp1),-1,mpi_enreg_fft,nfft,1,ngfft,0) 1008 call zerosym(work,cplex,ngfft(1),ngfft(2),ngfft(3),comm_fft=my_comm_fft,distribfft=my_distribfft) 1009 call fourdp(cplex,work,nhat12(:,:,isp1),+1,mpi_enreg_fft,nfft,1,ngfft,0) 1010 end do 1011 ABI_FREE(work) 1012 ! Destroy fake mpi_enreg 1013 call unset_mpi_enreg_fft(mpi_enreg_fft) 1014 if (.not.present(distribfft)) then 1015 call destroy_distribfft(my_distribfft) 1016 ABI_FREE(my_distribfft) 1017 end if 1018 end if 1019 1020 !Destroy atom table used for parallelism 1021 call free_my_atmtab(my_atmtab,my_atmtab_allocated) 1022 1023 DBG_EXIT("COLL") 1024 1025 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))]
SOURCE
1068 subroutine pawnhatfr(ider,idir,ipert,my_natom,natom,nspden,ntypat,& 1069 & pawang,pawfgrtab,pawrhoij,pawtab,rprimd, & 1070 & mpi_atmtab,comm_atom) ! optional arguments (parallelism) 1071 1072 implicit none 1073 1074 !Arguments ------------------------------------ 1075 !scalars 1076 integer,intent(in) :: ider,idir,ipert,my_natom,natom,nspden,ntypat 1077 integer,optional,intent(in) :: comm_atom 1078 type(pawang_type),intent(in) :: pawang 1079 !arrays 1080 integer,optional,target,intent(in) :: mpi_atmtab(:) 1081 real(dp),intent(in) :: rprimd(3,3) 1082 type(pawfgrtab_type),intent(inout) :: pawfgrtab(my_natom) 1083 type(pawrhoij_type),intent(in) :: pawrhoij(my_natom) 1084 type(pawtab_type),intent(in) :: pawtab(ntypat) 1085 1086 !Local variables------------------------------- 1087 !scalars 1088 integer :: iatom,iatom_tot,ic,ils,ilslm,irhoij,isel,ispden,istr,itypat,jrhoij 1089 integer :: klm,klmn,lm_size,lmn2_size,lm0,lmax,lmin,mua,mub,mm,mu,my_comm_atom,nfgd,nu,optgr0,optgr1,optgr2 1090 logical :: my_atmtab_allocated,my_pert,paral_atom 1091 real(dp) :: contrib,ro 1092 !arrays 1093 integer,parameter :: voigt(3,3)=reshape((/1,6,5,6,2,4,5,4,3/),(/3,3/)) 1094 integer,parameter :: alpha(9)=(/1,2,3,3,3,2,2,1,1/),beta(9)=(/1,2,3,2,1,1,3,3,2/) 1095 integer,pointer :: my_atmtab(:) 1096 real(dp),allocatable :: nhatfr_tmp(:,:),nhatfrgr_tmp(:,:,:) 1097 1098 ! ************************************************************************* 1099 1100 DBG_ENTER("COLL") 1101 1102 !Only relevant for atomic displacement and strain perturbation 1103 if (ipert>natom.and.ipert/=natom+3.and.ipert/=natom+4) return 1104 1105 !Set up parallelism over atoms 1106 paral_atom=(present(comm_atom).and.(my_natom/=natom)) 1107 nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab 1108 my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom 1109 call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,my_natom_ref=my_natom) 1110 1111 !Compatibility tests 1112 if (my_natom>0) then 1113 if ((pawfgrtab(1)%gylm_allocated==0.or.pawfgrtab(1)%gylmgr_allocated==0).and. & 1114 & pawfgrtab(1)%rfgd_allocated==0) then 1115 ABI_BUG('pawnhatfr: pawfgrtab()%rfgd array must be allocated!') 1116 end if 1117 if (pawrhoij(1)%qphase/=1) then 1118 ABI_BUG('pawnhatfr: not supposed to be called with qphase=2!') 1119 end if 1120 end if 1121 1122 my_pert = (ipert<=natom).or.ipert==natom+3.or.ipert==natom+4 1123 1124 !Get correct index of strain pertubation 1125 if (ipert==natom+3) istr = idir 1126 if (ipert==natom+4) istr = idir + 3 1127 1128 !Loops over atoms 1129 do iatom=1,my_natom 1130 iatom_tot=iatom;if (paral_atom) iatom_tot=my_atmtab(iatom) 1131 1132 ! Eventually allocate frozen nhat points 1133 if (my_pert) then 1134 if (pawfgrtab(iatom)%nhatfr_allocated==0) then 1135 if (allocated(pawfgrtab(iatom)%nhatfr)) then 1136 ABI_FREE(pawfgrtab(iatom)%nhatfr) 1137 end if 1138 ABI_MALLOC(pawfgrtab(iatom)%nhatfr,(pawfgrtab(iatom)%nfgd,nspden)) 1139 pawfgrtab(iatom)%nhatfr_allocated=1 1140 end if 1141 if (ider==1.and.pawfgrtab(iatom)%nhatfrgr_allocated==0) then 1142 if (allocated(pawfgrtab(iatom)%nhatfrgr)) then 1143 ABI_FREE(pawfgrtab(iatom)%nhatfrgr) 1144 end if 1145 ABI_MALLOC(pawfgrtab(iatom)%nhatfrgr,(3,pawfgrtab(iatom)%nfgd,nspden)) 1146 pawfgrtab(iatom)%nhatfrgr_allocated=1 1147 end if 1148 end if 1149 1150 ! Select if frozen part of nhat exists for the current perturbation 1151 if ((.not.my_pert).or.(pawfgrtab(iatom)%nhatfr_allocated==0)) cycle 1152 1153 ! Some atom-dependent quantities 1154 itypat=pawfgrtab(iatom)%itypat 1155 lm_size=pawfgrtab(iatom)%l_size**2 1156 lmn2_size=pawtab(itypat)%lmn2_size 1157 1158 ! Eventually compute g_l(r).Y_lm(r) factors for the current atom (if not already done) 1159 nfgd=pawfgrtab(iatom)%nfgd 1160 if ((pawfgrtab(iatom)%gylmgr_allocated==0).or. & 1161 & (pawfgrtab(iatom)%gylmgr2_allocated==0.and.ider==1)) then 1162 optgr0=0;optgr1=0;optgr2=0 1163 if(ipert==natom+3.or.ipert==natom+4)then 1164 if (pawfgrtab(iatom)%gylm_allocated==0) then 1165 if (allocated(pawfgrtab(iatom)%gylm)) then 1166 ABI_FREE(pawfgrtab(iatom)%gylm) 1167 end if 1168 ABI_MALLOC(pawfgrtab(iatom)%gylm,(nfgd,lm_size)) 1169 pawfgrtab(iatom)%gylm_allocated=2;optgr0=1 1170 end if 1171 end if 1172 if (pawfgrtab(iatom)%gylmgr_allocated==0) then 1173 if (allocated(pawfgrtab(iatom)%gylmgr)) then 1174 ABI_FREE(pawfgrtab(iatom)%gylmgr) 1175 end if 1176 ABI_MALLOC(pawfgrtab(iatom)%gylmgr,(3,nfgd,lm_size)) 1177 pawfgrtab(iatom)%gylmgr_allocated=2;optgr1=1 1178 end if 1179 if (ider==1.and.pawfgrtab(iatom)%gylmgr2_allocated==0) then 1180 if (allocated(pawfgrtab(iatom)%gylmgr2)) then 1181 ABI_FREE(pawfgrtab(iatom)%gylmgr2) 1182 end if 1183 ABI_MALLOC(pawfgrtab(iatom)%gylmgr2,(6,nfgd,lm_size)) 1184 pawfgrtab(iatom)%gylmgr2_allocated=2;optgr2=1 1185 end if 1186 if (optgr0+optgr1+optgr2>0) then 1187 call pawgylm(pawfgrtab(iatom)%gylm,pawfgrtab(iatom)%gylmgr,pawfgrtab(iatom)%gylmgr2,& 1188 & lm_size,nfgd,optgr0,optgr1,optgr2,pawtab(itypat),pawfgrtab(iatom)%rfgd) 1189 end if 1190 end if 1191 1192 1193 ! ============ Phonons ==================================== 1194 if (ipert<=natom) then 1195 1196 ! Loop over spin components 1197 do ispden=1,nspden 1198 1199 ABI_MALLOC(nhatfr_tmp,(3,nfgd)) 1200 nhatfr_tmp=zero 1201 if (ider==1) then 1202 ABI_MALLOC(nhatfrgr_tmp,(3,nfgd,3)) 1203 nhatfrgr_tmp=zero 1204 end if 1205 1206 jrhoij=1 1207 do irhoij=1,pawrhoij(iatom)%nrhoijsel 1208 klmn=pawrhoij(iatom)%rhoijselect(irhoij) 1209 klm =pawtab(itypat)%indklmn(1,klmn) 1210 lmin=pawtab(itypat)%indklmn(3,klmn) 1211 lmax=pawtab(itypat)%indklmn(4,klmn) 1212 1213 if (nspden/=2) then 1214 ro=pawrhoij(iatom)%rhoijp(jrhoij,ispden) 1215 else 1216 if (ispden==1) then 1217 ro=pawrhoij(iatom)%rhoijp(jrhoij,1)+pawrhoij(iatom)%rhoijp(jrhoij,2) 1218 else if (ispden==2) then 1219 ro=pawrhoij(iatom)%rhoijp(jrhoij,1) 1220 end if 1221 end if 1222 ro=pawtab(itypat)%dltij(klmn)*ro 1223 1224 do ils=lmin,lmax,2 1225 lm0=ils**2+ils+1 1226 do mm=-ils,ils 1227 ilslm=lm0+mm;isel=pawang%gntselect(lm0+mm,klm) 1228 if (isel>0) then 1229 do ic=1,nfgd 1230 do mu=1,3 1231 contrib=-ro*pawtab(itypat)%qijl(ilslm,klmn)& 1232 & *pawfgrtab(iatom)%gylmgr(mu,ic,ilslm) 1233 nhatfr_tmp(mu,ic)=nhatfr_tmp(mu,ic)+contrib 1234 end do 1235 end do 1236 if (ider==1) then 1237 do ic=1,nfgd 1238 do nu=1,3 1239 do mu=1,3 1240 contrib=-ro*pawtab(itypat)%qijl(ilslm,klmn) & 1241 & *pawfgrtab(iatom)%gylmgr2(voigt(mu,nu),ic,ilslm) 1242 nhatfrgr_tmp(mu,ic,nu)=nhatfrgr_tmp(mu,ic,nu)+contrib 1243 end do 1244 end do 1245 end do 1246 end if 1247 end if 1248 end do 1249 end do 1250 jrhoij=jrhoij+pawrhoij(iatom)%cplex_rhoij 1251 end do 1252 1253 ! Convert from cartesian to reduced coordinates 1254 do ic=1,nfgd 1255 pawfgrtab(iatom)%nhatfr(ic,ispden)= & 1256 & rprimd(1,idir)*nhatfr_tmp(1,ic) & 1257 & +rprimd(2,idir)*nhatfr_tmp(2,ic) & 1258 & +rprimd(3,idir)*nhatfr_tmp(3,ic) 1259 end do 1260 if (ider==1) then 1261 do nu=1,3 1262 do ic=1,nfgd 1263 pawfgrtab(iatom)%nhatfrgr(nu,ic,ispden)= & 1264 & rprimd(1,idir)*nhatfrgr_tmp(1,ic,nu) & 1265 & +rprimd(2,idir)*nhatfrgr_tmp(2,ic,nu) & 1266 & +rprimd(3,idir)*nhatfrgr_tmp(3,ic,nu) 1267 end do 1268 end do 1269 end if 1270 ABI_FREE(nhatfr_tmp) 1271 if (ider==1) then 1272 ABI_FREE(nhatfrgr_tmp) 1273 end if 1274 ! End loop over spin components 1275 end do ! ispden 1276 1277 1278 ! ============ Elastic tensor =============================== 1279 else if (ipert==natom+3.or.ipert==natom+4) then 1280 ! Loop over spin components 1281 pawfgrtab(iatom)%nhatfr(:,:) = zero 1282 do ispden=1,nspden 1283 jrhoij=1 1284 do irhoij=1,pawrhoij(iatom)%nrhoijsel 1285 klmn=pawrhoij(iatom)%rhoijselect(irhoij) 1286 klm =pawtab(itypat)%indklmn(1,klmn) 1287 lmin=pawtab(itypat)%indklmn(3,klmn) 1288 lmax=pawtab(itypat)%indklmn(4,klmn) 1289 if (nspden/=2) then 1290 ro=pawrhoij(iatom)%rhoijp(jrhoij,ispden) 1291 else 1292 if (ispden==1) then 1293 ro=pawrhoij(iatom)%rhoijp(jrhoij,1)+pawrhoij(iatom)%rhoijp(jrhoij,2) 1294 else if (ispden==2) then 1295 ro=pawrhoij(iatom)%rhoijp(jrhoij,1) 1296 end if 1297 end if 1298 ro=pawtab(itypat)%dltij(klmn)*ro 1299 do ils=lmin,lmax,2 1300 lm0=ils**2+ils+1 1301 do mm=-ils,ils 1302 ilslm=lm0+mm;isel=pawang%gntselect(lm0+mm,klm) 1303 if (isel>0) then 1304 ! Sum{[Q_ij_q^LM^(1)]} 1305 do ic=1,nfgd 1306 mua=alpha(istr);mub=beta(istr) 1307 pawfgrtab(iatom)%nhatfr(ic,ispden) = pawfgrtab(iatom)%nhatfr(ic,ispden)+& 1308 & ro*pawtab(itypat)%qijl(ilslm,klmn)*half*(& 1309 & pawfgrtab(iatom)%gylmgr(mua,ic,ilslm)*pawfgrtab(iatom)%rfgd(mub,ic)& 1310 & +pawfgrtab(iatom)%gylmgr(mub,ic,ilslm)*pawfgrtab(iatom)%rfgd(mua,ic)) 1311 end do 1312 ! Add volume contribution 1313 if(istr<=3)then 1314 do ic=1,nfgd 1315 pawfgrtab(iatom)%nhatfr(ic,ispden) = pawfgrtab(iatom)%nhatfr(ic,ispden)+& 1316 & ro*pawtab(itypat)%qijl(ilslm,klmn)*pawfgrtab(iatom)%gylm(ic,ilslm) 1317 end do 1318 end if 1319 if (ider==1) then 1320 ABI_ERROR("nhatgr not implemented for strain perturbationxs") 1321 ! do ic=1,nfgd 1322 ! do nu=1,6 1323 ! do mu=1,6 1324 ! contrib=-ro*pawtab(itypat)%qijl(ilslm,klmn) & 1325 !& *pawfgrtab(iatom)%gylmgr2(voigt(mu,nu),ic,ilslm) 1326 ! nhatfrgr_tmp(mu,ic,nu)=nhatfrgr_tmp(mu,ic,nu)+contrib 1327 ! end do 1328 ! end do 1329 ! end do 1330 end if 1331 end if 1332 end do 1333 end do 1334 jrhoij=jrhoij+pawrhoij(iatom)%cplex_rhoij 1335 end do 1336 end do ! ispden 1337 end if 1338 1339 ! Eventually free temporary space for g_l(r).Y_lm(r) gradients and exp(-i.q.r) 1340 if (pawfgrtab(iatom)%gylmgr_allocated==2) then 1341 ABI_FREE(pawfgrtab(iatom)%gylmgr) 1342 ABI_MALLOC(pawfgrtab(iatom)%gylmgr,(0,0,0)) 1343 pawfgrtab(iatom)%gylmgr_allocated=0 1344 end if 1345 if (pawfgrtab(iatom)%gylmgr2_allocated==2) then 1346 ABI_FREE(pawfgrtab(iatom)%gylmgr2) 1347 ABI_MALLOC(pawfgrtab(iatom)%gylmgr2,(0,0,0)) 1348 pawfgrtab(iatom)%gylmgr2_allocated=0 1349 end if 1350 1351 ! End loop on atoms 1352 end do 1353 1354 !Destroy atom table used for parallelism 1355 call free_my_atmtab(my_atmtab,my_atmtab_allocated) 1356 1357 DBG_EXIT("COLL") 1358 1359 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)
SOURCE
1417 subroutine pawsushat(atindx,cprj_k,gbound_diel,gylmg_diel,iband1,iband2,ispinor1,ispinor2,istwf_k,kg_diel,& 1418 & lmax_diel,mgfftdiel,natom,nband,ndiel4,ndiel5,ndiel6,& 1419 & ngfftdiel,npwdiel,nspinor,ntypat,optreal,& 1420 & pawang,pawtab,ph3d_diel,typat,wfprod,wfraug, & 1421 & mpi_atmtab,comm_atom,comm_fft,me_g0,paral_kgb,distribfft) ! optional arguments (parallelism) 1422 1423 implicit none 1424 1425 !Arguments --------------------------------------------- 1426 !scalars 1427 integer,intent(in) :: iband1,iband2,ispinor1,ispinor2,istwf_k,lmax_diel,mgfftdiel 1428 integer,intent(in) :: natom,nband,ndiel4,ndiel5,ndiel6,npwdiel,nspinor 1429 integer,intent(in) :: ntypat,optreal 1430 integer,optional,intent(in) :: me_g0,comm_atom,comm_fft,paral_kgb 1431 type(distribfft_type),optional,intent(in),target :: distribfft 1432 type(pawang_type),intent(in) :: pawang 1433 !arrays 1434 integer,intent(in) :: atindx(natom),gbound_diel(2*mgfftdiel+8,2) 1435 integer,intent(in) :: kg_diel(3,npwdiel),ngfftdiel(18),typat(natom) 1436 integer,optional,target,intent(in) :: mpi_atmtab(:) 1437 real(dp),intent(in) :: gylmg_diel(npwdiel,lmax_diel**2,ntypat) 1438 real(dp),intent(in) :: ph3d_diel(2,npwdiel,natom) 1439 real(dp),intent(inout) :: wfprod(2,npwdiel*(1-optreal)) 1440 real(dp),intent(inout) :: wfraug(2,ndiel4,ndiel5,ndiel6*optreal) 1441 type(pawcprj_type),intent(in) :: cprj_k(natom,nspinor*nband) 1442 type(pawtab_type),intent(in) :: pawtab(ntypat) 1443 1444 !Local variables --------------------------------------- 1445 !scalars 1446 integer :: cplex,iatm,iatom,iatom_tot,ibsp1,ibsp2,ierr,il,ilmn,ils,ilslm,ipw 1447 integer :: itypat,j0lmn,jlmn,klm,klmn,lmax,lmin,mm,my_comm_atom,my_comm_fft,my_natom,tim_fourwf 1448 real(dp) :: phil1,phil2,sgn,weight_dum,wf1,wf2 1449 logical :: my_atmtab_allocated,parity,paral_atom 1450 type(distribfft_type),pointer :: my_distribfft 1451 type(mpi_type) :: mpi_enreg_fft 1452 !arrays 1453 integer,pointer :: my_atmtab(:) 1454 real(dp) :: ro(2),ro_ql(2) 1455 real(dp),allocatable :: dummy(:,:),wfprod_paw(:,:),wfraug_paw(:,:,:,:) 1456 1457 ! ************************************************************************* 1458 1459 DBG_ENTER("COLL") 1460 1461 if (present(comm_fft)) then 1462 if ((.not.present(paral_kgb)).or.(.not.present(me_g0))) then 1463 ABI_BUG('Need paral_kgb and me_g0 with comm_fft !') 1464 end if 1465 end if 1466 1467 !Set up parallelism over atoms 1468 paral_atom=(present(comm_atom).and.(my_natom/=natom)) 1469 nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab 1470 my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom 1471 call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom) 1472 my_natom=natom;if (paral_atom) my_natom=size(my_atmtab) 1473 1474 cplex=1;if (istwf_k>1) cplex=2 1475 ABI_MALLOC(wfprod_paw,(2,npwdiel)) 1476 wfprod_paw(:,:)=zero 1477 ibsp1=(iband1-1)*nspinor+ispinor1 1478 ibsp2=(iband2-1)*nspinor+ispinor2 1479 1480 !------------------------------------------------------------------------ 1481 !----- Loop over atoms 1482 !------------------------------------------------------------------------ 1483 do iatom=1,my_natom 1484 iatom_tot=iatom;if (paral_atom) iatom_tot=my_atmtab(iatom) 1485 iatm=atindx(iatom_tot) 1486 itypat=typat(iatom_tot) 1487 1488 ! ------------------------------------------------------------------------ 1489 ! ----- Loop over ij channels (basis components) 1490 ! ------------------------------------------------------------------------ 1491 do jlmn=1,pawtab(itypat)%lmn_size 1492 j0lmn=jlmn*(jlmn-1)/2 1493 do ilmn=1,jlmn 1494 klmn=j0lmn+ilmn 1495 klm =pawtab(itypat)%indklmn(1,klmn) 1496 lmin=pawtab(itypat)%indklmn(3,klmn) 1497 lmax=pawtab(itypat)%indklmn(4,klmn) 1498 1499 ro(1)=cprj_k(iatm,ibsp1)%cp(1,ilmn)*cprj_k(iatm,ibsp2)%cp(1,jlmn) 1500 if (cplex==2) then 1501 ro(1)=ro(1)+cprj_k(iatm,ibsp1)%cp(2,ilmn)*cprj_k(iatm,ibsp2)%cp(2,jlmn) 1502 ro(2)=cprj_k(iatm,ibsp1)%cp(2,ilmn)*cprj_k(iatm,ibsp2)%cp(1,jlmn) & 1503 & -cprj_k(iatm,ibsp1)%cp(1,ilmn)*cprj_k(iatm,ibsp2)%cp(2,jlmn) 1504 end if 1505 ro(1:cplex)=ro(1:cplex)*pawtab(itypat)%dltij(klmn) 1506 1507 do ils=lmin,lmax,2 1508 il=mod(ils,4);parity=(mod(il,2)==0) 1509 sgn=one;if (il>1) sgn=-one 1510 1511 do mm=-ils,ils 1512 ilslm=ils*ils+ils+mm+1 1513 if (pawang%gntselect(ilslm,klm)>0) then 1514 1515 ro_ql(1:cplex)=pawtab(itypat)%qijl(ilslm,klmn)*ro(1:cplex) 1516 1517 ! Compute: Sum_{ijR} [ cpi* cpj qij^l (-i)^l g_l(g) S_lm(g) ] 1518 1519 if (cplex==1) then 1520 if (parity) then 1521 do ipw=1,npwdiel 1522 phil1= sgn*ph3d_diel(1,ipw,iatm) ! (i)^l.exp(i.g.R) 1523 phil2= sgn*ph3d_diel(2,ipw,iatm) 1524 wf1= phil1*ro_ql(1) ! cpi* cpj qij^l (-i)^l.exp(-i.g.R) 1525 wf2=-phil2*ro_ql(1) 1526 wfprod_paw(1,ipw)=wfprod_paw(1,ipw)+wf1*gylmg_diel(ipw,ilslm,itypat) 1527 wfprod_paw(2,ipw)=wfprod_paw(2,ipw)+wf2*gylmg_diel(ipw,ilslm,itypat) 1528 end do 1529 else 1530 do ipw=1,npwdiel 1531 phil1=-sgn*ph3d_diel(2,ipw,iatm) ! (i)^l.exp(i.g.R) 1532 phil2= sgn*ph3d_diel(1,ipw,iatm) 1533 wf1= phil1*ro_ql(1) ! cpi* cpj qij^l (-i)^l.exp(-i.g.R) 1534 wf2=-phil2*ro_ql(1) 1535 wfprod_paw(1,ipw)=wfprod_paw(1,ipw)+wf1*gylmg_diel(ipw,ilslm,itypat) 1536 wfprod_paw(2,ipw)=wfprod_paw(2,ipw)+wf2*gylmg_diel(ipw,ilslm,itypat) 1537 end do 1538 end if 1539 1540 else 1541 1542 if (parity) then 1543 do ipw=1,npwdiel 1544 phil1= sgn*ph3d_diel(1,ipw,iatm) ! (i)^l.exp(i.g.R) 1545 phil2= sgn*ph3d_diel(2,ipw,iatm) 1546 wf1=phil1*ro_ql(1)+phil2*ro_ql(2) ! cpi* cpj qij^l (-i)^l.exp(-i.g.R) 1547 wf2=phil1*ro_ql(2)-phil2*ro_ql(1) 1548 wfprod_paw(1,ipw)=wfprod_paw(1,ipw)+wf1*gylmg_diel(ipw,ilslm,itypat) 1549 wfprod_paw(2,ipw)=wfprod_paw(2,ipw)+wf2*gylmg_diel(ipw,ilslm,itypat) 1550 end do 1551 else 1552 do ipw=1,npwdiel 1553 phil1=-sgn*ph3d_diel(2,ipw,iatm) ! (i)^l.exp(i.g.R) 1554 phil2= sgn*ph3d_diel(1,ipw,iatm) 1555 wf1=phil1*ro_ql(1)+phil2*ro_ql(2) ! cpi* cpj qij^l (-i)^l.exp(-i.g.R) 1556 wf2=phil1*ro_ql(2)-phil2*ro_ql(1) 1557 wfprod_paw(1,ipw)=wfprod_paw(1,ipw)+wf1*gylmg_diel(ipw,ilslm,itypat) 1558 wfprod_paw(2,ipw)=wfprod_paw(2,ipw)+wf2*gylmg_diel(ipw,ilslm,itypat) 1559 end do 1560 end if 1561 1562 end if 1563 end if 1564 end do 1565 end do 1566 1567 ! ----- End loop over ij channels 1568 end do 1569 end do 1570 1571 ! ----- End loop over atoms 1572 end do 1573 1574 !Reduction in case of parallelism over atoms 1575 if (paral_atom) then 1576 call xmpi_sum(wfprod_paw,my_comm_atom,ierr) 1577 end if 1578 1579 if (optreal==0) then 1580 1581 ! === Output in reciprocal space 1582 wfprod(:,:)=wfprod(:,:)+wfprod_paw(:,:) 1583 1584 else 1585 ! === Output in reciprocal space 1586 tim_fourwf=17;weight_dum=0 1587 ! Create fake mpi_enreg to wrap fourdp 1588 if (present(distribfft)) then 1589 my_distribfft => distribfft 1590 else 1591 ABI_MALLOC(my_distribfft,) 1592 call init_distribfft_seq(my_distribfft,'c',ngfftdiel(2),ngfftdiel(3),'fourwf') 1593 end if 1594 call initmpi_seq(mpi_enreg_fft) 1595 ABI_FREE(mpi_enreg_fft%distribfft) 1596 if (present(comm_fft)) then 1597 call set_mpi_enreg_fft(mpi_enreg_fft,comm_fft,my_distribfft,me_g0,paral_kgb) 1598 my_comm_fft=comm_fft 1599 mpi_enreg_fft%paral_kgb = paral_kgb 1600 else 1601 my_comm_fft=xmpi_comm_self 1602 mpi_enreg_fft%paral_kgb = 0 1603 mpi_enreg_fft%distribfft => my_distribfft 1604 end if 1605 ! do FFT 1606 ABI_MALLOC(wfraug_paw,(2,ndiel4,ndiel5,ndiel6)) 1607 call fourwf(1,dummy,wfprod_paw,dummy,wfraug_paw,gbound_diel,gbound_diel,& 1608 & istwf_k,kg_diel,kg_diel,mgfftdiel,mpi_enreg_fft,1,ngfftdiel,1,npwdiel,& 1609 & ndiel4,ndiel5,ndiel6,0,tim_fourwf,weight_dum,weight_dum) 1610 wfraug(:,:,:,:)=wfraug(:,:,:,:)+wfraug_paw(:,:,:,:) 1611 ABI_FREE(wfraug_paw) 1612 call unset_mpi_enreg_fft(mpi_enreg_fft) 1613 if (.not.present(distribfft)) then 1614 call destroy_distribfft(my_distribfft) 1615 ABI_FREE(my_distribfft) 1616 end if 1617 end if 1618 1619 ABI_FREE(wfprod_paw) 1620 1621 !Destroy atom table used for parallelism 1622 call free_my_atmtab(my_atmtab,my_atmtab_allocated) 1623 1624 DBG_EXIT("COLL") 1625 1626 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
SOURCE
1917 subroutine wvl_nhatgrid(atindx1,geocode,h,i3s,natom,natom_tot,& 1918 & nattyp,ntypat,n1,n1i,n2,n2i,n3,n3pi,optcut,optgr0,optgr1,optgr2,optrad,& 1919 & pawfgrtab,pawtab,psppar,rprimd,shift,xred) 1920 1921 implicit none 1922 1923 !Arguments --------------------------------------------- 1924 !scalars 1925 integer,intent(in) :: i3s,natom,natom_tot,ntypat,optcut,optgr0,optgr1,optgr2,optrad 1926 integer,intent(in) :: n1,n2,n3,n1i,n2i,n3pi,shift 1927 real(dp),intent(in) :: h(3) 1928 character(1),intent(in) :: geocode 1929 !integer,intent(in),optional :: mpi_comm_wvl 1930 !arrays 1931 integer,intent(in) :: atindx1(natom),nattyp(ntypat) 1932 real(dp),intent(in) :: psppar(0:4,0:6,ntypat),rprimd(3,3) 1933 real(dp),intent(inout) :: xred(3,natom) 1934 type(pawfgrtab_type),intent(inout) :: pawfgrtab(natom) 1935 type(pawtab_type),intent(in) :: pawtab(ntypat) 1936 1937 !Local variables ------------------------------ 1938 !scalars 1939 !buffer to be added at the end of the last dimension of an array to control bounds_check 1940 integer :: iat,iatm,iatom,iatom_tot,itypat,lm_size,nfgd 1941 real(dp) :: rloc,rshp,xcart(3,natom) 1942 !arrays 1943 integer,allocatable :: ifftsph_tmp(:) 1944 real(dp) :: hh(3) !fine grid spacing for wavelets 1945 real(dp) :: tsec(2) 1946 real(dp),allocatable :: rfgd_tmp(:,:) 1947 1948 ! ************************************************************************* 1949 1950 DBG_ENTER("COLL") 1951 1952 #if !defined HAVE_BIGDFT 1953 BIGDFT_NOTENABLED_ERROR() 1954 #endif 1955 1956 call timab(559,1,tsec) 1957 1958 !Set up parallelism for wvl 1959 !for debug: use me_wvl=xmpi_comm_rank(MPI_COMM_WORLD) 1960 !if (present(mpi_comm_wvl)) then 1961 !me_wvl=xmpi_comm_rank(mpi_comm_wvl) 1962 !nproc_wvl=xmpi_comm_size(mpi_comm_wvl) 1963 !else 1964 !me_wvl=0;nproc_wvl=1 1965 !end if 1966 !Pending: parallelism over atoms: see nhatgrid 1967 1968 if (natom_tot<natom) then ! This test has to be remove when natom_tot is used 1969 ABI_BUG(' natom_tot<natom !') 1970 end if 1971 1972 !Fine grid 1973 hh(:)=0.5d0*h(:) 1974 1975 !Compute xcart from xred 1976 call xred2xcart(natom,rprimd,xcart,xred) 1977 1978 !Loop over types of atom 1979 iatm=0 1980 do itypat=1,ntypat 1981 1982 rloc=psppar(0,0,itypat) 1983 if (optcut==1) then 1984 rshp=pawtab(itypat)%rpaw 1985 else 1986 rshp=pawtab(itypat)%rshp 1987 end if 1988 1989 ! Loop over atoms 1990 do iat=1,nattyp(itypat) 1991 iatm=iatm+1;iatom=atindx1(iatm) 1992 iatom_tot=iatom; !if (paral_atom) iatom_tot=my_atmtab(iatom) 1993 lm_size=pawfgrtab(iatom)%l_size**2 1994 1995 ! Determine FFT points and r-R vectors around the atom 1996 call pawrfgd_wvl(geocode,hh,ifftsph_tmp,i3s,n1,n1i,n2,n2i,n3,n3pi,nfgd,rshp,rloc,& 1997 & rfgd_tmp,shift,xcart(:,iatom_tot)) 1998 1999 ! Allocate arrays defining sphere (and related data) around current atom 2000 if (allocated(pawfgrtab(iatom)%ifftsph)) then 2001 ABI_FREE(pawfgrtab(iatom)%ifftsph) 2002 end if 2003 ABI_MALLOC(pawfgrtab(iatom)%ifftsph,(nfgd)) 2004 pawfgrtab(iatom)%nfgd=nfgd 2005 pawfgrtab(iatom)%ifftsph(1:nfgd)=ifftsph_tmp(1:nfgd) 2006 2007 if (optrad==1) then 2008 if (allocated(pawfgrtab(iatom)%rfgd)) then 2009 ABI_FREE(pawfgrtab(iatom)%rfgd) 2010 end if 2011 ABI_MALLOC(pawfgrtab(iatom)%rfgd,(3,nfgd)) 2012 pawfgrtab(iatom)%rfgd_allocated=1 2013 pawfgrtab(iatom)%rfgd(1:3,1:nfgd)=rfgd_tmp(1:3,1:nfgd) 2014 end if 2015 2016 if (optgr0==1) then 2017 if (allocated(pawfgrtab(iatom)%gylm)) then 2018 ABI_FREE(pawfgrtab(iatom)%gylm) 2019 end if 2020 ABI_MALLOC(pawfgrtab(iatom)%gylm,(nfgd,lm_size)) 2021 pawfgrtab(iatom)%gylm_allocated=1 2022 end if 2023 2024 if (optgr1==1) then 2025 if (allocated(pawfgrtab(iatom)%gylmgr)) then 2026 ABI_FREE(pawfgrtab(iatom)%gylmgr) 2027 end if 2028 ABI_MALLOC(pawfgrtab(iatom)%gylmgr,(3,nfgd,lm_size)) 2029 pawfgrtab(iatom)%gylmgr_allocated=1 2030 end if 2031 2032 if (optgr2==1) then 2033 if (allocated(pawfgrtab(iatom)%gylmgr2)) then 2034 ABI_FREE(pawfgrtab(iatom)%gylmgr2) 2035 end if 2036 ABI_MALLOC(pawfgrtab(iatom)%gylmgr2,(6,nfgd,lm_size)) 2037 pawfgrtab(iatom)%gylmgr2_allocated=1 2038 end if 2039 2040 ! Calculate g_l(r-R)*Y_lm(r-R) for each r around the atom R 2041 if (optgr0+optgr1+optgr2>0) then 2042 call pawgylm(pawfgrtab(iatom)%gylm,pawfgrtab(iatom)%gylmgr,pawfgrtab(iatom)%gylmgr2,& 2043 & lm_size,nfgd,optgr0,optgr1,optgr2,pawtab(itypat),rfgd_tmp(:,1:nfgd)) 2044 end if 2045 2046 ! End loops over types/atoms 2047 ABI_FREE(ifftsph_tmp) 2048 ABI_FREE(rfgd_tmp) 2049 end do 2050 end do 2051 2052 call timab(559,2,tsec) 2053 2054 DBG_EXIT("COLL") 2055 2056 end subroutine wvl_nhatgrid