TABLE OF CONTENTS
- ABINIT/m_paw_denpot
- m_paw_denpot/paw_mknewh0
- m_paw_denpot/pawaccenergy
- m_paw_denpot/pawaccenergy_nospin
- m_paw_denpot/pawdenpot
- m_paw_denpot/pawdensities
- m_paw_denpot/pawkindensities
ABINIT/m_paw_denpot [ Modules ]
NAME
m_paw_denpot
FUNCTION
This module contains routines related to PAW on-site densities and on-site potentials.
COPYRIGHT
Copyright (C) 2018-2024 ABINIT group (FJ, MT) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
SOURCE
16 #if defined HAVE_CONFIG_H 17 #include "config.h" 18 #endif 19 20 #include "abi_common.h" 21 22 MODULE m_paw_denpot 23 24 use defs_basis 25 use m_abicore 26 use m_errors 27 use m_xmpi 28 use m_time, only : timab 29 30 use m_pawang, only : pawang_type 31 use m_pawrad, only : pawrad_type,pawrad_deducer0,poisson,simp_gen 32 use m_pawtab, only : pawtab_type 33 use m_paw_an, only : paw_an_type 34 use m_paw_ij, only : paw_ij_type 35 use m_pawfgrtab, only : pawfgrtab_type 36 use m_pawrhoij, only : pawrhoij_type 37 use m_pawdij, only : pawdijhartree,pawdiju_euijkl,pawdijnd,pawdijso,pawxpot,pawdijfock,symdij,symdij_all 38 use m_pawxc, only : pawxc,pawxc_dfpt,pawxcm,pawxcm_dfpt,pawxcpositron,pawxcmpositron, & 39 & pawxc_get_usekden,pawxc_is_tb09 40 use m_paw_finegrid, only : pawgylm 41 use m_paral_atom, only : get_my_atmtab,free_my_atmtab 42 use m_paw_correlations, only : pawuenergy,pawxenergy,setnoccmmp 43 use m_paral_atom, only : get_my_atmtab,free_my_atmtab 44 45 use m_crystal, only : crystal_t 46 use m_electronpositron, only : electronpositron_type,electronpositron_calctype 47 48 #ifdef HAVE_FC_ISO_C_BINDING 49 use, intrinsic :: iso_c_binding, only : c_ptr,c_loc,c_f_pointer 50 #endif 51 52 implicit none 53 54 private 55 56 !public procedures. 57 public :: pawdenpot ! Compute different (PAW) energies, densities and potentials inside PAW spheres 58 public :: pawdensities ! Compute PAW on-site densities (all-electron, pseudo and compensation) 59 public :: pawkindensities ! Compute PAW on-site kinetic energy densities (all-electron, pseudo) 60 public :: pawaccenergy ! Accumulate the atomic contribution of a PAW on-site energy 61 public :: pawaccenergy_nospin ! As pawaccenergy, but with no spin polarization 62 public :: paw_mknewh0 ! Compute bare PAW on-site Hamiltonian (-> GW calculations) 63 64 CONTAINS !========================================================================================
m_paw_denpot/paw_mknewh0 [ Functions ]
[ Top ] [ m_paw_denpot ] [ Functions ]
NAME
paw_mknewh0
FUNCTION
Calculates the new bare PAW Hamiltonian in the case of quasi-particle self-consistent GW calculations.
INPUTS
mpi_atmtab(:)=--optional-- indexes of the atoms treated by current proc comm_atom=--optional-- MPI communicator over atoms my_natom=number of atoms treated by current processor nsppol=1 for unpolarized, 2 for spin-polarized nspden=number of spin-density components nfftf=(effective) number of FFT grid points (for this proc) for the "fine" grid pawspnorb=flag: 1 if spin-orbit coupling is activated pawprtvol=control print volume and debugging output for PAW Cryst<crystal_t>=Info on unit cell and its symmetries Pawtab(ntypat*usepaw)<type(pawtab_type)>=paw tabulated starting data Paw_an(natom) <type(paw_an_type)>=paw arrays given on angular mesh Pawang<type(pawang_type)>=paw angular mesh and related data Pawfgrtab(natom) <type(pawfgrtab_type)>=atomic data given on fine rectangular grid vxc(nfftf,nspden)=exchange-correlation potential vxc_val(nfftf,nspden)=valence only exchange-correlation potential vtrial(nfftf,nspden)=potential (Hartree+XC+loc)
SIDE EFFECTS
Paw_ij(natom*usepaw)<Paw_ij_type)>=paw arrays given on (i,j) channels
At output: new value for Paw_ij()%dij
SOURCE
2076 subroutine paw_mknewh0(my_natom,nsppol,nspden,nfftf,pawspnorb,pawprtvol,Cryst,& 2077 & Pawtab,Paw_an,Paw_ij,Pawang,Pawfgrtab,vxc,vxc_val,vtrial,& 2078 & mpi_atmtab,comm_atom) ! optional arguments (parallelism) 2079 2080 !Arguments ------------------------------------ 2081 !scalars 2082 integer,intent(in) :: my_natom,nsppol,nspden,nfftf,pawprtvol,pawspnorb 2083 integer,optional,intent(in) :: comm_atom 2084 !arrays 2085 integer,optional,target,intent(in) :: mpi_atmtab(:) 2086 real(dp),intent(in) :: vxc(nfftf,nspden),vxc_val(nfftf,nspden),vtrial(nfftf,nspden) 2087 type(crystal_t),intent(in) :: Cryst 2088 type(Pawang_type),intent(in) :: Pawang 2089 type(Pawtab_type),target,intent(in) :: Pawtab(Cryst%ntypat) 2090 type(Paw_an_type),intent(in) :: Paw_an(my_natom) 2091 type(Paw_ij_type),intent(inout) :: Paw_ij(my_natom) 2092 type(Pawfgrtab_type),intent(inout) :: Pawfgrtab(my_natom) 2093 2094 !Local variables------------------------------- 2095 !scalars 2096 integer,parameter :: ipert0=0 2097 integer :: iat,iat_tot,idij,ndij,option_dij 2098 integer :: itypat,lmn_size,j0lmn,jlmn,ilmn,klmn,klmn1,klm 2099 integer :: lmin,lmax,mm,isel,lm_size,lmn2_size,my_comm_atom,cplex_dij 2100 integer :: ils,ilslm,ic,lm0 2101 integer :: nsploop,is2fft,qphase 2102 real(dp) :: gylm,qijl 2103 logical :: ltest,my_atmtab_allocated,paral_atom 2104 character(len=500) :: msg 2105 !arrays 2106 integer, pointer :: indklmn_(:,:) 2107 integer,pointer :: my_atmtab(:) 2108 real(dp) :: rdum(1),rdum2(1) 2109 real(dp),allocatable :: prod_hloc(:,:),prodhxc_core(:,:) 2110 real(dp),allocatable :: dijhl_hat(:,:),dijhmxc_val(:,:) 2111 2112 ! ************************************************************************* 2113 2114 DBG_ENTER("COLL") 2115 2116 call wrtout(std_out,'Assembling PAW strengths for the bare Hamiltonian','COLL') 2117 2118 !== Set up parallelism over atoms === 2119 paral_atom=(present(comm_atom).and.(my_natom/=Cryst%natom)) 2120 nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab 2121 my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom 2122 call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,Cryst%natom,my_natom_ref=my_natom) 2123 2124 if (my_natom>0) then 2125 2126 ! === Test if required pointers in paw_ij are allocated === 2127 ltest = (allocated(Paw_ij(1)%dijxc).and.allocated(Paw_ij(1)%dijxc_val) ) 2128 ABI_CHECK(ltest,'dijxc or dijxc_val not calculated') 2129 2130 ltest=(allocated(Paw_ij(1)%dijhat)) !.and.Paw_ij(1)%has_dijhat==2) 2131 ABI_CHECK(ltest,'dijhat not calculated') 2132 2133 ltest=(allocated(Paw_ij(1)%dijhartree)) !.and.Paw_ij(1)%has_dijhartree==2) 2134 ABI_CHECK(ltest,'dijhartree not calculated') 2135 2136 if (ANY(Pawtab(:)%usepawu/=0)) then 2137 do iat=1,my_natom 2138 iat_tot=iat;if (paral_atom) iat_tot=my_atmtab(iat) 2139 itypat=Cryst%typat(iat_tot) 2140 if (Pawtab(itypat)%usepawu/=0) then 2141 ltest=(allocated(Paw_ij(iat)%dijU) ) !.and.Paw_ij(iat)%has_dijU==2) 2142 write(msg,'(a,i3,a)')" For atom no. ",iat," %dijU(iat) has not been calculated." 2143 ABI_CHECK(ltest,msg) 2144 end if 2145 end do 2146 end if 2147 2148 if (pawspnorb>0) then 2149 do iat=1,my_natom 2150 ltest=(allocated(Paw_ij(iat)%dijso) ) !.and.Paw_ij(iat)%has_dijso==2) 2151 write(msg,'(a,i3,a)')" For atom no. ",iat," %dijso(iat) has not been calculated." 2152 ABI_CHECK(ltest,msg) 2153 end do 2154 end if 2155 end if ! my_natom>0 2156 2157 !== Construct the new PAW H0 Hamiltonian === 2158 do iat=1,my_natom 2159 iat_tot=iat;if (paral_atom) iat_tot=my_atmtab(iat) 2160 2161 itypat = Cryst%typat(iat_tot) 2162 lmn_size = Pawtab(itypat)%lmn_size 2163 lmn2_size = Pawtab(itypat)%lmn2_size 2164 lm_size = Paw_an(iat)%lm_size 2165 cplex_dij = Paw_ij(iat)%cplex_dij 2166 qphase = Paw_ij(iat)%qphase 2167 ndij = Paw_ij(iat)%ndij 2168 2169 ABI_CHECK(cplex_dij==1,'cplex_dij/=1 not implemented') 2170 ABI_CHECK(qphase==1,'qphase/=1 not implemented') 2171 ! 2172 ! Eventually compute g_l(r).Y_lm(r) factors for the current atom (if not already done) 2173 if (Pawfgrtab(iat)%gylm_allocated==0) then 2174 if (allocated(Pawfgrtab(iat)%gylm)) then 2175 ABI_FREE(Pawfgrtab(iat)%gylm) 2176 end if 2177 ABI_MALLOC(Pawfgrtab(iat)%gylm,(Pawfgrtab(iat)%nfgd,lm_size)) 2178 Pawfgrtab(iat)%gylm_allocated=2 2179 2180 call pawgylm(Pawfgrtab(iat)%gylm,rdum,rdum2,lm_size,& 2181 & Pawfgrtab(iat)%nfgd,1,0,0,Pawtab(itypat),Pawfgrtab(iat)%rfgd) 2182 end if 2183 2184 ! === Calculate LM contribution to dijhmxc_val for this atom === 2185 ! * Dijxc contains also the Hat term on the FFT mesh while Dijxc_val does not 2186 ! contain neither the hat term nor the LM sum of onsite terms (they should cancel each other) 2187 ! FIXME change paw_dij, otherwise I miss tnc in vxc 2188 ! * prodhxc_core is used to assemble $\int g_l Ylm (vtrial - vxc_val[tn+nhat] dr$ on the FFT mesh === 2189 ! * The following quantities do not depend on ij 2190 ABI_MALLOC(prod_hloc ,(lm_size,ndij)) 2191 ABI_MALLOC(prodhxc_core,(lm_size,ndij)) 2192 prod_hloc =zero 2193 prodhxc_core=zero 2194 do idij=1,ndij 2195 do ilslm=1,lm_size 2196 do ic=1,Pawfgrtab(iat)%nfgd 2197 is2fft=Pawfgrtab(iat)%ifftsph(ic) 2198 gylm=Pawfgrtab(iat)%gylm(ic,ilslm) 2199 prod_hloc (ilslm,idij)=prod_hloc (ilslm,idij) + (vtrial(is2fft,idij)-vxc(is2fft,idij))*gylm 2200 ! prodhxc_core(ilslm,idij)=prodhxc_core(ilslm,idij) + (vxc_val(is2fft,idij))*gylm 2201 prodhxc_core(ilslm,idij)=prodhxc_core(ilslm,idij) + (vtrial(is2fft,idij)-vxc_val(is2fft,idij))*gylm 2202 end do 2203 end do 2204 end do !idij 2205 2206 ! === Assembly the "Hat" contribution for this atom ==== 2207 ABI_MALLOC(dijhl_hat ,(cplex_dij*lmn2_size,ndij)) 2208 ABI_MALLOC(dijhmxc_val,(cplex_dij*lmn2_size,ndij)) 2209 dijhl_hat =zero 2210 dijhmxc_val=zero 2211 indklmn_ => Pawtab(itypat)%indklmn(1:6,1:lmn2_size) 2212 2213 do idij=1,ndij 2214 do klmn=1,lmn2_size 2215 klm =indklmn_(1,klmn) 2216 lmin=indklmn_(3,klmn) 2217 lmax=indklmn_(4,klmn) 2218 2219 ! === $\sum_lm q_ij^l prod* for each idij$ === 2220 do ils=lmin,lmax,2 2221 lm0=ils**2+ils+1 2222 do mm=-ils,ils 2223 ilslm=lm0+mm 2224 isel=Pawang%gntselect(lm0+mm,klm) 2225 if (isel>0) then 2226 qijl=Pawtab(itypat)%qijl(ilslm,klmn) 2227 dijhl_hat (klmn,idij)=dijhl_hat (klmn,idij) + prod_hloc (ilslm,idij)*qijl 2228 dijhmxc_val(klmn,idij)=dijhmxc_val(klmn,idij) +prodhxc_core(ilslm,idij)*qijl 2229 end if 2230 end do 2231 end do 2232 end do 2233 end do 2234 2235 ABI_FREE(prod_hloc) 2236 ABI_FREE(prodhxc_core) 2237 2238 ! * Normalization factor due to integration on the FFT mesh 2239 dijhl_hat = dijhl_hat *Cryst%ucvol/DBLE(nfftf) 2240 dijhmxc_val= dijhmxc_val*Cryst%ucvol/DBLE(nfftf) 2241 2242 ! === Now assembly the bare Hamiltonian === 2243 ! * Loop over density components overwriting %dij 2244 nsploop=nsppol; if (Paw_ij(iat)%ndij==4) nsploop=4 2245 2246 do idij=1,nsploop 2247 klmn1=1 2248 2249 do jlmn=1,lmn_size 2250 j0lmn=jlmn*(jlmn-1)/2 2251 do ilmn=1,jlmn 2252 klmn=j0lmn+ilmn 2253 2254 ! The following gives back the input dij. 2255 ! since dijxc contains the hat term done on the FFT mesh 2256 if (.FALSE.) then 2257 Paw_ij(iat)%dij(klmn,idij) = & 2258 & Pawtab(itypat)%dij0 (klmn) & 2259 & +Paw_ij(iat)%dijhartree(klmn) & 2260 & +Paw_ij(iat)%dijxc (klmn,idij) & 2261 & +dijhl_hat (klmn,idij) 2262 2263 else 2264 ! === Make nonlocal part of h0 removing the valence contribution === 2265 ! Remeber that XC contains already the Hat contribution 2266 Paw_ij(iat)%dij(klmn,idij) = & 2267 & Pawtab(itypat)%dij0 (klmn) & 2268 & +Paw_ij(iat)%dijhartree(klmn) & 2269 & +Paw_ij(iat)%dijxc (klmn,idij) & ! 2 lines to get the d1-dt1 XC core contribution + XC hat (core+val) 2270 & -Paw_ij(iat)%dijxc_val (klmn,idij) & ! I suppose that the "hat" term on the FFT mesh in included in both. 2271 & +dijhmxc_val(klmn,idij) ! Local + Hartree - XC val contribution to the "hat" term. 2272 2273 ! Add the U contribution to the 2274 ! if (.FALSE. .and. Pawtab(itypat)%usepawu/=0) then 2275 if (.TRUE. .and. Pawtab(itypat)%usepawu/=0) then 2276 Paw_ij(iat)%dij(klmn,idij) = Paw_ij(iat)%dij(klmn,idij) + Paw_ij(iat)%dijU(klmn,idij) 2277 end if 2278 end if 2279 ! TODO dijso, dijU, vpawx? 2280 ! Just to be consistent, update some values. 2281 !$Paw_ij(iat)%dijhat(klmn,idij)=Paw_ij(iat)%dijhat(klmn,idij)-dijhmxc_val(klmn,idij) 2282 2283 end do !ilmn 2284 end do !jlmn 2285 end do !idij 2286 2287 ! this is to be consistent? 2288 ! deallocate(Paw_ij(iat)%dijvxc_val) 2289 ABI_FREE(dijhl_hat) 2290 ABI_FREE(dijhmxc_val) 2291 end do !iat 2292 2293 !=== Symmetrize total Dij === 2294 option_dij=0 ! For total Dij. 2295 #if 0 2296 if (paral_atom) then 2297 call symdij(Cryst%gprimd,Cryst%indsym,ipert0,my_natom,Cryst%natom,Cryst%nsym,Cryst%ntypat,option_dij,& 2298 & Paw_ij,Pawang,pawprtvol,Pawtab,Cryst%rprimd,Cryst%symafm,Cryst%symrec,& 2299 & comm_atom=my_comm_atom,mpi_atmtab=my_atmtab) 2300 else 2301 call symdij(Cryst%gprimd,,Cryst%indsym,ipert0,my_natom,Cryst%natom,Cryst%nsym,Cryst%ntypat,option_dij,& 2302 & Paw_ij,Pawang,pawprtvol,Pawtab,Cryst%rprimd,Cryst%symafm,Cryst%symrec) 2303 end if 2304 #else 2305 if (paral_atom) then 2306 call symdij_all(Cryst%gprimd,Cryst%indsym,ipert0,my_natom,Cryst%natom,Cryst%nsym,Cryst%ntypat,& 2307 & Paw_ij,Pawang,pawprtvol,Pawtab,Cryst%rprimd,Cryst%symafm,Cryst%symrec,& 2308 & comm_atom=my_comm_atom,mpi_atmtab=my_atmtab) 2309 else 2310 call symdij_all(Cryst%gprimd,Cryst%indsym,ipert0,my_natom,Cryst%natom,Cryst%nsym,Cryst%ntypat,& 2311 & Paw_ij,Pawang,pawprtvol,Pawtab,Cryst%rprimd,Cryst%symafm,Cryst%symrec) 2312 end if 2313 #endif 2314 2315 !Destroy atom table used for parallelism 2316 call free_my_atmtab(my_atmtab,my_atmtab_allocated) 2317 2318 DBG_EXIT("COLL") 2319 2320 end subroutine paw_mknewh0
m_paw_denpot/pawaccenergy [ Functions ]
[ Top ] [ m_paw_denpot ] [ Functions ]
NAME
pawaccenergy
FUNCTION
Accumulate an on-site PAW energy by adding the contribution of the current atom. This contribution has the form: Sum_ij[Rhoij.Dij]
INPUTS
pawrhoij<type(pawrhoij_type)>= datastructure containing Rho_ij values dij(cplex_dij*qphase_dij*lmn2_size,nspden_dij)= array containing D_ij values cplex_dij= 2 if dij is COMPLEX (as in the spin-orbit case), 1 if dij is REAL qphase_dij= 2 if dij has a exp(iqR) phase, 1 if not nspden_dij= number of spin components for dij pawtab<type(pawtab_type)>=paw tabulated starting data
OUTPUT
SIDE EFFECTS
epaw= PAW on-site energy. At output, the contribution of the current atom has been added to epaw. [epaw_im]= imaginary part of PAW on-site energy. At output, the contribution of the current atom has been added to epaw. This imaginary p rt only exists in a few cases (f.i. non-stationnary expression of 2nd-order energy)
NOTES
* The general form for Dij is:
D^{s1,s2}_ij = D1^{s1,s2}_ij.cos(qr) + i.D2^{s1,s2}_ij.sin(qr)
= [D1re^{s1,s2}_ij + i.D1im^{s1,s2}_ij).cos(qr)]
+ i.[D2re^{s1,s2}_ij + i.D2im^{s1,s2}_ij).sin(qr)]
where
ij are the partial waves channels
s1,s2 are spin/spinor components
q is the wave vector of the phase
D1^{s1,s2}_ij.cos(qr) is stored in the the first half of paw_ij%dij and corresponds to iq=1
D2^{s1,s2}_ij.sin(qr) is stored in the the 2nd half of paw_ij%dij and corresponds to iq=2
D1^{s1,s2}_ij.cos(qr) and D2^{s1,s2}_ij.sin(qr) are complex if cplex_dij=2
* The same for Rho_ij
* The contribution to the PAW on-site energy is:
Sum_ij_s1s2[Rho^{s2,s1}_ij * D^{s1,s2}_ij]
Note the order of s1/s2 indices, especially for Rho_ij.
The present implementation follows eq(15) in Hobbs et al, PRB 62, 11556(2000)
rho^{s1,s2}^_ij = Sum[<Psi^s2|pi><pj|Psi^s1] (s1 and s2 exponents inverted)
SOURCE
1829 subroutine pawaccenergy(epaw,pawrhoij,dij,cplex_dij,qphase_dij,nspden_dij,pawtab,epaw_im) 1830 1831 !Arguments --------------------------------------------- 1832 !scalars 1833 integer,intent(in) :: cplex_dij,qphase_dij,nspden_dij 1834 real(dp),intent(inout) :: epaw 1835 real(dp),intent(inout),optional :: epaw_im 1836 type(pawrhoij_type),intent(in),target :: pawrhoij 1837 type(pawtab_type),intent(in) :: pawtab 1838 !arrays 1839 real(dp),intent(in) :: dij(cplex_dij*qphase_dij*pawtab%lmn2_size,nspden_dij) 1840 1841 !Local variables --------------------------------------- 1842 !scalars 1843 integer :: cplex_rhoij,iq,iq0_dij,iq0_rhoij,irhoij,isp_dij,isp_rhoij,jrhoij 1844 integer :: klmn,kklmn,krhoij,lmn2_size,nspden_rhoij,nsploop 1845 logical :: add_imaginary 1846 real(dp) :: etmp 1847 character(len=500) :: msg 1848 !arrays 1849 real(dp),pointer :: rhoij(:,:) 1850 1851 ! ************************************************************************* 1852 1853 DBG_ENTER("COLL") 1854 1855 !Compatibility tests 1856 if (pawrhoij%qphase/=qphase_dij) then 1857 msg='pawaccenergy: pawrhoij%qphase/=qphase_dij!' 1858 ABI_BUG(msg) 1859 end if 1860 if (pawrhoij%nspden>nspden_dij.and.nspden_dij/=1) then 1861 msg='pawaccenergy: pawrhoij%nspden>nspden_dij!' 1862 ABI_BUG(msg) 1863 end if 1864 1865 !Useful data 1866 nspden_rhoij=pawrhoij%nspden 1867 lmn2_size=pawtab%lmn2_size 1868 1869 !Special treatment for nspden 1870 nsploop=nspden_rhoij 1871 if (nspden_dij==1.and.nspden_rhoij==4) nsploop=1 1872 1873 !Non-collinear case: need a temporary rhoij 1874 if (nspden_rhoij==4.and.nspden_dij==4) then 1875 cplex_rhoij=2 1876 ABI_MALLOC(rhoij,(2*lmn2_size,4)) 1877 else 1878 cplex_rhoij=pawrhoij%cplex_rhoij 1879 rhoij => pawrhoij%rhoijp 1880 end if 1881 1882 add_imaginary=(cplex_dij==2.and.cplex_rhoij==2) 1883 1884 !Loop over qphase components 1885 do iq=1,qphase_dij 1886 iq0_rhoij=(iq-1)*lmn2_size*cplex_rhoij 1887 iq0_dij =(iq-1)*lmn2_size*cplex_dij 1888 1889 ! Non-collinear case 1890 if (nspden_rhoij==4.and.nspden_dij==4) then 1891 rhoij(:,:)=zero 1892 jrhoij=(iq-1)*lmn2_size*pawrhoij%cplex_rhoij+1 ; krhoij=1 1893 do irhoij=1,pawrhoij%nrhoijsel 1894 klmn=pawrhoij%rhoijselect(irhoij) 1895 rhoij(krhoij ,1)= half*(pawrhoij%rhoijp(jrhoij,1)+pawrhoij%rhoijp(jrhoij,4)) 1896 rhoij(krhoij ,2)= half*(pawrhoij%rhoijp(jrhoij,1)-pawrhoij%rhoijp(jrhoij,4)) 1897 !Be careful we store rhoij^21 in rhoij(:,3) and rhoij^12 in rhoij(:,4) 1898 !because of the inversion of spins in rhoij definition 1899 rhoij(krhoij ,3)= half*pawrhoij%rhoijp(jrhoij,2) 1900 rhoij(krhoij+1,3)= half*pawrhoij%rhoijp(jrhoij,3) 1901 rhoij(krhoij ,4)= half*pawrhoij%rhoijp(jrhoij,2) 1902 rhoij(krhoij+1,4)=-half*pawrhoij%rhoijp(jrhoij,3) 1903 if (pawrhoij%cplex_rhoij==2) then 1904 rhoij(krhoij+1,1)= half*(pawrhoij%rhoijp(jrhoij+1,1)+pawrhoij%rhoijp(jrhoij+1,4)) 1905 rhoij(krhoij+1,2)= half*(pawrhoij%rhoijp(jrhoij+1,1)-pawrhoij%rhoijp(jrhoij+1,4)) 1906 !Be careful we store rhoij^21 in rhoij(:,3) and rhoij^12 in rhoij(:,4) 1907 !because of the inversion of spins in rhoij definition 1908 rhoij(krhoij ,3)= rhoij(krhoij ,3)-half*pawrhoij%rhoijp(jrhoij+1,3) 1909 rhoij(krhoij+1,3)= rhoij(krhoij+1,3)+half*pawrhoij%rhoijp(jrhoij+1,2) 1910 rhoij(krhoij ,4)= rhoij(krhoij ,4)+half*pawrhoij%rhoijp(jrhoij+1,3) 1911 rhoij(krhoij+1,4)= rhoij(krhoij+1,4)+half*pawrhoij%rhoijp(jrhoij+1,2) 1912 end if 1913 jrhoij=jrhoij+pawrhoij%cplex_rhoij ; krhoij=krhoij+2 1914 end do 1915 iq0_rhoij=0 1916 end if 1917 1918 ! Contribution to on-site energy (real part) 1919 do isp_rhoij=1,nsploop 1920 isp_dij=min(isp_rhoij,nspden_dij) 1921 jrhoij=iq0_rhoij+1 1922 do irhoij=1,pawrhoij%nrhoijsel 1923 klmn=pawrhoij%rhoijselect(irhoij) 1924 kklmn=iq0_dij+cplex_dij*(klmn-1)+1 1925 etmp=rhoij(jrhoij,isp_rhoij)*dij(kklmn,isp_dij) 1926 if (add_imaginary) etmp=etmp-rhoij(jrhoij+1,isp_rhoij)*dij(kklmn+1,isp_dij) 1927 epaw=epaw+etmp*pawtab%dltij(klmn) 1928 jrhoij=jrhoij+cplex_rhoij 1929 end do 1930 end do ! nsploop 1931 1932 ! Contribution to on-site energy (imaginary part) 1933 if (present(epaw_im).and.qphase_dij==2) then 1934 do isp_rhoij=1,nsploop 1935 isp_dij=min(isp_rhoij,nspden_dij) 1936 jrhoij=iq0_rhoij+1 1937 do irhoij=1,pawrhoij%nrhoijsel 1938 klmn=pawrhoij%rhoijselect(irhoij) 1939 if (iq==1) then 1940 kklmn=lmn2_size*cplex_dij+cplex_dij*(klmn-1)+1 1941 etmp=-rhoij(jrhoij,isp_rhoij)*dij(kklmn,isp_dij) 1942 if (add_imaginary) etmp=etmp+rhoij(jrhoij+1,isp_rhoij)*dij(kklmn+1,isp_dij) 1943 end if 1944 if (iq==2) then 1945 kklmn=cplex_dij*(klmn-1)+1 1946 etmp=rhoij(jrhoij,isp_rhoij)*dij(kklmn,isp_dij) 1947 if (add_imaginary) etmp=etmp-rhoij(jrhoij+1,isp_rhoij)*dij(kklmn+1,isp_dij) 1948 end if 1949 epaw_im=epaw_im+etmp*pawtab%dltij(klmn) 1950 jrhoij=jrhoij+cplex_rhoij 1951 end do 1952 end do ! nsploop 1953 end if 1954 1955 end do ! qphase 1956 1957 if (nspden_rhoij==4.and.nspden_dij==4) then 1958 ABI_FREE(rhoij) 1959 end if 1960 1961 DBG_EXIT("COLL") 1962 1963 end subroutine pawaccenergy
m_paw_denpot/pawaccenergy_nospin [ Functions ]
[ Top ] [ m_paw_denpot ] [ Functions ]
NAME
pawaccenergy_nospin
FUNCTION
Accumulate an on-site PAW energy by adding the contribution of the current atom. This contribution has the form: Sum_ij[Rhoij.Dij] Applies only for Dij without spin components (as f.i. Dij^Hartree). This routine is a wrapper to pawaccenergy.
INPUTS
pawrhoij<type(pawrhoij_type)>= datastructure containing Rho_ij values dij(cplex_dij*qphase_dij*lmn2_size)= array containing D_ij values cplex_dij= 2 if dij is COMPLEX (as in the spin-orbit case), 1 if dij is REAL qphase_dij= 2 if dij has a exp(iqR) phase, 1 if not pawtab<type(pawtab_type)>=paw tabulated starting data
OUTPUT
SIDE EFFECTS
epaw= PAW on-site energy. At output, the contribution of the current atom has been added to epaw. [epaw_im]= imaginary part of PAW on-site energy. At output, the contribution of the current atom has been added to epaw. This imaginary part only exists in a few cases (f.i. non-stationnary expression of 2nd-order energy)
SOURCE
1997 subroutine pawaccenergy_nospin(epaw,pawrhoij,dij,cplex_dij,qphase_dij,pawtab,epaw_im) 1998 1999 !Arguments --------------------------------------------- 2000 !scalars 2001 integer,intent(in) :: cplex_dij,qphase_dij 2002 real(dp),intent(inout) :: epaw 2003 real(dp),intent(inout),optional :: epaw_im 2004 type(pawrhoij_type),intent(in),target :: pawrhoij 2005 type(pawtab_type),intent(in) :: pawtab 2006 !arrays 2007 real(dp),intent(in),target :: dij(cplex_dij*qphase_dij*pawtab%lmn2_size) 2008 2009 !Local variables --------------------------------------- 2010 !scalars 2011 integer :: size_dij 2012 #ifdef HAVE_FC_ISO_C_BINDING 2013 type(C_PTR) :: cptr 2014 #endif 2015 !arrays 2016 real(dp), ABI_CONTIGUOUS pointer :: dij_2D(:,:) 2017 2018 ! ************************************************************************* 2019 2020 size_dij=size(dij) 2021 2022 #ifdef HAVE_FC_ISO_C_BINDING 2023 cptr=c_loc(dij(1)) 2024 call c_f_pointer(cptr,dij_2D,shape=[size_dij,1]) 2025 #else 2026 ABI_MALLOC(dij_2D,(size_dij,1)) 2027 dij_2D=reshape(dij,[size_dij,1]) 2028 #endif 2029 2030 if (present(epaw_im)) then 2031 call pawaccenergy(epaw,pawrhoij,dij_2D,cplex_dij,qphase_dij,1,pawtab,epaw_im=epaw_im) 2032 else 2033 call pawaccenergy(epaw,pawrhoij,dij_2D,cplex_dij,qphase_dij,1,pawtab) 2034 end if 2035 2036 #ifndef HAVE_FC_ISO_C_BINDING 2037 ABI_FREE(dij_2D) 2038 #endif 2039 2040 end subroutine pawaccenergy_nospin
m_paw_denpot/pawdenpot [ Functions ]
[ Top ] [ m_paw_denpot ] [ Functions ]
NAME
pawdenpot
FUNCTION
Compute different (PAW) energies, densities and potentials (or potential-like quantities) inside PAW spheres Can also compute first-order densities potentials and second-order energies (RF calculations).
INPUTS
electronpositron <type(electronpositron_type)>=quantities for the electron-positron annihilation (optional argument) [hyb_mixing, hyb_mixing_sr]= -- optional-- mixing factors for the global (resp. screened) XC hybrid functional ipert=index of perturbation (used only for RF calculation ; set ipert<=0 for GS calculations. ixc= choice of exchange-correlation scheme (see above, and below) mpi_atmtab(:)=--optional-- indexes of the atoms treated by current proc comm_atom=--optional-- MPI communicator over atoms my_natom=number of atoms treated by current processor natom=total number of atoms in cell nspden=number of spin-density components ntypat=number of types of atoms in unit cell. nucdipmom(3,natom) nuclear dipole moments nzlmopt= if -1, compute all LM-moments of densities initialize "lmselect" (index of non-zero LM-moments of densities) if 0, compute all LM-moments of densities force "lmselect" to .true. (index of non-zero LM-moments of densities) if 1, compute only non-zero LM-moments of densities (stored before) option=0: compute both energies and potentials 1: compute only potentials 2: compute only energies paw_an(my_natom) <type(paw_an_type)>=paw arrays given on angular mesh paw_an0(my_natom) <type(paw_an_type)>=paw arrays given on angular mesh for Ground-State paw_ij(my_natom) <type(paw_ij_type)>=paw arrays given on (i,j) channels pawang <type(pawang_type)>=paw angular mesh and related data pawprtvol=control print volume and debugging output for PAW pawrad(ntypat) <type(pawrad_type)>=paw radial mesh and related data pawrhoij(my_natom) <type(pawrhoij_type)>= paw rhoij occupancies and related data pawspnorb=flag: 1 if spin-orbit coupling is activated pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data pawxcdev=Choice of XC development (0=no dev. (use of angular mesh) ; 1 or 2=dev. on moments) ucvol=unit cell volume (bohr^3) xclevel= XC functional level xc_denpos= lowest allowed density (usually for the computation of the XC functionals) xc_taupos= lowest allowed kinetic energy density (for mGGA XC functionals) znucl(ntypat)=gives the nuclear charge for all types of atoms
OUTPUT
paw_ij(my_natom)%dijhartree(qphase*lmn2_size)=Hartree contribution to dij; Enters into calculation of hartree energy ==== if option=0 or 2 epaw=contribution to total energy from the PAW "on-site" part epawdc=contribution to total double-counting energy from the PAW "on-site" part ==== if option=0 or 2 and ipert<=0 compch_sph=compensation charge integral inside spheres computed over spherical meshes ==== if (option=0 or 1) and paw_an(:)%has_vxc=1 paw_an(my_natom)%vxc1(cplex*mesh_size,:,nspden)=XC potential calculated from "on-site" density paw_an(my_natom)%vxct1(cplex*mesh_size,:,nspden)=XC potential calculated from "on-site" pseudo density ==== if (option=0 or 1) and paw_an(:)%has_vxctau=1 paw_an(my_natom)%vxctau1(cplex*mesh_size,:,nspden)=1st deriv. of XC energy wrt to kinetic energy density (all electron) paw_an(my_natom)%vxcttau1(cplex*mesh_size,:,nspden)=1st deriv. of XC energy wrt to kinetic energy density (pseudo) ==== if paw_an(iatom_tot)%has_vxcval==1 compute also XC potentials neglecting core charge paw_an(my_natom)%vxc1_val(cplex*mesh_size,:nspden)=XC potential calculated from spherical valence density paw_an(my_natom)%vxct1_val(cplex*mesh_size,:nspden)=XC potential calculated from spherical valence pseudo density ==== if nzlmopt==-1, paw_an(iatom_tot)%lnmselect(lm_size,nspden)=select the non-zero LM-moments of rho1 and trho1 ==== if paw_an(:)%has_vhartree=1 paw_an(my_natom)%vh1(cplex*mesh_size,1,1)=Hartree total potential calculated from "on-site" density ==== if pawspnorb>0 paw_ij(my_natom)%dijso(qphase*cplex_dij*lmn2_size,nspden)=spin-orbit contribution to dij
NOTES
Response function calculations:
In order to compute first- or second-order quantities, paw_an (resp. paw_ij) datastructures
must contain first-order quantities, namely paw_an1 (resp. paw_ij1).
SOURCE
146 subroutine pawdenpot(compch_sph,epaw,epawdc,ipert,ixc,& 147 & my_natom,natom,nspden,ntypat,nucdipmom,nzlmopt,option,paw_an,paw_an0,& 148 & paw_ij,pawang,pawprtvol,pawrad,pawrhoij,pawspnorb,pawtab,pawxcdev,spnorbscl,& 149 & xclevel,xc_denpos,xc_taupos,ucvol,znucl,& 150 & electronpositron,mpi_atmtab,comm_atom,vpotzero,hyb_mixing,hyb_mixing_sr) ! optional arguments 151 152 !Arguments --------------------------------------------- 153 !scalars 154 integer,intent(in) :: ipert,ixc,my_natom,natom,nspden,ntypat,nzlmopt,option,pawprtvol 155 integer,intent(in) :: pawspnorb,pawxcdev,xclevel 156 integer,optional,intent(in) :: comm_atom 157 real(dp), intent(in) :: spnorbscl,xc_denpos,xc_taupos,ucvol 158 real(dp),intent(in),optional :: hyb_mixing,hyb_mixing_sr 159 real(dp),intent(out) :: compch_sph,epaw,epawdc 160 type(electronpositron_type),pointer,optional :: electronpositron 161 type(pawang_type),intent(in) :: pawang 162 !arrays 163 integer,optional,target,intent(in) :: mpi_atmtab(:) 164 real(dp),intent(in) :: nucdipmom(3,natom),znucl(ntypat) 165 real(dp),intent(out),optional :: vpotzero(2) 166 type(paw_an_type),intent(inout) :: paw_an(my_natom) 167 type(paw_an_type), intent(in) :: paw_an0(my_natom) 168 type(paw_ij_type),intent(inout) :: paw_ij(my_natom) 169 type(pawrad_type),intent(in) :: pawrad(ntypat) 170 type(pawrhoij_type),intent(inout) :: pawrhoij(my_natom) 171 type(pawtab_type),intent(in) :: pawtab(ntypat) 172 173 !Local variables --------------------------------------- 174 !scalars 175 integer, parameter :: PAWU_ALGO_1=1,PAWU_ALGO_2=2 176 integer, parameter :: PAWU_FLL=1,PAWU_AMF=2 177 integer :: cplex,cplex_dij,cplex_rhoij,has_kxc,has_k3xc,has_vxctau 178 integer :: iatom,iatom_tot,idum,ierr,ii,ipositron,iq,iq0_dij,iq0_rhoij 179 integer :: itypat,itypat0,lm_size,lmn2_size,mesh_size 180 integer :: my_comm_atom,ndij,nkxc1,nk3xc1,nsppol,opt_compch,pawu_algo,pawu_dblec 181 integer :: qphase,usecore,usekden,usetcore,usepawu,usexcnhat,usenhat,usefock 182 logical :: keep_vhartree,my_atmtab_allocated,need_kxc,need_k3xc,need_vxctau 183 logical :: xc_is_tb09,non_magnetic_xc,paral_atom,temp_vxc 184 real(dp) :: e1t10,e1xc,e1xcdc,efock,efockdc,eexc,eexcdc,eexdctemp 185 real(dp) :: eexc_val,eexcdc_val,eexex,eexexdc,eextemp,eh2 186 real(dp) :: edftumdc,edftumdcdc,edftufll,enucdip,etmp,espnorb,etild1xc,etild1xcdc 187 real(dp) :: exccore,exchmix,hyb_mixing_,hyb_mixing_sr_,rdum 188 character(len=3) :: pertstrg 189 character(len=500) :: msg 190 !arrays 191 integer :: idum1(0),idum3(0,0,0) 192 integer,pointer :: my_atmtab(:) 193 logical,allocatable :: lmselect_cur(:),lmselect_cur_ep(:),lmselect_ep(:),lmselect_tmp(:) 194 real(dp) :: mpiarr(7),tsec(2) 195 real(dp),allocatable :: dij_ep(:),dijfock_vv(:,:),dijfock_cv(:,:) 196 real(dp),allocatable :: one_over_rad2(:),kxc_tmp(:,:,:),k3xc_tmp(:,:,:) 197 real(dp),allocatable :: nhat1(:,:,:),nhat1_ep(:,:,:) 198 real(dp) :: rdum2(0,0),rdum3(0,0,0),rdum3a(0,0,0),rdum4(0,0,0,0) 199 real(dp),allocatable :: rho(:),rho1(:,:,:),rho1_ep(:,:,:),rho1xx(:,:,:) 200 real(dp),allocatable :: tau1(:,:,:),ttau1(:,:,:), trho1(:,:,:),trho1_ep(:,:,:) 201 real(dp),allocatable :: vh(:),vxc_tmp(:,:,:),vxctau_tmp(:,:,:) 202 203 ! ************************************************************************* 204 205 DBG_ENTER("COLL") 206 207 call timab(560,1,tsec) 208 209 !Various inits 210 hyb_mixing_ =zero ; if(present(hyb_mixing)) hyb_mixing_ =hyb_mixing 211 hyb_mixing_sr_=zero ; if(present(hyb_mixing_sr)) hyb_mixing_sr_=hyb_mixing_sr 212 usefock=0;if (abs(hyb_mixing_)>tol8.or.abs(hyb_mixing_sr_)>tol8) usefock=1 213 usexcnhat=maxval(pawtab(1:ntypat)%usexcnhat) 214 usekden=pawxc_get_usekden(ixc) 215 xc_is_tb09=pawxc_is_tb09(ixc) 216 usenhat = usexcnhat 217 keep_vhartree=(maxval(paw_an(:)%has_vhartree)>0) 218 if (keep_vhartree) usenhat = 1 219 compch_sph=-1.d5 220 opt_compch=0;if (option/=1.and.ipert<=0) opt_compch=1 221 if (opt_compch==1) compch_sph=zero 222 nsppol=1;if (my_natom>0) nsppol=pawrhoij(1)%nsppol 223 pertstrg=" ";if (ipert>0) pertstrg="(1)" 224 225 !Various checks 226 if(nzlmopt/=0.and.nzlmopt/=1.and.nzlmopt/=-1) then 227 msg='invalid value for variable "nzlmopt"!' 228 ABI_BUG(msg) 229 end if 230 if (my_natom>0) then 231 if(paw_ij(1)%has_dijhartree==0.and.ipert/=natom+1.and.ipert/=natom+10) then 232 msg='dijhartree must be allocated!' 233 ABI_BUG(msg) 234 end if 235 if(paw_ij(1)%has_dijU==0.and.pawtab(1)%usepawu/=0.and. & 236 & ((ipert>0.and.ipert/=natom+1.and.ipert/=natom+10).or.pawtab(1)%usepawu<0)) then 237 msg='dijU must be allocated!' 238 ABI_BUG(msg) 239 end if 240 if (pawrhoij(1)%qphase<paw_an(1)%cplex) then 241 msg='pawrhoij()%qphase must be >=paw_an()%cplex!' 242 ABI_BUG(msg) 243 end if 244 if (ipert>0.and.(ipert<=natom.or.ipert==natom+2).and.paw_an0(1)%has_kxc/=2) then 245 msg='XC kernels for ground state must be in memory!' 246 ABI_BUG(msg) 247 end if 248 if(paw_an(1)%has_vxc==0.and.(option==0.or.option==1).and. & 249 & .not.(ipert==natom+1.or.ipert==natom+10)) then 250 msg='vxc1 and vxct1 must be allocated!' 251 ABI_BUG(msg) 252 end if 253 if(paw_an(1)%has_vxctau==0.and.(option==0.or.option==1).and.usekden==1) then 254 msg='vxctau1 and vxcttau1 must be allocated!' 255 ABI_BUG(msg) 256 end if 257 if (ipert>0.and.paw_an(1)%has_vxctau==1.and.usekden==1) then 258 if (ipert .NE. natom+1) then 259 msg='computation of vxctau not compatible with RF (ipert>0)!' 260 ABI_BUG(msg) 261 end if 262 end if 263 if (ipert>0.and.paw_an(1)%has_vhartree==1) then 264 msg='computation of vhartree not compatible with RF (ipert>0)!' 265 ABI_BUG(msg) 266 end if 267 if (ipert>0.and.paw_an(1)%has_vxcval==1.and.(option==0.or.option==1)) then 268 msg='computation of vxc_val not compatible with RF (ipert>0)!' 269 ABI_BUG(msg) 270 end if 271 end if 272 273 ipositron=0 274 if (present(electronpositron)) then 275 ipositron=electronpositron_calctype(electronpositron) 276 if (ipositron==1.and.pawtab(1)%has_kij/=2) then 277 msg='kij must be in memory for electronpositron%calctype=1!' 278 ABI_BUG(msg) 279 end if 280 if (ipert>0) then 281 msg='electron-positron calculation not available for ipert>0!' 282 ABI_ERROR(msg) 283 end if 284 end if 285 286 !Set up parallelism over atoms 287 paral_atom=(present(comm_atom).and.(my_natom/=natom)) 288 nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab 289 my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom 290 call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,& 291 & my_natom_ref=my_natom) 292 293 !For some perturbations, nothing to do 294 if (ipert==natom+1.or.ipert==natom+10) then 295 if (option/=1) then 296 epaw=zero;epawdc=zero 297 end if 298 return 299 end if 300 301 !Init energies 302 if (option/=1) then 303 e1xc=zero ; e1xcdc=zero 304 etild1xc=zero ; etild1xcdc=zero 305 exccore=zero ; eh2=zero ; e1t10=zero 306 edftumdc=zero ; edftumdcdc=zero ; edftufll=zero 307 eexex=zero ; eexexdc=zero 308 eextemp=zero ; eexdctemp=zero 309 espnorb=zero ; enucdip=zero 310 efock=zero ; efockdc=zero 311 if (ipositron/=0) then 312 electronpositron%e_paw =zero 313 electronpositron%e_pawdc=zero 314 end if 315 end if 316 !vpotzero is needed for both the energy and the potential 317 if (present(vpotzero)) vpotzero(:)=zero 318 319 !Select PAW+U algo, different for DFT and DFPT 320 usepawu=maxval(pawtab(1:ntypat)%usepawu) 321 ii=minval(pawtab(1:ntypat)%usepawu);if (ii<0) usepawu=ii 322 non_magnetic_xc=(mod(abs(usepawu),10)==4) 323 324 !if PAW+U, compute noccmmp^{\sigma}_{m,m'} occupation matrix 325 if (usepawu/=0.and.ipert<=0.and.ipositron/=1) then 326 if (paral_atom) then 327 call setnoccmmp(1,0,rdum4,0,0,idum3,my_natom,natom,0,1,nsppol,0,ntypat,& 328 & paw_ij,pawang,pawprtvol,pawrhoij,pawtab,rdum2,idum1,idum1,0,usepawu,& 329 & comm_atom=my_comm_atom,mpi_atmtab=mpi_atmtab) 330 else 331 call setnoccmmp(1,0,rdum4,0,0,idum3,my_natom,natom,0,1,nsppol,0,ntypat,& 332 & paw_ij,pawang,pawprtvol,pawrhoij,pawtab,rdum2,idum1,idum1,0,usepawu) 333 end if 334 end if 335 336 !Print some titles 337 if (abs(pawprtvol)>=2) then 338 if (nzlmopt<1) write(msg, '(6a)') ch10,' PAW TEST:',ch10,& 339 ' ====== Moments of (n1-tn1)',trim(pertstrg),' =========' 340 if (nzlmopt==1) write(msg, '(6a)') ch10,' PAW TEST:',ch10,& 341 ' ==== Non-zero Moments of (n1-tn1)',trim(pertstrg),' ====' 342 call wrtout(std_out,msg,'COLL') 343 if (usexcnhat/=0) then 344 write(msg, '(6a)')' The moments of (n1-tn1-nhat1)',trim(pertstrg),' must be very small...' 345 call wrtout(std_out,msg,'COLL') 346 end if 347 end if 348 349 !================ Big loop on atoms ======================= 350 !========================================================== 351 352 do iatom=1,my_natom 353 iatom_tot=iatom;if (paral_atom) iatom_tot=my_atmtab(iatom) 354 itypat=pawrhoij(iatom)%itypat 355 exchmix=pawtab(itypat)%exchmix 356 lmn2_size=paw_ij(iatom)%lmn2_size 357 lm_size=paw_an(iatom)%lm_size 358 mesh_size=pawtab(itypat)%mesh_size 359 usecore=1;usetcore =pawtab(itypat)%usetcore 360 if (ipert/=0) usecore=0 ! This is true for phonons and Efield pert. 361 if (ipert/=0) usetcore=0 ! This is true for phonons and Efield pert. 362 has_kxc =paw_an(iatom)%has_kxc ;need_kxc =(has_kxc ==1) 363 has_k3xc=paw_an(iatom)%has_k3xc;need_k3xc=(has_k3xc==1) 364 has_vxctau=paw_an(iatom)%has_vxctau ;need_vxctau =(has_vxctau>=1.and.usekden==1) 365 cplex=paw_an(iatom)%cplex 366 cplex_dij=paw_ij(iatom)%cplex_dij 367 cplex_rhoij=pawrhoij(iatom)%cplex_rhoij 368 qphase=pawrhoij(iatom)%qphase 369 ndij=paw_ij(iatom)%ndij 370 iq0_rhoij=cplex_rhoij*lmn2_size 371 iq0_dij=cplex_dij*lmn2_size 372 usepawu=pawtab(itypat)%usepawu 373 pawu_algo=merge(PAWU_ALGO_1,PAWU_ALGO_2,ipert<=0.and.usepawu>=0) 374 pawu_dblec=merge(PAWU_FLL,PAWU_AMF,abs(usepawu)==1.or.abs(usepawu)==4) 375 376 ! Allocations of "on-site" densities 377 ABI_MALLOC(rho1 ,(cplex*mesh_size,lm_size,nspden)) 378 ABI_MALLOC(trho1,(cplex*mesh_size,lm_size,nspden)) 379 ABI_MALLOC(nhat1,(cplex*mesh_size,lm_size,nspden*usenhat)) 380 rho1(:,:,:)=zero;trho1(:,:,:)=zero;nhat1(:,:,:)=zero 381 if (usekden==1) then 382 ABI_MALLOC(tau1 ,(cplex*mesh_size,lm_size,nspden)) 383 ABI_MALLOC(ttau1,(cplex*mesh_size,lm_size,nspden)) 384 end if 385 if (ipositron/=0) then ! Additional allocation for the electron-positron case 386 ABI_MALLOC(rho1_ep ,(cplex*mesh_size,lm_size,nspden)) 387 ABI_MALLOC(trho1_ep,(cplex*mesh_size,lm_size,nspden)) 388 ABI_MALLOC(nhat1_ep,(cplex*mesh_size,lm_size,nspden*usenhat)) 389 end if 390 ABI_MALLOC(lmselect_cur,(lm_size)) 391 lmselect_cur(:)=.true. 392 if (nzlmopt==1) lmselect_cur(:)=paw_an(iatom)%lmselect(:) 393 394 ! Store some usefull quantities 395 itypat0=0;if (iatom>1) itypat0=pawrhoij(iatom-1)%itypat 396 if (itypat/=itypat0) then 397 ABI_MALLOC(one_over_rad2,(mesh_size)) 398 one_over_rad2(1)=zero 399 one_over_rad2(2:mesh_size)=one/pawrad(itypat)%rad(2:mesh_size)**2 400 end if 401 402 ! Need to allocate vxc1 in particular cases 403 if (pawspnorb>0.and.ipert==0.and.option==2.and.ipositron/=1.and. & 404 & cplex_rhoij==2.and.paw_an(iatom)%has_vxc==0) then 405 ! These should already be allocated in paw_an_init! 406 if (allocated(paw_an(iatom)%vxc1)) then 407 ABI_FREE(paw_an(iatom)%vxc1) 408 end if 409 if (pawxcdev==0)then 410 ABI_MALLOC(paw_an(iatom)%vxc1,(cplex*mesh_size,paw_an(iatom)%angl_size,nspden)) 411 else 412 ABI_MALLOC(paw_an(iatom)%vxc1,(cplex*mesh_size,lm_size,nspden)) 413 end if 414 paw_an(iatom)%has_vxc=1 415 temp_vxc=.true. 416 else 417 temp_vxc=.false. 418 end if 419 420 ! ===== Compute "on-site" densities (n1, ntild1, nhat1) ===== 421 ! ========================================================== 422 423 call pawdensities(compch_sph,cplex,iatom_tot,lmselect_cur,paw_an(iatom)%lmselect,lm_size,& 424 & nhat1,nspden,nzlmopt,opt_compch,1-usenhat,-1,1,pawang,pawprtvol,pawrad(itypat),& 425 & pawrhoij(iatom),pawtab(itypat),rho1,trho1,one_over_rad2=one_over_rad2) 426 427 if (usekden==1) then 428 ABI_MALLOC(lmselect_tmp,(lm_size)) 429 lmselect_tmp(:)=.true. 430 call pawkindensities(cplex,lmselect_tmp,lm_size,nspden,-1,1,-1,& 431 & pawang,pawrad(itypat),pawrhoij(iatom),pawtab(itypat),tau1,ttau1,& 432 & one_over_rad2=one_over_rad2) 433 ABI_FREE(lmselect_tmp) 434 end if 435 436 if (ipositron/=0) then 437 ! Electron-positron calculation: need additional on-site densities: 438 ! if ipositron==1, need electronic on-site densities 439 ! if ipositron==2, need positronic on-site densities 440 ABI_MALLOC(lmselect_ep,(lm_size)) 441 ABI_MALLOC(lmselect_cur_ep,(lm_size)) 442 lmselect_cur_ep(:)=.true. 443 if (nzlmopt==1) lmselect_cur_ep(:)=electronpositron%lmselect_ep(1:lm_size,iatom) 444 call pawdensities(rdum,cplex,iatom_tot,lmselect_cur_ep,lmselect_ep,& 445 & lm_size,nhat1_ep,nspden,nzlmopt,0,1-usenhat,-1,0,pawang,0,pawrad(itypat),& 446 & electronpositron%pawrhoij_ep(iatom),pawtab(itypat),& 447 & rho1_ep,trho1_ep,one_over_rad2=one_over_rad2) 448 if (nzlmopt<1) electronpositron%lmselect_ep(1:lm_size,iatom)=lmselect_ep(1:lm_size) 449 ABI_FREE(lmselect_ep) 450 ABI_FREE(lmselect_cur_ep) 451 end if 452 453 ! =========== Compute XC potentials and energies =========== 454 ! ========================================================== 455 456 ! Temporary storage 457 nkxc1 =0;if (paw_an(iatom)%has_kxc /=0) nkxc1 =paw_an(iatom)%nkxc1 458 nk3xc1=0;if (paw_an(iatom)%has_k3xc/=0.and.pawxcdev==0) nk3xc1=paw_an(iatom)%nk3xc1 459 if (pawxcdev/=0) then 460 ABI_MALLOC(vxc_tmp,(cplex*mesh_size,lm_size,nspden)) 461 if (need_kxc) then 462 ABI_MALLOC(kxc_tmp,(mesh_size,lm_size,nkxc1)) 463 end if 464 if (need_k3xc) then 465 msg = 'Computation of k3xc with pawxcdev/=0 is not implemented yet!' 466 ABI_BUG(msg) 467 end if 468 end if 469 if (pawxcdev==0) then 470 ABI_MALLOC(vxc_tmp,(cplex*mesh_size,pawang%angl_size,nspden)) 471 vxc_tmp(:,:,:)=zero 472 if (need_kxc) then 473 ABI_MALLOC(kxc_tmp,(mesh_size,pawang%angl_size,nkxc1)) 474 end if 475 if (need_k3xc) then 476 ABI_MALLOC(k3xc_tmp,(mesh_size,pawang%angl_size,nk3xc1)) 477 end if 478 if (need_vxctau) then 479 ABI_MALLOC(vxctau_tmp,(cplex*mesh_size,pawang%angl_size,nspden)) 480 vxctau_tmp(:,:,:)=zero 481 end if 482 end if 483 idum=0 484 if (.not.allocated(vxc_tmp)) then 485 ABI_MALLOC(vxc_tmp,(0,0,0)) 486 end if 487 if (.not.allocated(kxc_tmp)) then 488 ABI_MALLOC(kxc_tmp,(0,0,0)) 489 end if 490 if (.not.allocated(k3xc_tmp)) then 491 ABI_MALLOC(k3xc_tmp,(0,0,0)) 492 end if 493 if (.not.allocated(vxctau_tmp)) then 494 ABI_MALLOC(vxctau_tmp,(0,0,0)) 495 end if 496 497 ! ===== Vxc1 term ===== 498 if (ipositron/=1) then 499 if (pawxcdev/=0) then 500 if (ipert==0) then 501 call pawxcm(pawtab(itypat)%coredens,eexc,eexcdc,idum,hyb_mixing_,ixc,kxc_tmp,lm_size,& 502 & paw_an(iatom)%lmselect,nhat1,nkxc1,non_magnetic_xc,mesh_size,nspden,option,& 503 & pawang,pawrad(itypat),pawxcdev,rho1,usecore,0,vxc_tmp,xclevel,xc_denpos) 504 else 505 call pawxcm_dfpt(pawtab(itypat)%coredens,cplex,cplex,eexc,ixc,paw_an0(iatom)%kxc1,lm_size,& 506 & paw_an(iatom)%lmselect,nhat1,paw_an0(iatom)%nkxc1,non_magnetic_xc,mesh_size,nspden,option,& 507 & pawang,pawrad(itypat),rho1,usecore,0,vxc_tmp,xclevel) 508 eexcdc=zero 509 end if 510 else 511 if (ipert==0) then 512 call pawxc(pawtab(itypat)%coredens,eexc,eexcdc,hyb_mixing_,ixc,kxc_tmp,k3xc_tmp,lm_size,& 513 & paw_an(iatom)%lmselect,nhat1,nkxc1,nk3xc1,non_magnetic_xc,mesh_size,nspden,option,& 514 & pawang,pawrad(itypat),rho1,usecore,0,vxc_tmp,xclevel,xc_denpos,& 515 & coretau=pawtab(itypat)%coretau,taur=tau1,vxctau=vxctau_tmp,xc_taupos=xc_taupos) 516 else 517 call pawxc_dfpt(pawtab(itypat)%coredens,cplex,cplex,eexc,ixc,paw_an0(iatom)%kxc1,lm_size,& 518 & paw_an(iatom)%lmselect,nhat1,paw_an0(iatom)%nkxc1,non_magnetic_xc,mesh_size,nspden,option,& 519 & pawang,pawrad(itypat),rho1,usecore,0,paw_an0(iatom)%vxc1,vxc_tmp,xclevel) 520 eexcdc=zero 521 end if 522 end if 523 if (option/=1) then 524 e1xc=e1xc+eexc 525 e1xcdc=e1xcdc+eexcdc 526 end if 527 if (option<2.or.temp_vxc) paw_an(iatom)%vxc1(:,:,:)=vxc_tmp(:,:,:) 528 if (option<2.and.need_vxctau) paw_an(iatom)%vxctau1(:,:,:)=vxctau_tmp(:,:,:) 529 if (need_kxc .and.nkxc1>0 ) paw_an(iatom)%kxc1(:,:,:) =kxc_tmp(:,:,:) 530 if (need_k3xc.and.nk3xc1>0) paw_an(iatom)%k3xc1(:,:,:)=k3xc_tmp(:,:,:) 531 else ! ipositron==1 532 if (option<2.or.temp_vxc) paw_an(iatom)%vxc1(:,:,:)=zero 533 if (need_kxc.and.nkxc1>0) paw_an(iatom)%kxc1(:,:,:)=zero 534 end if 535 536 ! Additional electron-positron XC term (if ipositron/=0) 537 if (ipositron/=0) then 538 if (pawxcdev/=0) then 539 call pawxcmpositron(ipositron,pawtab(itypat)%coredens,eexc,eexcdc,electronpositron%ixcpositron,& 540 & lm_size,paw_an(iatom)%lmselect,electronpositron%lmselect_ep(1:lm_size,iatom),& 541 & nhat1,nhat1_ep,mesh_size,nspden,option,pawang,pawrad(itypat),pawxcdev,& 542 & electronpositron%posdensity0_limit,rho1,rho1_ep,usecore,0,vxc_tmp,xc_denpos) 543 else 544 call pawxcpositron(ipositron,pawtab(itypat)%coredens,eexc,eexcdc,electronpositron%ixcpositron,& 545 & lm_size,paw_an(iatom)%lmselect,electronpositron%lmselect_ep(1:lm_size,iatom),& 546 & nhat1,nhat1_ep,mesh_size,nspden,option,pawang,pawrad(itypat),& 547 & electronpositron%posdensity0_limit,rho1,rho1_ep,usecore,0,vxc_tmp,xc_denpos) 548 end if 549 if (option/=1) then 550 electronpositron%e_paw =electronpositron%e_paw +eexc 551 electronpositron%e_pawdc=electronpositron%e_pawdc+eexcdc 552 end if 553 if (option<2.or.temp_vxc) paw_an(iatom)%vxc1(:,:,:)=paw_an(iatom)%vxc1(:,:,:)+vxc_tmp(:,:,:) 554 if (need_kxc.and.nkxc1>0) paw_an(iatom)%kxc1(:,:,:)=paw_an(iatom)%kxc1(:,:,:)+kxc_tmp(:,:,:) 555 end if 556 557 ! ===== tVxc1 term ===== 558 if (ipositron/=1) then 559 if (pawxcdev/=0) then 560 if (ipert==0) then 561 call pawxcm(pawtab(itypat)%tcoredens(:,1),& 562 & eexc,eexcdc,idum,hyb_mixing_,ixc,kxc_tmp,lm_size,& 563 & paw_an(iatom)%lmselect,nhat1,nkxc1,non_magnetic_xc,mesh_size,nspden,option,& 564 & pawang,pawrad(itypat),pawxcdev,trho1,usetcore,2*usexcnhat,vxc_tmp,xclevel,xc_denpos) 565 else 566 call pawxcm_dfpt(pawtab(itypat)%tcoredens(:,1),& 567 & cplex,cplex,eexc,ixc,paw_an0(iatom)%kxct1,lm_size,& 568 & paw_an(iatom)%lmselect,nhat1,paw_an0(iatom)%nkxc1,non_magnetic_xc,mesh_size,nspden,option,& 569 & pawang,pawrad(itypat),trho1,usetcore,2*usexcnhat,vxc_tmp,xclevel) 570 eexcdc=zero 571 end if 572 else 573 if (ipert==0) then 574 call pawxc(pawtab(itypat)%tcoredens(:,1),& 575 & eexc,eexcdc,hyb_mixing_,ixc,kxc_tmp,k3xc_tmp,lm_size,& 576 & paw_an(iatom)%lmselect,nhat1,nkxc1,nk3xc1,non_magnetic_xc,mesh_size,nspden,option,& 577 & pawang,pawrad(itypat),trho1,usetcore,2*usexcnhat,vxc_tmp,xclevel,xc_denpos,& 578 & coretau=pawtab(itypat)%tcoretau,taur=ttau1,vxctau=vxctau_tmp,xc_taupos=xc_taupos) 579 else 580 call pawxc_dfpt(pawtab(itypat)%tcoredens(:,1),& 581 & cplex,cplex,eexc,ixc,paw_an0(iatom)%kxct1,lm_size,& 582 & paw_an(iatom)%lmselect,nhat1,paw_an0(iatom)%nkxc1,non_magnetic_xc,mesh_size,nspden,option,& 583 & pawang,pawrad(itypat),trho1,usetcore,2*usexcnhat,paw_an0(iatom)%vxct1,vxc_tmp,xclevel) 584 eexcdc=zero 585 end if 586 end if 587 if (option/=1) then 588 etild1xc=etild1xc+eexc 589 etild1xcdc=etild1xcdc+eexcdc 590 end if 591 if (option<2) paw_an(iatom)%vxct1(:,:,:)=vxc_tmp(:,:,:) 592 if (option<2.and.need_vxctau) paw_an(iatom)%vxcttau1(:,:,:)=vxctau_tmp(:,:,:) 593 if (need_kxc.and. nkxc1>0 ) paw_an(iatom)%kxct1(:,:,:) =kxc_tmp(:,:,:) 594 if (need_k3xc.and.nk3xc1>0) paw_an(iatom)%k3xct1(:,:,:)=k3xc_tmp(:,:,:) 595 else ! ipositron==1 596 if (option<2) paw_an(iatom)%vxct1(:,:,:)=zero 597 if (need_kxc.and.nkxc1>0) paw_an(iatom)%kxct1(:,:,:)=zero 598 end if 599 600 ! Additional electron-positron XC term (if ipositron/=0) 601 if (ipositron/=0) then 602 if (pawxcdev/=0) then 603 call pawxcmpositron(ipositron,pawtab(itypat)%tcoredens(:,1),& 604 & eexc,eexcdc,electronpositron%ixcpositron,& 605 & lm_size,paw_an(iatom)%lmselect,electronpositron%lmselect_ep(1:lm_size,iatom),& 606 & nhat1,nhat1_ep,mesh_size,nspden,option,pawang,pawrad(itypat),pawxcdev,& 607 & electronpositron%posdensity0_limit,trho1,trho1_ep,usetcore,2*usexcnhat,vxc_tmp,xc_denpos) 608 else 609 call pawxcpositron(ipositron,pawtab(itypat)%tcoredens(:,1),& 610 & eexc,eexcdc,electronpositron%ixcpositron,& 611 & lm_size,paw_an(iatom)%lmselect,electronpositron%lmselect_ep(1:lm_size,iatom),& 612 & nhat1,nhat1_ep,mesh_size,nspden,option,pawang,pawrad(itypat),& 613 & electronpositron%posdensity0_limit,trho1,trho1_ep,usetcore,2*usexcnhat,vxc_tmp,xc_denpos) 614 end if 615 if (option/=1) then 616 electronpositron%e_paw =electronpositron%e_paw -eexc 617 electronpositron%e_pawdc=electronpositron%e_pawdc-eexcdc 618 end if 619 if (option<2) paw_an(iatom)%vxct1(:,:,:)=paw_an(iatom)%vxct1(:,:,:)+vxc_tmp(:,:,:) 620 if (need_kxc.and.nkxc1>0) paw_an(iatom)%kxct1(:,:,:)=paw_an(iatom)%kxct1(:,:,:)+kxc_tmp(:,:,:) 621 end if 622 623 ! Update flags defining the state of vxc and kxc 624 if (option<2) paw_an(iatom)%has_vxc=2 625 if (option<2.and.need_vxctau) paw_an(iatom)%has_vxctau=2 626 if (need_kxc.and.nkxc1>0) paw_an(iatom)%has_kxc=2 627 628 ! Update core XC contribution to energy 629 if (option/=1.and.ipositron/=1) exccore=exccore+pawtab(itypat)%exccore 630 631 ! =========== Compute valence-only XC potentials =========== 632 ! ========================================================== 633 if (ipert==0.and.paw_an(iatom)%has_vxcval==1.and.(option==0.or.option==1)) then 634 if (.not.allocated(paw_an(iatom)%vxc1_val).or..not.allocated(paw_an(iatom)%vxct1_val)) then 635 msg=' vxc1_val and vxct1_val must be associated' 636 ABI_BUG(msg) 637 end if 638 ! ===== Vxc1_val term, vxc[n1] ===== 639 if (pawxcdev/=0) then 640 write(msg,'(4a,es16.6)')ch10,& 641 & ' pawdenpot : Computing valence-only v_xc[n1] using moments ',ch10,& 642 & ' Min density rho1 = ',MINVAL(rho1) 643 call wrtout(std_out,msg,'COLL') 644 call pawxcm(pawtab(itypat)%coredens,eexc_val,eexcdc_val,idum,hyb_mixing_,ixc,kxc_tmp,lm_size,& 645 & paw_an(iatom)%lmselect,nhat1,nkxc1,non_magnetic_xc,mesh_size,nspden,option,& 646 & pawang,pawrad(itypat),pawxcdev,rho1,0,0,vxc_tmp,xclevel,xc_denpos) 647 else 648 write(msg,'(2a)')ch10,' pawdenpot : Computing valence-only v_xc[n1] using angular mesh ' 649 call wrtout(std_out,msg,'COLL') 650 651 call pawxc(pawtab(itypat)%coredens,eexc_val,eexcdc_val,hyb_mixing_,ixc,kxc_tmp,k3xc_tmp,lm_size,& 652 & paw_an(iatom)%lmselect,nhat1,nkxc1,nk3xc1,non_magnetic_xc,mesh_size,nspden,option,& 653 & pawang,pawrad(itypat),rho1,0,0,vxc_tmp,xclevel,xc_denpos) 654 end if 655 if (option<2) paw_an(iatom)%vxc1_val(:,:,:)=vxc_tmp(:,:,:) 656 657 ! ===== tVxc1_val term ===== 658 if (pawxcdev/=0) then 659 if (usexcnhat/=0) then 660 write(msg,'(4a,e16.6,2a,es16.6)')ch10,& 661 & ' pawdenpot : Computing valence-only v_xc[tn1+nhat] using moments ',ch10,& 662 & ' Min density trho1 = ',MINVAL(trho1),ch10,& 663 & ' Min density trho1 + nhat = ',MINVAL(trho1+nhat1) 664 else 665 write(msg,'(4a,e16.6)')ch10,& 666 & ' pawdenpot : Computing valence-only v_xc[tn1] using moments ',ch10,& 667 & ' Min density trho1 = ',MINVAL(trho1) 668 end if 669 call wrtout(std_out,msg,'COLL') 670 call pawxcm(pawtab(itypat)%tcoredens(:,1),& 671 & eexc_val,eexcdc_val,idum,hyb_mixing_,ixc,kxc_tmp,lm_size,& 672 & paw_an(iatom)%lmselect,nhat1,nkxc1,non_magnetic_xc,mesh_size,nspden,option,& 673 & pawang,pawrad(itypat),pawxcdev,trho1,0,2*usexcnhat,vxc_tmp,xclevel,xc_denpos) 674 else 675 write(msg,'(2a)')ch10,' pawdenpot : Computing valence-only v_xc[tn1+nhat] using angular mesh' 676 call wrtout(std_out,msg,'COLL') 677 call pawxc(pawtab(itypat)%tcoredens(:,1),& 678 & eexc_val,eexcdc_val,hyb_mixing_,ixc,kxc_tmp,k3xc_tmp,lm_size,& 679 & paw_an(iatom)%lmselect,nhat1,nkxc1,nk3xc1,non_magnetic_xc,mesh_size,nspden,option,& 680 & pawang,pawrad(itypat),trho1,0,2*usexcnhat,vxc_tmp,xclevel,xc_denpos) 681 end if 682 if (option<2) then 683 paw_an(iatom)%vxct1_val(:,:,:)=vxc_tmp(:,:,:) 684 paw_an(iatom)%has_vxcval=2 685 end if 686 end if ! valence-only XC potentials 687 688 ABI_FREE(vxc_tmp) 689 if (usekden==1) then 690 ABI_FREE(vxctau_tmp) 691 end if 692 ABI_FREE(kxc_tmp) 693 ABI_FREE(k3xc_tmp) 694 695 ! ===== Compute first part of local exact-exchange energy term ===== 696 ! ===== Also compute corresponding potential ===== 697 ! ================================================================== 698 699 if (pawtab(itypat)%useexexch/=0.and.ipert==0.and.ipositron/=1) then 700 701 ! ===== Re-compute a partial "on-site" density n1 (only l=lexexch contrib.) 702 ABI_MALLOC(rho1xx,(mesh_size,lm_size,nspden)) 703 ABI_MALLOC(lmselect_tmp,(lm_size)) 704 lmselect_tmp(:)=lmselect_cur(:) 705 call pawdensities(rdum,cplex,iatom_tot,lmselect_cur,lmselect_tmp,lm_size,rdum3,nspden,& 706 & 1,0,2,pawtab(itypat)%lexexch,0,pawang,pawprtvol,pawrad(itypat),& 707 & pawrhoij(iatom),pawtab(itypat),rho1xx,rdum3a,one_over_rad2=one_over_rad2) 708 ABI_FREE(lmselect_tmp) 709 ! ===== Re-compute Exc1 and Vxc1; for local exact-exchange, this is done in GGA only 710 ABI_MALLOC(vxc_tmp,(mesh_size,lm_size,nspden)) 711 ABI_MALLOC(kxc_tmp,(mesh_size,lm_size,nkxc1)) 712 call pawxcm(pawtab(itypat)%coredens,eextemp,eexdctemp,pawtab(itypat)%useexexch,hyb_mixing_,ixc,kxc_tmp,lm_size,& 713 & paw_an(iatom)%lmselect,nhat1,nkxc1,non_magnetic_xc,mesh_size,nspden,option,pawang,pawrad(itypat),pawxcdev,& 714 & rho1xx,0,0,vxc_tmp,xclevel,xc_denpos) 715 if (option/=1) then 716 e1xc=e1xc-eextemp*exchmix 717 e1xcdc=e1xcdc-eexdctemp*exchmix 718 end if 719 if (option<2) then 720 paw_an(iatom)%vxc_ex(:,:,:)=vxc_tmp(:,:,:) 721 paw_an(iatom)%has_vxc_ex=2 722 end if 723 ABI_FREE(rho1xx) 724 ABI_FREE(vxc_tmp) 725 ABI_FREE(kxc_tmp) 726 727 end if ! useexexch 728 729 itypat0=0;if (iatom<my_natom) itypat0=pawrhoij(iatom+1)%itypat 730 if (itypat/=itypat0) then 731 ABI_FREE(one_over_rad2) 732 end if 733 734 ABI_FREE(lmselect_cur) 735 736 ! ==== Compute Hartree potential terms and some energy terms ==== 737 ! =============================================================== 738 739 ! Hartree Dij computation 740 if (ipositron/=1) then 741 call pawdijhartree(paw_ij(iatom)%dijhartree,cplex,nspden,pawrhoij(iatom),pawtab(itypat)) 742 else 743 paw_ij(iatom)%dijhartree(:)=zero 744 end if 745 paw_ij(iatom)%has_dijhartree=2 746 747 ! Hartree energy computation 748 if (option/=1) then 749 call pawaccenergy_nospin(eh2,pawrhoij(iatom),paw_ij(iatom)%dijhartree,1,qphase,pawtab(itypat)) 750 end if 751 752 ! Electron-positron calculation: 753 ! - Compute Dij due to fixed particles (elec. or pos. depending on calctype) 754 ! - Compute contribution to energy 755 ! - Add electron and positron 756 if (ipositron/=0) then 757 ABI_CHECK(qphase==1,'qphase should be 1 for electron-positron!') 758 ABI_MALLOC(dij_ep,(qphase*lmn2_size)) 759 call pawdijhartree(dij_ep,qphase,nspden,electronpositron%pawrhoij_ep(iatom),pawtab(itypat)) 760 if (option/=1) then 761 etmp=zero 762 call pawaccenergy_nospin(etmp,pawrhoij(iatom),dij_ep,1,1,pawtab(itypat)) 763 electronpositron%e_paw =electronpositron%e_paw -etmp 764 electronpositron%e_pawdc=electronpositron%e_pawdc-etmp 765 end if 766 paw_ij(iatom)%dijhartree(:)=paw_ij(iatom)%dijhartree(:)-dij_ep(:) 767 ABI_FREE(dij_ep) 768 end if 769 770 ! Compute 1st moment of total Hartree potential VH(n_Z+n_core+n1) 771 ! Equation 10 (density) and up to 43 (Hartree potential of density) 772 ! of Kresse and Joubert PRB 59 1758 (1999) [[cite:Kresse1999]] 773 keep_vhartree=(paw_an(iatom)%has_vhartree>0) 774 if ((pawspnorb>0.and.ipert==0.and.ipositron/=1).or.keep_vhartree) then 775 776 !In the first clause case, would it not be simpler just to turn on has_vhartree? 777 if (.not. allocated(paw_an(iatom)%vh1)) then 778 ABI_MALLOC(paw_an(iatom)%vh1,(cplex*mesh_size,1,1)) 779 end if 780 if (.not. allocated(paw_an(iatom)%vht1)) then 781 ABI_MALLOC(paw_an(iatom)%vht1,(cplex*mesh_size,1,1)) 782 end if 783 ABI_MALLOC(rho,(mesh_size)) 784 ABI_MALLOC(vh,(mesh_size)) 785 786 ! Construct vh1 and tvh1 787 do iq=1,cplex 788 !Construct vh1 789 ! The sqrt(4pi) factor comes from the fact we are calculating the spherical moments, 790 ! and for the 00 channel the prefactor of Y_00 = 2 sqrt(pi) 791 rho(1:mesh_size)=rho1(iq:cplex*mesh_size:cplex,1,1)*four_pi*pawrad(itypat)%rad(1:mesh_size)**2 792 if (usecore==1) then 793 rho(1:mesh_size)=rho(1:mesh_size)+sqrt(four_pi)*pawtab(itypat)%coredens(1:mesh_size) & 794 & *four_pi*pawrad(itypat)%rad(1:mesh_size)**2 795 end if 796 call poisson(rho,0,pawrad(itypat),vh) 797 vh(2:mesh_size)=(vh(2:mesh_size)-sqrt(four_pi)*znucl(itypat))/pawrad(itypat)%rad(2:mesh_size) 798 call pawrad_deducer0(vh,mesh_size,pawrad(itypat)) 799 paw_an(iatom)%vh1(iq:cplex*mesh_size:cplex,1,1)=vh(1:mesh_size) 800 ! TODO: check this is equivalent to the previous version (commented) which explicitly recalculated VH(coredens) 801 ! DONE: numerically there are residual differences on abiref (7th digit). 802 ! paw_an(iatom)%vh1(2:mesh_size,1,1)=paw_an(iatom)%vh1(2:mesh_size,1,1)/pawrad(itypat)%rad(2:mesh_size) & 803 !& +sqrt(four_pi) * pawtab(itypat)%VHnZC(2:mesh_size) 804 805 !Same for vht1 806 rho(1:mesh_size)=trho1(iq:cplex*mesh_size:cplex,1,1)*four_pi*pawrad(itypat)%rad(1:mesh_size)**2 807 if (usenhat/=0) then 808 rho(1:mesh_size)=rho(1:mesh_size)+nhat1(iq:cplex*mesh_size:cplex,1,1) & 809 & *four_pi*pawrad(itypat)%rad(1:mesh_size)**2 810 end if 811 if (usetcore==1) then 812 rho(1:mesh_size)=rho(1:mesh_size)+sqrt(four_pi)*pawtab(itypat)%tcoredens(1:mesh_size,1) & 813 & *four_pi*pawrad(itypat)%rad(1:mesh_size)**2 814 end if 815 call poisson(rho,0,pawrad(itypat),vh) 816 vh(2:mesh_size)=(vh(2:mesh_size)-sqrt(four_pi)*znucl(itypat))/pawrad(itypat)%rad(2:mesh_size) 817 call pawrad_deducer0(vh,mesh_size,pawrad(itypat)) 818 paw_an(iatom)%vht1(iq:cplex*mesh_size:cplex,1,1)=vh(1:mesh_size) 819 820 end do ! cplex phase 821 822 paw_an(iatom)%has_vhartree=2 823 ABI_FREE(rho) 824 ABI_FREE(vh) 825 end if 826 827 ! ========= Compute PAW+U and energy contribution ========= 828 ! ========================================================== 829 830 if (usepawu/=0.and.usepawu<10.and.ipositron/=1.and.option/=1) then 831 832 if (pawu_algo==PAWU_ALGO_1) then 833 834 ! PAW+U energy computation from nocc_m_mp 835 call pawuenergy(iatom_tot,edftumdc,edftumdcdc,paw_ij(iatom)%noccmmp, & 836 & paw_ij(iatom)%nocctot,pawprtvol,pawtab(itypat)) 837 else 838 839 ! PAW+U energy computation from eU_ijkl 840 !First, compute DijU 841 call pawdiju_euijkl(paw_ij(iatom)%dijU,cplex_dij,cplex,ndij, & 842 & pawrhoij(iatom),pawtab(itypat)) 843 paw_ij(iatom)%has_dijU=2 844 !Then, compute energy 845 if (option/=1) then 846 etmp=zero 847 call pawaccenergy(etmp,pawrhoij(iatom),paw_ij(iatom)%dijU,cplex_dij,qphase,ndij,pawtab(itypat)) 848 edftumdc=edftumdc+half*etmp ; edftumdcdc=edftumdcdc-half*etmp 849 !Add FLL double-counting part 850 if (pawu_dblec==PAWU_FLL.and.ipert==0) then 851 ABI_CHECK(qphase==1,'BUG in pawdenpot: qphase should be 1 for Dble-C FLL term!') 852 call pawaccenergy_nospin(edftufll,pawrhoij(iatom),pawtab(itypat)%euij_fll,1,1,pawtab(itypat)) 853 end if 854 end if 855 856 end if ! DFT+U algo 857 end if ! Dij Hartree 858 859 ! ========= Compute nuclear dipole moment energy contribution ======== 860 ! ===================================================================== 861 862 if (any(abs(nucdipmom(:,iatom))>tol8).and.ipert==0.and.ipositron/=1) then 863 864 ABI_CHECK(cplex_rhoij==2,'BUG in pawdenpot: rhoij must be complex for ND moments!') 865 ABI_CHECK(qphase==1,'BUG in pawdenpot: qphase should be 1 for ND moments!') 866 867 ! Compute nuclear dipole contribution to Dij if necessary 868 if (paw_ij(iatom)%has_dijnd/=2) then 869 call pawdijnd(paw_ij(iatom)%dijnd,cplex_dij,ndij,nucdipmom(:,iatom),pawrad(itypat),pawtab(itypat)) 870 paw_ij(iatom)%has_dijnd=2 871 end if 872 873 ! Compute nuclear dipole contribution to energy 874 if (option/=1) then 875 call pawaccenergy_nospin(enucdip,pawrhoij(iatom),paw_ij(iatom)%dijnd,cplex_dij,1,pawtab(itypat)) 876 end if 877 878 end if 879 880 ! ========= Compute spin-orbit energy contribution ======== 881 ! ========================================================== 882 883 if (pawspnorb>0.and.ipert==0.and.ipositron/=1) then 884 885 ! Compute spin-orbit contribution to Dij 886 if (option/=2.or.cplex_rhoij==2) then 887 call pawdijso(paw_ij(iatom)%dijso,cplex_dij,cplex,ndij,nspden,pawang,pawrad(itypat),pawtab(itypat), & 888 & pawxcdev,spnorbscl,paw_an(iatom)%vh1,paw_an(iatom)%vxc1) 889 paw_ij(iatom)%has_dijso=2 890 end if 891 892 ! Compute spin-orbit contribution to on-site energy 893 if (option/=1.and.cplex_rhoij==2) then 894 call pawaccenergy(espnorb,pawrhoij(iatom),paw_ij(iatom)%dijso,cplex_dij,qphase,ndij,pawtab(itypat)) 895 end if 896 897 end if 898 899 ! === Compute 2nd part of local exact-exchange energy and potential === 900 ! ====================================================================== 901 902 if (pawtab(itypat)%useexexch/=0.and.ipert==0.and.ipositron/=1) then 903 904 ABI_CHECK(paw_ij(iatom)%nspden/=4,'BUG in pawdenpot: Local ex-exch. not implemented for nspden=4!') 905 if (option<2) then 906 call pawxpot(ndij,pawprtvol,pawrhoij(iatom),pawtab(itypat),paw_ij(iatom)%vpawx) 907 paw_ij(iatom)%has_exexch_pot=2 908 end if 909 if (option/=1) then 910 if (abs(pawprtvol)>=2) then 911 write(msg, '(2a)' )ch10,'======= PAW local exact exchange terms (in Hartree) ====' 912 call wrtout(std_out, msg,'COLL') 913 write(msg, '(2a,i4)' )ch10,' For Atom',iatom_tot 914 call wrtout(std_out, msg,'COLL') 915 end if 916 call pawxenergy(eexex,pawprtvol,pawrhoij(iatom),pawtab(itypat)) 917 end if 918 919 end if ! useexexch 920 921 ! ==== Compute Fock Dij term and Fock energy terms ==== 922 ! ===================================================== 923 924 if (usefock==1) then 925 926 if (ipositron/=1) then 927 928 ! Fock contribution to Dij 929 ABI_MALLOC(dijfock_vv,(cplex_dij*qphase*lmn2_size,ndij)) 930 ABI_MALLOC(dijfock_cv,(cplex_dij*qphase*lmn2_size,ndij)) 931 call pawdijfock(dijfock_vv,dijfock_cv,cplex_dij,cplex,hyb_mixing_,hyb_mixing_sr_, & 932 & ndij,pawrhoij(iatom),pawtab(itypat)) 933 paw_ij(iatom)%dijfock(:,:)=dijfock_vv(:,:)+dijfock_cv(:,:) 934 paw_ij(iatom)%has_dijfock=2 935 936 ! Fock contribution to energy 937 if (option/=1) then 938 dijfock_vv(:,:)=half*dijfock_vv(:,:) ; dijfock_cv(:,:)=dijfock_vv(:,:)+dijfock_cv(:,:) 939 call pawaccenergy(efock ,pawrhoij(iatom),dijfock_cv,cplex_dij,qphase,ndij,pawtab(itypat)) 940 call pawaccenergy(efockdc,pawrhoij(iatom),dijfock_vv,cplex_dij,qphase,ndij,pawtab(itypat)) 941 end if 942 943 ABI_FREE(dijfock_vv) 944 ABI_FREE(dijfock_cv) 945 end if 946 947 ! Special case for positron 948 if (ipositron==1) then 949 paw_ij(iatom)%dijfock(:,:)=zero 950 paw_ij(iatom)%has_dijfock=2 951 end if 952 953 end if 954 955 ! === Compute the zero of the potentials if requested ================== 956 ! ====================================================================== 957 958 if (pawtab(itypat)%usepotzero==1.and.present(vpotzero).and.ipert<=0) then 959 960 !Term 1 : beta 961 vpotzero(1)=vpotzero(1)-pawtab(itypat)%beta/ucvol 962 963 !Term 2 : \sum_ij rho_ij gamma_ij 964 etmp=zero 965 call pawaccenergy_nospin(etmp,pawrhoij(iatom),pawtab(itypat)%gammaij,1,1,pawtab(itypat)) 966 vpotzero(2)=vpotzero(2)-etmp/ucvol 967 968 end if 969 970 ! ======= Compute atomic contribution to the energy (Dij0) =========== 971 ! ====================================================================== 972 973 if (option/=1.and.ipert<=0) then 974 call pawaccenergy_nospin(e1t10,pawrhoij(iatom),pawtab(itypat)%dij0,1,1,pawtab(itypat)) 975 976 ! Positron special case (dij0 is opposite, except for kinetic term) 977 if (ipositron==1) then 978 ABI_MALLOC(dij_ep,(lmn2_size)) 979 dij_ep(:)=two*(pawtab(itypat)%kij(:)-pawtab(itypat)%dij0(:)) 980 call pawaccenergy_nospin(e1t10,pawrhoij(iatom),dij_ep,1,1,pawtab(itypat)) 981 ABI_FREE(dij_ep) 982 end if 983 984 end if 985 986 ! ========================================================== 987 ! No more need of some densities/potentials 988 989 ! Deallocate densities 990 ABI_FREE(rho1) 991 ABI_FREE(trho1) 992 ABI_FREE(nhat1) 993 if (usekden==1) then 994 ABI_FREE(tau1) 995 ABI_FREE(ttau1) 996 end if 997 if (ipositron/=0) then 998 ABI_FREE(rho1_ep) 999 ABI_FREE(trho1_ep) 1000 ABI_FREE(nhat1_ep) 1001 end if 1002 1003 ! Deallocate potentials 1004 if (.not.keep_vhartree) then 1005 paw_an(iatom)%has_vhartree=0 1006 if (allocated(paw_an(iatom)%vh1)) then 1007 ABI_FREE(paw_an(iatom)%vh1) 1008 end if 1009 end if 1010 if (temp_vxc) then 1011 paw_an(iatom)%has_vxc=0 1012 if (allocated(paw_an(iatom)%vxc1)) then 1013 ABI_FREE(paw_an(iatom)%vxc1) 1014 end if 1015 end if 1016 1017 ! =========== End loop on atoms ============================ 1018 ! ========================================================== 1019 1020 end do 1021 1022 !========== Assemble "on-site" energy terms =============== 1023 !========================================================== 1024 1025 if (option/=1) then 1026 if (ipert==0) then 1027 epaw=e1xc+half*eh2+e1t10-exccore-etild1xc+edftumdc+edftufll+eexex+espnorb+efock+enucdip 1028 epawdc=e1xc-e1xcdc-half*eh2-exccore-etild1xc+etild1xcdc+edftumdcdc-eexex-efockdc 1029 else 1030 epaw=e1xc-etild1xc+eh2+two*edftumdc 1031 epawdc=zero 1032 end if 1033 end if 1034 1035 !========== Reduction in case of parallelism ============== 1036 !========================================================== 1037 1038 if (paral_atom) then 1039 if (option/=1) then 1040 call timab(48,1,tsec) 1041 mpiarr=zero 1042 mpiarr(1)=compch_sph;mpiarr(2)=epaw;mpiarr(3)=epawdc 1043 if (ipositron/=0) then 1044 mpiarr(4)=electronpositron%e_paw 1045 mpiarr(5)=electronpositron%e_pawdc 1046 end if 1047 if (present(vpotzero)) then 1048 mpiarr(6)=vpotzero(1) 1049 mpiarr(7)=vpotzero(2) 1050 end if 1051 call xmpi_sum(mpiarr,my_comm_atom,ierr) 1052 compch_sph=mpiarr(1);epaw=mpiarr(2);epawdc=mpiarr(3) 1053 if (ipositron/=0) then 1054 electronpositron%e_paw=mpiarr(4) 1055 electronpositron%e_pawdc=mpiarr(5) 1056 end if 1057 if (present(vpotzero)) then 1058 vpotzero(1)=mpiarr(6) 1059 vpotzero(2)=mpiarr(7) 1060 end if 1061 call timab(48,2,tsec) 1062 end if 1063 end if 1064 1065 !Destroy atom table used for parallelism 1066 call free_my_atmtab(my_atmtab,my_atmtab_allocated) 1067 1068 call timab(560,2,tsec) 1069 1070 DBG_EXIT("COLL") 1071 1072 end subroutine pawdenpot
m_paw_denpot/pawdensities [ Functions ]
[ Top ] [ m_paw_denpot ] [ Functions ]
NAME
pawdensities
FUNCTION
Compute PAW on-site densities (all-electron, pseudo and compensation) for a given atom
INPUTS
cplex: if 1, on-site densities are REAL, if 2, COMPLEX (response function only)
iatom=index of current atom (note: this is the absolute index, not the index on current proc)
lm_size=number of (l,m) moments
lmselectin(lm_size)=flags selecting the non-zero LM-moments of on-site densities
(value of these flags at input; must be .TRUE. for nzlmopt/=1)
nspden=number of spin-density components
nzlmopt=if -1, compute all LM-moments of densities (lmselectin=.true. forced)
initialize "lmselectout" (index of non-zero LM-moments of densities)
if 0, compute all LM-moments of densities (lmselectin=.true. forced)
force "lmselectout" to .true. (index of non-zero LM-moments of densities)
if 1, compute only non-zero LM-moments of densities (stored before in "lmselectin")
one_over_rad2(mesh_size)= contains 1/r**2 for each point of the radial grid -optional argument-
opt_compch=flag controlling the accumulation of compensation charge density moments
inside PAW spheres (compch_sph)
opt_dens=flag controlling which on-site density(ies) is (are) computed
0: all on-site densities (all-electron, pseudo and compensation)
1: all-electron and pseudo densities (no compensation)
2: only all-electron density
opt_l=controls which l-moment(s) contribute to the density:
<0 : all l contribute
>=0: only l=opt_l contributes
Note: opt_l>=0 is only compatible with opt_dens=2
opt_print=1 if the densities moments have to be printed out (only if pawprtvol>=2)
pawang <type(pawang_type)>=paw angular mesh and related data
pawprtvol=control print volume and debugging output for PAW
pawrad <type(pawrad_type)>=paw radial mesh and related data (for the current atom type)
pawrhoij <type(pawrhoij_type)>= paw rhoij occupancies and related data (for the current atom)
pawtab <type(pawtab_type)>=paw tabulated starting data (for the current atom type)
OUTPUT
nhat1(cplex*mesh_size,lm_size,nspden)= compensation charge on-site density for current atom
rho1(cplex*mesh_size,lm_size,nspden)= all electron on-site density for current atom
trho1(cplex*mesh_size,lm_size,nspden)= pseudo on-site density for current atom
==== if nzlmopt/=1
lmselectout(lm_size)=flags selecting the non-zero LM-moments of on-site densities
(value of these flags at output if updated, i.e. if nzlmopt<1)
SIDE EFFECTS
==== if opt_compch==1
compch_sph=compensation charge integral inside spheres computed over spherical meshes
updated with the contribution of current atom
SOURCE
1129 subroutine pawdensities(compch_sph,cplex,iatom,lmselectin,lmselectout,lm_size,nhat1,nspden,nzlmopt,& 1130 & opt_compch,opt_dens,opt_l,opt_print,pawang,pawprtvol,pawrad,pawrhoij,pawtab,rho1,trho1,& 1131 & one_over_rad2) ! optional 1132 1133 !Arguments --------------------------------------------- 1134 !scalars 1135 integer,intent(in) :: cplex,iatom,lm_size,nspden,nzlmopt,opt_compch,opt_dens,opt_l,opt_print,pawprtvol 1136 real(dp),intent(inout) :: compch_sph 1137 type(pawang_type),intent(in) :: pawang 1138 type(pawrad_type),intent(in) :: pawrad 1139 type(pawrhoij_type),intent(in) :: pawrhoij 1140 type(pawtab_type),intent(in) :: pawtab 1141 !arrays 1142 logical,intent(in) :: lmselectin(lm_size) 1143 logical,intent(inout) :: lmselectout(lm_size) 1144 real(dp),intent(in),target,optional :: one_over_rad2(pawtab%mesh_size) 1145 real(dp),intent(out) :: nhat1(cplex*pawtab%mesh_size,lm_size,nspden*(1-((opt_dens+1)/2))) 1146 real(dp),intent(out) :: rho1(cplex*pawtab%mesh_size,lm_size,nspden) 1147 real(dp),intent(out) :: trho1(cplex*pawtab%mesh_size,lm_size,nspden*(1-(opt_dens/2))) 1148 !Local variables --------------------------------------- 1149 !scalars 1150 integer :: dplex,ii,ilm,iplex,iq0,ir,irhoij,isel,ispden,jrhoij 1151 integer :: klm,klmn,kln,ll,lmax,lmin,mesh_size 1152 real(dp) :: m1,mt1,rdum 1153 character(len=500) :: msg 1154 !arrays 1155 real(dp) :: compchspha(cplex),compchsphb(cplex),ro(cplex),ro_ql(cplex),ro_rg(cplex) 1156 real(dp),allocatable :: aa(:),bb(:) 1157 real(dp),pointer :: one_over_rad2_(:) 1158 1159 ! ************************************************************************* 1160 1161 DBG_ENTER("COLL") 1162 1163 !Compatibility tests 1164 if (opt_dens/=2.and.opt_l>=0) then 1165 msg='opt_dens/=2 incompatible with opt_l>=0!' 1166 ABI_BUG(msg) 1167 end if 1168 if(nzlmopt/=0.and.nzlmopt/=1.and.nzlmopt/=-1) then 1169 msg='invalid value for variable "nzlmopt"!' 1170 ABI_BUG(msg) 1171 end if 1172 if(nspden>pawrhoij%nspden) then 1173 msg='nspden must be <= pawrhoij%nspden!' 1174 ABI_BUG(msg) 1175 end if 1176 if (cplex>pawrhoij%qphase) then 1177 msg='cplex must be <= pawrhoij%qphase!' 1178 ABI_BUG(msg) 1179 end if 1180 if (nzlmopt/=1) then 1181 if (any(.not.lmselectin(1:lm_size))) then 1182 msg='With nzlmopt/=1, lmselectin must be true!' 1183 ABI_BUG(msg) 1184 end if 1185 end if 1186 if (pawang%gnt_option==0) then 1187 msg='pawang%gnt_option=0!' 1188 ABI_BUG(msg) 1189 end if 1190 1191 !Various inits 1192 rho1=zero 1193 if (opt_dens<2) trho1=zero 1194 if (opt_dens==0) nhat1=zero 1195 mesh_size=pawtab%mesh_size;dplex=cplex-1 1196 iq0=pawrhoij%cplex_rhoij*pawrhoij%lmn2_size 1197 if (nzlmopt<1) lmselectout(1:lm_size)=.true. 1198 if (present(one_over_rad2)) then 1199 one_over_rad2_ => one_over_rad2 1200 else 1201 ABI_MALLOC(one_over_rad2_,(mesh_size)) 1202 one_over_rad2_(1)=zero 1203 one_over_rad2_(2:mesh_size)=one/pawrad%rad(2:mesh_size)**2 1204 end if 1205 1206 !===== Compute "on-site" densities (n1, ntild1, nhat1) ===== 1207 !=========================================================== 1208 1209 do ispden=1,nspden 1210 1211 ! -- Loop over ij channels (basis components) 1212 jrhoij=1 1213 do irhoij=1,pawrhoij%nrhoijsel 1214 klmn=pawrhoij%rhoijselect(irhoij) 1215 klm =pawtab%indklmn(1,klmn) 1216 kln =pawtab%indklmn(2,klmn) 1217 lmin=pawtab%indklmn(3,klmn) 1218 lmax=pawtab%indklmn(4,klmn) 1219 1220 ! Retrieve rhoij 1221 if (pawrhoij%nspden/=2) then 1222 ro(1)=pawrhoij%rhoijp(jrhoij,ispden) 1223 if (cplex==2) ro(2)=pawrhoij%rhoijp(iq0+jrhoij,ispden) 1224 else 1225 if (ispden==1) then 1226 ro(1)=pawrhoij%rhoijp(jrhoij,1)+pawrhoij%rhoijp(jrhoij,2) 1227 if (cplex==2) ro(2)=pawrhoij%rhoijp(iq0+jrhoij,1)+pawrhoij%rhoijp(iq0+jrhoij,2) 1228 else if (ispden==2) then 1229 ro(1)=pawrhoij%rhoijp(jrhoij,1) 1230 if (cplex==2) ro(2)=pawrhoij%rhoijp(iq0+jrhoij,1) 1231 end if 1232 end if 1233 ro(1:cplex)=pawtab%dltij(klmn)*ro(1:cplex) 1234 1235 ! First option: all on-site densities are computed (opt_dens==0) 1236 ! -------------------------------------------------------------- 1237 if (opt_dens==0) then 1238 do ll=lmin,lmax,2 1239 do ilm=ll**2+1,(ll+1)**2 1240 if (lmselectin(ilm)) then 1241 isel=pawang%gntselect(ilm,klm) 1242 if (isel>0) then 1243 ro_ql(1:cplex)=ro(1:cplex)*pawtab%qijl(ilm,klmn) 1244 ro_rg(1:cplex)=ro(1:cplex)*pawang%realgnt(isel) 1245 ! == nhat1(r=0) 1246 nhat1(1:cplex,ilm,ispden)=nhat1(1:cplex,ilm,ispden) & 1247 & +ro_ql(1:cplex)*pawtab%shapefunc(1,ll+1) 1248 ! == rho1(r>0), trho1(r>0), nhat1(r>0) 1249 do ir=2,mesh_size 1250 rho1(cplex*ir-dplex:ir*cplex,ilm,ispden) =rho1(cplex*ir-dplex:ir*cplex,ilm,ispden)& 1251 & +ro_rg(1:cplex)*pawtab%phiphj(ir,kln)*one_over_rad2_(ir) 1252 trho1(cplex*ir-dplex:ir*cplex,ilm,ispden)=trho1(cplex*ir-dplex:ir*cplex,ilm,ispden)& 1253 & +ro_rg(1:cplex)*pawtab%tphitphj(ir,kln)*one_over_rad2_(ir) 1254 nhat1(cplex*ir-dplex:ir*cplex,ilm,ispden)=nhat1(cplex*ir-dplex:ir*cplex,ilm,ispden)& 1255 & +ro_ql(1:cplex)*pawtab%shapefunc(ir,ll+1) 1256 end do 1257 end if 1258 end if 1259 end do ! End loops over ll,lm 1260 end do 1261 1262 ! 2nd option: AE and pseudo densities are computed (opt_dens==1) 1263 ! -------------------------------------------------------------- 1264 else if (opt_dens==1) then 1265 do ll=lmin,lmax,2 1266 do ilm=ll**2+1,(ll+1)**2 1267 if (lmselectin(ilm)) then 1268 isel=pawang%gntselect(ilm,klm) 1269 if (isel>0) then 1270 ro_rg(1:cplex)=ro(1:cplex)*pawang%realgnt(isel) 1271 ! == rho1(r>0), trho1(r>0) 1272 do ir=2,mesh_size 1273 rho1(cplex*ir-dplex:ir*cplex,ilm,ispden) =rho1(cplex*ir-dplex:ir*cplex,ilm,ispden)& 1274 & +ro_rg(1:cplex)*pawtab%phiphj (ir,kln)*one_over_rad2_(ir) 1275 trho1(cplex*ir-dplex:ir*cplex,ilm,ispden)=trho1(cplex*ir-dplex:ir*cplex,ilm,ispden)& 1276 & +ro_rg(1:cplex)*pawtab%tphitphj(ir,kln)*one_over_rad2_(ir) 1277 end do 1278 end if 1279 end if 1280 end do ! End loops over ll,lm 1281 end do 1282 1283 ! 3rd option: only all-electron on-site density is computed (opt_dens==2) 1284 ! ----------------------------------------------------------------------- 1285 else if (opt_dens==2) then 1286 if (opt_l<0.or.(pawtab%indklmn(3,klmn)==0.and.pawtab%indklmn(4,klmn)==2*opt_l)) then 1287 do ll=lmin,lmax,2 1288 do ilm=ll**2+1,(ll+1)**2 1289 if (lmselectin(ilm)) then 1290 isel=pawang%gntselect(ilm,klm) 1291 if (isel>0) then 1292 ro_rg(1:cplex)=ro(1:cplex)*pawang%realgnt(isel) 1293 ! == rho1(r>0) 1294 do ir=2,mesh_size 1295 rho1(cplex*ir-dplex:ir*cplex,ilm,ispden) =rho1(cplex*ir-dplex:ir*cplex,ilm,ispden)& 1296 & +ro_rg(1:cplex)*pawtab%phiphj(ir,kln)*one_over_rad2_(ir) 1297 end do 1298 end if 1299 end if 1300 end do ! End loops over ll, lm 1301 end do 1302 end if 1303 end if 1304 1305 1306 ! -- End loop over ij channels 1307 jrhoij=jrhoij+pawrhoij%cplex_rhoij 1308 end do 1309 1310 ! Compute rho1(r=0) and trho1(r=0) 1311 if (cplex==2) then 1312 ABI_MALLOC(aa,(5)) 1313 ABI_MALLOC(bb,(5)) 1314 end if 1315 if (opt_dens==0.or.opt_dens==1) then 1316 do ll=0,pawtab%lcut_size-1 1317 do ilm=ll**2+1,(ll+1)**2 1318 if (lmselectin(ilm)) then 1319 if (cplex==1) then 1320 call pawrad_deducer0(rho1 (:,ilm,ispden),mesh_size,pawrad) 1321 call pawrad_deducer0(trho1(:,ilm,ispden),mesh_size,pawrad) 1322 else 1323 do ii=0,1 1324 do ir=2,5 1325 aa(ir)=rho1 (2*ir-ii,ilm,ispden) 1326 bb(ir)=trho1(2*ir-ii,ilm,ispden) 1327 end do 1328 call pawrad_deducer0(aa,5,pawrad) 1329 call pawrad_deducer0(bb,5,pawrad) 1330 rho1 (2-ii,ilm,ispden)=aa(1) 1331 trho1(2-ii,ilm,ispden)=bb(1) 1332 end do 1333 end if 1334 end if 1335 end do 1336 end do 1337 else 1338 do ll=0,pawtab%lcut_size-1 1339 do ilm=ll**2+1,(ll+1)**2 1340 if (lmselectin(ilm)) then 1341 if (cplex==1) then 1342 call pawrad_deducer0(rho1(:,ilm,ispden),mesh_size,pawrad) 1343 else 1344 do ii=0,1 1345 do ir=2,5 1346 aa(ir)=rho1 (2*ir-ii,ilm,ispden) 1347 end do 1348 call pawrad_deducer0(aa,5,pawrad) 1349 rho1(2-ii,ilm,ispden)=aa(1) 1350 end do 1351 end if 1352 end if 1353 end do 1354 end do 1355 end if 1356 if (cplex==2) then 1357 ABI_FREE(aa) 1358 ABI_FREE(bb) 1359 end if 1360 1361 ! -- Test moments of densities and store non-zero ones 1362 if (nzlmopt==-1) then 1363 do ll=0,pawtab%lcut_size-1 1364 do ilm=ll**2+1,(ll+1)**2 1365 m1=zero;mt1=zero 1366 if (cplex==1) then 1367 m1=maxval(abs(rho1 (1:mesh_size,ilm,ispden))) 1368 if (opt_dens<2) mt1=maxval(abs(trho1(1:mesh_size,ilm,ispden))) 1369 else 1370 do ir=1,mesh_size 1371 rdum=sqrt(rho1(2*ir-1,ilm,ispden)**2+rho1(2*ir,ilm,ispden)**2) 1372 m1=max(m1,rdum) 1373 end do 1374 if (opt_dens<2) then 1375 do ir=1,mesh_size 1376 rdum=sqrt(trho1(2*ir-1,ilm,ispden)**2+trho1(2*ir,ilm,ispden)**2) 1377 mt1=max(mt1,rdum) 1378 end do 1379 end if 1380 end if 1381 if (ispden==1) then 1382 if ((ilm>1).and.(m1<tol16).and.(mt1<tol16)) then 1383 lmselectout(ilm)=.false. 1384 end if 1385 else if (.not.(lmselectout(ilm))) then 1386 lmselectout(ilm)=((m1>=tol16).or.(mt1>=tol16)) 1387 end if 1388 end do 1389 end do 1390 end if 1391 1392 ! -- Compute integral of (n1-tn1) inside spheres 1393 if (opt_compch==1.and.ispden==1.and.opt_dens<2) then 1394 ABI_MALLOC(aa,(mesh_size)) 1395 aa(1:mesh_size)=(rho1(1:mesh_size,1,1)-trho1(1:mesh_size,1,1)) & 1396 & *pawrad%rad(1:mesh_size)**2 1397 call simp_gen(compchspha(1),aa,pawrad) 1398 compch_sph=compch_sph+compchspha(1)*sqrt(four_pi) 1399 ABI_FREE(aa) 1400 end if 1401 1402 ! -- Print out moments of densities (if requested) 1403 if (abs(pawprtvol)>=2.and.opt_print==1.and.opt_dens<2) then 1404 ABI_MALLOC(aa,(cplex*mesh_size)) 1405 ABI_MALLOC(bb,(cplex*mesh_size)) 1406 if (opt_dens==0) then 1407 write(msg,'(2a,i3,a,i1,3a)') ch10, & 1408 & ' Atom ',iatom,' (ispden=',ispden,'):',ch10,& 1409 & ' ******* Moment of (n1-tn1) ** Moment of (n1-tn1-nhat1)' 1410 else 1411 write(msg,'(2a,i3,a,i1,3a)') ch10, & 1412 & ' Atom ',iatom,' (ispden=',ispden,'):',ch10,& 1413 & ' ******* Moment of (n1-tn1)' 1414 end if 1415 call wrtout(std_out,msg,'PERS') 1416 do ll=0,pawtab%lcut_size-1 1417 do ilm=ll**2+1,(ll+1)**2 1418 if (lmselectin(ilm)) then 1419 do iplex=1,cplex 1420 if (opt_dens==0) then 1421 do ir=1,mesh_size 1422 ii=cplex*(ir-1)+iplex 1423 ro(1)=pawrad%rad(ir)**(2+ll) 1424 aa(ir)=ro(1)*(rho1(ii,ilm,ispden)-trho1(ii,ilm,ispden)) 1425 bb(ir)=ro(1)*nhat1(ii,ilm,ispden) 1426 end do 1427 call simp_gen(compchspha(iplex),aa,pawrad) 1428 call simp_gen(compchsphb(iplex),bb,pawrad) 1429 else 1430 do ir=1,mesh_size 1431 ii=cplex*(ir-1)+iplex 1432 ro(1)=pawrad%rad(ir)**(2+ll) 1433 aa(ir)=ro(1)*(rho1(ii,ilm,ispden)-trho1(ii,ilm,ispden)) 1434 end do 1435 call simp_gen(compchspha(iplex),aa,pawrad) 1436 end if 1437 end do 1438 if (opt_dens==0) then 1439 if (cplex==1) then 1440 write(msg,'(3x,a,2i2,2(a,es14.7))') & 1441 & 'l,m=',ll,ilm-(ll**2+ll+1),': M=',compchspha(1),& 1442 & ' ** M=',compchspha(1)-compchsphb(1) 1443 else 1444 write(msg,'(3x,a,2i2,2(a,2es14.7))') & 1445 & 'l,m=',ll,ilm-(ll**2+ll+1),': M=',compchspha(1:2),& 1446 & ' ** M=',compchspha(1:2)-compchsphb(1:2) 1447 end if 1448 else 1449 if (cplex==1) then 1450 write(msg,'(3x,a,2i2,a,es14.7)') & 1451 & 'l,m=',ll,ilm-(ll**2+ll+1),': M=',compchspha(1) 1452 else 1453 write(msg,'(3x,a,2i2,a,2es14.7)') & 1454 & 'l,m=',ll,ilm-(ll**2+ll+1),': M=',compchspha(1:2) 1455 end if 1456 end if 1457 call wrtout(std_out,msg,'PERS') 1458 end if 1459 end do 1460 end do 1461 ABI_FREE(aa) 1462 ABI_FREE(bb) 1463 end if 1464 1465 ! ----- End loop over spin components 1466 end do 1467 1468 if (.not.present(one_over_rad2)) then 1469 ABI_FREE(one_over_rad2_) 1470 end if 1471 1472 DBG_EXIT("COLL") 1473 1474 end subroutine pawdensities
m_paw_denpot/pawkindensities [ Functions ]
[ Top ] [ m_paw_denpot ] [ Functions ]
NAME
pawkindensities
FUNCTION
Compute PAW on-site kinetic energy densities (all-electron, pseudo) for a given atom
INPUTS
cplex: if 1, on-site densities are REAL, if 2, COMPLEX (response function only) lm_size=number of (l,m) moments lmselectin(lm_size)=flags selecting the non-zero LM-moments of on-site kinetic energy densities (value of these flags at input; must be .TRUE. for nzlmopt/=1) nspden=number of spin-density components nzlmopt=if -1, compute all LM-moments of densities (lmselectin=.true. forced) initialize "lmselectout" (index of non-zero LM-moments of densities) if 0, compute all LM-moments of densities (lmselectin=.true. forced) force "lmselectout" to .true. (index of non-zero LM-moments of densities) if 1, compute only non-zero LM-moments of densities (stored before in "lmselectin") one_over_rad2(mesh_size)= contains 1/r**2 for each point of the radial grid -optional argument- opt_dens=flag controlling which on-site kinetic energy density(ies) is (are) computed 0,1: all-electron and pseudo on-site kinetic energy densities 2: only all-electron density opt_l=controls which l-moment(s) contribute to the kinetic energy density: <0 : all l contribute >=0: only l=opt_l contributes Note: opt_l>=0 is only compatible with opt_dens=2 pawang <type(pawang_type)>=paw angular mesh and related data pawrad <type(pawrad_type)>=paw radial mesh and related data (for the current atom type) pawrhoij <type(pawrhoij_type)>= paw rhoij occupancies and related data (for the current atom) pawtab <type(pawtab_type)>=paw tabulated starting data (for the current atom type)
OUTPUT
tau1(cplex*mesh_size,lm_size,nspden)= on site kinetic energy density ttau1(cplex*mesh_size,lm_size,nspden)]= pseudo on site kinetic energy density
SOURCE
1516 subroutine pawkindensities(cplex,lmselectin,lm_size,nspden,nzlmopt,& 1517 & opt_dens,opt_l,pawang,pawrad,pawrhoij,pawtab,tau1,ttau1,& 1518 & one_over_rad2) ! optional 1519 1520 !Arguments --------------------------------------------- 1521 !scalars 1522 integer,intent(in) :: cplex,lm_size,nspden,nzlmopt,opt_dens,opt_l 1523 type(pawang_type),intent(in) :: pawang 1524 type(pawrad_type),intent(in) :: pawrad 1525 type(pawrhoij_type),intent(in) :: pawrhoij 1526 type(pawtab_type),intent(in) :: pawtab 1527 !arrays 1528 logical,intent(in) :: lmselectin(lm_size) 1529 real(dp),intent(in),target,optional :: one_over_rad2(pawtab%mesh_size) 1530 real(dp),intent(out),optional :: tau1(cplex*pawtab%mesh_size,lm_size,nspden*(1-(opt_dens/2))) 1531 real(dp),intent(out),optional :: ttau1(cplex*pawtab%mesh_size,lm_size,nspden*(1-(opt_dens/2))) 1532 !Local variables --------------------------------------- 1533 !scalars 1534 integer :: dplex,ii,iq0,ir,irhoij,isel,ispden,jrhoij 1535 integer :: ilmn,ilm,ilm1,iln,jlmn,jlm1,jln,klm,klmn,ll,lmax,lmin,mesh_size 1536 real(dp) :: phiphj,tphitphj 1537 character(len=500) :: msg 1538 !arrays 1539 real(dp) :: ro(cplex),ro_rg(cplex) 1540 real(dp),allocatable :: aa(:),bb(:) 1541 real(dp),pointer :: one_over_rad2_(:) 1542 1543 ! ************************************************************************* 1544 1545 DBG_ENTER("COLL") 1546 1547 !Compatibility tests 1548 if (nzlmopt/=-1) then 1549 msg='nzlmopt/=-1 has not not been tested (might be wrong)!' 1550 ABI_BUG(msg) 1551 end if 1552 if (opt_dens/=2.and.opt_l>=0) then 1553 msg='opt_dens/=2 incompatible with opt_l>=0!' 1554 ABI_BUG(msg) 1555 end if 1556 if(nzlmopt/=0.and.nzlmopt/=1.and.nzlmopt/=-1) then 1557 msg='invalid value for variable "nzlmopt"!' 1558 ABI_BUG(msg) 1559 end if 1560 if(nspden>pawrhoij%nspden) then 1561 msg='nspden must be <= pawrhoij%nspden!' 1562 ABI_BUG(msg) 1563 end if 1564 if (cplex>pawrhoij%qphase) then 1565 msg='cplex must be <= pawrhoij%qphase!' 1566 ABI_BUG(msg) 1567 end if 1568 if (nzlmopt/=1) then 1569 if (any(.not.lmselectin(1:lm_size))) then 1570 msg='With nzlmopt/=1, lmselectin must be true!' 1571 ABI_BUG(msg) 1572 end if 1573 end if 1574 if (pawang%gnt_option==0) then 1575 msg='pawang%gnt_option=0!' 1576 ABI_BUG(msg) 1577 end if 1578 if (pawang%nabgnt_option==0) then 1579 msg='pawang%nabgnt_option=0!' 1580 ABI_BUG(msg) 1581 end if 1582 1583 !Various inits 1584 tau1=zero 1585 if (opt_dens<2) ttau1=zero 1586 mesh_size=pawtab%mesh_size;dplex=cplex-1 1587 iq0=pawrhoij%cplex_rhoij*pawrhoij%lmn2_size 1588 if (present(one_over_rad2)) then 1589 one_over_rad2_ => one_over_rad2 1590 else 1591 ABI_MALLOC(one_over_rad2_,(mesh_size)) 1592 one_over_rad2_(1)=zero 1593 one_over_rad2_(2:mesh_size)=one/pawrad%rad(2:mesh_size)**2 1594 end if 1595 1596 !=== Compute "on-site" kin. energy densities (n1, ntild1) ===== 1597 !============================================================== 1598 1599 do ispden=1,nspden 1600 1601 ! -- Loop over ij channels (basis components) 1602 jrhoij=1 1603 do irhoij=1,pawrhoij%nrhoijsel 1604 klmn=pawrhoij%rhoijselect(irhoij) 1605 klm =pawtab%indklmn(1,klmn) 1606 lmin=pawtab%indklmn(3,klmn) 1607 lmax=pawtab%indklmn(4,klmn) 1608 ilmn=pawtab%indklmn(7,klmn) ! (l,m,n) orbital 1 1609 jlmn=pawtab%indklmn(8,klmn) ! (l,m,n) orbital 2 1610 ilm1=pawtab%indklmn(5,klmn) ! (l,m) orbital 1 1611 jlm1=pawtab%indklmn(6,klmn) ! (l,m) orbital 2 1612 iln=pawtab%indlmn(5,ilmn) ! (l,n) orbital 1 1613 jln=pawtab%indlmn(5,jlmn) ! (l,n) orbital 2 1614 1615 ! Retrieve rhoij 1616 if (pawrhoij%nspden/=2) then 1617 ro(1)=pawrhoij%rhoijp(jrhoij,ispden) 1618 if (cplex==2) ro(2)=pawrhoij%rhoijp(iq0+jrhoij,ispden) 1619 else 1620 if (ispden==1) then 1621 ro(1)=pawrhoij%rhoijp(jrhoij,1)+pawrhoij%rhoijp(jrhoij,2) 1622 if (cplex==2) ro(2)=pawrhoij%rhoijp(iq0+jrhoij,1)+pawrhoij%rhoijp(iq0+jrhoij,2) 1623 else if (ispden==2) then 1624 ro(1)=pawrhoij%rhoijp(jrhoij,1) 1625 if (cplex==2) ro(2)=pawrhoij%rhoijp(iq0+jrhoij,1) 1626 end if 1627 end if 1628 ! Apply factor 1/2 (because tau=1/2 * Sum_ij[rhoij.Nabla_phi_i*Nabla_phi_j]) 1629 ro(1:cplex)=half*pawtab%dltij(klmn)*ro(1:cplex) 1630 1631 ! First option: AE and PS on-site kin. energy densities (opt_dens==0 or 1) 1632 ! ------------------------------------------------------------------------ 1633 if (opt_dens==0.or.opt_dens==1) then 1634 1635 ! Compute part of tau_lm depending on gaunt coefficients 1636 do ll=lmin,lmax,2 1637 do ilm=ll**2+1,(ll+1)**2 1638 if (lmselectin(ilm)) then 1639 isel=pawang%gntselect(ilm,klm) 1640 if (isel>0) then 1641 ro_rg(1:cplex)=ro(1:cplex)*pawang%realgnt(isel) 1642 do ir=2,mesh_size 1643 phiphj=pawtab%nablaphi(ir,iln)*pawtab%nablaphi(ir,jln) 1644 tphitphj=pawtab%tnablaphi(ir,iln)*pawtab%tnablaphi(ir,jln) 1645 tau1(cplex*ir-dplex:ir*cplex,ilm,ispden)=tau1(cplex*ir-dplex:ir*cplex,ilm,ispden)& 1646 & +ro_rg(1:cplex)*phiphj*one_over_rad2_(ir) 1647 ttau1(cplex*ir-dplex:ir*cplex,ilm,ispden)=ttau1(cplex*ir-dplex:ir*cplex,ilm,ispden)& 1648 & +ro_rg(1:cplex)*tphitphj*one_over_rad2_(ir) 1649 end do 1650 end if 1651 end if 1652 end do ! End loops over ll,lm 1653 end do 1654 1655 ! Compute the part of tau_lm depending on nablagaunt coefficients 1656 do ilm=1,(1+lmax)**2 1657 if (lmselectin(ilm)) then 1658 isel=pawang%nablagntselect(ilm,ilm1,jlm1) 1659 if (isel>0) then 1660 ro_rg(1:cplex)=ro(1:cplex)*pawang%nablarealgnt(isel) 1661 do ir=2,mesh_size 1662 phiphj=pawtab%phi(ir,iln)*pawtab%phi(ir,jln) 1663 tphitphj=pawtab%tphi(ir,iln)*pawtab%tphi(ir,jln) 1664 tau1(cplex*ir-dplex:ir*cplex,ilm,ispden)=tau1(cplex*ir-dplex:ir*cplex,ilm,ispden)& 1665 & +ro_rg(1:cplex)*phiphj*one_over_rad2_(ir)**2 1666 ttau1(cplex*ir-dplex:ir*cplex,ilm,ispden)=ttau1(cplex*ir-dplex:ir*cplex,ilm,ispden)& 1667 & +ro_rg(1:cplex)*tphitphj*one_over_rad2_(ir)**2 1668 end do 1669 end if 1670 end if 1671 end do 1672 1673 ! 2nd option: AE on-site kinetic energy density only (opt_dens==2) 1674 ! ---------------------------------------------------------------- 1675 else if (opt_dens==2) then 1676 1677 ! Compute part of tau_lm depending on gaunt coefficients 1678 do ll=lmin,lmax,2 1679 do ilm=ll**2+1,(ll+1)**2 1680 if (lmselectin(ilm)) then 1681 isel=pawang%gntselect(ilm,klm) 1682 if (isel>0) then 1683 ro_rg(1:cplex)=ro(1:cplex)*pawang%realgnt(isel) 1684 do ir=2,mesh_size 1685 phiphj=pawtab%nablaphi(ir,iln)*pawtab%nablaphi(ir,jln) 1686 tau1(cplex*ir-dplex:ir*cplex,ilm,ispden)=tau1(cplex*ir-dplex:ir*cplex,ilm,ispden)& 1687 & +ro_rg(1:cplex)*phiphj*one_over_rad2_(ir) 1688 end do 1689 end if 1690 end if 1691 end do ! End loops over ll,lm 1692 end do 1693 1694 ! Compute the part of tau_lm depending on nablagaunt coefficients 1695 do ilm=1,(1+lmax)**2 1696 if (lmselectin(ilm)) then 1697 isel=pawang%nablagntselect(ilm,ilm1,jlm1) 1698 if (isel>0) then 1699 ro_rg(1:cplex)=ro(1:cplex)*pawang%nablarealgnt(isel) 1700 do ir=2,mesh_size 1701 phiphj=pawtab%phi(ir,iln)*pawtab%phi(ir,jln) 1702 tau1(cplex*ir-dplex:ir*cplex,ilm,ispden)=tau1(cplex*ir-dplex:ir*cplex,ilm,ispden)& 1703 & +ro_rg(1:cplex)*phiphj*one_over_rad2_(ir)**2 1704 end do 1705 end if 1706 end if 1707 end do 1708 1709 end if 1710 1711 ! -- End loop over ij channels 1712 jrhoij=jrhoij+pawrhoij%cplex_rhoij 1713 end do 1714 1715 ! Compute tau1(r=0) and ttau1(r=0) 1716 if (cplex==2) then 1717 ABI_MALLOC(aa,(5)) 1718 ABI_MALLOC(bb,(5)) 1719 end if 1720 if (opt_dens==0.or.opt_dens==1) then 1721 do ll=0,pawtab%lcut_size-1 1722 do ilm=ll**2+1,(ll+1)**2 1723 if (lmselectin(ilm)) then 1724 if (cplex==1) then 1725 call pawrad_deducer0(tau1 (:,ilm,ispden),mesh_size,pawrad) 1726 call pawrad_deducer0(ttau1(:,ilm,ispden),mesh_size,pawrad) 1727 else 1728 do ii=0,1 1729 do ir=2,5 1730 aa(ir)=tau1 (2*ir-ii,ilm,ispden) 1731 bb(ir)=ttau1(2*ir-ii,ilm,ispden) 1732 end do 1733 call pawrad_deducer0(aa,5,pawrad) 1734 call pawrad_deducer0(bb,5,pawrad) 1735 tau1 (2-ii,ilm,ispden)=aa(1) 1736 ttau1(2-ii,ilm,ispden)=bb(1) 1737 end do 1738 end if 1739 end if 1740 end do 1741 end do 1742 else if (opt_dens==2) then 1743 do ll=0,pawtab%lcut_size-1 1744 do ilm=ll**2+1,(ll+1)**2 1745 if (lmselectin(ilm)) then 1746 if (cplex==1) then 1747 call pawrad_deducer0(tau1(:,ilm,ispden),mesh_size,pawrad) 1748 else 1749 do ii=0,1 1750 do ir=2,5 1751 aa(ir)=tau1(2*ir-ii,ilm,ispden) 1752 end do 1753 call pawrad_deducer0(aa,5,pawrad) 1754 tau1(2-ii,ilm,ispden)=aa(1) 1755 end do 1756 end if 1757 end if 1758 end do 1759 end do 1760 end if 1761 if (cplex==2) then 1762 ABI_FREE(aa) 1763 ABI_FREE(bb) 1764 end if 1765 1766 ! ----- End loop over spin components 1767 end do 1768 1769 if (.not.present(one_over_rad2)) then 1770 ABI_FREE(one_over_rad2_) 1771 end if 1772 1773 DBG_EXIT("COLL") 1774 1775 end subroutine pawkindensities