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