Line data Source code
1 : !!****m* ABINIT/m_paw_denpot
2 : !! NAME
3 : !! m_paw_denpot
4 : !!
5 : !! FUNCTION
6 : !! This module contains routines related to PAW on-site densities and on-site potentials.
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 2018-2026 ABINIT group (FJ, MT)
10 : !! This file is distributed under the terms of the
11 : !! GNU General Public License, see ~abinit/COPYING
12 : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : !!
14 : !! SOURCE
15 :
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 defs_abitypes
26 : use m_dtset
27 : use defs_datatypes, only : pseudopotential_type
28 : use m_paw_init
29 : use m_abicore
30 : use m_errors
31 : use m_xmpi
32 : use m_time, only : timab
33 :
34 : use m_rcpaw
35 : use m_paw_atom_solve
36 : use m_extfpmd, only : extfpmd_type
37 : use m_paw_atomorb, only : orb_relaxed_core
38 : use m_pawang, only : pawang_type
39 : use m_pawrad, only : pawrad_type,pawrad_deducer0,poisson,simp_gen
40 : use m_pawtab, only : pawtab_type
41 : use m_paw_an, only : paw_an_type
42 : use m_paw_ij, only : paw_ij_type
43 : use m_pawfgrtab, only : pawfgrtab_type
44 : use m_pawrhoij, only : pawrhoij_type
45 : use m_pawdij, only : pawdijhartree,pawdiju_euijkl,pawdijnd,pawdijso,pawxpot,pawdijfock,symdij,symdij_all
46 : use m_pawxc, only : pawxc,pawxc_dfpt,pawxcm,pawxcm_dfpt,pawxcpositron,pawxcmpositron, &
47 : & pawxc_get_usekden
48 : use m_paw_energies, only : paw_energies_type,paw_energies_setzero,paw_energies_to_array,n_paw_energies
49 : use m_paw_finegrid, only : pawgylm
50 : use m_paral_atom, only : get_my_atmtab,free_my_atmtab
51 : use m_paw_correlations, only : pawuenergy,pawxenergy,setnoccmmp
52 : use m_paral_atom, only : get_my_atmtab,free_my_atmtab
53 :
54 : use m_crystal, only : crystal_t
55 : use m_electronpositron, only : electronpositron_type,electronpositron_calctype
56 :
57 : #ifdef HAVE_FC_ISO_C_BINDING
58 : use, intrinsic :: iso_c_binding, only : c_ptr,c_loc,c_f_pointer
59 : #endif
60 :
61 : implicit none
62 :
63 : private
64 :
65 : !public procedures.
66 : public :: pawdenpot ! Compute different (PAW) energies, densities and potentials inside PAW spheres
67 : public :: pawdensities ! Compute PAW on-site densities (all-electron, pseudo and compensation)
68 : public :: pawkindensities ! Compute PAW on-site kinetic energy densities (all-electron, pseudo)
69 : public :: pawaccenergy ! Accumulate the atomic contribution of a PAW on-site energy
70 : public :: pawaccenergy_nospin ! As pawaccenergy, but with no spin polarization
71 : public :: paw_mknewh0 ! Compute bare PAW on-site Hamiltonian (-> GW calculations)
72 : public :: paw_relax_core ! Relax PAW core
73 : !public :: paw_energies_setzero ! Set all energies in a paw_energies datastructure to zero
74 :
75 : CONTAINS !========================================================================================
76 : !!***
77 :
78 : !----------------------------------------------------------------------
79 :
80 : !!****f* m_paw_denpot/pawdenpot
81 : !! NAME
82 : !! pawdenpot
83 : !!
84 : !! FUNCTION
85 : !! Compute different (PAW) energies, densities and potentials (or potential-like quantities)
86 : !! inside PAW spheres
87 : !! Can also compute first-order densities potentials and second-order energies (RF calculations).
88 : !!
89 : !! INPUTS
90 : !! electronpositron <type(electronpositron_type)>=quantities for the electron-positron annihilation (optional argument)
91 : !! extfpmd <type(extfpmd_type)>=--optional--extended first-principles molecular dynamics type
92 : !! [hyb_mixing, hyb_mixing_sr]= -- optional-- mixing factors for the global (resp. screened) XC hybrid functional
93 : !! el_temp=electronic temperature (hartree)
94 : !! gprimd=reciprocal lattice translations
95 : !! ipert=index of perturbation (used only for RF calculation ; set ipert<=0 for GS calculations.
96 : !! ixc= choice of exchange-correlation scheme (see above, and below)
97 : !! mpi_atmtab(:)=--optional-- indexes of the atoms treated by current proc
98 : !! comm_atom=--optional-- MPI communicator over atoms
99 : !! my_natom=number of atoms treated by current processor
100 : !! natom=total number of atoms in cell
101 : !! nspden=number of spin-density components
102 : !! ntypat=number of types of atoms in unit cell.
103 : !! nucdipmom(3,natom) nuclear dipole moments
104 : !! nzlmopt= if -1, compute all LM-moments of densities
105 : !! initialize "lmselect" (index of non-zero LM-moments of densities)
106 : !! if 0, compute all LM-moments of densities
107 : !! force "lmselect" to .true. (index of non-zero LM-moments of densities)
108 : !! if 1, compute only non-zero LM-moments of densities (stored before)
109 : !! option=0: compute both energies and potentials
110 : !! 1: compute only potentials
111 : !! 2: compute only energies
112 : !! paw_an(my_natom) <type(paw_an_type)>=paw arrays given on angular mesh
113 : !! paw_an0(my_natom) <type(paw_an_type)>=paw arrays given on angular mesh for Ground-State
114 : ! used only if ipert>0; must be set equal to paw_an for GS calc.
115 : !! paw_ij(my_natom) <type(paw_ij_type)>=paw arrays given on (i,j) channels
116 : !! pawang <type(pawang_type)>=paw angular mesh and related data
117 : !! pawprtvol=control print volume and debugging output for PAW
118 : !! pawrad(ntypat) <type(pawrad_type)>=paw radial mesh and related data
119 : !! pawrhoij(my_natom) <type(pawrhoij_type)>= paw rhoij occupancies and related data
120 : !! pawspnorb=flag: 1 if spin-orbit coupling is activated
121 : !! pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data
122 : !! pawxcdev=Choice of XC development (0=no dev. (use of angular mesh) ; 1 or 2=dev. on moments)
123 : !! ucvol=unit cell volume (bohr^3)
124 : !! xclevel= XC functional level
125 : !! xc_denpos= lowest allowed density (usually for the computation of the XC functionals)
126 : !! xc_taupos= lowest allowed kinetic energy density (for mGGA XC functionals)
127 : !! xred=atom positions in reduced coordinates
128 : !! znucl(ntypat)=gives the nuclear charge for all types of atoms
129 : !!
130 : !! OUTPUT
131 : !! paw_ij(my_natom)%dijhartree(qphase*lmn2_size)=Hartree contribution to dij;
132 : !! Enters into calculation of hartree energy
133 : !! ==== if option=0 or 2
134 : !! paw_energies <type(pawang_type)>=several contributions to on-site PAW energies
135 : !! %epaw= total on-site PAW energy (direct scheme)
136 : !! %epaw_dc= total on-site PAW energy (double counting scheme)
137 : !! %epaw_core= core contribution to PAW energy (direct scheme)
138 : !! %epaw_core_dc= core contribution to PAW energy (double counting scheme)
139 : !! %epaw_xc= exchange-correlation on-site contribution to PAW energy
140 : !! %entropy_paw= on-site PAW contribution to total entropy
141 : !! ==== if option=0 or 2 and ipert<=0
142 : !! compch_sph=compensation charge integral inside spheres computed over spherical meshes
143 : !! ==== if (option=0 or 1) and paw_an(:)%has_vxc=1
144 : !! paw_an(my_natom)%vxc1(cplex*mesh_size,:,nspden)=XC potential calculated from "on-site" density
145 : !! paw_an(my_natom)%vxct1(cplex*mesh_size,:,nspden)=XC potential calculated from "on-site" pseudo density
146 : !! ==== if (option=0 or 1) and paw_an(:)%has_vxctau=1
147 : !! paw_an(my_natom)%vxctau1(cplex*mesh_size,:,nspden)=1st deriv. of XC energy wrt to kinetic energy density (all electron)
148 : !! paw_an(my_natom)%vxcttau1(cplex*mesh_size,:,nspden)=1st deriv. of XC energy wrt to kinetic energy density (pseudo)
149 : !! ==== if paw_an(iatom_tot)%has_vxcval==1 compute also XC potentials neglecting core charge
150 : !! paw_an(my_natom)%vxc1_val(cplex*mesh_size,:nspden)=XC potential calculated from spherical valence density
151 : !! paw_an(my_natom)%vxct1_val(cplex*mesh_size,:nspden)=XC potential calculated from spherical valence pseudo density
152 : !! ==== if nzlmopt==-1,
153 : !! paw_an(iatom_tot)%lnmselect(lm_size,nspden)=select the non-zero LM-moments of rho1 and trho1
154 : !! ==== if paw_an(:)%has_vhartree=1
155 : !! paw_an(my_natom)%vh1(cplex*mesh_size,1,1)=Hartree total potential calculated from "on-site" density
156 : !! ==== if pawspnorb>0
157 : !! paw_ij(my_natom)%dijso(qphase*cplex_dij*lmn2_size,nspden)=spin-orbit contribution to dij
158 : !!
159 : !! NOTES
160 : !! Response function calculations:
161 : !! In order to compute first- or second-order quantities, paw_an (resp. paw_ij) datastructures
162 : !! must contain first-order quantities, namely paw_an1 (resp. paw_ij1).
163 : !!
164 : !! SOURCE
165 :
166 49030 : subroutine pawdenpot(compch_sph,el_temp,gprimd,ipert,ixc,my_natom,natom,nspden,ntypat,nucdipmom,&
167 49030 : & nzlmopt,option,paw_an,paw_an0,paw_energies,paw_ij,pawang,pawprtvol,pawrad,pawrhoij,&
168 49030 : & pawspnorb,pawtab,pawxcdev,spnorbscl,xclevel,xc_denpos,xc_taupos,xred,ucvol,znucl,spinaxis,&
169 24515 : & electronpositron,mpi_atmtab,comm_atom,vpotzero,hyb_mixing,hyb_mixing_sr,rcpaw,extfpmd) ! optional arguments
170 :
171 : !Arguments ---------------------------------------------
172 : !scalars
173 : integer,intent(in) :: ipert,ixc,my_natom,natom,nspden,ntypat,nzlmopt,option,pawprtvol
174 : integer,intent(in) :: pawspnorb,pawxcdev,xclevel
175 : integer,optional,intent(in) :: comm_atom
176 : real(dp),intent(in) :: spnorbscl,xc_denpos,xc_taupos,ucvol,el_temp
177 : real(dp),intent(in),optional :: hyb_mixing,hyb_mixing_sr
178 : real(dp),intent(out) :: compch_sph
179 : type(electronpositron_type),pointer,optional :: electronpositron
180 : type(pawang_type),intent(in) :: pawang
181 : type(paw_energies_type),intent(out) :: paw_energies
182 : type(rcpaw_type),pointer,intent(inout),optional :: rcpaw
183 : type(extfpmd_type),pointer,intent(in),optional :: extfpmd
184 : !arrays
185 : integer,optional,target,intent(in) :: mpi_atmtab(:)
186 : real(dp),intent(in) :: gprimd(3,3),nucdipmom(3,natom),xred(3,natom),znucl(ntypat),spinaxis(3)
187 : real(dp),intent(out),optional :: vpotzero(2)
188 : type(paw_an_type),intent(inout) :: paw_an(my_natom)
189 : type(paw_an_type), intent(in) :: paw_an0(my_natom)
190 : type(paw_ij_type),intent(inout) :: paw_ij(my_natom)
191 : type(pawrad_type),intent(in) :: pawrad(ntypat)
192 : type(pawrhoij_type),intent(inout) :: pawrhoij(my_natom)
193 : type(pawtab_type),intent(in) :: pawtab(ntypat)
194 :
195 : !Local variables ---------------------------------------
196 : !scalars
197 : integer, parameter :: PAWU_ALGO_1=1,PAWU_ALGO_2=2
198 : integer, parameter :: PAWU_FLL=1,PAWU_AMF=2
199 : integer :: add_core_energy,cplex,cplex_dij,cplex_rhoij,has_kxc,has_k3xc,has_vxctau
200 : integer :: iatom,iatom_tot,idum,ierr,ii,ipositron,iq,iq0_dij,iq0_rhoij
201 : integer :: itypat,itypat0,lm_size,lmn2_size,mesh_size
202 : integer :: my_comm_atom,ndij,nkxc1,nk3xc1,nsppol,opt_compch,pawu_algo,pawu_dblec
203 : integer :: ilmn,ilm,iln,j0lmn,jlm,jlmn,jln,klmn,ispden
204 : integer :: qphase,usecore,usekden,usetcore,usepawu,usexcnhat,usenhat,usefock
205 : logical :: keep_vhartree,my_atmtab_allocated,need_kxc,need_k3xc,need_vxctau,extfpmd_pawsph
206 : logical :: non_magnetic_xc,paral_atom,temp_vxc,eijkl_is_sym,rcpaw_has_valdens,usercpaw
207 : real(dp) :: e1t10,e1xc,e1xcdc,efock,efockdc,eexc,ssxc,eexcdc,eexdctemp
208 : real(dp) :: eexc_val,ssxc_val,eexcdc_val,eexex,eexexdc,eextemp,ssxtemp,eh2
209 : real(dp) :: edftumdc,edftumdcdc,edftufll,ehnzc,ekincore,enucdip,etmp,espnorb,etild1xc,etild1xcdc
210 : real(dp) :: s1xc,stild1xc,sxccore,extfpmd_rho
211 : real(dp) :: exccore,exchmix,hyb_mixing_,hyb_mixing_sr_,rdum
212 : real(dp) :: intvh,intg,eshift,eh2dc,ehpw
213 : character(len=3) :: pertstrg
214 : character(len=500) :: msg
215 : !arrays
216 : integer :: idum1(0),idum3(0,0,0)
217 24515 : integer,pointer :: my_atmtab(:)
218 24515 : logical,allocatable :: lmselect_cur(:),lmselect_cur_ep(:),lmselect_ep(:),lmselect_tmp(:)
219 : real(dp) :: tsec(2)
220 24515 : real(dp),allocatable :: dij_ep(:),dijfock_vv(:,:),dijfock_cv(:,:),ff(:)
221 24515 : real(dp),allocatable :: one_over_rad2(:),kxc_tmp(:,:,:),k3xc_tmp(:,:,:)
222 24515 : real(dp),allocatable :: mpiarr(:),nhat1(:,:,:),nhat1_ep(:,:,:)
223 : real(dp) :: rdum2(0,0),rdum3(0,0,0),rdum3a(0,0,0),rdum4(0,0,0,0)
224 24515 : real(dp),allocatable :: rho(:),rho1(:,:,:),rho1_ep(:,:,:),rho1xx(:,:,:)
225 24515 : real(dp),allocatable :: tau1(:,:,:),ttau1(:,:,:), trho1(:,:,:),trho1_ep(:,:,:)
226 24515 : real(dp),allocatable :: vh(:),vxc_tmp(:,:,:),vxctau_tmp(:,:,:)
227 :
228 : ! *************************************************************************
229 :
230 : DBG_ENTER("COLL")
231 :
232 24515 : call timab(560,1,tsec)
233 :
234 : !Various inits
235 24515 : hyb_mixing_ =zero ; if(present(hyb_mixing)) hyb_mixing_ =hyb_mixing
236 24515 : hyb_mixing_sr_=zero ; if(present(hyb_mixing_sr)) hyb_mixing_sr_=hyb_mixing_sr
237 24515 : usercpaw=.false. ; if (present(rcpaw)) usercpaw=associated(rcpaw)
238 24515 : usefock=0;if (abs(hyb_mixing_)>tol8.or.abs(hyb_mixing_sr_)>tol8) usefock=1
239 63110 : usexcnhat=maxval(pawtab(1:ntypat)%usexcnhat)
240 63110 : add_core_energy=minval(pawtab(1:ntypat)%add_core_energy)
241 24515 : usekden=pawxc_get_usekden(ixc)
242 24515 : usenhat = usexcnhat
243 77044 : keep_vhartree=(maxval(paw_an(:)%has_vhartree)>0)
244 24515 : if (keep_vhartree) usenhat = 1
245 24515 : compch_sph=-1.d5
246 24515 : opt_compch=0;if (option/=1.and.ipert<=0) opt_compch=1
247 24515 : if (opt_compch==1) compch_sph=zero
248 24515 : nsppol=1;if (my_natom>0) nsppol=pawrhoij(1)%nsppol
249 24515 : pertstrg=" ";if (ipert>0) pertstrg="(1)"
250 :
251 : !Various checks
252 24515 : if(nzlmopt/=0.and.nzlmopt/=1.and.nzlmopt/=-1) then
253 0 : msg='invalid value for variable "nzlmopt"!'
254 0 : ABI_BUG(msg)
255 : end if
256 24515 : if (my_natom>0) then
257 22731 : if(paw_ij(1)%has_dijhartree==0.and.ipert/=natom+1.and.ipert/=natom+10) then
258 0 : msg='dijhartree must be allocated!'
259 0 : ABI_BUG(msg)
260 : end if
261 22731 : if(paw_ij(1)%has_dijU==0.and.pawtab(1)%usepawu/=0.and. &
262 : & ((ipert>0.and.ipert/=natom+1.and.ipert/=natom+10).or.pawtab(1)%usepawu<0)) then
263 0 : msg='dijU must be allocated!'
264 0 : ABI_BUG(msg)
265 : end if
266 22731 : if (pawrhoij(1)%qphase<paw_an(1)%cplex) then
267 0 : msg='pawrhoij()%qphase must be >=paw_an()%cplex!'
268 0 : ABI_BUG(msg)
269 : end if
270 22731 : if (ipert>0.and.(ipert<=natom.or.ipert==natom+2).and.paw_an0(1)%has_kxc/=2) then
271 0 : msg='XC kernels for ground state must be in memory!'
272 0 : ABI_BUG(msg)
273 : end if
274 22731 : if(paw_an(1)%has_vxc==0.and.(option==0.or.option==1).and. &
275 : & .not.(ipert==natom+1.or.ipert==natom+10)) then
276 0 : msg='vxc1 and vxct1 must be allocated!'
277 0 : ABI_BUG(msg)
278 : end if
279 22731 : if(paw_an(1)%has_vxctau==0.and.(option==0.or.option==1).and.usekden==1) then
280 0 : msg='vxctau1 and vxcttau1 must be allocated!'
281 0 : ABI_BUG(msg)
282 : end if
283 22731 : if (ipert>0.and.paw_an(1)%has_vxctau==1.and.usekden==1) then
284 84 : if (ipert .NE. natom+1) then
285 0 : msg='computation of vxctau not compatible with RF (ipert>0)!'
286 0 : ABI_BUG(msg)
287 : end if
288 : end if
289 22731 : if (ipert>0.and.paw_an(1)%has_vhartree==1) then
290 0 : msg='computation of vhartree not compatible with RF (ipert>0)!'
291 0 : ABI_BUG(msg)
292 : end if
293 22731 : if (ipert>0.and.paw_an(1)%has_vxcval==1.and.(option==0.or.option==1)) then
294 0 : msg='computation of vxc_val not compatible with RF (ipert>0)!'
295 0 : ABI_BUG(msg)
296 : end if
297 : end if
298 :
299 24515 : ipositron=0
300 24515 : if (present(electronpositron)) then
301 15988 : ipositron=electronpositron_calctype(electronpositron)
302 15988 : if (ipositron==1.and.pawtab(1)%has_kij/=2) then
303 0 : msg='kij must be in memory for electronpositron%calctype=1!'
304 0 : ABI_BUG(msg)
305 : end if
306 15988 : if (ipert>0) then
307 0 : msg='electron-positron calculation not available for ipert>0!'
308 0 : ABI_ERROR(msg)
309 : end if
310 : end if
311 :
312 : !Set up parallelism over atoms
313 24515 : paral_atom=(present(comm_atom).and.(my_natom/=natom))
314 24515 : nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab
315 24515 : my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom
316 : call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,&
317 24515 : & my_natom_ref=my_natom)
318 :
319 : !Init energies
320 24515 : if (option/=1) then
321 20452 : call paw_energies_setzero(paw_energies)
322 20452 : e1xc=zero ; e1xcdc=zero
323 20452 : etild1xc=zero ; etild1xcdc=zero
324 20452 : s1xc=zero ; stild1xc=zero ; sxccore=zero
325 20452 : exccore=zero ; eh2=zero ; e1t10=zero
326 20452 : eh2dc=zero ; ehpw = zero
327 20452 : edftumdc=zero ; edftumdcdc=zero ; edftufll=zero
328 20452 : eexex=zero ; eexexdc=zero
329 20452 : eextemp=zero ; eexdctemp=zero
330 20452 : espnorb=zero ; enucdip=zero
331 20452 : efock=zero ; efockdc=zero
332 20452 : ekincore=zero ; ehnzc=zero
333 20452 : if (ipositron/=0) then
334 488 : electronpositron%e_paw =zero
335 488 : electronpositron%e_pawdc=zero
336 : end if
337 : end if
338 :
339 : !For some perturbations, nothing to do
340 24515 : if (ipert==natom+1.or.ipert==natom+10) return
341 :
342 : !vpotzero is needed for both the energy and the potential
343 22877 : if (present(vpotzero)) vpotzero(:)=zero
344 :
345 : !Select PAW+U algo, different for DFT and DFPT
346 58714 : usepawu=maxval(pawtab(1:ntypat)%usepawu)
347 58714 : ii=minval(pawtab(1:ntypat)%usepawu);if (ii<0) usepawu=ii
348 22877 : non_magnetic_xc=(mod(abs(usepawu),10)==4)
349 :
350 : !if PAW+U, compute noccmmp^{\sigma}_{m,m'} occupation matrix
351 22877 : if (usepawu/=0.and.ipert<=0.and.ipositron/=1) then
352 2283 : if (paral_atom) then
353 : call setnoccmmp(1,0,rdum4,0,0,idum3,my_natom,natom,0,1,nsppol,0,ntypat,&
354 : & paw_ij,pawang,pawprtvol,pawrhoij,pawtab,rdum2,idum1,idum1,0,usepawu,&
355 160 : & comm_atom=my_comm_atom,mpi_atmtab=mpi_atmtab)
356 : else
357 : call setnoccmmp(1,0,rdum4,0,0,idum3,my_natom,natom,0,1,nsppol,0,ntypat,&
358 2123 : & paw_ij,pawang,pawprtvol,pawrhoij,pawtab,rdum2,idum1,idum1,0,usepawu)
359 : end if
360 : end if
361 :
362 : !Print some titles
363 22877 : if (abs(pawprtvol)>=2) then
364 1810 : if (nzlmopt<1) write(msg, '(6a)') ch10,' PAW TEST:',ch10,&
365 392 : ' ====== Moments of (n1-tn1)',trim(pertstrg),' ========='
366 3032 : if (nzlmopt==1) write(msg, '(6a)') ch10,' PAW TEST:',ch10,&
367 2836 : ' ==== Non-zero Moments of (n1-tn1)',trim(pertstrg),' ===='
368 1614 : call wrtout(std_out,msg,'COLL')
369 1614 : if (usexcnhat/=0) then
370 1311 : write(msg, '(6a)')' The moments of (n1-tn1-nhat1)',trim(pertstrg),' must be very small...'
371 1311 : call wrtout(std_out,msg,'COLL')
372 : end if
373 : end if
374 :
375 22877 : extfpmd_rho=zero
376 22877 : extfpmd_pawsph=.false.
377 22877 : if(present(extfpmd)) then
378 15129 : if(associated(extfpmd)) then
379 90 : extfpmd%eshift_paw=zero
380 90 : if(extfpmd%pawsph>=1) then
381 7 : extfpmd_pawsph=.true.
382 7 : extfpmd_rho=extfpmd%nelect/ucvol
383 7 : usenhat=1
384 : endif
385 : endif
386 : endif
387 :
388 :
389 : !================ Big loop on atoms =======================
390 : !==========================================================
391 :
392 72298 : do iatom=1,my_natom
393 49421 : iatom_tot=iatom;if (paral_atom) iatom_tot=my_atmtab(iatom)
394 49421 : itypat=pawrhoij(iatom)%itypat
395 49421 : exchmix=pawtab(itypat)%exchmix
396 49421 : lmn2_size=paw_ij(iatom)%lmn2_size
397 49421 : lm_size=paw_an(iatom)%lm_size
398 49421 : mesh_size=pawtab(itypat)%mesh_size
399 49421 : usecore=1;usetcore =pawtab(itypat)%usetcore
400 49421 : if (ipert/=0) usecore=0 ! This is true for phonons and Efield pert.
401 49421 : if (ipert/=0) usetcore=0 ! This is true for phonons and Efield pert.
402 49421 : has_kxc =paw_an(iatom)%has_kxc ;need_kxc =(has_kxc ==1)
403 49421 : has_k3xc=paw_an(iatom)%has_k3xc;need_k3xc=(has_k3xc==1)
404 49421 : has_vxctau=paw_an(iatom)%has_vxctau ;need_vxctau =(has_vxctau>=1.and.usekden==1)
405 49421 : cplex=paw_an(iatom)%cplex
406 49421 : cplex_dij=paw_ij(iatom)%cplex_dij
407 49421 : cplex_rhoij=pawrhoij(iatom)%cplex_rhoij
408 49421 : qphase=pawrhoij(iatom)%qphase
409 49421 : ndij=paw_ij(iatom)%ndij
410 49421 : iq0_rhoij=cplex_rhoij*lmn2_size
411 49421 : iq0_dij=cplex_dij*lmn2_size
412 49421 : usepawu=pawtab(itypat)%usepawu
413 49421 : pawu_algo=merge(PAWU_ALGO_1,PAWU_ALGO_2,ipert<=0.and.usepawu>=0)
414 49421 : pawu_dblec=merge(PAWU_FLL,PAWU_AMF,abs(usepawu)==1.or.abs(usepawu)==4)
415 :
416 : ! Allocations of "on-site" densities
417 247105 : ABI_MALLOC(rho1 ,(cplex*mesh_size,lm_size,nspden))
418 197684 : ABI_MALLOC(trho1,(cplex*mesh_size,lm_size,nspden))
419 247105 : ABI_MALLOC(nhat1,(cplex*mesh_size,lm_size,nspden*usenhat))
420 1671198938 : rho1(:,:,:)=zero;trho1(:,:,:)=zero;nhat1(:,:,:)=zero
421 49421 : if (usekden==1) then
422 1936 : ABI_MALLOC(tau1 ,(cplex*mesh_size,lm_size,nspden))
423 1936 : ABI_MALLOC(ttau1,(cplex*mesh_size,lm_size,nspden))
424 : end if
425 49421 : if (ipositron/=0) then ! Additional allocation for the electron-positron case
426 23504 : ABI_MALLOC(rho1_ep ,(cplex*mesh_size,lm_size,nspden))
427 23504 : ABI_MALLOC(trho1_ep,(cplex*mesh_size,lm_size,nspden))
428 23504 : ABI_MALLOC(nhat1_ep,(cplex*mesh_size,lm_size,nspden*usenhat))
429 : end if
430 148263 : ABI_MALLOC(lmselect_cur,(lm_size))
431 566723 : lmselect_cur(:)=.true.
432 438677 : if (nzlmopt==1) lmselect_cur(:)=paw_an(iatom)%lmselect(:)
433 :
434 : ! Store some usefull quantities
435 49421 : itypat0=0;if (iatom>1) itypat0=pawrhoij(iatom-1)%itypat
436 49421 : if (itypat/=itypat0) then
437 86691 : ABI_MALLOC(one_over_rad2,(mesh_size))
438 28897 : one_over_rad2(1)=zero
439 32467877 : one_over_rad2(2:mesh_size)=one/pawrad(itypat)%rad(2:mesh_size)**2
440 : end if
441 :
442 : ! Need to allocate vxc1 in particular cases
443 : if (pawspnorb>0.and.ipert==0.and.option==2.and.ipositron/=1.and. &
444 49421 : & cplex_rhoij==2.and.paw_an(iatom)%has_vxc==0) then
445 : ! These should already be allocated in paw_an_init!
446 0 : if (allocated(paw_an(iatom)%vxc1)) then
447 0 : ABI_FREE(paw_an(iatom)%vxc1)
448 : end if
449 0 : if (pawxcdev==0)then
450 0 : ABI_MALLOC(paw_an(iatom)%vxc1,(cplex*mesh_size,paw_an(iatom)%angl_size,nspden))
451 : else
452 0 : ABI_MALLOC(paw_an(iatom)%vxc1,(cplex*mesh_size,lm_size,nspden))
453 : end if
454 0 : paw_an(iatom)%has_vxc=1
455 0 : temp_vxc=.true.
456 : else
457 : temp_vxc=.false.
458 : end if
459 :
460 : ! ===== Compute "on-site" densities (n1, ntild1, nhat1) =====
461 : ! ==========================================================
462 :
463 49421 : rcpaw_has_valdens=.false.
464 49421 : if(usercpaw) then
465 35 : rcpaw_has_valdens=rcpaw%val(iatom)%has_dens
466 : endif
467 35 : if(rcpaw_has_valdens) then
468 1870548 : rho1=rcpaw%val(iatom)%rho1
469 1870548 : trho1=rcpaw%val(iatom)%trho1
470 827958 : if(usenhat==1) nhat1=rcpaw%val(iatom)%nhat1
471 30 : if(opt_compch==1) compch_sph=compch_sph+rcpaw%val(iatom)%compch_sph
472 : else
473 : call pawdensities(compch_sph,cplex,iatom_tot,lmselect_cur,paw_an(iatom)%lmselect,lm_size,&
474 : & nhat1,nspden,nzlmopt,opt_compch,1-usenhat,-1,1,pawang,pawprtvol,pawrad(itypat),&
475 49391 : & pawrhoij(iatom),pawtab(itypat),rho1,trho1,extfpmd_rho=extfpmd_rho/nspden,one_over_rad2=one_over_rad2)
476 : endif
477 49421 : if (usekden==1) then
478 1452 : ABI_MALLOC(lmselect_tmp,(lm_size))
479 9960 : lmselect_tmp(:)=.true.
480 : call pawkindensities(cplex,lmselect_tmp,lm_size,nspden,-1,1,-1,&
481 : & pawang,pawrad(itypat),pawrhoij(iatom),pawtab(itypat),tau1,ttau1,&
482 484 : & one_over_rad2=one_over_rad2)
483 484 : ABI_FREE(lmselect_tmp)
484 : end if
485 :
486 49421 : if (ipositron/=0) then
487 : ! Electron-positron calculation: need additional on-site densities:
488 : ! if ipositron==1, need electronic on-site densities
489 : ! if ipositron==2, need positronic on-site densities
490 17628 : ABI_MALLOC(lmselect_ep,(lm_size))
491 11752 : ABI_MALLOC(lmselect_cur_ep,(lm_size))
492 58760 : lmselect_cur_ep(:)=.true.
493 41408 : if (nzlmopt==1) lmselect_cur_ep(:)=electronpositron%lmselect_ep(1:lm_size,iatom)
494 : call pawdensities(rdum,cplex,iatom_tot,lmselect_cur_ep,lmselect_ep,&
495 : & lm_size,nhat1_ep,nspden,nzlmopt,0,1-usenhat,-1,0,pawang,0,pawrad(itypat),&
496 : & electronpositron%pawrhoij_ep(iatom),pawtab(itypat),&
497 5876 : & rho1_ep,trho1_ep,one_over_rad2=one_over_rad2)
498 23228 : if (nzlmopt<1) electronpositron%lmselect_ep(1:lm_size,iatom)=lmselect_ep(1:lm_size)
499 5876 : ABI_FREE(lmselect_ep)
500 5876 : ABI_FREE(lmselect_cur_ep)
501 : end if
502 :
503 : ! =========== Compute XC potentials and energies ===========
504 : ! ==========================================================
505 :
506 : ! Temporary storage
507 49421 : nkxc1 =0;if (paw_an(iatom)%has_kxc /=0) nkxc1 =paw_an(iatom)%nkxc1
508 49421 : nk3xc1=0;if (paw_an(iatom)%has_k3xc/=0.and.pawxcdev==0) nk3xc1=paw_an(iatom)%nk3xc1
509 49421 : if (pawxcdev/=0) then
510 226835 : ABI_MALLOC(vxc_tmp,(cplex*mesh_size,lm_size,nspden))
511 45367 : if (need_kxc) then
512 315 : ABI_MALLOC(kxc_tmp,(mesh_size,lm_size,nkxc1))
513 : end if
514 45367 : if (need_k3xc) then
515 0 : msg = 'Computation of k3xc with pawxcdev/=0 is not implemented yet!'
516 0 : ABI_BUG(msg)
517 : end if
518 : end if
519 49421 : if (pawxcdev==0) then
520 20270 : ABI_MALLOC(vxc_tmp,(cplex*mesh_size,pawang%angl_size,nspden))
521 1354292671 : vxc_tmp(:,:,:)=zero
522 4054 : if (need_kxc) then
523 160 : ABI_MALLOC(kxc_tmp,(mesh_size,pawang%angl_size,nkxc1))
524 : end if
525 4054 : if (need_k3xc) then
526 40 : ABI_MALLOC(k3xc_tmp,(mesh_size,pawang%angl_size,nk3xc1))
527 : end if
528 4054 : if (need_vxctau) then
529 1936 : ABI_MALLOC(vxctau_tmp,(cplex*mesh_size,pawang%angl_size,nspden))
530 108267776 : vxctau_tmp(:,:,:)=zero
531 : end if
532 : end if
533 49421 : idum=0
534 49421 : if (.not.allocated(vxc_tmp)) then
535 0 : ABI_MALLOC(vxc_tmp,(0,0,0))
536 : end if
537 49421 : if (.not.allocated(kxc_tmp)) then
538 49326 : ABI_MALLOC(kxc_tmp,(0,0,0))
539 : end if
540 49421 : if (.not.allocated(k3xc_tmp)) then
541 49413 : ABI_MALLOC(k3xc_tmp,(0,0,0))
542 : end if
543 49421 : if (.not.allocated(vxctau_tmp)) then
544 48937 : ABI_MALLOC(vxctau_tmp,(0,0,0))
545 : end if
546 :
547 : ! ===== Vxc1 term =====
548 49421 : if (ipositron/=1) then
549 47030 : if (pawxcdev/=0) then
550 42976 : if (ipert==0) then
551 : call pawxcm(pawtab(itypat)%coredens,eexc,eexcdc,ssxc,idum,hyb_mixing_,ixc,kxc_tmp,lm_size,&
552 : & paw_an(iatom)%lmselect,nhat1,nkxc1,non_magnetic_xc,mesh_size,nspden,option,&
553 34662 : & pawang,pawrad(itypat),pawxcdev,rho1,usecore,0,vxc_tmp,xclevel,xc_denpos,el_temp)
554 : else
555 : call pawxcm_dfpt(pawtab(itypat)%coredens,cplex,cplex,eexc,ixc,paw_an0(iatom)%kxc1,lm_size,&
556 : & paw_an(iatom)%lmselect,nhat1,paw_an0(iatom)%nkxc1,non_magnetic_xc,mesh_size,nspden,option,&
557 8314 : & pawang,pawrad(itypat),rho1,usecore,0,vxc_tmp,xclevel)
558 8314 : eexcdc=zero
559 8314 : ssxc=zero ! Forcing entropy to zero when using dfpt for now
560 : end if
561 : else
562 4054 : if (ipert==0) then
563 : call pawxc(pawtab(itypat)%coredens,eexc,eexcdc,ssxc,hyb_mixing_,ixc,kxc_tmp,k3xc_tmp,lm_size,&
564 : & paw_an(iatom)%lmselect,nhat1,nkxc1,nk3xc1,non_magnetic_xc,mesh_size,nspden,option,&
565 : & pawang,pawrad(itypat),rho1,usecore,0,vxc_tmp,xclevel,xc_denpos,el_temp,&
566 974 : & coretau=pawtab(itypat)%coretau,taur=tau1,vxctau=vxctau_tmp,xc_taupos=xc_taupos)
567 : else
568 : call pawxc_dfpt(pawtab(itypat)%coredens,cplex,cplex,eexc,ixc,paw_an0(iatom)%kxc1,lm_size,&
569 : & paw_an(iatom)%lmselect,nhat1,paw_an0(iatom)%nkxc1,non_magnetic_xc,mesh_size,nspden,option,&
570 3080 : & pawang,pawrad(itypat),rho1,usecore,0,paw_an0(iatom)%vxc1,vxc_tmp,xclevel)
571 3080 : eexcdc=zero
572 3080 : ssxc=zero ! Forcing entropy to zero when using dfpt for now
573 : end if
574 : end if
575 47030 : if (option/=1) then
576 39658 : e1xc=e1xc+eexc
577 39658 : e1xcdc=e1xcdc+eexcdc
578 39658 : s1xc=s1xc+ssxc
579 : end if
580 1413778647 : if (option<2.or.temp_vxc) paw_an(iatom)%vxc1(:,:,:)=vxc_tmp(:,:,:)
581 96162490 : if (option<2.and.need_vxctau) paw_an(iatom)%vxctau1(:,:,:)=vxctau_tmp(:,:,:)
582 22938684 : if (need_kxc .and.nkxc1>0 ) paw_an(iatom)%kxc1(:,:,:) =kxc_tmp(:,:,:)
583 5039992 : if (need_k3xc.and.nk3xc1>0) paw_an(iatom)%k3xc1(:,:,:)=k3xc_tmp(:,:,:)
584 47030 : if(extfpmd_pawsph) then
585 42 : ABI_MALLOC(ff,(mesh_size))
586 19712 : ff=zero
587 28 : do ispden=1,nspden
588 19712 : ff(1:mesh_size)=paw_an(iatom)%vxc1(1:mesh_size,1,ispden)*sqrt(four_pi)*pawrad(itypat)%rad(1:mesh_size)**2
589 14 : call simp_gen(eshift,ff,pawrad(itypat))
590 28 : extfpmd%eshift_paw=extfpmd%eshift_paw+eshift/ucvol/nspden
591 : enddo
592 14 : ABI_FREE(ff)
593 : endif
594 : else ! ipositron==1
595 28175007 : if (option<2.or.temp_vxc) paw_an(iatom)%vxc1(:,:,:)=zero
596 2391 : if (need_kxc.and.nkxc1>0) paw_an(iatom)%kxc1(:,:,:)=zero
597 : end if
598 :
599 : ! Additional electron-positron XC term (if ipositron/=0)
600 49421 : if (ipositron/=0) then
601 5876 : if (pawxcdev/=0) then
602 : call pawxcmpositron(ipositron,pawtab(itypat)%coredens,eexc,eexcdc,electronpositron%ixcpositron,&
603 : & lm_size,paw_an(iatom)%lmselect,electronpositron%lmselect_ep(1:lm_size,iatom),&
604 : & nhat1,nhat1_ep,mesh_size,nspden,option,pawang,pawrad(itypat),pawxcdev,&
605 5876 : & electronpositron%posdensity0_limit,rho1,rho1_ep,usecore,0,vxc_tmp,xc_denpos)
606 : else
607 : call pawxcpositron(ipositron,pawtab(itypat)%coredens,eexc,eexcdc,electronpositron%ixcpositron,&
608 : & lm_size,paw_an(iatom)%lmselect,electronpositron%lmselect_ep(1:lm_size,iatom),&
609 : & nhat1,nhat1_ep,mesh_size,nspden,option,pawang,pawrad(itypat),&
610 0 : & electronpositron%posdensity0_limit,rho1,rho1_ep,usecore,0,vxc_tmp,xc_denpos)
611 : end if
612 5876 : if (option/=1) then
613 5876 : electronpositron%e_paw =electronpositron%e_paw +eexc
614 5876 : electronpositron%e_pawdc=electronpositron%e_pawdc+eexcdc
615 : end if
616 70820332 : if (option<2.or.temp_vxc) paw_an(iatom)%vxc1(:,:,:)=paw_an(iatom)%vxc1(:,:,:)+vxc_tmp(:,:,:)
617 5876 : if (need_kxc.and.nkxc1>0) paw_an(iatom)%kxc1(:,:,:)=paw_an(iatom)%kxc1(:,:,:)+kxc_tmp(:,:,:)
618 : end if
619 :
620 : ! ===== tVxc1 term =====
621 49421 : if (ipositron/=1) then
622 47030 : if (pawxcdev/=0) then
623 42976 : if (ipert==0) then
624 : call pawxcm(pawtab(itypat)%tcoredens(:,1),&
625 : & eexc,eexcdc,ssxc,idum,hyb_mixing_,ixc,kxc_tmp,lm_size,&
626 : & paw_an(iatom)%lmselect,nhat1,nkxc1,non_magnetic_xc,mesh_size,nspden,option,&
627 34662 : & pawang,pawrad(itypat),pawxcdev,trho1,usetcore,2*usexcnhat,vxc_tmp,xclevel,xc_denpos,el_temp)
628 : else
629 : call pawxcm_dfpt(pawtab(itypat)%tcoredens(:,1),&
630 : & cplex,cplex,eexc,ixc,paw_an0(iatom)%kxct1,lm_size,&
631 : & paw_an(iatom)%lmselect,nhat1,paw_an0(iatom)%nkxc1,non_magnetic_xc,mesh_size,nspden,option,&
632 8314 : & pawang,pawrad(itypat),trho1,usetcore,2*usexcnhat,vxc_tmp,xclevel)
633 8314 : eexcdc=zero
634 8314 : ssxc=zero ! Forcing entropy to zero when using dfpt for now
635 : end if
636 : else
637 4054 : if (ipert==0) then
638 : call pawxc(pawtab(itypat)%tcoredens(:,1),&
639 : & eexc,eexcdc,ssxc,hyb_mixing_,ixc,kxc_tmp,k3xc_tmp,lm_size,&
640 : & paw_an(iatom)%lmselect,nhat1,nkxc1,nk3xc1,non_magnetic_xc,mesh_size,nspden,option,&
641 : & pawang,pawrad(itypat),trho1,usetcore,2*usexcnhat,vxc_tmp,xclevel,xc_denpos,el_temp,&
642 974 : & coretau=pawtab(itypat)%tcoretau,taur=ttau1,vxctau=vxctau_tmp,xc_taupos=xc_taupos)
643 : else
644 : call pawxc_dfpt(pawtab(itypat)%tcoredens(:,1),&
645 : & cplex,cplex,eexc,ixc,paw_an0(iatom)%kxct1,lm_size,&
646 : & paw_an(iatom)%lmselect,nhat1,paw_an0(iatom)%nkxc1,non_magnetic_xc,mesh_size,nspden,option,&
647 3080 : & pawang,pawrad(itypat),trho1,usetcore,2*usexcnhat,paw_an0(iatom)%vxct1,vxc_tmp,xclevel)
648 3080 : eexcdc=zero
649 3080 : ssxc=zero ! Forcing entropy to zero when using dfpt for now
650 : end if
651 : end if
652 47030 : if (option/=1) then
653 39658 : etild1xc=etild1xc+eexc
654 39658 : etild1xcdc=etild1xcdc+eexcdc
655 39658 : stild1xc=stild1xc+ssxc
656 : end if
657 1413778647 : if (option<2) paw_an(iatom)%vxct1(:,:,:)=vxc_tmp(:,:,:)
658 96162490 : if (option<2.and.need_vxctau) paw_an(iatom)%vxcttau1(:,:,:)=vxctau_tmp(:,:,:)
659 22938684 : if (need_kxc.and. nkxc1>0 ) paw_an(iatom)%kxct1(:,:,:) =kxc_tmp(:,:,:)
660 5039992 : if (need_k3xc.and.nk3xc1>0) paw_an(iatom)%k3xct1(:,:,:)=k3xc_tmp(:,:,:)
661 47030 : if(extfpmd_pawsph) then
662 42 : ABI_MALLOC(ff,(mesh_size))
663 19712 : ff=zero
664 28 : do ispden=1,nspden
665 19712 : ff(1:mesh_size)=paw_an(iatom)%vxct1(1:mesh_size,1,ispden)*sqrt(four_pi)*pawrad(itypat)%rad(1:mesh_size)**2
666 14 : call simp_gen(eshift,ff,pawrad(itypat))
667 28 : extfpmd%eshift_paw=extfpmd%eshift_paw-eshift/ucvol/nspden
668 : enddo
669 14 : ABI_FREE(ff)
670 : endif
671 : else ! ipositron==1
672 28175007 : if (option<2) paw_an(iatom)%vxct1(:,:,:)=zero
673 2391 : if (need_kxc.and.nkxc1>0) paw_an(iatom)%kxct1(:,:,:)=zero
674 : end if
675 :
676 : ! Additional electron-positron XC term (if ipositron/=0)
677 49421 : if (ipositron/=0) then
678 5876 : if (pawxcdev/=0) then
679 : call pawxcmpositron(ipositron,pawtab(itypat)%tcoredens(:,1),&
680 : & eexc,eexcdc,electronpositron%ixcpositron,&
681 : & lm_size,paw_an(iatom)%lmselect,electronpositron%lmselect_ep(1:lm_size,iatom),&
682 : & nhat1,nhat1_ep,mesh_size,nspden,option,pawang,pawrad(itypat),pawxcdev,&
683 5876 : & electronpositron%posdensity0_limit,trho1,trho1_ep,usetcore,2*usexcnhat,vxc_tmp,xc_denpos)
684 : else
685 : call pawxcpositron(ipositron,pawtab(itypat)%tcoredens(:,1),&
686 : & eexc,eexcdc,electronpositron%ixcpositron,&
687 : & lm_size,paw_an(iatom)%lmselect,electronpositron%lmselect_ep(1:lm_size,iatom),&
688 : & nhat1,nhat1_ep,mesh_size,nspden,option,pawang,pawrad(itypat),&
689 0 : & electronpositron%posdensity0_limit,trho1,trho1_ep,usetcore,2*usexcnhat,vxc_tmp,xc_denpos)
690 : end if
691 5876 : if (option/=1) then
692 5876 : electronpositron%e_paw =electronpositron%e_paw -eexc
693 5876 : electronpositron%e_pawdc=electronpositron%e_pawdc-eexcdc
694 : end if
695 70820332 : if (option<2) paw_an(iatom)%vxct1(:,:,:)=paw_an(iatom)%vxct1(:,:,:)+vxc_tmp(:,:,:)
696 5876 : if (need_kxc.and.nkxc1>0) paw_an(iatom)%kxct1(:,:,:)=paw_an(iatom)%kxct1(:,:,:)+kxc_tmp(:,:,:)
697 : end if
698 :
699 : ! Update flags defining the state of vxc and kxc
700 49421 : if (option<2) paw_an(iatom)%has_vxc=2
701 49421 : if (option<2.and.need_vxctau) paw_an(iatom)%has_vxctau=2
702 49421 : if (need_kxc.and.nkxc1>0) paw_an(iatom)%has_kxc=2
703 :
704 : ! Update core XC contribution to energy
705 49421 : if (option/=1.and.ipositron/=1) then
706 39658 : exccore=exccore+pawtab(itypat)%exccore
707 39658 : sxccore=sxccore+pawtab(itypat)%sxccore
708 : end if
709 :
710 : ! =========== Compute valence-only XC potentials ===========
711 : ! ==========================================================
712 49421 : if (ipert==0.and.paw_an(iatom)%has_vxcval==1.and.(option==0.or.option==1)) then
713 31 : if (.not.allocated(paw_an(iatom)%vxc1_val).or..not.allocated(paw_an(iatom)%vxct1_val)) then
714 0 : msg=' vxc1_val and vxct1_val must be associated'
715 0 : ABI_BUG(msg)
716 : end if
717 : ! ===== Vxc1_val term, vxc[n1] =====
718 31 : if (pawxcdev/=0) then
719 31 : write(msg,'(4a,es16.6)')ch10,&
720 31 : & ' pawdenpot : Computing valence-only v_xc[n1] using moments ',ch10,&
721 394086 : & ' Min density rho1 = ',MINVAL(rho1)
722 31 : call wrtout(std_out,msg,'COLL')
723 : call pawxcm(pawtab(itypat)%coredens,eexc_val,eexcdc_val,ssxc_val,idum,hyb_mixing_,ixc,kxc_tmp,lm_size,&
724 : & paw_an(iatom)%lmselect,nhat1,nkxc1,non_magnetic_xc,mesh_size,nspden,option,&
725 31 : & pawang,pawrad(itypat),pawxcdev,rho1,0,0,vxc_tmp,xclevel,xc_denpos,el_temp)
726 : else
727 0 : write(msg,'(2a)')ch10,' pawdenpot : Computing valence-only v_xc[n1] using angular mesh '
728 0 : call wrtout(std_out,msg,'COLL')
729 :
730 : call pawxc(pawtab(itypat)%coredens,eexc_val,eexcdc_val,ssxc_val,hyb_mixing_,ixc,kxc_tmp,k3xc_tmp,lm_size,&
731 : & paw_an(iatom)%lmselect,nhat1,nkxc1,nk3xc1,non_magnetic_xc,mesh_size,nspden,option,&
732 0 : & pawang,pawrad(itypat),rho1,0,0,vxc_tmp,xclevel,xc_denpos,el_temp)
733 : end if
734 394055 : if (option<2) paw_an(iatom)%vxc1_val(:,:,:)=vxc_tmp(:,:,:)
735 :
736 : ! ===== tVxc1_val term =====
737 31 : if (pawxcdev/=0) then
738 31 : if (usexcnhat/=0) then
739 31 : write(msg,'(4a,e16.6,2a,es16.6)')ch10,&
740 31 : & ' pawdenpot : Computing valence-only v_xc[tn1+nhat] using moments ',ch10,&
741 394055 : & ' Min density trho1 = ',MINVAL(trho1),ch10,&
742 394086 : & ' Min density trho1 + nhat = ',MINVAL(trho1+nhat1)
743 : else
744 0 : write(msg,'(4a,e16.6)')ch10,&
745 0 : & ' pawdenpot : Computing valence-only v_xc[tn1] using moments ',ch10,&
746 0 : & ' Min density trho1 = ',MINVAL(trho1)
747 : end if
748 31 : call wrtout(std_out,msg,'COLL')
749 : call pawxcm(pawtab(itypat)%tcoredens(:,1),&
750 : & eexc_val,eexcdc_val,ssxc_val,idum,hyb_mixing_,ixc,kxc_tmp,lm_size,&
751 : & paw_an(iatom)%lmselect,nhat1,nkxc1,non_magnetic_xc,mesh_size,nspden,option,&
752 31 : & pawang,pawrad(itypat),pawxcdev,trho1,0,2*usexcnhat,vxc_tmp,xclevel,xc_denpos,el_temp)
753 : else
754 0 : write(msg,'(2a)')ch10,' pawdenpot : Computing valence-only v_xc[tn1+nhat] using angular mesh'
755 0 : call wrtout(std_out,msg,'COLL')
756 : call pawxc(pawtab(itypat)%tcoredens(:,1),&
757 : & eexc_val,eexcdc_val,ssxc_val,hyb_mixing_,ixc,kxc_tmp,k3xc_tmp,lm_size,&
758 : & paw_an(iatom)%lmselect,nhat1,nkxc1,nk3xc1,non_magnetic_xc,mesh_size,nspden,option,&
759 0 : & pawang,pawrad(itypat),trho1,0,2*usexcnhat,vxc_tmp,xclevel,xc_denpos,el_temp)
760 : end if
761 : if (option<2) then
762 394055 : paw_an(iatom)%vxct1_val(:,:,:)=vxc_tmp(:,:,:)
763 31 : paw_an(iatom)%has_vxcval=2
764 : end if
765 : end if ! valence-only XC potentials
766 :
767 49421 : ABI_FREE(vxc_tmp)
768 49421 : ABI_FREE(vxctau_tmp)
769 49421 : ABI_FREE(kxc_tmp)
770 49421 : ABI_FREE(k3xc_tmp)
771 :
772 : ! ===== Compute first part of local exact-exchange energy term =====
773 : ! ===== Also compute corresponding potential =====
774 : ! ==================================================================
775 :
776 49421 : if (pawtab(itypat)%useexexch/=0.and.ipert==0.and.ipositron/=1) then
777 :
778 : ! ===== Re-compute a partial "on-site" density n1 (only l=lexexch contrib.)
779 210 : ABI_MALLOC(rho1xx,(mesh_size,lm_size,nspden))
780 126 : ABI_MALLOC(lmselect_tmp,(lm_size))
781 1092 : lmselect_tmp(:)=lmselect_cur(:)
782 : call pawdensities(rdum,cplex,iatom_tot,lmselect_cur,lmselect_tmp,lm_size,rdum3,nspden,&
783 : & 1,0,2,pawtab(itypat)%lexexch,0,pawang,pawprtvol,pawrad(itypat),&
784 42 : & pawrhoij(iatom),pawtab(itypat),rho1xx,rdum3a,one_over_rad2=one_over_rad2)
785 42 : ABI_FREE(lmselect_tmp)
786 : ! ===== Re-compute Exc1 and Vxc1; for local exact-exchange, this is done in GGA only
787 210 : ABI_MALLOC(vxc_tmp,(mesh_size,lm_size,nspden))
788 210 : ABI_MALLOC(kxc_tmp,(mesh_size,lm_size,nkxc1))
789 : call pawxcm(pawtab(itypat)%coredens,eextemp,eexdctemp,ssxtemp,pawtab(itypat)%useexexch,hyb_mixing_,ixc,kxc_tmp,lm_size,&
790 : & paw_an(iatom)%lmselect,nhat1,nkxc1,non_magnetic_xc,mesh_size,nspden,option,pawang,pawrad(itypat),pawxcdev,&
791 42 : & rho1xx,0,0,vxc_tmp,xclevel,xc_denpos,el_temp)
792 42 : if (option/=1) then
793 42 : e1xc=e1xc-eextemp*exchmix
794 42 : e1xcdc=e1xcdc-eexdctemp*exchmix
795 42 : s1xc=s1xc-ssxtemp*exchmix
796 : end if
797 42 : if (option<2) then
798 2270226 : paw_an(iatom)%vxc_ex(:,:,:)=vxc_tmp(:,:,:)
799 42 : paw_an(iatom)%has_vxc_ex=2
800 : end if
801 42 : ABI_FREE(rho1xx)
802 42 : ABI_FREE(vxc_tmp)
803 42 : ABI_FREE(kxc_tmp)
804 :
805 : end if ! useexexch
806 :
807 49421 : itypat0=0;if (iatom<my_natom) itypat0=pawrhoij(iatom+1)%itypat
808 49421 : if (itypat/=itypat0) then
809 28897 : ABI_FREE(one_over_rad2)
810 : end if
811 :
812 49421 : ABI_FREE(lmselect_cur)
813 :
814 : ! ==== Compute Hartree potential terms and some energy terms ====
815 : ! ===============================================================
816 :
817 : ! Hartree Dij computation
818 49421 : if (ipositron/=1) then
819 47030 : eijkl_is_sym=.true.
820 47030 : if(usercpaw) then
821 35 : eijkl_is_sym=rcpaw%eijkl_is_sym(itypat)
822 : endif
823 : call pawdijhartree(paw_ij(iatom)%dijhartree,cplex,nspden,pawrhoij(iatom),pawtab(itypat),&
824 47030 : & is_sym=eijkl_is_sym)
825 47030 : if (option/=1) then
826 39658 : call pawaccenergy_nospin(eh2dc,pawrhoij(iatom),paw_ij(iatom)%dijhartree,1,qphase,pawtab(itypat))
827 : end if
828 47030 : if(extfpmd_pawsph) then
829 42 : ABI_MALLOC(vh,(mesh_size))
830 28 : ABI_MALLOC(rho,(mesh_size))
831 28 : ABI_MALLOC(ff,(mesh_size))
832 19712 : ff=zero
833 19712 : vh=zero
834 : ! vhnzc-vhtnzc
835 : ff(2:mesh_size)=(pawtab(itypat)%vhnzc(2:mesh_size)-pawtab(itypat)%vhtnzc(2:mesh_size))*&
836 19698 : & four_pi*pawrad(itypat)%rad(2:mesh_size)**2
837 19712 : do ii=1,mesh_size
838 19712 : if(pawrad(itypat)%rad(ii)>pawtab(itypat)%rpaw) then
839 84 : ff(ii)=zero
840 : endif
841 : enddo
842 14 : call simp_gen(eshift,ff,pawrad(itypat))
843 14 : ehpw=ehpw+eshift*extfpmd_rho
844 14 : extfpmd%eshift_paw=extfpmd%eshift_paw+eshift/ucvol
845 19712 : rho(1:mesh_size)=extfpmd_rho*four_pi*pawrad(itypat)%rad(1:mesh_size)**2
846 19712 : do ii=1,mesh_size
847 19712 : if(pawrad(itypat)%rad(ii)>pawtab(itypat)%rpaw) then
848 84 : rho(ii)=zero
849 : endif
850 : enddo
851 14 : call poisson(rho,0,pawrad(itypat),vh)
852 19698 : do ii=2,mesh_size
853 19698 : vh(ii)=vh(ii)/pawrad(itypat)%rad(ii)
854 : enddo
855 14 : call pawrad_deducer0(vh,mesh_size,pawrad(itypat))
856 686 : do jlmn=1,pawtab(itypat)%lmn_size
857 672 : j0lmn=jlmn*(jlmn-1)/2
858 672 : jlm=pawtab(itypat)%indlmn(4,jlmn);jln=pawtab(itypat)%indlmn(5,jlmn)
859 17150 : do ilmn=1,jlmn
860 16464 : klmn=j0lmn+ilmn
861 16464 : ilm=pawtab(itypat)%indlmn(4,ilmn);iln=pawtab(itypat)%indlmn(5,ilmn)
862 17136 : if (jlm==ilm) then
863 : ff(1:mesh_size)=vh(1:mesh_size)*(pawtab(itypat)%phi(1:mesh_size,iln)*pawtab(itypat)%phi(1:mesh_size,jln)-&
864 1892352 : & pawtab(itypat)%tphi(1:mesh_size,iln)*pawtab(itypat)%tphi(1:mesh_size,jln))
865 1344 : call simp_gen(intvh,ff,pawrad(itypat))
866 1344 : paw_ij(iatom)%dijhartree(klmn)=paw_ij(iatom)%dijhartree(klmn)+intvh
867 1892352 : ff(1:mesh_size)=vh(1:mesh_size)*pawtab(itypat)%shapefunc(1:mesh_size,1)*pawrad(itypat)%rad(1:mesh_size)**2
868 1344 : call simp_gen(intvh,ff,pawrad(itypat))
869 1344 : intg=pawtab(itypat)%qijl(1,klmn)
870 1344 : paw_ij(iatom)%dijhartree(klmn)=paw_ij(iatom)%dijhartree(klmn)-intvh*intg*sqrt(four_pi)
871 : end if
872 : end do
873 : end do
874 19712 : rho=zero
875 28 : do ispden=1,nspden
876 19712 : rho(1:mesh_size)=sqrt(four_pi)*(rho1(1:mesh_size,1,ispden)-trho1(1:mesh_size,1,ispden)-nhat1(1:mesh_size,1,ispden))*pawrad(itypat)%rad(1:mesh_size)**2
877 14 : call poisson(rho,0,pawrad(itypat),vh)
878 19698 : do ii=2,mesh_size
879 19698 : vh(ii)=vh(ii)/pawrad(itypat)%rad(ii)
880 : enddo
881 14 : call pawrad_deducer0(vh,mesh_size,pawrad(itypat))
882 19712 : vh(1:mesh_size)=vh(1:mesh_size)*four_pi*pawrad(itypat)%rad(1:mesh_size)**2
883 19712 : do ii=1,mesh_size
884 19712 : if(pawrad(itypat)%rad(ii)>pawtab(itypat)%rpaw) then
885 84 : vh(ii)=zero
886 : endif
887 : enddo
888 14 : call simp_gen(eshift,vh,pawrad(itypat))
889 28 : extfpmd%eshift_paw=extfpmd%eshift_paw+eshift/ucvol/nspden
890 : enddo
891 14 : ABI_FREE(vh)
892 14 : ABI_FREE(ff)
893 28 : ABI_FREE(rho)
894 : endif
895 : else
896 89811 : paw_ij(iatom)%dijhartree(:)=zero
897 : end if
898 49421 : paw_ij(iatom)%has_dijhartree=2
899 :
900 : ! Hartree energy computation
901 49421 : if (option/=1) then
902 42049 : call pawaccenergy_nospin(eh2,pawrhoij(iatom),paw_ij(iatom)%dijhartree,1,qphase,pawtab(itypat))
903 42049 : if(extfpmd_pawsph) then
904 14 : eh2dc=eh2dc+two*eshift*extfpmd_rho
905 14 : eh2=eh2+eshift*extfpmd_rho
906 : endif
907 : end if
908 :
909 : ! Core + nucleus Hartree energy accumulation
910 : ! Core kinetic energy accumulation
911 49421 : if (option/=1) then
912 42049 : ehnzc=ehnzc+pawtab(itypat)%ehnzc
913 42049 : ekincore=ekincore+pawtab(itypat)%ekincore
914 : end if
915 :
916 : ! Electron-positron calculation:
917 : ! - Compute Dij due to fixed particles (elec. or pos. depending on calctype)
918 : ! - Compute contribution to energy
919 : ! - Add electron and positron
920 49421 : if (ipositron/=0) then
921 5876 : ABI_CHECK(qphase==1,'qphase should be 1 for electron-positron!')
922 17628 : ABI_MALLOC(dij_ep,(qphase*lmn2_size))
923 5876 : call pawdijhartree(dij_ep,qphase,nspden,electronpositron%pawrhoij_ep(iatom),pawtab(itypat))
924 5876 : if (option/=1) then
925 5876 : etmp=zero
926 5876 : call pawaccenergy_nospin(etmp,pawrhoij(iatom),dij_ep,1,1,pawtab(itypat))
927 5876 : electronpositron%e_paw =electronpositron%e_paw -etmp
928 5876 : electronpositron%e_pawdc=electronpositron%e_pawdc-etmp
929 : end if
930 218756 : paw_ij(iatom)%dijhartree(:)=paw_ij(iatom)%dijhartree(:)-dij_ep(:)
931 5876 : ABI_FREE(dij_ep)
932 : end if
933 :
934 : ! Compute 1st moment of total Hartree potential VH(n_Z+n_core+n1)
935 : ! Equation 10 (density) and up to 43 (Hartree potential of density)
936 : ! of Kresse and Joubert PRB 59 1758 (1999) [[cite:Kresse1999]]
937 49421 : keep_vhartree=(paw_an(iatom)%has_vhartree>0)
938 196501 : if ((pawspnorb>0.and.ipert==0.and.ipositron/=1).or.keep_vhartree.or.&
939 : (any(abs(nucdipmom(:,iatom))>tol8))) then
940 :
941 : !In the first clause case, would it not be simpler just to turn on has_vhartree?
942 2619 : if (.not. allocated(paw_an(iatom)%vh1)) then
943 7185 : ABI_MALLOC(paw_an(iatom)%vh1,(cplex*mesh_size,1,1))
944 : end if
945 2619 : if (.not. allocated(paw_an(iatom)%vht1)) then
946 522 : ABI_MALLOC(paw_an(iatom)%vht1,(cplex*mesh_size,1,1))
947 : end if
948 7857 : ABI_MALLOC(rho,(mesh_size))
949 5238 : ABI_MALLOC(vh,(mesh_size))
950 :
951 : ! Construct vh1 and tvh1
952 5238 : do iq=1,cplex
953 : !Construct vh1
954 : ! The sqrt(4pi) factor comes from the fact we are calculating the spherical moments,
955 : ! and for the 00 channel the prefactor of Y_00 = 2 sqrt(pi)
956 2955491 : rho(1:mesh_size)=rho1(iq:cplex*mesh_size:cplex,1,1)*four_pi*pawrad(itypat)%rad(1:mesh_size)**2
957 2619 : if (usecore==1) then
958 : rho(1:mesh_size)=rho(1:mesh_size)+sqrt(four_pi)*pawtab(itypat)%coredens(1:mesh_size) &
959 2955491 : & *four_pi*pawrad(itypat)%rad(1:mesh_size)**2
960 : end if
961 2619 : call poisson(rho,0,pawrad(itypat),vh)
962 2952872 : vh(2:mesh_size)=(vh(2:mesh_size)-sqrt(four_pi)*znucl(itypat))/pawrad(itypat)%rad(2:mesh_size)
963 2619 : call pawrad_deducer0(vh,mesh_size,pawrad(itypat))
964 2955491 : paw_an(iatom)%vh1(iq:cplex*mesh_size:cplex,1,1)=vh(1:mesh_size)
965 : ! TODO: check this is equivalent to the previous version (commented) which explicitly recalculated VH(coredens)
966 : ! DONE: numerically there are residual differences on abiref (7th digit).
967 : ! paw_an(iatom)%vh1(2:mesh_size,1,1)=paw_an(iatom)%vh1(2:mesh_size,1,1)/pawrad(itypat)%rad(2:mesh_size) &
968 : !& +sqrt(four_pi) * pawtab(itypat)%VHnZC(2:mesh_size)
969 :
970 : !Same for vht1
971 2955491 : rho(1:mesh_size)=trho1(iq:cplex*mesh_size:cplex,1,1)*four_pi*pawrad(itypat)%rad(1:mesh_size)**2
972 2619 : if (usenhat/=0) then
973 : rho(1:mesh_size)=rho(1:mesh_size)+nhat1(iq:cplex*mesh_size:cplex,1,1) &
974 937391 : & *four_pi*pawrad(itypat)%rad(1:mesh_size)**2
975 : end if
976 2619 : if (usetcore==1) then
977 : rho(1:mesh_size)=rho(1:mesh_size)+sqrt(four_pi)*pawtab(itypat)%tcoredens(1:mesh_size,1) &
978 2042243 : & *four_pi*pawrad(itypat)%rad(1:mesh_size)**2
979 : end if
980 2619 : call poisson(rho,0,pawrad(itypat),vh)
981 2952872 : vh(2:mesh_size)=(vh(2:mesh_size)-sqrt(four_pi)*znucl(itypat))/pawrad(itypat)%rad(2:mesh_size)
982 2619 : call pawrad_deducer0(vh,mesh_size,pawrad(itypat))
983 2958110 : paw_an(iatom)%vht1(iq:cplex*mesh_size:cplex,1,1)=vh(1:mesh_size)
984 :
985 : end do ! cplex phase
986 :
987 2619 : paw_an(iatom)%has_vhartree=2
988 2619 : ABI_FREE(rho)
989 2619 : ABI_FREE(vh)
990 : end if
991 :
992 : ! ========= Compute PAW+U and energy contribution =========
993 : ! ==========================================================
994 :
995 49421 : if (usepawu/=0.and.usepawu<10.and.ipositron/=1.and.option/=1) then
996 :
997 2685 : if (pawu_algo==PAWU_ALGO_1) then
998 :
999 : ! PAW+U energy computation from nocc_m_mp
1000 : call pawuenergy(iatom_tot,edftumdc,edftumdcdc,paw_ij(iatom)%noccmmp, &
1001 2463 : & paw_ij(iatom)%nocctot,pawprtvol,pawtab(itypat))
1002 : else
1003 :
1004 : ! PAW+U energy computation from eU_ijkl
1005 : !First, compute DijU
1006 : call pawdiju_euijkl(paw_ij(iatom)%dijU,cplex_dij,cplex,ndij, &
1007 222 : & pawrhoij(iatom),pawtab(itypat))
1008 222 : paw_ij(iatom)%has_dijU=2
1009 : !Then, compute energy
1010 : if (option/=1) then
1011 222 : etmp=zero
1012 222 : call pawaccenergy(etmp,pawrhoij(iatom),paw_ij(iatom)%dijU,cplex_dij,qphase,ndij,pawtab(itypat))
1013 222 : edftumdc=edftumdc+half*etmp ; edftumdcdc=edftumdcdc-half*etmp
1014 : !Add FLL double-counting part
1015 222 : if (pawu_dblec==PAWU_FLL.and.ipert==0) then
1016 192 : ABI_CHECK(qphase==1,'BUG in pawdenpot: qphase should be 1 for Dble-C FLL term!')
1017 192 : call pawaccenergy_nospin(edftufll,pawrhoij(iatom),pawtab(itypat)%euij_fll,1,1,pawtab(itypat))
1018 : end if
1019 : end if
1020 :
1021 : end if ! DFT+U algo
1022 : end if ! Dij Hartree
1023 :
1024 : ! ========= Compute nuclear dipole moment energy contribution ========
1025 : ! =====================================================================
1026 :
1027 196501 : if (any(abs(nucdipmom(:,iatom))>tol8).and.ipert==0.and.ipositron/=1) then
1028 :
1029 535 : ABI_CHECK(cplex_rhoij==2,'BUG in pawdenpot: rhoij must be complex for ND moments!')
1030 535 : ABI_CHECK(qphase==1,'BUG in pawdenpot: qphase should be 1 for ND moments!')
1031 :
1032 : ! Compute nuclear dipole contribution to Dij if necessary
1033 535 : if (paw_ij(iatom)%has_dijnd/=2) then
1034 : call pawdijnd(paw_ij(iatom)%dijnd,cplex_dij,gprimd,iatom,natom,ndij,nspden,&
1035 : & nucdipmom(1:3,1:natom),pawang,pawrad(itypat),pawtab(itypat),pawxcdev,qphase,&
1036 535 : & paw_an(iatom)%vh1,paw_an(iatom)%vxc1,xred,znucl(itypat),paw_ij(iatom)%zora)
1037 535 : paw_ij(iatom)%has_dijnd=2
1038 : end if
1039 :
1040 : ! Compute nuclear dipole contribution to energy
1041 535 : if (option/=1) then
1042 : call pawaccenergy_nospin(enucdip,pawrhoij(iatom),paw_ij(iatom)%dijnd,&
1043 520 : & cplex_dij,1,pawtab(itypat))
1044 : end if
1045 :
1046 : end if
1047 :
1048 : ! ========= Compute spin-orbit energy contribution ========
1049 : ! ==========================================================
1050 :
1051 49421 : if (pawspnorb>0.and.ipert==0.and.ipositron/=1) then
1052 :
1053 : ! Compute spin-orbit contribution to Dij
1054 2220 : if (option/=2.or.cplex_rhoij==2) then
1055 : call pawdijso(paw_ij(iatom)%dijso,cplex_dij,cplex,ndij,nspden,pawang,&
1056 : & pawrad(itypat),pawtab(itypat),pawxcdev,spnorbscl,paw_an(iatom)%vh1,&
1057 : & paw_an(iatom)%vxc1,znucl(itypat),paw_ij(iatom)%zora,&
1058 2220 : & nucdipmom=nucdipmom(1:3,iatom),spinaxis=spinaxis)
1059 2220 : paw_ij(iatom)%has_dijso=2
1060 : end if
1061 :
1062 : ! Compute spin-orbit contribution to on-site energy
1063 2220 : if (option/=1.and.cplex_rhoij==2) then
1064 : call pawaccenergy(espnorb,pawrhoij(iatom),paw_ij(iatom)%dijso,&
1065 2190 : & cplex_dij,qphase,ndij,pawtab(itypat))
1066 : end if
1067 :
1068 : end if
1069 :
1070 : ! === Compute 2nd part of local exact-exchange energy and potential ===
1071 : ! ======================================================================
1072 :
1073 49421 : if (pawtab(itypat)%useexexch/=0.and.ipert==0.and.ipositron/=1) then
1074 :
1075 42 : ABI_CHECK(paw_ij(iatom)%nspden/=4,'BUG in pawdenpot: Local ex-exch. not implemented for nspden=4!')
1076 42 : if (option<2) then
1077 42 : call pawxpot(ndij,pawprtvol,pawrhoij(iatom),pawtab(itypat),paw_ij(iatom)%vpawx)
1078 42 : paw_ij(iatom)%has_exexch_pot=2
1079 : end if
1080 42 : if (option/=1) then
1081 42 : if (abs(pawprtvol)>=2) then
1082 0 : write(msg, '(2a)' )ch10,'======= PAW local exact exchange terms (in Hartree) ===='
1083 0 : call wrtout(std_out, msg,'COLL')
1084 0 : write(msg, '(2a,i4)' )ch10,' For Atom',iatom_tot
1085 0 : call wrtout(std_out, msg,'COLL')
1086 : end if
1087 42 : call pawxenergy(eexex,pawprtvol,pawrhoij(iatom),pawtab(itypat))
1088 : end if
1089 :
1090 : end if ! useexexch
1091 :
1092 : ! ==== Compute Fock Dij term and Fock energy terms ====
1093 : ! =====================================================
1094 :
1095 49421 : if (usefock==1) then
1096 :
1097 620 : if (ipositron/=1) then
1098 :
1099 : ! Fock contribution to Dij
1100 2480 : ABI_MALLOC(dijfock_vv,(cplex_dij*qphase*lmn2_size,ndij))
1101 1860 : ABI_MALLOC(dijfock_cv,(cplex_dij*qphase*lmn2_size,ndij))
1102 : call pawdijfock(dijfock_vv,dijfock_cv,cplex_dij,cplex,hyb_mixing_,hyb_mixing_sr_, &
1103 620 : & ndij,pawrhoij(iatom),pawtab(itypat))
1104 23385 : paw_ij(iatom)%dijfock(:,:)=dijfock_vv(:,:)+dijfock_cv(:,:)
1105 620 : paw_ij(iatom)%has_dijfock=2
1106 :
1107 : ! Fock contribution to energy
1108 620 : if (option/=1) then
1109 37000 : dijfock_vv(:,:)=half*dijfock_vv(:,:) ; dijfock_cv(:,:)=dijfock_vv(:,:)+dijfock_cv(:,:)
1110 498 : call pawaccenergy(efock ,pawrhoij(iatom),dijfock_cv,cplex_dij,qphase,ndij,pawtab(itypat))
1111 498 : call pawaccenergy(efockdc,pawrhoij(iatom),dijfock_vv,cplex_dij,qphase,ndij,pawtab(itypat))
1112 : end if
1113 :
1114 620 : ABI_FREE(dijfock_vv)
1115 620 : ABI_FREE(dijfock_cv)
1116 : end if
1117 :
1118 : ! Special case for positron
1119 620 : if (ipositron==1) then
1120 0 : paw_ij(iatom)%dijfock(:,:)=zero
1121 0 : paw_ij(iatom)%has_dijfock=2
1122 : end if
1123 :
1124 : end if
1125 :
1126 : ! === Compute the zero of the potentials if requested ==================
1127 : ! ======================================================================
1128 :
1129 49421 : if (pawtab(itypat)%usepotzero==1.and.present(vpotzero).and.ipert<=0) then
1130 :
1131 : !Term 1 : beta
1132 63 : vpotzero(1)=vpotzero(1)-pawtab(itypat)%beta/ucvol
1133 :
1134 : !Term 2 : \sum_ij rho_ij gamma_ij
1135 63 : etmp=zero
1136 63 : call pawaccenergy_nospin(etmp,pawrhoij(iatom),pawtab(itypat)%gammaij,1,1,pawtab(itypat))
1137 63 : vpotzero(2)=vpotzero(2)-etmp/ucvol
1138 :
1139 : end if
1140 :
1141 : ! ======= Compute atomic contribution to the energy (Dij0) ===========
1142 : ! ======================================================================
1143 :
1144 49421 : if (option/=1.and.ipert<=0) then
1145 36598 : call pawaccenergy_nospin(e1t10,pawrhoij(iatom),pawtab(itypat)%dij0,1,1,pawtab(itypat))
1146 : ! Positron special case (dij0 is opposite, except for kinetic term)
1147 36598 : if (ipositron==1) then
1148 7173 : ABI_MALLOC(dij_ep,(lmn2_size))
1149 89811 : dij_ep(:)=two*(pawtab(itypat)%kij(:)-pawtab(itypat)%dij0(:))
1150 2391 : call pawaccenergy_nospin(e1t10,pawrhoij(iatom),dij_ep,1,1,pawtab(itypat))
1151 2391 : ABI_FREE(dij_ep)
1152 : end if
1153 :
1154 : end if
1155 :
1156 : ! ==========================================================
1157 : ! No more need of some densities/potentials
1158 :
1159 : ! Deallocate densities
1160 49421 : ABI_FREE(rho1)
1161 49421 : ABI_FREE(trho1)
1162 49421 : ABI_FREE(nhat1)
1163 49421 : if (usekden==1) then
1164 484 : ABI_FREE(tau1)
1165 484 : ABI_FREE(ttau1)
1166 : end if
1167 49421 : if (ipositron/=0) then
1168 5876 : ABI_FREE(rho1_ep)
1169 5876 : ABI_FREE(trho1_ep)
1170 5876 : ABI_FREE(nhat1_ep)
1171 : end if
1172 :
1173 : ! Deallocate potentials
1174 49421 : if (.not.keep_vhartree) then
1175 49197 : paw_an(iatom)%has_vhartree=0
1176 49197 : if (allocated(paw_an(iatom)%vh1)) then
1177 2395 : ABI_FREE(paw_an(iatom)%vh1)
1178 : end if
1179 : end if
1180 72298 : if (temp_vxc) then
1181 0 : paw_an(iatom)%has_vxc=0
1182 0 : if (allocated(paw_an(iatom)%vxc1)) then
1183 0 : ABI_FREE(paw_an(iatom)%vxc1)
1184 : end if
1185 : end if
1186 :
1187 : ! =========== End loop on atoms ============================
1188 : ! ==========================================================
1189 :
1190 : end do
1191 :
1192 : !========== Assemble "on-site" energy terms ===============
1193 : !==========================================================
1194 :
1195 22877 : if (option/=1) then
1196 18979 : if (ipert==0) then
1197 16302 : paw_energies%epaw=e1xc+half*eh2+e1t10-exccore-etild1xc+edftumdc+edftufll+eexex+espnorb+efock+enucdip+ehpw
1198 16302 : paw_energies%epaw_dc=e1xc-e1xcdc-half*eh2dc-exccore-etild1xc+etild1xcdc+edftumdcdc-eexex-efockdc
1199 16302 : paw_energies%epaw_xc=e1xc-etild1xc !+eexex
1200 16302 : paw_energies%entropy_paw=s1xc-sxccore-stild1xc ! PAW entropy coming from finite-temperature xc functionals
1201 16302 : paw_energies%epaw_core=zero ; paw_energies%epaw_core_dc=zero
1202 16302 : if (abs(ekincore)>tiny(zero).and.ipositron/=1) then
1203 6295 : paw_energies%epaw_core=exccore+ekincore+ehnzc
1204 6295 : paw_energies%epaw_core_dc=paw_energies%epaw_core
1205 : end if
1206 16302 : if (add_core_energy==1.and.(.not.usercpaw)) then
1207 5 : paw_energies%epaw=paw_energies%epaw+paw_energies%epaw_core
1208 5 : paw_energies%epaw_dc=paw_energies%epaw_dc+paw_energies%epaw_core_dc
1209 : endif
1210 : else
1211 2677 : paw_energies%epaw=e1xc-etild1xc+eh2+two*edftumdc
1212 2677 : paw_energies%epaw_dc=zero
1213 2677 : paw_energies%epaw_core=zero
1214 2677 : paw_energies%epaw_core_dc=zero
1215 2677 : paw_energies%entropy_paw=zero ! Force PAW entropy contribution to zero when using RF.
1216 : end if
1217 : end if
1218 :
1219 : !========== Reduction in case of parallelism ==============
1220 : !==========================================================
1221 :
1222 22877 : if (paral_atom) then
1223 7766 : if (option/=1) then
1224 6296 : call timab(48,1,tsec)
1225 6296 : ABI_MALLOC(mpiarr,(5+n_paw_energies))
1226 75552 : mpiarr=zero ; ii=0
1227 6296 : call paw_energies_to_array(paw_energies,mpiarr(ii+1:ii+n_paw_energies),1)
1228 6296 : ii=ii+n_paw_energies
1229 6296 : mpiarr(ii+1)=compch_sph ; ii=ii+1
1230 6296 : if (ipositron/=0) then
1231 0 : mpiarr(ii+1)=electronpositron%e_paw
1232 0 : mpiarr(ii+2)=electronpositron%e_pawdc
1233 0 : ii=ii+2
1234 : end if
1235 6296 : if (present(vpotzero)) then
1236 14586 : mpiarr(ii+1:ii+2)=vpotzero(1:2)
1237 6296 : ii=ii+2
1238 : end if
1239 6296 : call xmpi_sum(mpiarr,my_comm_atom,ierr)
1240 6296 : ii=0
1241 6296 : call paw_energies_to_array(paw_energies,mpiarr(ii+1:ii+n_paw_energies),-1)
1242 6296 : ii=ii+n_paw_energies
1243 6296 : compch_sph=mpiarr(ii+1) ; ii=ii+1
1244 6296 : if (ipositron/=0) then
1245 0 : electronpositron%e_paw=mpiarr(ii+1)
1246 0 : electronpositron%e_pawdc=mpiarr(ii+2)
1247 0 : ii=ii+2
1248 : end if
1249 6296 : if (present(vpotzero)) then
1250 14586 : vpotzero(1:2)=mpiarr(ii+1:ii+2)
1251 6296 : ii=ii+2
1252 : end if
1253 6296 : ABI_FREE(mpiarr)
1254 12592 : call timab(48,2,tsec)
1255 : end if
1256 7766 : if(extfpmd_pawsph) then
1257 0 : call xmpi_sum(extfpmd%eshift_paw,my_comm_atom,ierr)
1258 : endif
1259 : end if
1260 :
1261 : !In case we have an entropy associated with PAW contribution
1262 : !(e.g. using finite-temperature exchange-correlation functionals),
1263 : !we retrieve exchange-correlation internal energies e_paw, e_pawdc
1264 : !using entropy spaw.
1265 22877 : if(option/=1.and.abs(paw_energies%entropy_paw)>tiny(zero)) then
1266 14 : paw_energies%epaw=paw_energies%epaw+el_temp*paw_energies%entropy_paw
1267 14 : paw_energies%epaw_dc=paw_energies%epaw_dc+el_temp*paw_energies%entropy_paw
1268 14 : paw_energies%epaw_xc=paw_energies%epaw_xc+el_temp*paw_energies%entropy_paw
1269 : end if
1270 :
1271 : !Destroy atom table used for parallelism
1272 22877 : call free_my_atmtab(my_atmtab,my_atmtab_allocated)
1273 :
1274 22877 : call timab(560,2,tsec)
1275 :
1276 : DBG_EXIT("COLL")
1277 :
1278 49030 : end subroutine pawdenpot
1279 : !!***
1280 :
1281 : !----------------------------------------------------------------------
1282 :
1283 : !!****f* m_paw_denpot/pawdensities
1284 : !! NAME
1285 : !! pawdensities
1286 : !!
1287 : !! FUNCTION
1288 : !! Compute PAW on-site densities (all-electron, pseudo and compensation) for a given atom
1289 : !!
1290 : !! INPUTS
1291 : !! cplex: if 1, on-site densities are REAL, if 2, COMPLEX (response function only)
1292 : !! iatom=index of current atom (note: this is the absolute index, not the index on current proc)
1293 : !! lm_size=number of (l,m) moments
1294 : !! lmselectin(lm_size)=flags selecting the non-zero LM-moments of on-site densities
1295 : !! (value of these flags at input; must be .TRUE. for nzlmopt/=1)
1296 : !! nspden=number of spin-density components
1297 : !! nzlmopt=if -1, compute all LM-moments of densities (lmselectin=.true. forced)
1298 : !! initialize "lmselectout" (index of non-zero LM-moments of densities)
1299 : !! if 0, compute all LM-moments of densities (lmselectin=.true. forced)
1300 : !! force "lmselectout" to .true. (index of non-zero LM-moments of densities)
1301 : !! if 1, compute only non-zero LM-moments of densities (stored before in "lmselectin")
1302 : !! one_over_rad2(mesh_size)= contains 1/r**2 for each point of the radial grid -optional argument-
1303 : !! opt_compch=flag controlling the accumulation of compensation charge density moments
1304 : !! inside PAW spheres (compch_sph)
1305 : !! opt_dens=flag controlling which on-site density(ies) is (are) computed
1306 : !! 0: all on-site densities (all-electron, pseudo and compensation)
1307 : !! 1: all-electron and pseudo densities (no compensation)
1308 : !! 2: only all-electron density
1309 : !! opt_l=controls which l-moment(s) contribute to the density:
1310 : !! <0 : all l contribute
1311 : !! >=0: only l=opt_l contributes
1312 : !! Note: opt_l>=0 is only compatible with opt_dens=2
1313 : !! opt_print=1 if the densities moments have to be printed out (only if pawprtvol>=2)
1314 : !! pawang <type(pawang_type)>=paw angular mesh and related data
1315 : !! pawprtvol=control print volume and debugging output for PAW
1316 : !! pawrad <type(pawrad_type)>=paw radial mesh and related data (for the current atom type)
1317 : !! pawrhoij <type(pawrhoij_type)>= paw rhoij occupancies and related data (for the current atom)
1318 : !! pawtab <type(pawtab_type)>=paw tabulated starting data (for the current atom type)
1319 : !!
1320 : !! OUTPUT
1321 : !! nhat1(cplex*mesh_size,lm_size,nspden)= compensation charge on-site density for current atom
1322 : !! rho1(cplex*mesh_size,lm_size,nspden)= all electron on-site density for current atom
1323 : !! trho1(cplex*mesh_size,lm_size,nspden)= pseudo on-site density for current atom
1324 : !! ==== if nzlmopt/=1
1325 : !! lmselectout(lm_size)=flags selecting the non-zero LM-moments of on-site densities
1326 : !! (value of these flags at output if updated, i.e. if nzlmopt<1)
1327 : !!
1328 : !! SIDE EFFECTS
1329 : !! ==== if opt_compch==1
1330 : !! compch_sph=compensation charge integral inside spheres computed over spherical meshes
1331 : !! updated with the contribution of current atom
1332 : !!
1333 : !! SOURCE
1334 :
1335 69247 : subroutine pawdensities(compch_sph,cplex,iatom,lmselectin,lmselectout,lm_size,nhat1,nspden,nzlmopt,&
1336 69247 : & opt_compch,opt_dens,opt_l,opt_print,pawang,pawprtvol,pawrad,pawrhoij,pawtab,rho1,trho1,&
1337 : & extfpmd_rho,one_over_rad2) ! optional
1338 :
1339 : !Arguments ---------------------------------------------
1340 : !scalars
1341 : integer,intent(in) :: cplex,iatom,lm_size,nspden,nzlmopt,opt_compch,opt_dens,opt_l,opt_print,pawprtvol
1342 : real(dp),intent(inout) :: compch_sph
1343 : real(dp),intent(in),optional :: extfpmd_rho
1344 : type(pawang_type),intent(in) :: pawang
1345 : type(pawrad_type),intent(in) :: pawrad
1346 : type(pawrhoij_type),intent(in) :: pawrhoij
1347 : type(pawtab_type),intent(in) :: pawtab
1348 : !arrays
1349 : logical,intent(in) :: lmselectin(lm_size)
1350 : logical,intent(inout) :: lmselectout(lm_size)
1351 : real(dp),intent(in),target,optional :: one_over_rad2(pawtab%mesh_size)
1352 : real(dp),intent(out) :: nhat1(cplex*pawtab%mesh_size,lm_size,nspden*(1-((opt_dens+1)/2)))
1353 : real(dp),intent(out) :: rho1(cplex*pawtab%mesh_size,lm_size,nspden)
1354 : real(dp),intent(out) :: trho1(cplex*pawtab%mesh_size,lm_size,nspden*(1-(opt_dens/2)))
1355 : !Local variables ---------------------------------------
1356 : !scalars
1357 : integer :: dplex,ii,ilm,iplex,iq0,ir,irhoij,isel,ispden,jrhoij
1358 : integer :: klm,klmn,kln,ll,lmax,lmin,mesh_size
1359 : real(dp) :: m1,mt1,rdum
1360 : character(len=500) :: msg
1361 : !arrays
1362 138494 : real(dp) :: compchspha(cplex),compchsphb(cplex),ro(cplex),ro_ql(cplex),ro_rg(cplex)
1363 69247 : real(dp),allocatable :: aa(:),bb(:)
1364 69247 : real(dp),pointer :: one_over_rad2_(:)
1365 :
1366 : ! *************************************************************************
1367 :
1368 : DBG_ENTER("COLL")
1369 :
1370 : !Compatibility tests
1371 69196 : if (opt_dens/=2.and.opt_l>=0) then
1372 0 : msg='opt_dens/=2 incompatible with opt_l>=0!'
1373 0 : ABI_BUG(msg)
1374 : end if
1375 69247 : if(nzlmopt/=0.and.nzlmopt/=1.and.nzlmopt/=-1) then
1376 0 : msg='invalid value for variable "nzlmopt"!'
1377 0 : ABI_BUG(msg)
1378 : end if
1379 69247 : if(nspden>pawrhoij%nspden) then
1380 0 : msg='nspden must be <= pawrhoij%nspden!'
1381 0 : ABI_BUG(msg)
1382 : end if
1383 69247 : if (cplex>pawrhoij%qphase) then
1384 0 : msg='cplex must be <= pawrhoij%qphase!'
1385 0 : ABI_BUG(msg)
1386 : end if
1387 69247 : if (nzlmopt/=1) then
1388 268206 : if (any(.not.lmselectin(1:lm_size))) then
1389 0 : msg='With nzlmopt/=1, lmselectin must be true!'
1390 0 : ABI_BUG(msg)
1391 : end if
1392 : end if
1393 69247 : if (pawang%gnt_option==0) then
1394 0 : msg='pawang%gnt_option=0!'
1395 0 : ABI_BUG(msg)
1396 : end if
1397 :
1398 : !Various inits
1399 895125056 : rho1=zero
1400 892747376 : if (opt_dens<2) trho1=zero
1401 312336162 : if (opt_dens==0) nhat1=zero
1402 69247 : mesh_size=pawtab%mesh_size;dplex=cplex-1
1403 69247 : iq0=pawrhoij%cplex_rhoij*pawrhoij%lmn2_size
1404 312309 : if (nzlmopt<1) lmselectout(1:lm_size)=.true.
1405 69247 : if (present(one_over_rad2)) then
1406 55309 : one_over_rad2_ => one_over_rad2
1407 : else
1408 41814 : ABI_MALLOC(one_over_rad2_,(mesh_size))
1409 13938 : one_over_rad2_(1)=zero
1410 11065792 : one_over_rad2_(2:mesh_size)=one/pawrad%rad(2:mesh_size)**2
1411 : end if
1412 :
1413 : !===== Compute "on-site" densities (n1, ntild1, nhat1) =====
1414 : !===========================================================
1415 :
1416 152280 : do ispden=1,nspden
1417 :
1418 : ! -- Loop over ij channels (basis components)
1419 83033 : jrhoij=1
1420 3101635 : do irhoij=1,pawrhoij%nrhoijsel
1421 3018602 : klmn=pawrhoij%rhoijselect(irhoij)
1422 3018602 : klm =pawtab%indklmn(1,klmn)
1423 3018602 : kln =pawtab%indklmn(2,klmn)
1424 3018602 : lmin=pawtab%indklmn(3,klmn)
1425 3018602 : lmax=pawtab%indklmn(4,klmn)
1426 :
1427 : ! Retrieve rhoij
1428 3018602 : if (pawrhoij%nspden/=2) then
1429 2211724 : ro(1)=pawrhoij%rhoijp(jrhoij,ispden)
1430 2211724 : if (cplex==2) ro(2)=pawrhoij%rhoijp(iq0+jrhoij,ispden)
1431 : else
1432 806878 : if (ispden==1) then
1433 403439 : ro(1)=pawrhoij%rhoijp(jrhoij,1)+pawrhoij%rhoijp(jrhoij,2)
1434 403439 : if (cplex==2) ro(2)=pawrhoij%rhoijp(iq0+jrhoij,1)+pawrhoij%rhoijp(iq0+jrhoij,2)
1435 403439 : else if (ispden==2) then
1436 403439 : ro(1)=pawrhoij%rhoijp(jrhoij,1)
1437 403439 : if (cplex==2) ro(2)=pawrhoij%rhoijp(iq0+jrhoij,1)
1438 : end if
1439 : end if
1440 6094771 : ro(1:cplex)=pawtab%dltij(klmn)*ro(1:cplex)
1441 :
1442 : ! First option: all on-site densities are computed (opt_dens==0)
1443 : ! --------------------------------------------------------------
1444 3018602 : if (opt_dens==0) then
1445 4406976 : do ll=lmin,lmax,2
1446 17230950 : do ilm=ll**2+1,(ll+1)**2
1447 15780862 : if (lmselectin(ilm)) then
1448 11044031 : isel=pawang%gntselect(ilm,klm)
1449 11044031 : if (isel>0) then
1450 5971011 : ro_ql(1:cplex)=ro(1:cplex)*pawtab%qijl(ilm,klmn)
1451 5971011 : ro_rg(1:cplex)=ro(1:cplex)*pawang%realgnt(isel)
1452 : ! == nhat1(r=0)
1453 : nhat1(1:cplex,ilm,ispden)=nhat1(1:cplex,ilm,ispden) &
1454 5971011 : & +ro_ql(1:cplex)*pawtab%shapefunc(1,ll+1)
1455 : ! == rho1(r>0), trho1(r>0), nhat1(r>0)
1456 2108324153 : do ir=2,mesh_size
1457 : rho1(cplex*ir-dplex:ir*cplex,ilm,ispden) =rho1(cplex*ir-dplex:ir*cplex,ilm,ispden)&
1458 4236216710 : & +ro_rg(1:cplex)*pawtab%phiphj(ir,kln)*one_over_rad2_(ir)
1459 : trho1(cplex*ir-dplex:ir*cplex,ilm,ispden)=trho1(cplex*ir-dplex:ir*cplex,ilm,ispden)&
1460 4236216710 : & +ro_rg(1:cplex)*pawtab%tphitphj(ir,kln)*one_over_rad2_(ir)
1461 : nhat1(cplex*ir-dplex:ir*cplex,ilm,ispden)=nhat1(cplex*ir-dplex:ir*cplex,ilm,ispden)&
1462 4239177224 : & +ro_ql(1:cplex)*pawtab%shapefunc(ir,ll+1)
1463 : end do
1464 : end if
1465 : end if
1466 : end do ! End loops over ll,lm
1467 : end do
1468 :
1469 : ! 2nd option: AE and pseudo densities are computed (opt_dens==1)
1470 : ! --------------------------------------------------------------
1471 1568514 : else if (opt_dens==1) then
1472 4366202 : do ll=lmin,lmax,2
1473 14460602 : do ilm=ll**2+1,(ll+1)**2
1474 12899154 : if (lmselectin(ilm)) then
1475 9252110 : isel=pawang%gntselect(ilm,klm)
1476 9252110 : if (isel>0) then
1477 5367745 : ro_rg(1:cplex)=ro(1:cplex)*pawang%realgnt(isel)
1478 : ! == rho1(r>0), trho1(r>0)
1479 3184217243 : do ir=2,mesh_size
1480 : rho1(cplex*ir-dplex:ir*cplex,ilm,ispden) =rho1(cplex*ir-dplex:ir*cplex,ilm,ispden)&
1481 6407528575 : & +ro_rg(1:cplex)*pawtab%phiphj (ir,kln)*one_over_rad2_(ir)
1482 : trho1(cplex*ir-dplex:ir*cplex,ilm,ispden)=trho1(cplex*ir-dplex:ir*cplex,ilm,ispden)&
1483 6410193539 : & +ro_rg(1:cplex)*pawtab%tphitphj(ir,kln)*one_over_rad2_(ir)
1484 : end do
1485 : end if
1486 : end if
1487 : end do ! End loops over ll,lm
1488 : end do
1489 :
1490 : ! 3rd option: only all-electron on-site density is computed (opt_dens==2)
1491 : ! -----------------------------------------------------------------------
1492 7066 : else if (opt_dens==2) then
1493 7066 : if (opt_l<0.or.(pawtab%indklmn(3,klmn)==0.and.pawtab%indklmn(4,klmn)==2*opt_l)) then
1494 16241 : do ll=lmin,lmax,2
1495 78346 : do ilm=ll**2+1,(ll+1)**2
1496 71280 : if (lmselectin(ilm)) then
1497 55649 : isel=pawang%gntselect(ilm,klm)
1498 55649 : if (isel>0) then
1499 24414 : ro_rg(1:cplex)=ro(1:cplex)*pawang%realgnt(isel)
1500 : ! == rho1(r>0)
1501 13245262 : do ir=2,mesh_size
1502 : rho1(cplex*ir-dplex:ir*cplex,ilm,ispden) =rho1(cplex*ir-dplex:ir*cplex,ilm,ispden)&
1503 26478317 : & +ro_rg(1:cplex)*pawtab%phiphj(ir,kln)*one_over_rad2_(ir)
1504 : end do
1505 : end if
1506 : end if
1507 : end do ! End loops over ll, lm
1508 : end do
1509 : end if
1510 : end if
1511 :
1512 :
1513 : ! -- End loop over ij channels
1514 3101635 : jrhoij=jrhoij+pawrhoij%cplex_rhoij
1515 : end do
1516 :
1517 : ! Compute rho1(r=0) and trho1(r=0)
1518 83033 : if (cplex==2) then
1519 1520 : ABI_MALLOC(aa,(5))
1520 1520 : ABI_MALLOC(bb,(5))
1521 : end if
1522 83033 : if (opt_dens==0.or.opt_dens==1) then
1523 349118 : do ll=0,pawtab%lcut_size-1
1524 1305876 : do ilm=ll**2+1,(ll+1)**2
1525 1222936 : if (lmselectin(ilm)) then
1526 754781 : if (cplex==1) then
1527 740477 : call pawrad_deducer0(rho1 (:,ilm,ispden),mesh_size,pawrad)
1528 740477 : call pawrad_deducer0(trho1(:,ilm,ispden),mesh_size,pawrad)
1529 : else
1530 42912 : do ii=0,1
1531 143040 : do ir=2,5
1532 114432 : aa(ir)=rho1 (2*ir-ii,ilm,ispden)
1533 143040 : bb(ir)=trho1(2*ir-ii,ilm,ispden)
1534 : end do
1535 28608 : call pawrad_deducer0(aa,5,pawrad)
1536 28608 : call pawrad_deducer0(bb,5,pawrad)
1537 28608 : rho1 (2-ii,ilm,ispden)=aa(1)
1538 42912 : trho1(2-ii,ilm,ispden)=bb(1)
1539 : end do
1540 : end if
1541 : end if
1542 : end do
1543 : end do
1544 : else
1545 542 : do ll=0,pawtab%lcut_size-1
1546 2739 : do ilm=ll**2+1,(ll+1)**2
1547 2646 : if (lmselectin(ilm)) then
1548 1361 : if (cplex==1) then
1549 1361 : call pawrad_deducer0(rho1(:,ilm,ispden),mesh_size,pawrad)
1550 : else
1551 0 : do ii=0,1
1552 0 : do ir=2,5
1553 0 : aa(ir)=rho1 (2*ir-ii,ilm,ispden)
1554 : end do
1555 0 : call pawrad_deducer0(aa,5,pawrad)
1556 0 : rho1(2-ii,ilm,ispden)=aa(1)
1557 : end do
1558 : end if
1559 : end if
1560 : end do
1561 : end do
1562 : end if
1563 83033 : if (cplex==2) then
1564 1520 : ABI_FREE(aa)
1565 1520 : ABI_FREE(bb)
1566 : end if
1567 :
1568 : ! Add extfpmd contribution
1569 83033 : if(present(extfpmd_rho)) then
1570 100087805 : rho1(:,1,:)=rho1(:,1,:)+sqrt(four_pi)*extfpmd_rho
1571 61256 : if(opt_dens<=1) then
1572 100087805 : trho1(:,1,:)=trho1(:,1,:)+sqrt(four_pi)*extfpmd_rho
1573 : endif
1574 : endif
1575 :
1576 : ! -- Test moments of densities and store non-zero ones
1577 83033 : if (nzlmopt==-1) then
1578 65150 : do ll=0,pawtab%lcut_size-1
1579 223966 : do ilm=ll**2+1,(ll+1)**2
1580 158816 : m1=zero;mt1=zero
1581 158816 : if (cplex==1) then
1582 123326233 : m1=maxval(abs(rho1 (1:mesh_size,ilm,ispden)))
1583 123218746 : if (opt_dens<2) mt1=maxval(abs(trho1(1:mesh_size,ilm,ispden)))
1584 : else
1585 3367052 : do ir=1,mesh_size
1586 3362596 : rdum=sqrt(rho1(2*ir-1,ilm,ispden)**2+rho1(2*ir,ilm,ispden)**2)
1587 3367052 : m1=max(m1,rdum)
1588 : end do
1589 4456 : if (opt_dens<2) then
1590 3367052 : do ir=1,mesh_size
1591 3362596 : rdum=sqrt(trho1(2*ir-1,ilm,ispden)**2+trho1(2*ir,ilm,ispden)**2)
1592 3367052 : mt1=max(mt1,rdum)
1593 : end do
1594 : end if
1595 : end if
1596 207908 : if (ispden==1) then
1597 141876 : if ((ilm>1).and.(m1<tol16).and.(mt1<tol16)) then
1598 64711 : lmselectout(ilm)=.false.
1599 : end if
1600 16940 : else if (.not.(lmselectout(ilm))) then
1601 6090 : lmselectout(ilm)=((m1>=tol16).or.(mt1>=tol16))
1602 : end if
1603 : end do
1604 : end do
1605 : end if
1606 :
1607 : ! -- Compute integral of (n1-tn1) inside spheres
1608 83033 : if (opt_compch==1.and.ispden==1.and.opt_dens<2) then
1609 109797 : ABI_MALLOC(aa,(mesh_size))
1610 : aa(1:mesh_size)=(rho1(1:mesh_size,1,1)-trho1(1:mesh_size,1,1)) &
1611 48498472 : & *pawrad%rad(1:mesh_size)**2
1612 36599 : call simp_gen(compchspha(1),aa,pawrad)
1613 36599 : compch_sph=compch_sph+compchspha(1)*sqrt(four_pi)
1614 36599 : ABI_FREE(aa)
1615 : end if
1616 :
1617 : ! -- Print out moments of densities (if requested)
1618 152280 : if (abs(pawprtvol)>=2.and.opt_print==1.and.opt_dens<2) then
1619 15642 : ABI_MALLOC(aa,(cplex*mesh_size))
1620 10428 : ABI_MALLOC(bb,(cplex*mesh_size))
1621 5214 : if (opt_dens==0) then
1622 4288 : write(msg,'(2a,i3,a,i1,3a)') ch10, &
1623 4288 : & ' Atom ',iatom,' (ispden=',ispden,'):',ch10,&
1624 8576 : & ' ******* Moment of (n1-tn1) ** Moment of (n1-tn1-nhat1)'
1625 : else
1626 926 : write(msg,'(2a,i3,a,i1,3a)') ch10, &
1627 926 : & ' Atom ',iatom,' (ispden=',ispden,'):',ch10,&
1628 1852 : & ' ******* Moment of (n1-tn1)'
1629 : end if
1630 5214 : call wrtout(std_out,msg,'PERS')
1631 31544 : do ll=0,pawtab%lcut_size-1
1632 166694 : do ilm=ll**2+1,(ll+1)**2
1633 161480 : if (lmselectin(ilm)) then
1634 200240 : do iplex=1,cplex
1635 200240 : if (opt_dens==0) then
1636 65325578 : do ir=1,mesh_size
1637 65243848 : ii=cplex*(ir-1)+iplex
1638 65243848 : ro(1)=pawrad%rad(ir)**(2+ll)
1639 65243848 : aa(ir)=ro(1)*(rho1(ii,ilm,ispden)-trho1(ii,ilm,ispden))
1640 65325578 : bb(ir)=ro(1)*nhat1(ii,ilm,ispden)
1641 : end do
1642 81730 : call simp_gen(compchspha(iplex),aa,pawrad)
1643 81730 : call simp_gen(compchsphb(iplex),bb,pawrad)
1644 : else
1645 18220270 : do ir=1,mesh_size
1646 18201880 : ii=cplex*(ir-1)+iplex
1647 18201880 : ro(1)=pawrad%rad(ir)**(2+ll)
1648 18220270 : aa(ir)=ro(1)*(rho1(ii,ilm,ispden)-trho1(ii,ilm,ispden))
1649 : end do
1650 18390 : call simp_gen(compchspha(iplex),aa,pawrad)
1651 : end if
1652 : end do
1653 100120 : if (opt_dens==0) then
1654 81730 : if (cplex==1) then
1655 : write(msg,'(3x,a,2i2,2(a,es14.7))') &
1656 81730 : & 'l,m=',ll,ilm-(ll**2+ll+1),': M=',compchspha(1),&
1657 163460 : & ' ** M=',compchspha(1)-compchsphb(1)
1658 : else
1659 : write(msg,'(3x,a,2i2,2(a,2es14.7))') &
1660 0 : & 'l,m=',ll,ilm-(ll**2+ll+1),': M=',compchspha(1:2),&
1661 0 : & ' ** M=',compchspha(1:2)-compchsphb(1:2)
1662 : end if
1663 : else
1664 18390 : if (cplex==1) then
1665 : write(msg,'(3x,a,2i2,a,es14.7)') &
1666 18390 : & 'l,m=',ll,ilm-(ll**2+ll+1),': M=',compchspha(1)
1667 : else
1668 : write(msg,'(3x,a,2i2,a,2es14.7)') &
1669 0 : & 'l,m=',ll,ilm-(ll**2+ll+1),': M=',compchspha(1:2)
1670 : end if
1671 : end if
1672 100120 : call wrtout(std_out,msg,'PERS')
1673 : end if
1674 : end do
1675 : end do
1676 5214 : ABI_FREE(aa)
1677 5214 : ABI_FREE(bb)
1678 : end if
1679 :
1680 : ! ----- End loop over spin components
1681 : end do
1682 :
1683 69247 : if (.not.present(one_over_rad2)) then
1684 13938 : ABI_FREE(one_over_rad2_)
1685 : end if
1686 :
1687 : DBG_EXIT("COLL")
1688 :
1689 138494 : end subroutine pawdensities
1690 : !!***
1691 :
1692 : !----------------------------------------------------------------------
1693 :
1694 : !!****f* m_paw_denpot/pawkindensities
1695 : !! NAME
1696 : !! pawkindensities
1697 : !!
1698 : !! FUNCTION
1699 : !! Compute PAW on-site kinetic energy densities (all-electron, pseudo) for a given atom
1700 : !!
1701 : !! INPUTS
1702 : !! cplex: if 1, on-site densities are REAL, if 2, COMPLEX (response function only)
1703 : !! lm_size=number of (l,m) moments
1704 : !! lmselectin(lm_size)=flags selecting the non-zero LM-moments of on-site kinetic energy densities
1705 : !! (value of these flags at input; must be .TRUE. for nzlmopt/=1)
1706 : !! nspden=number of spin-density components
1707 : !! nzlmopt=if -1, compute all LM-moments of densities (lmselectin=.true. forced)
1708 : !! initialize "lmselectout" (index of non-zero LM-moments of densities)
1709 : !! if 0, compute all LM-moments of densities (lmselectin=.true. forced)
1710 : !! force "lmselectout" to .true. (index of non-zero LM-moments of densities)
1711 : !! if 1, compute only non-zero LM-moments of densities (stored before in "lmselectin")
1712 : !! one_over_rad2(mesh_size)= contains 1/r**2 for each point of the radial grid -optional argument-
1713 : !! opt_dens=flag controlling which on-site kinetic energy density(ies) is (are) computed
1714 : !! 0,1: all-electron and pseudo on-site kinetic energy densities
1715 : !! 2: only all-electron density
1716 : !! opt_l=controls which l-moment(s) contribute to the kinetic energy density:
1717 : !! <0 : all l contribute
1718 : !! >=0: only l=opt_l contributes
1719 : !! Note: opt_l>=0 is only compatible with opt_dens=2
1720 : !! pawang <type(pawang_type)>=paw angular mesh and related data
1721 : !! pawrad <type(pawrad_type)>=paw radial mesh and related data (for the current atom type)
1722 : !! pawrhoij <type(pawrhoij_type)>= paw rhoij occupancies and related data (for the current atom)
1723 : !! pawtab <type(pawtab_type)>=paw tabulated starting data (for the current atom type)
1724 : !!
1725 : !! OUTPUT
1726 : !! tau1(cplex*mesh_size,lm_size,nspden)= on site kinetic energy density
1727 : !! ttau1(cplex*mesh_size,lm_size,nspden)]= pseudo on site kinetic energy density
1728 : !!
1729 : !! SOURCE
1730 :
1731 484 : subroutine pawkindensities(cplex,lmselectin,lm_size,nspden,nzlmopt,&
1732 968 : & opt_dens,opt_l,pawang,pawrad,pawrhoij,pawtab,tau1,ttau1,&
1733 : & one_over_rad2) ! optional
1734 :
1735 : !Arguments ---------------------------------------------
1736 : !scalars
1737 : integer,intent(in) :: cplex,lm_size,nspden,nzlmopt,opt_dens,opt_l
1738 : type(pawang_type),intent(in) :: pawang
1739 : type(pawrad_type),intent(in) :: pawrad
1740 : type(pawrhoij_type),intent(in) :: pawrhoij
1741 : type(pawtab_type),intent(in) :: pawtab
1742 : !arrays
1743 : logical,intent(in) :: lmselectin(lm_size)
1744 : real(dp),intent(in),target,optional :: one_over_rad2(pawtab%mesh_size)
1745 : real(dp),intent(out),optional :: tau1(cplex*pawtab%mesh_size,lm_size,nspden*(1-(opt_dens/2)))
1746 : real(dp),intent(out),optional :: ttau1(cplex*pawtab%mesh_size,lm_size,nspden*(1-(opt_dens/2)))
1747 : !Local variables ---------------------------------------
1748 : !scalars
1749 : integer :: dplex,ii,iq0,ir,irhoij,isel,ispden,jrhoij
1750 : integer :: ilmn,ilm,ilm1,iln,jlmn,jlm1,jln,klm,klmn,ll,lmax,lmin,mesh_size
1751 : real(dp) :: phiphj,tphitphj
1752 : character(len=500) :: msg
1753 : !arrays
1754 968 : real(dp) :: ro(cplex),ro_rg(cplex)
1755 484 : real(dp),allocatable :: aa(:),bb(:)
1756 484 : real(dp),pointer :: one_over_rad2_(:)
1757 :
1758 : ! *************************************************************************
1759 :
1760 : DBG_ENTER("COLL")
1761 :
1762 : !Compatibility tests
1763 484 : if (nzlmopt/=-1) then
1764 0 : msg='nzlmopt/=-1 has not not been tested (might be wrong)!'
1765 0 : ABI_BUG(msg)
1766 : end if
1767 484 : if (opt_dens/=2.and.opt_l>=0) then
1768 0 : msg='opt_dens/=2 incompatible with opt_l>=0!'
1769 0 : ABI_BUG(msg)
1770 : end if
1771 484 : if(nzlmopt/=0.and.nzlmopt/=1.and.nzlmopt/=-1) then
1772 0 : msg='invalid value for variable "nzlmopt"!'
1773 0 : ABI_BUG(msg)
1774 : end if
1775 484 : if(nspden>pawrhoij%nspden) then
1776 0 : msg='nspden must be <= pawrhoij%nspden!'
1777 0 : ABI_BUG(msg)
1778 : end if
1779 484 : if (cplex>pawrhoij%qphase) then
1780 0 : msg='cplex must be <= pawrhoij%qphase!'
1781 0 : ABI_BUG(msg)
1782 : end if
1783 484 : if (nzlmopt/=1) then
1784 9960 : if (any(.not.lmselectin(1:lm_size))) then
1785 0 : msg='With nzlmopt/=1, lmselectin must be true!'
1786 0 : ABI_BUG(msg)
1787 : end if
1788 : end if
1789 484 : if (pawang%gnt_option==0) then
1790 0 : msg='pawang%gnt_option=0!'
1791 0 : ABI_BUG(msg)
1792 : end if
1793 484 : if (pawang%nabgnt_option==0) then
1794 0 : msg='pawang%nabgnt_option=0!'
1795 0 : ABI_BUG(msg)
1796 : end if
1797 :
1798 : !Various inits
1799 13701850 : tau1=zero
1800 13702334 : if (opt_dens<2) ttau1=zero
1801 484 : mesh_size=pawtab%mesh_size;dplex=cplex-1
1802 484 : iq0=pawrhoij%cplex_rhoij*pawrhoij%lmn2_size
1803 484 : if (present(one_over_rad2)) then
1804 484 : one_over_rad2_ => one_over_rad2
1805 : else
1806 0 : ABI_MALLOC(one_over_rad2_,(mesh_size))
1807 0 : one_over_rad2_(1)=zero
1808 0 : one_over_rad2_(2:mesh_size)=one/pawrad%rad(2:mesh_size)**2
1809 : end if
1810 :
1811 : !=== Compute "on-site" kin. energy densities (n1, ntild1) =====
1812 : !==============================================================
1813 :
1814 968 : do ispden=1,nspden
1815 :
1816 : ! -- Loop over ij channels (basis components)
1817 484 : jrhoij=1
1818 48594 : do irhoij=1,pawrhoij%nrhoijsel
1819 48110 : klmn=pawrhoij%rhoijselect(irhoij)
1820 48110 : klm =pawtab%indklmn(1,klmn)
1821 48110 : lmin=pawtab%indklmn(3,klmn)
1822 48110 : lmax=pawtab%indklmn(4,klmn)
1823 48110 : ilmn=pawtab%indklmn(7,klmn) ! (l,m,n) orbital 1
1824 48110 : jlmn=pawtab%indklmn(8,klmn) ! (l,m,n) orbital 2
1825 48110 : ilm1=pawtab%indklmn(5,klmn) ! (l,m) orbital 1
1826 48110 : jlm1=pawtab%indklmn(6,klmn) ! (l,m) orbital 2
1827 48110 : iln=pawtab%indlmn(5,ilmn) ! (l,n) orbital 1
1828 48110 : jln=pawtab%indlmn(5,jlmn) ! (l,n) orbital 2
1829 :
1830 : ! Retrieve rhoij
1831 48110 : if (pawrhoij%nspden/=2) then
1832 48110 : ro(1)=pawrhoij%rhoijp(jrhoij,ispden)
1833 48110 : if (cplex==2) ro(2)=pawrhoij%rhoijp(iq0+jrhoij,ispden)
1834 : else
1835 0 : if (ispden==1) then
1836 0 : ro(1)=pawrhoij%rhoijp(jrhoij,1)+pawrhoij%rhoijp(jrhoij,2)
1837 0 : if (cplex==2) ro(2)=pawrhoij%rhoijp(iq0+jrhoij,1)+pawrhoij%rhoijp(iq0+jrhoij,2)
1838 0 : else if (ispden==2) then
1839 0 : ro(1)=pawrhoij%rhoijp(jrhoij,1)
1840 0 : if (cplex==2) ro(2)=pawrhoij%rhoijp(iq0+jrhoij,1)
1841 : end if
1842 : end if
1843 : ! Apply factor 1/2 (because tau=1/2 * Sum_ij[rhoij.Nabla_phi_i*Nabla_phi_j])
1844 96220 : ro(1:cplex)=half*pawtab%dltij(klmn)*ro(1:cplex)
1845 :
1846 : ! First option: AE and PS on-site kin. energy densities (opt_dens==0 or 1)
1847 : ! ------------------------------------------------------------------------
1848 48110 : if (opt_dens==0.or.opt_dens==1) then
1849 :
1850 : ! Compute part of tau_lm depending on gaunt coefficients
1851 149425 : do ll=lmin,lmax,2
1852 609178 : do ilm=ll**2+1,(ll+1)**2
1853 561068 : if (lmselectin(ilm)) then
1854 459753 : isel=pawang%gntselect(ilm,klm)
1855 459753 : if (isel>0) then
1856 208628 : ro_rg(1:cplex)=ro(1:cplex)*pawang%realgnt(isel)
1857 151430465 : do ir=2,mesh_size
1858 151326151 : phiphj=pawtab%nablaphi(ir,iln)*pawtab%nablaphi(ir,jln)
1859 151326151 : tphitphj=pawtab%tnablaphi(ir,iln)*pawtab%tnablaphi(ir,jln)
1860 : tau1(cplex*ir-dplex:ir*cplex,ilm,ispden)=tau1(cplex*ir-dplex:ir*cplex,ilm,ispden)&
1861 302652302 : & +ro_rg(1:cplex)*phiphj*one_over_rad2_(ir)
1862 : ttau1(cplex*ir-dplex:ir*cplex,ilm,ispden)=ttau1(cplex*ir-dplex:ir*cplex,ilm,ispden)&
1863 302756616 : & +ro_rg(1:cplex)*tphitphj*one_over_rad2_(ir)
1864 : end do
1865 : end if
1866 : end if
1867 : end do ! End loops over ll,lm
1868 : end do
1869 :
1870 : ! Compute the part of tau_lm depending on nablagaunt coefficients
1871 792359 : do ilm=1,(1+lmax)**2
1872 792359 : if (lmselectin(ilm)) then
1873 744249 : isel=pawang%nablagntselect(ilm,ilm1,jlm1)
1874 744249 : if (isel>0) then
1875 190586 : ro_rg(1:cplex)=ro(1:cplex)*pawang%nablarealgnt(isel)
1876 138367015 : do ir=2,mesh_size
1877 138271722 : phiphj=pawtab%phi(ir,iln)*pawtab%phi(ir,jln)
1878 138271722 : tphitphj=pawtab%tphi(ir,iln)*pawtab%tphi(ir,jln)
1879 : tau1(cplex*ir-dplex:ir*cplex,ilm,ispden)=tau1(cplex*ir-dplex:ir*cplex,ilm,ispden)&
1880 276543444 : & +ro_rg(1:cplex)*phiphj*one_over_rad2_(ir)**2
1881 : ttau1(cplex*ir-dplex:ir*cplex,ilm,ispden)=ttau1(cplex*ir-dplex:ir*cplex,ilm,ispden)&
1882 276638737 : & +ro_rg(1:cplex)*tphitphj*one_over_rad2_(ir)**2
1883 : end do
1884 : end if
1885 : end if
1886 : end do
1887 :
1888 : ! 2nd option: AE on-site kinetic energy density only (opt_dens==2)
1889 : ! ----------------------------------------------------------------
1890 0 : else if (opt_dens==2) then
1891 :
1892 : ! Compute part of tau_lm depending on gaunt coefficients
1893 0 : do ll=lmin,lmax,2
1894 0 : do ilm=ll**2+1,(ll+1)**2
1895 0 : if (lmselectin(ilm)) then
1896 0 : isel=pawang%gntselect(ilm,klm)
1897 0 : if (isel>0) then
1898 0 : ro_rg(1:cplex)=ro(1:cplex)*pawang%realgnt(isel)
1899 0 : do ir=2,mesh_size
1900 0 : phiphj=pawtab%nablaphi(ir,iln)*pawtab%nablaphi(ir,jln)
1901 : tau1(cplex*ir-dplex:ir*cplex,ilm,ispden)=tau1(cplex*ir-dplex:ir*cplex,ilm,ispden)&
1902 0 : & +ro_rg(1:cplex)*phiphj*one_over_rad2_(ir)
1903 : end do
1904 : end if
1905 : end if
1906 : end do ! End loops over ll,lm
1907 : end do
1908 :
1909 : ! Compute the part of tau_lm depending on nablagaunt coefficients
1910 0 : do ilm=1,(1+lmax)**2
1911 0 : if (lmselectin(ilm)) then
1912 0 : isel=pawang%nablagntselect(ilm,ilm1,jlm1)
1913 0 : if (isel>0) then
1914 0 : ro_rg(1:cplex)=ro(1:cplex)*pawang%nablarealgnt(isel)
1915 0 : do ir=2,mesh_size
1916 0 : phiphj=pawtab%phi(ir,iln)*pawtab%phi(ir,jln)
1917 : tau1(cplex*ir-dplex:ir*cplex,ilm,ispden)=tau1(cplex*ir-dplex:ir*cplex,ilm,ispden)&
1918 0 : & +ro_rg(1:cplex)*phiphj*one_over_rad2_(ir)**2
1919 : end do
1920 : end if
1921 : end if
1922 : end do
1923 :
1924 : end if
1925 :
1926 : ! -- End loop over ij channels
1927 48594 : jrhoij=jrhoij+pawrhoij%cplex_rhoij
1928 : end do
1929 :
1930 : ! Compute tau1(r=0) and ttau1(r=0)
1931 484 : if (cplex==2) then
1932 0 : ABI_MALLOC(aa,(5))
1933 0 : ABI_MALLOC(bb,(5))
1934 : end if
1935 484 : if (opt_dens==0.or.opt_dens==1) then
1936 2576 : do ll=0,pawtab%lcut_size-1
1937 12052 : do ilm=ll**2+1,(ll+1)**2
1938 11568 : if (lmselectin(ilm)) then
1939 9476 : if (cplex==1) then
1940 9476 : call pawrad_deducer0(tau1 (:,ilm,ispden),mesh_size,pawrad)
1941 9476 : call pawrad_deducer0(ttau1(:,ilm,ispden),mesh_size,pawrad)
1942 : else
1943 0 : do ii=0,1
1944 0 : do ir=2,5
1945 0 : aa(ir)=tau1 (2*ir-ii,ilm,ispden)
1946 0 : bb(ir)=ttau1(2*ir-ii,ilm,ispden)
1947 : end do
1948 0 : call pawrad_deducer0(aa,5,pawrad)
1949 0 : call pawrad_deducer0(bb,5,pawrad)
1950 0 : tau1 (2-ii,ilm,ispden)=aa(1)
1951 0 : ttau1(2-ii,ilm,ispden)=bb(1)
1952 : end do
1953 : end if
1954 : end if
1955 : end do
1956 : end do
1957 0 : else if (opt_dens==2) then
1958 0 : do ll=0,pawtab%lcut_size-1
1959 0 : do ilm=ll**2+1,(ll+1)**2
1960 0 : if (lmselectin(ilm)) then
1961 0 : if (cplex==1) then
1962 0 : call pawrad_deducer0(tau1(:,ilm,ispden),mesh_size,pawrad)
1963 : else
1964 0 : do ii=0,1
1965 0 : do ir=2,5
1966 0 : aa(ir)=tau1(2*ir-ii,ilm,ispden)
1967 : end do
1968 0 : call pawrad_deducer0(aa,5,pawrad)
1969 0 : tau1(2-ii,ilm,ispden)=aa(1)
1970 : end do
1971 : end if
1972 : end if
1973 : end do
1974 : end do
1975 : end if
1976 968 : if (cplex==2) then
1977 0 : ABI_FREE(aa)
1978 0 : ABI_FREE(bb)
1979 : end if
1980 :
1981 : ! ----- End loop over spin components
1982 : end do
1983 :
1984 484 : if (.not.present(one_over_rad2)) then
1985 0 : ABI_FREE(one_over_rad2_)
1986 : end if
1987 :
1988 : DBG_EXIT("COLL")
1989 :
1990 1936 : end subroutine pawkindensities
1991 : !!***
1992 :
1993 : !----------------------------------------------------------------------
1994 :
1995 : !!****f* m_paw_denpot/pawaccenergy
1996 : !! NAME
1997 : !! pawaccenergy
1998 : !!
1999 : !! FUNCTION
2000 : !! Accumulate an on-site PAW energy by adding the contribution of the current atom.
2001 : !! This contribution has the form: Sum_ij[Rhoij.Dij]
2002 : !!
2003 : !! INPUTS
2004 : !! pawrhoij<type(pawrhoij_type)>= datastructure containing Rho_ij values
2005 : !! dij(cplex_dij*qphase_dij*lmn2_size,nspden_dij)= array containing D_ij values
2006 : !! cplex_dij= 2 if dij is COMPLEX (as in the spin-orbit case), 1 if dij is REAL
2007 : !! qphase_dij= 2 if dij has a exp(iqR) phase, 1 if not
2008 : !! nspden_dij= number of spin components for dij
2009 : !! pawtab<type(pawtab_type)>=paw tabulated starting data
2010 : !!
2011 : !! OUTPUT
2012 : !!
2013 : !! SIDE EFFECTS
2014 : !! epaw= PAW on-site energy. At output, the contribution of the current atom
2015 : !! has been added to epaw.
2016 : !! [epaw_im]= imaginary part of PAW on-site energy. At output, the contribution
2017 : !! of the current atom has been added to epaw.
2018 : !! This imaginary p rt only exists in a few cases (f.i. non-stationnary
2019 : !! expression of 2nd-order energy)
2020 : !!
2021 : !! NOTES
2022 : !! * The general form for Dij is:
2023 : !! D^{s1,s2}_ij = D1^{s1,s2}_ij.cos(qr) + i.D2^{s1,s2}_ij.sin(qr)
2024 : !! = [D1re^{s1,s2}_ij + i.D1im^{s1,s2}_ij).cos(qr)]
2025 : !! + i.[D2re^{s1,s2}_ij + i.D2im^{s1,s2}_ij).sin(qr)]
2026 : !! where
2027 : !! ij are the partial waves channels
2028 : !! s1,s2 are spin/spinor components
2029 : !! q is the wave vector of the phase
2030 : !! D1^{s1,s2}_ij.cos(qr) is stored in the the first half of paw_ij%dij and corresponds to iq=1
2031 : !! D2^{s1,s2}_ij.sin(qr) is stored in the the 2nd half of paw_ij%dij and corresponds to iq=2
2032 : !! D1^{s1,s2}_ij.cos(qr) and D2^{s1,s2}_ij.sin(qr) are complex if cplex_dij=2
2033 : !!
2034 : !! * The same for Rho_ij
2035 : !!
2036 : !! * The contribution to the PAW on-site energy is:
2037 : !! Sum_ij_s1s2[Rho^{s2,s1}_ij * D^{s1,s2}_ij]
2038 : !! Note the order of s1/s2 indices, especially for Rho_ij.
2039 : !! The present implementation follows eq(15) in Hobbs et al, PRB 62, 11556(2000)
2040 : !! rho^{s1,s2}^_ij = Sum[<Psi^s2|pi><pj|Psi^s1] (s1 and s2 exponents inverted)
2041 : !!
2042 : !! SOURCE
2043 :
2044 139863 : subroutine pawaccenergy(epaw,pawrhoij,dij,cplex_dij,qphase_dij,nspden_dij,pawtab,epaw_im)
2045 :
2046 : !Arguments ---------------------------------------------
2047 : !scalars
2048 : integer,intent(in) :: cplex_dij,qphase_dij,nspden_dij
2049 : real(dp),intent(inout) :: epaw
2050 : real(dp),intent(inout),optional :: epaw_im
2051 : type(pawrhoij_type),intent(in),target :: pawrhoij
2052 : type(pawtab_type),intent(in) :: pawtab
2053 : !arrays
2054 : real(dp),intent(in) :: dij(cplex_dij*qphase_dij*pawtab%lmn2_size,nspden_dij)
2055 :
2056 : !Local variables ---------------------------------------
2057 : !scalars
2058 : integer :: cplex_rhoij,iq,iq0_dij,iq0_rhoij,irhoij,isp_dij,isp_rhoij,jrhoij
2059 : integer :: klmn,kklmn,krhoij,lmn2_size,nspden_rhoij,nsploop
2060 : logical :: add_imaginary
2061 : real(dp) :: etmp
2062 : character(len=500) :: msg
2063 : !arrays
2064 139863 : real(dp),pointer :: rhoij(:,:)
2065 :
2066 : ! *************************************************************************
2067 :
2068 : DBG_ENTER("COLL")
2069 :
2070 : !Compatibility tests
2071 139863 : if (pawrhoij%qphase/=qphase_dij) then
2072 0 : msg='pawaccenergy: pawrhoij%qphase/=qphase_dij!'
2073 0 : ABI_BUG(msg)
2074 : end if
2075 139863 : if (pawrhoij%nspden>nspden_dij.and.nspden_dij/=1) then
2076 0 : msg='pawaccenergy: pawrhoij%nspden>nspden_dij!'
2077 0 : ABI_BUG(msg)
2078 : end if
2079 :
2080 : !Useful data
2081 139863 : nspden_rhoij=pawrhoij%nspden
2082 139863 : lmn2_size=pawtab%lmn2_size
2083 :
2084 : !Special treatment for nspden
2085 139863 : nsploop=nspden_rhoij
2086 139863 : if (nspden_dij==1.and.nspden_rhoij==4) nsploop=1
2087 :
2088 : !Non-collinear case: need a temporary rhoij
2089 11421 : if (nspden_rhoij==4.and.nspden_dij==4) then
2090 2322 : cplex_rhoij=2
2091 9288 : ABI_MALLOC(rhoij,(2*lmn2_size,4))
2092 : else
2093 137541 : cplex_rhoij=pawrhoij%cplex_rhoij
2094 137541 : rhoij => pawrhoij%rhoijp
2095 : end if
2096 :
2097 139863 : add_imaginary=(cplex_dij==2.and.cplex_rhoij==2)
2098 :
2099 : !Loop over qphase components
2100 281154 : do iq=1,qphase_dij
2101 141291 : iq0_rhoij=(iq-1)*lmn2_size*cplex_rhoij
2102 141291 : iq0_dij =(iq-1)*lmn2_size*cplex_dij
2103 :
2104 : ! Non-collinear case
2105 141291 : if (nspden_rhoij==4.and.nspden_dij==4) then
2106 2002562 : rhoij(:,:)=zero
2107 2322 : jrhoij=(iq-1)*lmn2_size*pawrhoij%cplex_rhoij+1 ; krhoij=1
2108 228708 : do irhoij=1,pawrhoij%nrhoijsel
2109 226386 : klmn=pawrhoij%rhoijselect(irhoij)
2110 226386 : rhoij(krhoij ,1)= half*(pawrhoij%rhoijp(jrhoij,1)+pawrhoij%rhoijp(jrhoij,4))
2111 226386 : rhoij(krhoij ,2)= half*(pawrhoij%rhoijp(jrhoij,1)-pawrhoij%rhoijp(jrhoij,4))
2112 : !Be careful we store rhoij^21 in rhoij(:,3) and rhoij^12 in rhoij(:,4)
2113 : !because of the inversion of spins in rhoij definition
2114 226386 : rhoij(krhoij ,3)= half*pawrhoij%rhoijp(jrhoij,2)
2115 226386 : rhoij(krhoij+1,3)= half*pawrhoij%rhoijp(jrhoij,3)
2116 226386 : rhoij(krhoij ,4)= half*pawrhoij%rhoijp(jrhoij,2)
2117 226386 : rhoij(krhoij+1,4)=-half*pawrhoij%rhoijp(jrhoij,3)
2118 226386 : if (pawrhoij%cplex_rhoij==2) then
2119 226386 : rhoij(krhoij+1,1)= half*(pawrhoij%rhoijp(jrhoij+1,1)+pawrhoij%rhoijp(jrhoij+1,4))
2120 226386 : rhoij(krhoij+1,2)= half*(pawrhoij%rhoijp(jrhoij+1,1)-pawrhoij%rhoijp(jrhoij+1,4))
2121 : !Be careful we store rhoij^21 in rhoij(:,3) and rhoij^12 in rhoij(:,4)
2122 : !because of the inversion of spins in rhoij definition
2123 226386 : rhoij(krhoij ,3)= rhoij(krhoij ,3)-half*pawrhoij%rhoijp(jrhoij+1,3)
2124 226386 : rhoij(krhoij+1,3)= rhoij(krhoij+1,3)+half*pawrhoij%rhoijp(jrhoij+1,2)
2125 226386 : rhoij(krhoij ,4)= rhoij(krhoij ,4)+half*pawrhoij%rhoijp(jrhoij+1,3)
2126 226386 : rhoij(krhoij+1,4)= rhoij(krhoij+1,4)+half*pawrhoij%rhoijp(jrhoij+1,2)
2127 : end if
2128 228708 : jrhoij=jrhoij+pawrhoij%cplex_rhoij ; krhoij=krhoij+2
2129 : end do
2130 : iq0_rhoij=0
2131 : end if
2132 :
2133 : ! Contribution to on-site energy (real part)
2134 306346 : do isp_rhoij=1,nsploop
2135 165055 : isp_dij=min(isp_rhoij,nspden_dij)
2136 165055 : jrhoij=iq0_rhoij+1
2137 6248065 : do irhoij=1,pawrhoij%nrhoijsel
2138 5941719 : klmn=pawrhoij%rhoijselect(irhoij)
2139 5941719 : kklmn=iq0_dij+cplex_dij*(klmn-1)+1
2140 5941719 : etmp=rhoij(jrhoij,isp_rhoij)*dij(kklmn,isp_dij)
2141 5941719 : if (add_imaginary) etmp=etmp-rhoij(jrhoij+1,isp_rhoij)*dij(kklmn+1,isp_dij)
2142 5941719 : epaw=epaw+etmp*pawtab%dltij(klmn)
2143 6106774 : jrhoij=jrhoij+cplex_rhoij
2144 : end do
2145 : end do ! nsploop
2146 :
2147 : ! Contribution to on-site energy (imaginary part)
2148 281154 : if (present(epaw_im).and.qphase_dij==2) then
2149 1224 : do isp_rhoij=1,nsploop
2150 624 : isp_dij=min(isp_rhoij,nspden_dij)
2151 624 : jrhoij=iq0_rhoij+1
2152 26068 : do irhoij=1,pawrhoij%nrhoijsel
2153 24844 : klmn=pawrhoij%rhoijselect(irhoij)
2154 24844 : if (iq==1) then
2155 12422 : kklmn=lmn2_size*cplex_dij+cplex_dij*(klmn-1)+1
2156 12422 : etmp=-rhoij(jrhoij,isp_rhoij)*dij(kklmn,isp_dij)
2157 12422 : if (add_imaginary) etmp=etmp+rhoij(jrhoij+1,isp_rhoij)*dij(kklmn+1,isp_dij)
2158 : end if
2159 24844 : if (iq==2) then
2160 12422 : kklmn=cplex_dij*(klmn-1)+1
2161 12422 : etmp=rhoij(jrhoij,isp_rhoij)*dij(kklmn,isp_dij)
2162 12422 : if (add_imaginary) etmp=etmp-rhoij(jrhoij+1,isp_rhoij)*dij(kklmn+1,isp_dij)
2163 : end if
2164 24844 : epaw_im=epaw_im+etmp*pawtab%dltij(klmn)
2165 25468 : jrhoij=jrhoij+cplex_rhoij
2166 : end do
2167 : end do ! nsploop
2168 : end if
2169 :
2170 : end do ! qphase
2171 :
2172 139863 : if (nspden_rhoij==4.and.nspden_dij==4) then
2173 2322 : ABI_FREE(rhoij)
2174 : end if
2175 :
2176 : DBG_EXIT("COLL")
2177 :
2178 139863 : end subroutine pawaccenergy
2179 : !!***
2180 :
2181 : !----------------------------------------------------------------------
2182 :
2183 : !!****f* m_paw_denpot/pawaccenergy_nospin
2184 : !! NAME
2185 : !! pawaccenergy_nospin
2186 : !!
2187 : !! FUNCTION
2188 : !! Accumulate an on-site PAW energy by adding the contribution of the current atom.
2189 : !! This contribution has the form: Sum_ij[Rhoij.Dij]
2190 : !! Applies only for Dij without spin components (as f.i. Dij^Hartree).
2191 : !! This routine is a wrapper to pawaccenergy.
2192 : !!
2193 : !! INPUTS
2194 : !! pawrhoij<type(pawrhoij_type)>= datastructure containing Rho_ij values
2195 : !! dij(cplex_dij*qphase_dij*lmn2_size)= array containing D_ij values
2196 : !! cplex_dij= 2 if dij is COMPLEX (as in the spin-orbit case), 1 if dij is REAL
2197 : !! qphase_dij= 2 if dij has a exp(iqR) phase, 1 if not
2198 : !! pawtab<type(pawtab_type)>=paw tabulated starting data
2199 : !!
2200 : !! OUTPUT
2201 : !!
2202 : !! SIDE EFFECTS
2203 : !! epaw= PAW on-site energy. At output, the contribution of the current atom
2204 : !! has been added to epaw.
2205 : !! [epaw_im]= imaginary part of PAW on-site energy. At output, the contribution
2206 : !! of the current atom has been added to epaw.
2207 : !! This imaginary part only exists in a few cases (f.i. non-stationnary
2208 : !! expression of 2nd-order energy)
2209 : !!
2210 : !! SOURCE
2211 :
2212 136455 : subroutine pawaccenergy_nospin(epaw,pawrhoij,dij,cplex_dij,qphase_dij,pawtab,epaw_im)
2213 :
2214 : !Arguments ---------------------------------------------
2215 : !scalars
2216 : integer,intent(in) :: cplex_dij,qphase_dij
2217 : real(dp),intent(inout) :: epaw
2218 : real(dp),intent(inout),optional :: epaw_im
2219 : type(pawrhoij_type),intent(in),target :: pawrhoij
2220 : type(pawtab_type),intent(in) :: pawtab
2221 : !arrays
2222 : real(dp),intent(in),target :: dij(cplex_dij*qphase_dij*pawtab%lmn2_size)
2223 :
2224 : !Local variables ---------------------------------------
2225 : !scalars
2226 : integer :: size_dij
2227 : #ifdef HAVE_FC_ISO_C_BINDING
2228 : type(C_PTR) :: cptr
2229 : #endif
2230 : !arrays
2231 136455 : real(dp), ABI_CONTIGUOUS pointer :: dij_2D(:,:)
2232 :
2233 : ! *************************************************************************
2234 :
2235 136455 : size_dij=size(dij)
2236 :
2237 : #ifdef HAVE_FC_ISO_C_BINDING
2238 136455 : cptr=c_loc(dij(1))
2239 136455 : call c_f_pointer(cptr,dij_2D,shape=[size_dij,1])
2240 : #else
2241 : ABI_MALLOC(dij_2D,(size_dij,1))
2242 : dij_2D=reshape(dij,[size_dij,1])
2243 : #endif
2244 :
2245 136455 : if (present(epaw_im)) then
2246 9108 : call pawaccenergy(epaw,pawrhoij,dij_2D,cplex_dij,qphase_dij,1,pawtab,epaw_im=epaw_im)
2247 : else
2248 127347 : call pawaccenergy(epaw,pawrhoij,dij_2D,cplex_dij,qphase_dij,1,pawtab)
2249 : end if
2250 :
2251 : #ifndef HAVE_FC_ISO_C_BINDING
2252 : ABI_FREE(dij_2D)
2253 : #endif
2254 :
2255 136455 : end subroutine pawaccenergy_nospin
2256 : !!***
2257 :
2258 : !----------------------------------------------------------------------
2259 :
2260 : !!****f* m_paw_denpot/paw_mknewh0
2261 : !! NAME
2262 : !! paw_mknewh0
2263 : !!
2264 : !! FUNCTION
2265 : !! Calculates the new bare PAW Hamiltonian in the case of quasi-particle self-consistent GW calculations.
2266 : !!
2267 : !! INPUTS
2268 : !! mpi_atmtab(:)=--optional-- indexes of the atoms treated by current proc
2269 : !! comm_atom=--optional-- MPI communicator over atoms
2270 : !! my_natom=number of atoms treated by current processor
2271 : !! nsppol=1 for unpolarized, 2 for spin-polarized
2272 : !! nspden=number of spin-density components
2273 : !! nfftf=(effective) number of FFT grid points (for this proc) for the "fine" grid
2274 : !! pawspnorb=flag: 1 if spin-orbit coupling is activated
2275 : !! pawprtvol=control print volume and debugging output for PAW
2276 : !! Cryst<crystal_t>=Info on unit cell and its symmetries
2277 : !! Pawtab(ntypat*usepaw)<type(pawtab_type)>=paw tabulated starting data
2278 : !! Paw_an(natom) <type(paw_an_type)>=paw arrays given on angular mesh
2279 : !! Pawang<type(pawang_type)>=paw angular mesh and related data
2280 : !! Pawfgrtab(natom) <type(pawfgrtab_type)>=atomic data given on fine rectangular grid
2281 : !! vxc(nfftf,nspden)=exchange-correlation potential
2282 : !! vxc_val(nfftf,nspden)=valence only exchange-correlation potential
2283 : !! vtrial(nfftf,nspden)=potential (Hartree+XC+loc)
2284 : !!
2285 : !! SIDE EFFECTS
2286 : !! Paw_ij(natom*usepaw)<Paw_ij_type)>=paw arrays given on (i,j) channels
2287 : !! At output: new value for Paw_ij()%dij
2288 : !!
2289 : !! SOURCE
2290 :
2291 0 : subroutine paw_mknewh0(my_natom,nsppol,nspden,nfftf,pawspnorb,pawprtvol,Cryst,&
2292 0 : & Pawtab,Paw_an,Paw_ij,Pawang,Pawfgrtab,vxc,vxc_val,vtrial,&
2293 0 : & mpi_atmtab,comm_atom) ! optional arguments (parallelism)
2294 :
2295 : !Arguments ------------------------------------
2296 : !scalars
2297 : integer,intent(in) :: my_natom,nsppol,nspden,nfftf,pawprtvol,pawspnorb
2298 : integer,optional,intent(in) :: comm_atom
2299 : !arrays
2300 : integer,optional,target,intent(in) :: mpi_atmtab(:)
2301 : real(dp),intent(in) :: vxc(nfftf,nspden),vxc_val(nfftf,nspden),vtrial(nfftf,nspden)
2302 : type(crystal_t),intent(in) :: Cryst
2303 : type(Pawang_type),intent(in) :: Pawang
2304 : type(Pawtab_type),target,intent(in) :: Pawtab(Cryst%ntypat)
2305 : type(Paw_an_type),intent(in) :: Paw_an(my_natom)
2306 : type(Paw_ij_type),intent(inout) :: Paw_ij(my_natom)
2307 : type(Pawfgrtab_type),intent(inout) :: Pawfgrtab(my_natom)
2308 :
2309 : !Local variables-------------------------------
2310 : !scalars
2311 : integer,parameter :: ipert0=0
2312 : integer :: iat,iat_tot,idij,ndij,option_dij
2313 : integer :: itypat,lmn_size,j0lmn,jlmn,ilmn,klmn,klmn1,klm
2314 : integer :: lmin,lmax,mm,isel,lm_size,lmn2_size,my_comm_atom,cplex_dij
2315 : integer :: ils,ilslm,ic,lm0
2316 : integer :: nsploop,is2fft,qphase
2317 : real(dp) :: gylm,qijl
2318 : logical :: ltest,my_atmtab_allocated,paral_atom
2319 : character(len=500) :: msg
2320 : !arrays
2321 0 : integer, pointer :: indklmn_(:,:)
2322 0 : integer,pointer :: my_atmtab(:)
2323 : real(dp) :: rdum(1),rdum2(1)
2324 0 : real(dp),allocatable :: prod_hloc(:,:),prodhxc_core(:,:)
2325 0 : real(dp),allocatable :: dijhl_hat(:,:),dijhmxc_val(:,:)
2326 :
2327 : ! *************************************************************************
2328 :
2329 : DBG_ENTER("COLL")
2330 :
2331 0 : call wrtout(std_out,'Assembling PAW strengths for the bare Hamiltonian','COLL')
2332 :
2333 : !== Set up parallelism over atoms ===
2334 0 : paral_atom=(present(comm_atom).and.(my_natom/=Cryst%natom))
2335 0 : nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab
2336 0 : my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom
2337 0 : call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,Cryst%natom,my_natom_ref=my_natom)
2338 :
2339 0 : if (my_natom>0) then
2340 :
2341 : ! === Test if required pointers in paw_ij are allocated ===
2342 0 : ltest = (allocated(Paw_ij(1)%dijxc).and.allocated(Paw_ij(1)%dijxc_val) )
2343 0 : ABI_CHECK(ltest,'dijxc or dijxc_val not calculated')
2344 :
2345 0 : ltest=(allocated(Paw_ij(1)%dijhat)) !.and.Paw_ij(1)%has_dijhat==2)
2346 0 : ABI_CHECK(ltest,'dijhat not calculated')
2347 :
2348 0 : ltest=(allocated(Paw_ij(1)%dijhartree)) !.and.Paw_ij(1)%has_dijhartree==2)
2349 0 : ABI_CHECK(ltest,'dijhartree not calculated')
2350 :
2351 0 : if (ANY(Pawtab(:)%usepawu/=0)) then
2352 0 : do iat=1,my_natom
2353 0 : iat_tot=iat;if (paral_atom) iat_tot=my_atmtab(iat)
2354 0 : itypat=Cryst%typat(iat_tot)
2355 0 : if (Pawtab(itypat)%usepawu/=0) then
2356 0 : ltest=(allocated(Paw_ij(iat)%dijU) ) !.and.Paw_ij(iat)%has_dijU==2)
2357 0 : write(msg,'(a,i3,a)')" For atom no. ",iat," %dijU(iat) has not been calculated."
2358 0 : ABI_CHECK(ltest,msg)
2359 : end if
2360 : end do
2361 : end if
2362 :
2363 0 : if (pawspnorb>0) then
2364 0 : do iat=1,my_natom
2365 0 : ltest=(allocated(Paw_ij(iat)%dijso) ) !.and.Paw_ij(iat)%has_dijso==2)
2366 0 : write(msg,'(a,i3,a)')" For atom no. ",iat," %dijso(iat) has not been calculated."
2367 0 : ABI_CHECK(ltest,msg)
2368 : end do
2369 : end if
2370 : end if ! my_natom>0
2371 :
2372 : !== Construct the new PAW H0 Hamiltonian ===
2373 0 : do iat=1,my_natom
2374 0 : iat_tot=iat;if (paral_atom) iat_tot=my_atmtab(iat)
2375 :
2376 0 : itypat = Cryst%typat(iat_tot)
2377 0 : lmn_size = Pawtab(itypat)%lmn_size
2378 0 : lmn2_size = Pawtab(itypat)%lmn2_size
2379 0 : lm_size = Paw_an(iat)%lm_size
2380 0 : cplex_dij = Paw_ij(iat)%cplex_dij
2381 0 : qphase = Paw_ij(iat)%qphase
2382 0 : ndij = Paw_ij(iat)%ndij
2383 :
2384 0 : ABI_CHECK(cplex_dij==1,'cplex_dij/=1 not implemented')
2385 0 : ABI_CHECK(qphase==1,'qphase/=1 not implemented')
2386 : !
2387 : ! Eventually compute g_l(r).Y_lm(r) factors for the current atom (if not already done)
2388 0 : if (Pawfgrtab(iat)%gylm_allocated==0) then
2389 0 : if (allocated(Pawfgrtab(iat)%gylm)) then
2390 0 : ABI_FREE(Pawfgrtab(iat)%gylm)
2391 : end if
2392 0 : ABI_MALLOC(Pawfgrtab(iat)%gylm,(Pawfgrtab(iat)%nfgd,lm_size))
2393 0 : Pawfgrtab(iat)%gylm_allocated=2
2394 :
2395 : call pawgylm(Pawfgrtab(iat)%gylm,rdum,rdum2,lm_size,&
2396 0 : & Pawfgrtab(iat)%nfgd,1,0,0,Pawtab(itypat),Pawfgrtab(iat)%rfgd)
2397 : end if
2398 :
2399 : ! === Calculate LM contribution to dijhmxc_val for this atom ===
2400 : ! * Dijxc contains also the Hat term on the FFT mesh while Dijxc_val does not
2401 : ! contain neither the hat term nor the LM sum of onsite terms (they should cancel each other)
2402 : ! FIXME change paw_dij, otherwise I miss tnc in vxc
2403 : ! * prodhxc_core is used to assemble $\int g_l Ylm (vtrial - vxc_val[tn+nhat] dr$ on the FFT mesh ===
2404 : ! * The following quantities do not depend on ij
2405 0 : ABI_MALLOC(prod_hloc ,(lm_size,ndij))
2406 0 : ABI_MALLOC(prodhxc_core,(lm_size,ndij))
2407 0 : prod_hloc =zero
2408 0 : prodhxc_core=zero
2409 0 : do idij=1,ndij
2410 0 : do ilslm=1,lm_size
2411 0 : do ic=1,Pawfgrtab(iat)%nfgd
2412 0 : is2fft=Pawfgrtab(iat)%ifftsph(ic)
2413 0 : gylm=Pawfgrtab(iat)%gylm(ic,ilslm)
2414 0 : prod_hloc (ilslm,idij)=prod_hloc (ilslm,idij) + (vtrial(is2fft,idij)-vxc(is2fft,idij))*gylm
2415 : ! prodhxc_core(ilslm,idij)=prodhxc_core(ilslm,idij) + (vxc_val(is2fft,idij))*gylm
2416 0 : prodhxc_core(ilslm,idij)=prodhxc_core(ilslm,idij) + (vtrial(is2fft,idij)-vxc_val(is2fft,idij))*gylm
2417 : end do
2418 : end do
2419 : end do !idij
2420 :
2421 : ! === Assembly the "Hat" contribution for this atom ====
2422 0 : ABI_MALLOC(dijhl_hat ,(cplex_dij*lmn2_size,ndij))
2423 0 : ABI_MALLOC(dijhmxc_val,(cplex_dij*lmn2_size,ndij))
2424 0 : dijhl_hat =zero
2425 0 : dijhmxc_val=zero
2426 0 : indklmn_ => Pawtab(itypat)%indklmn(1:6,1:lmn2_size)
2427 :
2428 0 : do idij=1,ndij
2429 0 : do klmn=1,lmn2_size
2430 0 : klm =indklmn_(1,klmn)
2431 0 : lmin=indklmn_(3,klmn)
2432 0 : lmax=indklmn_(4,klmn)
2433 :
2434 : ! === $\sum_lm q_ij^l prod* for each idij$ ===
2435 0 : do ils=lmin,lmax,2
2436 0 : lm0=ils**2+ils+1
2437 0 : do mm=-ils,ils
2438 0 : ilslm=lm0+mm
2439 0 : isel=Pawang%gntselect(lm0+mm,klm)
2440 0 : if (isel>0) then
2441 0 : qijl=Pawtab(itypat)%qijl(ilslm,klmn)
2442 0 : dijhl_hat (klmn,idij)=dijhl_hat (klmn,idij) + prod_hloc (ilslm,idij)*qijl
2443 0 : dijhmxc_val(klmn,idij)=dijhmxc_val(klmn,idij) +prodhxc_core(ilslm,idij)*qijl
2444 : end if
2445 : end do
2446 : end do
2447 : end do
2448 : end do
2449 :
2450 0 : ABI_FREE(prod_hloc)
2451 0 : ABI_FREE(prodhxc_core)
2452 :
2453 : ! * Normalization factor due to integration on the FFT mesh
2454 0 : dijhl_hat = dijhl_hat *Cryst%ucvol/DBLE(nfftf)
2455 0 : dijhmxc_val= dijhmxc_val*Cryst%ucvol/DBLE(nfftf)
2456 :
2457 : ! === Now assembly the bare Hamiltonian ===
2458 : ! * Loop over density components overwriting %dij
2459 0 : nsploop=nsppol; if (Paw_ij(iat)%ndij==4) nsploop=4
2460 :
2461 0 : do idij=1,nsploop
2462 : klmn1=1
2463 :
2464 0 : do jlmn=1,lmn_size
2465 0 : j0lmn=jlmn*(jlmn-1)/2
2466 0 : do ilmn=1,jlmn
2467 0 : klmn=j0lmn+ilmn
2468 :
2469 : ! The following gives back the input dij.
2470 : ! since dijxc contains the hat term done on the FFT mesh
2471 0 : if (.FALSE.) then
2472 : Paw_ij(iat)%dij(klmn,idij) = &
2473 : & Pawtab(itypat)%dij0 (klmn) &
2474 : & +Paw_ij(iat)%dijhartree(klmn) &
2475 : & +Paw_ij(iat)%dijxc (klmn,idij) &
2476 : & +dijhl_hat (klmn,idij)
2477 :
2478 : else
2479 : ! === Make nonlocal part of h0 removing the valence contribution ===
2480 : ! Remeber that XC contains already the Hat contribution
2481 : Paw_ij(iat)%dij(klmn,idij) = &
2482 : & Pawtab(itypat)%dij0 (klmn) &
2483 : & +Paw_ij(iat)%dijhartree(klmn) &
2484 : & +Paw_ij(iat)%dijxc (klmn,idij) & ! 2 lines to get the d1-dt1 XC core contribution + XC hat (core+val)
2485 : & -Paw_ij(iat)%dijxc_val (klmn,idij) & ! I suppose that the "hat" term on the FFT mesh in included in both.
2486 0 : & +dijhmxc_val(klmn,idij) ! Local + Hartree - XC val contribution to the "hat" term.
2487 :
2488 : ! Add the U contribution to the
2489 : ! if (.FALSE. .and. Pawtab(itypat)%usepawu/=0) then
2490 0 : if (.TRUE. .and. Pawtab(itypat)%usepawu/=0) then
2491 0 : Paw_ij(iat)%dij(klmn,idij) = Paw_ij(iat)%dij(klmn,idij) + Paw_ij(iat)%dijU(klmn,idij)
2492 : end if
2493 : end if
2494 : ! TODO dijso, dijU, vpawx?
2495 : ! Just to be consistent, update some values.
2496 : !$Paw_ij(iat)%dijhat(klmn,idij)=Paw_ij(iat)%dijhat(klmn,idij)-dijhmxc_val(klmn,idij)
2497 :
2498 : end do !ilmn
2499 : end do !jlmn
2500 : end do !idij
2501 :
2502 : ! this is to be consistent?
2503 : ! deallocate(Paw_ij(iat)%dijvxc_val)
2504 0 : ABI_FREE(dijhl_hat)
2505 0 : ABI_FREE(dijhmxc_val)
2506 : end do !iat
2507 :
2508 : !=== Symmetrize total Dij ===
2509 0 : option_dij=0 ! For total Dij.
2510 : #if 0
2511 : if (paral_atom) then
2512 : call symdij(Cryst%gprimd,Cryst%indsym,ipert0,my_natom,Cryst%natom,Cryst%nsym,Cryst%ntypat,option_dij,&
2513 : & Paw_ij,Pawang,pawprtvol,Pawtab,Cryst%rprimd,Cryst%symafm,Cryst%symrec,&
2514 : & comm_atom=my_comm_atom,mpi_atmtab=my_atmtab)
2515 : else
2516 : call symdij(Cryst%gprimd,,Cryst%indsym,ipert0,my_natom,Cryst%natom,Cryst%nsym,Cryst%ntypat,option_dij,&
2517 : & Paw_ij,Pawang,pawprtvol,Pawtab,Cryst%rprimd,Cryst%symafm,Cryst%symrec)
2518 : end if
2519 : #else
2520 0 : if (paral_atom) then
2521 : call symdij_all(Cryst%gprimd,Cryst%indsym,ipert0,my_natom,Cryst%natom,Cryst%nsym,Cryst%ntypat,&
2522 : & Paw_ij,Pawang,pawprtvol,Pawtab,Cryst%rprimd,Cryst%symafm,Cryst%symrec,&
2523 0 : & comm_atom=my_comm_atom,mpi_atmtab=my_atmtab)
2524 : else
2525 : call symdij_all(Cryst%gprimd,Cryst%indsym,ipert0,my_natom,Cryst%natom,Cryst%nsym,Cryst%ntypat,&
2526 0 : & Paw_ij,Pawang,pawprtvol,Pawtab,Cryst%rprimd,Cryst%symafm,Cryst%symrec)
2527 : end if
2528 : #endif
2529 :
2530 : !Destroy atom table used for parallelism
2531 0 : call free_my_atmtab(my_atmtab,my_atmtab_allocated)
2532 :
2533 : DBG_EXIT("COLL")
2534 :
2535 0 : end subroutine paw_mknewh0
2536 : !!***
2537 :
2538 : !----------------------------------------------------------------------
2539 :
2540 : !!****f* m_paw_denpot/paw_relax_core
2541 : !! NAME
2542 : !! paw_relax_core
2543 : !!
2544 : !! FUNCTION
2545 : !! Relaxes the PAW cores
2546 : !!
2547 : !! INPUTS
2548 : !!
2549 : !! SIDE EFFECTS
2550 : !!
2551 : !! SOURCE
2552 :
2553 38 : subroutine paw_relax_core(pawtab,pawrad,pawang,pawrhoij,ntypat,rcpaw,psps,dtset,&
2554 19 : & cplex,nzlmopt,option,ucvol,paw_an,my_natom,&
2555 19 : & mpi_atmtab,comm_atom,extfpmd)
2556 :
2557 : !Arguments ------------------------------------
2558 : !scalars
2559 : integer, intent(in) :: ntypat,cplex,nzlmopt,option
2560 : integer,intent(in) :: my_natom
2561 : integer,optional,intent(in) :: comm_atom
2562 : real(dp), intent(in) :: ucvol
2563 : type(extfpmd_type), pointer,intent(in), optional :: extfpmd
2564 : type(rcpaw_type),pointer, intent(inout) :: rcpaw
2565 : type(dataset_type), intent(in) :: dtset
2566 : type(pseudopotential_type),intent(inout) :: psps
2567 : type(pawang_type),intent(in),target :: pawang
2568 : !arrays
2569 : integer,optional,target,intent(in) :: mpi_atmtab(:)
2570 : type(pawtab_type), target,intent(inout) :: pawtab(ntypat)
2571 : type(pawrad_type), intent(in) :: pawrad(ntypat)
2572 : type(pawrhoij_type),intent(in) :: pawrhoij(my_natom)
2573 : type(paw_an_type),intent(inout) :: paw_an(my_natom)
2574 :
2575 : !Local variables-------------------------------
2576 : !scalars
2577 : integer :: iatom,itypat
2578 : integer :: mesh_size
2579 : integer :: ispden
2580 : integer :: lm_size
2581 : integer :: my_comm_atom,iat,ierr
2582 : integer :: opt_compch
2583 : logical :: my_atmtab_allocated,paral_atom
2584 : real(dp) :: extfpmd_rho
2585 : type(pawang_type),pointer :: pawang_
2586 : !arrays
2587 19 : integer,pointer :: my_atmtab(:)
2588 19 : logical,allocatable :: lmselect_cur(:)
2589 19 : real(dp), allocatable :: nval(:),nval_tmp(:)
2590 19 : real(dp), allocatable :: tnval(:),tnval_tmp(:)
2591 :
2592 : !******************************************************************************************
2593 :
2594 : !Set up parallelism over atoms
2595 19 : paral_atom=(present(comm_atom).and.(my_natom/=dtset%natom))
2596 19 : nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab
2597 19 : my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom
2598 19 : call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,dtset%natom,my_natom_ref=my_natom)
2599 :
2600 : ! Some setup
2601 19 : if(cplex.ne.1) then
2602 0 : ABI_ERROR('RCPAW: cplex not 1')
2603 : endif
2604 19 : opt_compch=0;if (option/=1) opt_compch=1
2605 19 : pawang_=>pawang
2606 19 : extfpmd_rho=zero
2607 19 : if(present(extfpmd)) then
2608 18 : if(associated(extfpmd)) then
2609 6 : if(extfpmd%pawsph>=1) then
2610 6 : extfpmd_rho=extfpmd%nelect/ucvol
2611 : endif
2612 : endif
2613 : endif
2614 : ! loop over atoms
2615 38 : do itypat=1,dtset%ntypat
2616 19 : mesh_size=pawtab(itypat)%mesh_size
2617 57 : ABI_MALLOC(nval,(mesh_size))
2618 38 : ABI_MALLOC(nval_tmp,(mesh_size))
2619 26913 : nval=zero
2620 38 : ABI_MALLOC(tnval,(mesh_size))
2621 38 : ABI_MALLOC(tnval_tmp,(mesh_size))
2622 26913 : tnval=zero
2623 50 : do iat=1,my_natom
2624 31 : iatom=iat;if (paral_atom) iatom=my_atmtab(iat)
2625 50 : if(dtset%typat(iatom)==itypat) then ! Average over atoms of same typat
2626 43809 : nval_tmp=zero
2627 43809 : tnval_tmp=zero
2628 31 : mesh_size=pawtab(itypat)%mesh_size
2629 31 : lm_size=pawtab(itypat)%lcut_size**2
2630 93 : ABI_MALLOC(lmselect_cur,(lm_size))
2631 1382 : lmselect_cur(:)=.true.
2632 473 : if (nzlmopt==1) lmselect_cur(:)=paw_an(iat)%lmselect(:)
2633 31 : rcpaw%val(iat)%compch_sph=zero
2634 : call pawdensities(rcpaw%val(iat)%compch_sph,cplex,itypat,lmselect_cur,paw_an(iat)%lmselect,lm_size,&
2635 : & rcpaw%val(iat)%nhat1,dtset%nspden,nzlmopt,&
2636 : & opt_compch,0,-1,1,pawang_,dtset%pawprtvol,pawrad(itypat),pawrhoij(iat),pawtab(itypat),rcpaw%val(iat)%rho1,rcpaw%val(iat)%trho1,&
2637 31 : & extfpmd_rho/dtset%nspden)
2638 31 : rcpaw%val(iat)%has_dens=.true.
2639 62 : do ispden=1,dtset%nspden
2640 : nval_tmp(1:pawtab(itypat)%mesh_size)=nval_tmp(1:pawtab(itypat)%mesh_size)+rcpaw%val(iat)%rho1(1:pawtab(itypat)%mesh_size,1,ispden)*&
2641 43809 : & sqrt(four*pi)*pawrad(itypat)%rad(1:pawtab(itypat)%mesh_size)**2
2642 : tnval_tmp(1:pawtab(itypat)%mesh_size)=tnval_tmp(1:pawtab(itypat)%mesh_size)+rcpaw%val(iat)%trho1(1:pawtab(itypat)%mesh_size,1,ispden)*&
2643 43840 : & sqrt(four*pi)*pawrad(itypat)%rad(1:pawtab(itypat)%mesh_size)**2
2644 : end do
2645 31 : ABI_FREE(lmselect_cur)
2646 43840 : nval=nval+nval_tmp
2647 43840 : tnval=tnval+tnval_tmp
2648 : endif
2649 : enddo
2650 19 : ABI_FREE(nval_tmp)
2651 19 : ABI_FREE(tnval_tmp)
2652 : ! mpi reduction
2653 19 : if(paral_atom) then
2654 0 : call xmpi_sum(nval,my_comm_atom,ierr)
2655 0 : call xmpi_bcast(nval,0,my_comm_atom,ierr)
2656 0 : call xmpi_sum(tnval,my_comm_atom,ierr)
2657 0 : call xmpi_bcast(tnval,0,my_comm_atom,ierr)
2658 : endif
2659 26913 : nval=nval/rcpaw%atm(itypat)%mult ! Average over atoms of same typat
2660 26913 : tnval=tnval/rcpaw%atm(itypat)%mult ! Average over atoms of same typat
2661 19 : if(rcpaw%atm(itypat)%mode(1,1,1)==orb_relaxed_core) then ! Relax the core
2662 19 : write(std_out,*) 'RCPAW: core relaxation for typat',itypat,psps%ziontypat(itypat)
2663 19 : if((rcpaw%istep==rcpaw%updatetnc+1.and.rcpaw%updatetnc>0).and.rcpaw%atm(itypat)%zcore_orig>0) then
2664 0 : write(std_out,*) 'RCPAW: freezing tnc at istep = ',rcpaw%istep
2665 : endif
2666 : call atompaw_solve(rcpaw%atp(itypat),pawrad(itypat),pawtab(itypat),&
2667 : & nval,tnval,psps%mqgrid_vl,psps%qgrid_vl,psps%epsatm(itypat),psps%vlspl(:,:,itypat),&
2668 : & psps%ziontypat(itypat),&
2669 : &(rcpaw%istep>=rcpaw%updatepaw(1).and.rcpaw%istep<=rcpaw%updatepaw(2)),&
2670 : &((rcpaw%istep<=rcpaw%updatetnc.or.rcpaw%updatetnc==0).and.rcpaw%atm(itypat)%zcore_orig>0),&
2671 35 : & rcpaw%atm(itypat))
2672 : endif
2673 19 : ABI_FREE(nval)
2674 38 : ABI_FREE(tnval)
2675 : enddo
2676 :
2677 : ! Updae PAW transform related quantities
2678 19 : if(rcpaw%istep>=rcpaw%updatepaw(1).and.rcpaw%istep<=rcpaw%updatepaw(2)+1.and.rcpaw%updatepaw(2)>0) then
2679 4 : call pawinit(zero,0,zero,zero,dtset%pawlcutd,0,0,0,0,0,pawang_,pawrad,0,pawtab,0,0,0,rcpaw_update=.true.)
2680 8 : rcpaw%eijkl_is_sym=.false.
2681 : endif
2682 :
2683 : ! Update core energies
2684 19 : call rcpaw_core_energies(rcpaw,ntypat)
2685 :
2686 : !Destroy atom table used for parallelism
2687 19 : call free_my_atmtab(my_atmtab,my_atmtab_allocated)
2688 :
2689 38 : end subroutine paw_relax_core
2690 : !!***
2691 :
2692 : END MODULE m_paw_denpot
2693 : !!***
2694 :
|
