Branch data 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 : 51040 : subroutine pawdenpot(compch_sph,el_temp,gprimd,ipert,ixc,my_natom,natom,nspden,ntypat,nucdipmom,&
167 : 51040 : & nzlmopt,option,paw_an,paw_an0,paw_energies,paw_ij,pawang,pawprtvol,pawrad,pawrhoij,&
168 : 51040 : & pawspnorb,pawtab,pawxcdev,spnorbscl,xclevel,xc_denpos,xc_taupos,xred,ucvol,znucl,spinaxis,&
169 [ + + + + : 25520 : & 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 : 25520 : integer,pointer :: my_atmtab(:)
218 : 25520 : logical,allocatable :: lmselect_cur(:),lmselect_cur_ep(:),lmselect_ep(:),lmselect_tmp(:)
219 : : real(dp) :: tsec(2)
220 : 25520 : real(dp),allocatable :: dij_ep(:),dijfock_vv(:,:),dijfock_cv(:,:),ff(:)
221 : 25520 : real(dp),allocatable :: one_over_rad2(:),kxc_tmp(:,:,:),k3xc_tmp(:,:,:)
222 : 25520 : 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 : 25520 : real(dp),allocatable :: rho(:),rho1(:,:,:),rho1_ep(:,:,:),rho1xx(:,:,:)
225 : 25520 : real(dp),allocatable :: tau1(:,:,:),ttau1(:,:,:), trho1(:,:,:),trho1_ep(:,:,:)
226 : 25520 : real(dp),allocatable :: vh(:),vxc_tmp(:,:,:),vxctau_tmp(:,:,:)
227 : :
228 : : ! *************************************************************************
229 : :
230 : : DBG_ENTER("COLL")
231 : :
232 : 25520 : call timab(560,1,tsec)
233 : :
234 : : !Various inits
235 [ + + ]: 25520 : hyb_mixing_ =zero ; if(present(hyb_mixing)) hyb_mixing_ =hyb_mixing
236 [ + + ]: 25520 : hyb_mixing_sr_=zero ; if(present(hyb_mixing_sr)) hyb_mixing_sr_=hyb_mixing_sr
237 [ + + ]: 25520 : usercpaw=.false. ; if (present(rcpaw)) usercpaw=associated(rcpaw)
238 [ + + + + ]: 25520 : usefock=0;if (abs(hyb_mixing_)>tol8.or.abs(hyb_mixing_sr_)>tol8) usefock=1
239 [ + + ]: 65518 : usexcnhat=maxval(pawtab(1:ntypat)%usexcnhat)
240 [ + + ]: 65518 : add_core_energy=minval(pawtab(1:ntypat)%add_core_energy)
241 : 25520 : usekden=pawxc_get_usekden(ixc)
242 : 25520 : usenhat = usexcnhat
243 [ + + ]: 80066 : keep_vhartree=(maxval(paw_an(:)%has_vhartree)>0)
244 [ + + ]: 25520 : if (keep_vhartree) usenhat = 1
245 : 25520 : compch_sph=-1.d5
246 [ + + + + ]: 25520 : opt_compch=0;if (option/=1.and.ipert<=0) opt_compch=1
247 [ + + ]: 25520 : if (opt_compch==1) compch_sph=zero
248 [ + + ]: 25520 : nsppol=1;if (my_natom>0) nsppol=pawrhoij(1)%nsppol
249 [ + + ]: 25520 : pertstrg=" ";if (ipert>0) pertstrg="(1)"
250 : :
251 : : !Various checks
252 [ - + ]: 25520 : 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 [ + + ]: 25520 : if (my_natom>0) then
257 [ - + - - : 23736 : 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 [ + + + + : 23736 : 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 [ - + ]: 23736 : 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 [ + + + + : 23736 : 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 [ - + - - : 23736 : 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 [ + + + + : 23736 : 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 [ + + + + : 23736 : 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 [ + + - + ]: 23736 : 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 [ + + - + : 23736 : 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 : 25520 : ipositron=0
300 [ + + ]: 25520 : if (present(electronpositron)) then
301 : 16377 : ipositron=electronpositron_calctype(electronpositron)
302 [ + + - + ]: 16377 : 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 [ - + ]: 16377 : 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 [ + + + + ]: 25520 : paral_atom=(present(comm_atom).and.(my_natom/=natom))
314 [ + + ]: 25520 : nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab
315 [ + + ]: 25520 : 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 : 25520 : & my_natom_ref=my_natom)
318 : :
319 : : !Init energies
320 [ + + ]: 25520 : if (option/=1) then
321 : 21422 : call paw_energies_setzero(paw_energies)
322 : 21422 : e1xc=zero ; e1xcdc=zero
323 : 21422 : etild1xc=zero ; etild1xcdc=zero
324 : 21422 : s1xc=zero ; stild1xc=zero ; sxccore=zero
325 : 21422 : exccore=zero ; eh2=zero ; e1t10=zero
326 : 21422 : eh2dc=zero ; ehpw = zero
327 : 21422 : edftumdc=zero ; edftumdcdc=zero ; edftufll=zero
328 : 21422 : eexex=zero ; eexexdc=zero
329 : 21422 : eextemp=zero ; eexdctemp=zero
330 : 21422 : espnorb=zero ; enucdip=zero
331 : 21422 : efock=zero ; efockdc=zero
332 : 21422 : ekincore=zero ; ehnzc=zero
333 [ + + ]: 21422 : 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 [ + + + + ]: 25520 : if (ipert==natom+1.or.ipert==natom+10) return
341 : :
342 : : !vpotzero is needed for both the energy and the potential
343 [ + + ]: 23738 : if (present(vpotzero)) vpotzero(:)=zero
344 : :
345 : : !Select PAW+U algo, different for DFT and DFPT
346 [ + + ]: 60834 : usepawu=maxval(pawtab(1:ntypat)%usepawu)
347 [ + + + + ]: 60834 : ii=minval(pawtab(1:ntypat)%usepawu);if (ii<0) usepawu=ii
348 : 23738 : non_magnetic_xc=(mod(abs(usepawu),10)==4)
349 : :
350 : : !if PAW+U, compute noccmmp^{\sigma}_{m,m'} occupation matrix
351 [ + + + + : 23738 : if (usepawu/=0.and.ipert<=0.and.ipositron/=1) then
+ - ]
352 [ + + ]: 2326 : 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 : 2166 : & paw_ij,pawang,pawprtvol,pawrhoij,pawtab,rdum2,idum1,idum1,0,usepawu)
359 : : end if
360 : : end if
361 : :
362 : : !Print some titles
363 [ + + ]: 23738 : if (abs(pawprtvol)>=2) then
364 [ + + ]: 1816 : if (nzlmopt<1) write(msg, '(6a)') ch10,' PAW TEST:',ch10,&
365 : 400 : ' ====== Moments of (n1-tn1)',trim(pertstrg),' ========='
366 [ + + ]: 3032 : if (nzlmopt==1) write(msg, '(6a)') ch10,' PAW TEST:',ch10,&
367 : 2832 : ' ==== Non-zero Moments of (n1-tn1)',trim(pertstrg),' ===='
368 : 1616 : call wrtout(std_out,msg,'COLL')
369 [ + + ]: 1616 : if (usexcnhat/=0) then
370 : 1313 : write(msg, '(6a)')' The moments of (n1-tn1-nhat1)',trim(pertstrg),' must be very small...'
371 : 1313 : call wrtout(std_out,msg,'COLL')
372 : : end if
373 : : end if
374 : :
375 : 23738 : extfpmd_rho=zero
376 : 23738 : extfpmd_pawsph=.false.
377 [ + + ]: 23738 : if(present(extfpmd)) then
378 [ + + ]: 15528 : 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 [ + + ]: 75032 : do iatom=1,my_natom
393 [ + + ]: 51294 : iatom_tot=iatom;if (paral_atom) iatom_tot=my_atmtab(iatom)
394 : 51294 : itypat=pawrhoij(iatom)%itypat
395 : 51294 : exchmix=pawtab(itypat)%exchmix
396 : 51294 : lmn2_size=paw_ij(iatom)%lmn2_size
397 : 51294 : lm_size=paw_an(iatom)%lm_size
398 : 51294 : mesh_size=pawtab(itypat)%mesh_size
399 : 51294 : usecore=1;usetcore =pawtab(itypat)%usetcore
400 [ + + ]: 51294 : if (ipert/=0) usecore=0 ! This is true for phonons and Efield pert.
401 [ + + ]: 51294 : if (ipert/=0) usetcore=0 ! This is true for phonons and Efield pert.
402 : 51294 : has_kxc =paw_an(iatom)%has_kxc ;need_kxc =(has_kxc ==1)
403 : 51294 : has_k3xc=paw_an(iatom)%has_k3xc;need_k3xc=(has_k3xc==1)
404 : 51294 : has_vxctau=paw_an(iatom)%has_vxctau ;need_vxctau =(has_vxctau>=1.and.usekden==1)
405 : 51294 : cplex=paw_an(iatom)%cplex
406 : 51294 : cplex_dij=paw_ij(iatom)%cplex_dij
407 : 51294 : cplex_rhoij=pawrhoij(iatom)%cplex_rhoij
408 : 51294 : qphase=pawrhoij(iatom)%qphase
409 : 51294 : ndij=paw_ij(iatom)%ndij
410 : 51294 : iq0_rhoij=cplex_rhoij*lmn2_size
411 : 51294 : iq0_dij=cplex_dij*lmn2_size
412 : 51294 : usepawu=pawtab(itypat)%usepawu
413 [ + + + + ]: 51294 : pawu_algo=merge(PAWU_ALGO_1,PAWU_ALGO_2,ipert<=0.and.usepawu>=0)
414 [ + + ]: 51294 : pawu_dblec=merge(PAWU_FLL,PAWU_AMF,abs(usepawu)==1.or.abs(usepawu)==4)
415 : :
416 : : ! Allocations of "on-site" densities
417 [ - + + - : 256470 : ABI_MALLOC(rho1 ,(cplex*mesh_size,lm_size,nspden))
- + + - -
+ + - + -
+ - - + -
+ ]
418 [ - + + - : 205176 : ABI_MALLOC(trho1,(cplex*mesh_size,lm_size,nspden))
- + + - -
+ + - - +
- + ]
419 [ - + + - : 256470 : ABI_MALLOC(nhat1,(cplex*mesh_size,lm_size,nspden*usenhat))
- + + - +
+ + - + -
+ + - + -
+ ]
420 [ + + + + : 1722921609 : rho1(:,:,:)=zero;trho1(:,:,:)=zero;nhat1(:,:,:)=zero
+ + + + +
+ + + + +
+ + + + ]
421 [ + + ]: 51294 : 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 [ + + ]: 51294 : 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 [ - + + - : 153882 : ABI_MALLOC(lmselect_cur,(lm_size))
+ - + - -
+ - + ]
431 [ + + ]: 589677 : lmselect_cur(:)=.true.
432 [ + + + + ]: 452329 : if (nzlmopt==1) lmselect_cur(:)=paw_an(iatom)%lmselect(:)
433 : :
434 : : ! Store some usefull quantities
435 [ + + ]: 51294 : itypat0=0;if (iatom>1) itypat0=pawrhoij(iatom-1)%itypat
436 [ + + ]: 51294 : if (itypat/=itypat0) then
437 [ - + + - : 90468 : ABI_MALLOC(one_over_rad2,(mesh_size))
+ - + - -
+ - + -
+ ]
438 : 30156 : one_over_rad2(1)=zero
439 [ + + ]: 33933699 : 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 [ + + + - : 51294 : & 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 : 51294 : rcpaw_has_valdens=.false.
464 [ + + ]: 51294 : 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 : 51264 : & pawrhoij(iatom),pawtab(itypat),rho1,trho1,extfpmd_rho=extfpmd_rho/nspden,one_over_rad2=one_over_rad2)
476 : : endif
477 [ + + ]: 51294 : 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 [ + + ]: 51294 : 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 [ + + ]: 51294 : nkxc1 =0;if (paw_an(iatom)%has_kxc /=0) nkxc1 =paw_an(iatom)%nkxc1
508 [ + + + - ]: 51294 : nk3xc1=0;if (paw_an(iatom)%has_k3xc/=0.and.pawxcdev==0) nk3xc1=paw_an(iatom)%nk3xc1
509 [ + + ]: 51294 : if (pawxcdev/=0) then
510 [ - + + - : 236200 : ABI_MALLOC(vxc_tmp,(cplex*mesh_size,lm_size,nspden))
- + + - -
+ + - + -
+ - - + -
+ ]
511 [ + + ]: 47240 : if (need_kxc) then
512 [ - + + - : 335 : ABI_MALLOC(kxc_tmp,(mesh_size,lm_size,nkxc1))
- + + - -
+ + - + -
+ - - + -
+ ]
513 : : end if
514 [ - + ]: 47240 : 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 [ + + ]: 51294 : 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 : 51294 : idum=0
534 [ - + ]: 51294 : if (.not.allocated(vxc_tmp)) then
535 [ # # # # ]: 0 : ABI_MALLOC(vxc_tmp,(0,0,0))
536 : : end if
537 [ + + ]: 51294 : if (.not.allocated(kxc_tmp)) then
538 [ - + - + ]: 51195 : ABI_MALLOC(kxc_tmp,(0,0,0))
539 : : end if
540 [ + + ]: 51294 : if (.not.allocated(k3xc_tmp)) then
541 [ - + - + ]: 51286 : ABI_MALLOC(k3xc_tmp,(0,0,0))
542 : : end if
543 [ + + ]: 51294 : if (.not.allocated(vxctau_tmp)) then
544 [ - + - + ]: 50810 : ABI_MALLOC(vxctau_tmp,(0,0,0))
545 : : end if
546 : :
547 : : ! ===== Vxc1 term =====
548 [ + + ]: 51294 : if (ipositron/=1) then
549 [ + + ]: 48903 : if (pawxcdev/=0) then
550 [ + + ]: 44849 : 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 : 36215 : & 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 : 8634 : & pawang,pawrad(itypat),rho1,usecore,0,vxc_tmp,xclevel)
558 : 8634 : eexcdc=zero
559 : 8634 : 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 [ + + ]: 48903 : if (option/=1) then
576 : 41383 : e1xc=e1xc+eexc
577 : 41383 : e1xcdc=e1xcdc+eexcdc
578 : 41383 : s1xc=s1xc+ssxc
579 : : end if
580 [ + + - + : 1436009399 : if (option<2.or.temp_vxc) paw_an(iatom)%vxc1(:,:,:)=vxc_tmp(:,:,:)
+ + + + +
+ ]
581 [ + + + + : 96164363 : if (option<2.and.need_vxctau) paw_an(iatom)%vxctau1(:,:,:)=vxctau_tmp(:,:,:)
+ + + + +
+ ]
582 [ + + + - : 23012173 : if (need_kxc .and.nkxc1>0 ) paw_an(iatom)%kxc1(:,:,:) =kxc_tmp(:,:,:)
+ + + + +
+ ]
583 [ + + + - : 5041865 : if (need_k3xc.and.nk3xc1>0) paw_an(iatom)%k3xc1(:,:,:)=k3xc_tmp(:,:,:)
+ + + + +
+ ]
584 [ + + ]: 48903 : 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 [ + + ]: 51294 : 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 [ + + ]: 51294 : if (ipositron/=1) then
622 [ + + ]: 48903 : if (pawxcdev/=0) then
623 [ + + ]: 44849 : 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 : 36215 : & 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 : 8634 : & pawang,pawrad(itypat),trho1,usetcore,2*usexcnhat,vxc_tmp,xclevel)
633 : 8634 : eexcdc=zero
634 : 8634 : 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 [ + + ]: 48903 : if (option/=1) then
653 : 41383 : etild1xc=etild1xc+eexc
654 : 41383 : etild1xcdc=etild1xcdc+eexcdc
655 : 41383 : stild1xc=stild1xc+ssxc
656 : : end if
657 [ + + + + : 1436009399 : if (option<2) paw_an(iatom)%vxct1(:,:,:)=vxc_tmp(:,:,:)
+ + + + ]
658 [ + + + + : 96164363 : if (option<2.and.need_vxctau) paw_an(iatom)%vxcttau1(:,:,:)=vxctau_tmp(:,:,:)
+ + + + +
+ ]
659 [ + + + - : 23012173 : if (need_kxc.and. nkxc1>0 ) paw_an(iatom)%kxct1(:,:,:) =kxc_tmp(:,:,:)
+ + + + +
+ ]
660 [ + + + - : 5041865 : if (need_k3xc.and.nk3xc1>0) paw_an(iatom)%k3xct1(:,:,:)=k3xc_tmp(:,:,:)
+ + + + +
+ ]
661 [ + + ]: 48903 : 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 [ + + ]: 51294 : 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 [ + + ]: 51294 : if (option<2) paw_an(iatom)%has_vxc=2
701 [ + + + + ]: 51294 : if (option<2.and.need_vxctau) paw_an(iatom)%has_vxctau=2
702 [ + + + - ]: 51294 : if (need_kxc.and.nkxc1>0) paw_an(iatom)%has_kxc=2
703 : :
704 : : ! Update core XC contribution to energy
705 [ + + + + ]: 51294 : if (option/=1.and.ipositron/=1) then
706 : 41383 : exccore=exccore+pawtab(itypat)%exccore
707 : 41383 : sxccore=sxccore+pawtab(itypat)%sxccore
708 : : end if
709 : :
710 : : ! =========== Compute valence-only XC potentials ===========
711 : : ! ==========================================================
712 [ + + + + : 51294 : if (ipert==0.and.paw_an(iatom)%has_vxcval==1.and.(option==0.or.option==1)) then
+ - ]
713 [ + - - + ]: 35 : 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 [ + - ]: 35 : if (pawxcdev/=0) then
719 : 35 : write(msg,'(4a,es16.6)')ch10,&
720 : 35 : & ' pawdenpot : Computing valence-only v_xc[n1] using moments ',ch10,&
721 [ + + + + : 494114 : & ' Min density rho1 = ',MINVAL(rho1)
+ + + + +
- - + -
- ]
722 : 35 : 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 : 35 : & 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 [ + + + + : 494079 : if (option<2) paw_an(iatom)%vxc1_val(:,:,:)=vxc_tmp(:,:,:)
+ + ]
735 : :
736 : : ! ===== tVxc1_val term =====
737 [ + - ]: 35 : if (pawxcdev/=0) then
738 [ + + ]: 35 : 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 : 4 : write(msg,'(4a,e16.6)')ch10,&
745 : 4 : & ' pawdenpot : Computing valence-only v_xc[tn1] using moments ',ch10,&
746 [ + + + + : 100028 : & ' Min density trho1 = ',MINVAL(trho1)
+ + + + +
- - + -
- ]
747 : : end if
748 : 35 : 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 : 35 : & 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 [ + + + + : 494079 : paw_an(iatom)%vxct1_val(:,:,:)=vxc_tmp(:,:,:)
+ + ]
763 : 35 : paw_an(iatom)%has_vxcval=2
764 : : end if
765 : : end if ! valence-only XC potentials
766 : :
767 [ - + ]: 51294 : ABI_FREE(vxc_tmp)
768 [ - + ]: 51294 : ABI_FREE(vxctau_tmp)
769 [ - + ]: 51294 : ABI_FREE(kxc_tmp)
770 [ - + ]: 51294 : ABI_FREE(k3xc_tmp)
771 : :
772 : : ! ===== Compute first part of local exact-exchange energy term =====
773 : : ! ===== Also compute corresponding potential =====
774 : : ! ==================================================================
775 : :
776 [ + + + - : 51294 : 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 [ + + ]: 51294 : itypat0=0;if (iatom<my_natom) itypat0=pawrhoij(iatom+1)%itypat
808 [ + + ]: 51294 : if (itypat/=itypat0) then
809 [ - + ]: 30156 : ABI_FREE(one_over_rad2)
810 : : end if
811 : :
812 [ - + ]: 51294 : ABI_FREE(lmselect_cur)
813 : :
814 : : ! ==== Compute Hartree potential terms and some energy terms ====
815 : : ! ===============================================================
816 : :
817 : : ! Hartree Dij computation
818 [ + + ]: 51294 : if (ipositron/=1) then
819 : 48903 : eijkl_is_sym=.true.
820 [ + + ]: 48903 : 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 : 48903 : & is_sym=eijkl_is_sym)
825 [ + + ]: 48903 : if (option/=1) then
826 : 41383 : call pawaccenergy_nospin(eh2dc,pawrhoij(iatom),paw_ij(iatom)%dijhartree,1,qphase,pawtab(itypat))
827 : : end if
828 [ + + ]: 48903 : 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 : 51294 : paw_ij(iatom)%has_dijhartree=2
899 : :
900 : : ! Hartree energy computation
901 [ + + ]: 51294 : if (option/=1) then
902 : 43774 : call pawaccenergy_nospin(eh2,pawrhoij(iatom),paw_ij(iatom)%dijhartree,1,qphase,pawtab(itypat))
903 [ + + ]: 43774 : 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 [ + + ]: 51294 : if (option/=1) then
912 : 43774 : ehnzc=ehnzc+pawtab(itypat)%ehnzc
913 : 43774 : 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 [ + + ]: 51294 : 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 : 51294 : keep_vhartree=(paw_an(iatom)%has_vhartree>0)
938 [ + + + + : 203681 : 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 [ + + ]: 2780 : if (.not. allocated(paw_an(iatom)%vh1)) then
943 [ - + + - : 7620 : ABI_MALLOC(paw_an(iatom)%vh1,(cplex*mesh_size,1,1))
+ - + - -
+ - + -
+ ]
944 : : end if
945 [ + + ]: 2780 : if (.not. allocated(paw_an(iatom)%vht1)) then
946 [ - + + - : 600 : ABI_MALLOC(paw_an(iatom)%vht1,(cplex*mesh_size,1,1))
+ - + - -
+ - + -
+ ]
947 : : end if
948 [ - + + - : 8340 : ABI_MALLOC(rho,(mesh_size))
+ - + - -
+ - + -
+ ]
949 [ - + + - : 5560 : ABI_MALLOC(vh,(mesh_size))
- + - + -
+ ]
950 : :
951 : : ! Construct vh1 and tvh1
952 [ + + ]: 5560 : 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 [ + + ]: 3155308 : rho(1:mesh_size)=rho1(iq:cplex*mesh_size:cplex,1,1)*four_pi*pawrad(itypat)%rad(1:mesh_size)**2
957 [ + - ]: 2780 : if (usecore==1) then
958 : : rho(1:mesh_size)=rho(1:mesh_size)+sqrt(four_pi)*pawtab(itypat)%coredens(1:mesh_size) &
959 [ + + ]: 3155308 : & *four_pi*pawrad(itypat)%rad(1:mesh_size)**2
960 : : end if
961 : 2780 : call poisson(rho,0,pawrad(itypat),vh)
962 [ + + ]: 3152528 : vh(2:mesh_size)=(vh(2:mesh_size)-sqrt(four_pi)*znucl(itypat))/pawrad(itypat)%rad(2:mesh_size)
963 : 2780 : call pawrad_deducer0(vh,mesh_size,pawrad(itypat))
964 [ + + ]: 3155308 : 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 [ + + ]: 3155308 : rho(1:mesh_size)=trho1(iq:cplex*mesh_size:cplex,1,1)*four_pi*pawrad(itypat)%rad(1:mesh_size)**2
972 [ + + ]: 2780 : if (usenhat/=0) then
973 : : rho(1:mesh_size)=rho(1:mesh_size)+nhat1(iq:cplex*mesh_size:cplex,1,1) &
974 [ + + ]: 956737 : & *four_pi*pawrad(itypat)%rad(1:mesh_size)**2
975 : : end if
976 [ + + ]: 2780 : if (usetcore==1) then
977 : : rho(1:mesh_size)=rho(1:mesh_size)+sqrt(four_pi)*pawtab(itypat)%tcoredens(1:mesh_size,1) &
978 [ + + ]: 2089084 : & *four_pi*pawrad(itypat)%rad(1:mesh_size)**2
979 : : end if
980 : 2780 : call poisson(rho,0,pawrad(itypat),vh)
981 [ + + ]: 3152528 : vh(2:mesh_size)=(vh(2:mesh_size)-sqrt(four_pi)*znucl(itypat))/pawrad(itypat)%rad(2:mesh_size)
982 : 2780 : call pawrad_deducer0(vh,mesh_size,pawrad(itypat))
983 [ + + ]: 3158088 : paw_an(iatom)%vht1(iq:cplex*mesh_size:cplex,1,1)=vh(1:mesh_size)
984 : :
985 : : end do ! cplex phase
986 : :
987 : 2780 : paw_an(iatom)%has_vhartree=2
988 [ - + ]: 2780 : ABI_FREE(rho)
989 [ - + ]: 2780 : ABI_FREE(vh)
990 : : end if
991 : :
992 : : ! ========= Compute PAW+U and energy contribution =========
993 : : ! ==========================================================
994 : :
995 [ + + + - : 51294 : if (usepawu/=0.and.usepawu<10.and.ipositron/=1.and.option/=1) then
+ + ]
996 : :
997 [ + + ]: 2843 : 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 : 2541 : & 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 : 302 : & pawrhoij(iatom),pawtab(itypat))
1008 : 302 : paw_ij(iatom)%has_dijU=2
1009 : : !Then, compute energy
1010 : : if (option/=1) then
1011 : 302 : etmp=zero
1012 : 302 : call pawaccenergy(etmp,pawrhoij(iatom),paw_ij(iatom)%dijU,cplex_dij,qphase,ndij,pawtab(itypat))
1013 : 302 : edftumdc=edftumdc+half*etmp ; edftumdcdc=edftumdcdc-half*etmp
1014 : : !Add FLL double-counting part
1015 [ + + + + ]: 302 : 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 [ + + + + : 203681 : if (any(abs(nucdipmom(:,iatom))>tol8).and.ipert==0.and.ipositron/=1) then
+ + + - +
- ]
1028 : :
1029 [ - + ]: 639 : ABI_CHECK(cplex_rhoij==2,'BUG in pawdenpot: rhoij must be complex for ND moments!')
1030 [ - + ]: 639 : 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 [ + - ]: 639 : 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 : 639 : & paw_an(iatom)%vh1,paw_an(iatom)%vxc1,xred,znucl(itypat),paw_ij(iatom)%zora)
1037 : 639 : paw_ij(iatom)%has_dijnd=2
1038 : : end if
1039 : :
1040 : : ! Compute nuclear dipole contribution to energy
1041 [ + + ]: 639 : if (option/=1) then
1042 : : call pawaccenergy_nospin(enucdip,pawrhoij(iatom),paw_ij(iatom)%dijnd,&
1043 : 620 : & cplex_dij,1,pawtab(itypat))
1044 : : end if
1045 : :
1046 : : end if
1047 : :
1048 : : ! ========= Compute spin-orbit energy contribution ========
1049 : : ! ==========================================================
1050 : :
1051 [ + + + - : 51294 : if (pawspnorb>0.and.ipert==0.and.ipositron/=1) then
+ - ]
1052 : :
1053 : : ! Compute spin-orbit contribution to Dij
1054 [ - + - - ]: 2277 : 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 : 2277 : & nucdipmom=nucdipmom(1:3,iatom),spinaxis=spinaxis)
1059 : 2277 : paw_ij(iatom)%has_dijso=2
1060 : : end if
1061 : :
1062 : : ! Compute spin-orbit contribution to on-site energy
1063 [ + + + + ]: 2277 : if (option/=1.and.cplex_rhoij==2) then
1064 : : call pawaccenergy(espnorb,pawrhoij(iatom),paw_ij(iatom)%dijso,&
1065 : 2247 : & 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 [ + + + - : 51294 : 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 [ + + ]: 51294 : 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 [ + + + - : 51294 : 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 [ + + + + ]: 51294 : if (option/=1.and.ipert<=0) then
1145 : 38163 : 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 [ + + ]: 38163 : 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 [ - + ]: 51294 : ABI_FREE(rho1)
1161 [ - + ]: 51294 : ABI_FREE(trho1)
1162 [ - + ]: 51294 : ABI_FREE(nhat1)
1163 [ + + ]: 51294 : if (usekden==1) then
1164 [ - + ]: 484 : ABI_FREE(tau1)
1165 [ - + ]: 484 : ABI_FREE(ttau1)
1166 : : end if
1167 [ + + ]: 51294 : 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 [ + + ]: 51294 : if (.not.keep_vhartree) then
1175 : 51054 : paw_an(iatom)%has_vhartree=0
1176 [ + + ]: 51054 : if (allocated(paw_an(iatom)%vh1)) then
1177 : 2540 : ABI_FREE(paw_an(iatom)%vh1)
1178 : : end if
1179 : : end if
1180 [ - + ]: 75032 : 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 [ + + ]: 23738 : if (option/=1) then
1196 [ + + ]: 19805 : if (ipert==0) then
1197 : 17088 : paw_energies%epaw=e1xc+half*eh2+e1t10-exccore-etild1xc+edftumdc+edftufll+eexex+espnorb+efock+enucdip+ehpw
1198 : 17088 : paw_energies%epaw_dc=e1xc-e1xcdc-half*eh2dc-exccore-etild1xc+etild1xcdc+edftumdcdc-eexex-efockdc
1199 : 17088 : paw_energies%epaw_xc=e1xc-etild1xc !+eexex
1200 : 17088 : paw_energies%entropy_paw=s1xc-sxccore-stild1xc ! PAW entropy coming from finite-temperature xc functionals
1201 : 17088 : paw_energies%epaw_core=zero ; paw_energies%epaw_core_dc=zero
1202 [ + + + + ]: 17088 : if (abs(ekincore)>tiny(zero).and.ipositron/=1) then
1203 : 7041 : paw_energies%epaw_core=exccore+ekincore+ehnzc
1204 : 7041 : paw_energies%epaw_core_dc=paw_energies%epaw_core
1205 : : end if
1206 [ + + ]: 17088 : 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 : 2717 : paw_energies%epaw=e1xc-etild1xc+eh2+two*edftumdc
1212 : 2717 : paw_energies%epaw_dc=zero
1213 : 2717 : paw_energies%epaw_core=zero
1214 : 2717 : paw_energies%epaw_core_dc=zero
1215 : 2717 : 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 [ + + ]: 23738 : 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 [ + + + + ]: 23738 : 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 : 23738 : call free_my_atmtab(my_atmtab,my_atmtab_allocated)
1273 : :
1274 : 23738 : call timab(560,2,tsec)
1275 : :
1276 : : DBG_EXIT("COLL")
1277 : :
1278 [ - + - + : 51040 : 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 : 71224 : subroutine pawdensities(compch_sph,cplex,iatom,lmselectin,lmselectout,lm_size,nhat1,nspden,nzlmopt,&
1336 : 71224 : & 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 : 142448 : real(dp) :: compchspha(cplex),compchsphb(cplex),ro(cplex),ro_ql(cplex),ro_rg(cplex)
1363 : 71224 : real(dp),allocatable :: aa(:),bb(:)
1364 : 71224 : real(dp),pointer :: one_over_rad2_(:)
1365 : :
1366 : : ! *************************************************************************
1367 : :
1368 : : DBG_ENTER("COLL")
1369 : :
1370 : : !Compatibility tests
1371 [ - + ]: 71173 : 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 [ - + ]: 71224 : 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 [ - + ]: 71224 : if(nspden>pawrhoij%nspden) then
1380 : 0 : msg='nspden must be <= pawrhoij%nspden!'
1381 : 0 : ABI_BUG(msg)
1382 : : end if
1383 [ - + ]: 71224 : if (cplex>pawrhoij%qphase) then
1384 : 0 : msg='cplex must be <= pawrhoij%qphase!'
1385 : 0 : ABI_BUG(msg)
1386 : : end if
1387 [ + + ]: 71224 : if (nzlmopt/=1) then
1388 [ + + + - : 280154 : 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 [ - + ]: 71224 : 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 [ + + + + : 919316936 : rho1=zero
+ + ]
1400 [ + + + + : 916939256 : if (opt_dens<2) trho1=zero
+ + + + ]
1401 [ + + + + : 319403163 : if (opt_dens==0) nhat1=zero
+ + + + ]
1402 : 71224 : mesh_size=pawtab%mesh_size;dplex=cplex-1
1403 : 71224 : iq0=pawrhoij%cplex_rhoij*pawrhoij%lmn2_size
1404 [ + + + + ]: 325220 : if (nzlmopt<1) lmselectout(1:lm_size)=.true.
1405 [ + + ]: 71224 : if (present(one_over_rad2)) then
1406 : 57182 : one_over_rad2_ => one_over_rad2
1407 : : else
1408 [ - + + - : 42126 : ABI_MALLOC(one_over_rad2_,(mesh_size))
+ - + - -
+ - + ]
1409 : 14042 : one_over_rad2_(1)=zero
1410 [ + + ]: 11102192 : 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 [ + + ]: 156783 : do ispden=1,nspden
1417 : :
1418 : : ! -- Loop over ij channels (basis components)
1419 : 85559 : jrhoij=1
1420 [ + + ]: 2956087 : do irhoij=1,pawrhoij%nrhoijsel
1421 : 2870528 : klmn=pawrhoij%rhoijselect(irhoij)
1422 : 2870528 : klm =pawtab%indklmn(1,klmn)
1423 : 2870528 : kln =pawtab%indklmn(2,klmn)
1424 : 2870528 : lmin=pawtab%indklmn(3,klmn)
1425 : 2870528 : lmax=pawtab%indklmn(4,klmn)
1426 : :
1427 : : ! Retrieve rhoij
1428 [ + + ]: 2870528 : if (pawrhoij%nspden/=2) then
1429 : 2008298 : ro(1)=pawrhoij%rhoijp(jrhoij,ispden)
1430 [ + + ]: 2008298 : if (cplex==2) ro(2)=pawrhoij%rhoijp(iq0+jrhoij,ispden)
1431 : : else
1432 [ + + ]: 862230 : if (ispden==1) then
1433 : 431115 : ro(1)=pawrhoij%rhoijp(jrhoij,1)+pawrhoij%rhoijp(jrhoij,2)
1434 [ + + ]: 431115 : if (cplex==2) ro(2)=pawrhoij%rhoijp(iq0+jrhoij,1)+pawrhoij%rhoijp(iq0+jrhoij,2)
1435 [ + - ]: 431115 : else if (ispden==2) then
1436 : 431115 : ro(1)=pawrhoij%rhoijp(jrhoij,1)
1437 [ + + ]: 431115 : if (cplex==2) ro(2)=pawrhoij%rhoijp(iq0+jrhoij,1)
1438 : : end if
1439 : : end if
1440 [ + + ]: 5798623 : ro(1:cplex)=pawtab%dltij(klmn)*ro(1:cplex)
1441 : :
1442 : : ! First option: all on-site densities are computed (opt_dens==0)
1443 : : ! --------------------------------------------------------------
1444 [ + + ]: 2870528 : if (opt_dens==0) then
1445 [ - + + + ]: 4075509 : do ll=lmin,lmax,2
1446 [ + + ]: 15737537 : do ilm=ll**2+1,(ll+1)**2
1447 [ + + ]: 14404514 : if (lmselectin(ilm)) then
1448 : 10120599 : isel=pawang%gntselect(ilm,klm)
1449 [ + + ]: 10120599 : if (isel>0) then
1450 [ + + ]: 5524391 : ro_ql(1:cplex)=ro(1:cplex)*pawtab%qijl(ilm,klmn)
1451 [ + + ]: 5524391 : 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 [ + + ]: 5524391 : & +ro_ql(1:cplex)*pawtab%shapefunc(1,ll+1)
1455 : : ! == rho1(r>0), trho1(r>0), nhat1(r>0)
1456 [ + + ]: 1861175313 : do ir=2,mesh_size
1457 : : rho1(cplex*ir-dplex:ir*cplex,ilm,ispden) =rho1(cplex*ir-dplex:ir*cplex,ilm,ispden)&
1458 [ + + ]: 3742365650 : & +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 [ + + ]: 3742365650 : & +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 [ + + ]: 3745102854 : & +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 [ + + ]: 1537505 : else if (opt_dens==1) then
1472 [ - + + + ]: 4272115 : do ll=lmin,lmax,2
1473 [ + + ]: 13960223 : do ilm=ll**2+1,(ll+1)**2
1474 [ + + ]: 12429738 : if (lmselectin(ilm)) then
1475 : 8935091 : isel=pawang%gntselect(ilm,klm)
1476 [ + + ]: 8935091 : if (isel>0) then
1477 [ + + ]: 5260329 : ro_rg(1:cplex)=ro(1:cplex)*pawang%realgnt(isel)
1478 : : ! == rho1(r>0), trho1(r>0)
1479 [ + + ]: 3141439412 : do ir=2,mesh_size
1480 : : rho1(cplex*ir-dplex:ir*cplex,ilm,ispden) =rho1(cplex*ir-dplex:ir*cplex,ilm,ispden)&
1481 [ + + ]: 6322080329 : & +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 [ + + ]: 6324691585 : & +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 [ + - ]: 7020 : else if (opt_dens==2) then
1493 [ + + + + : 7020 : if (opt_l<0.or.(pawtab%indklmn(3,klmn)==0.and.pawtab%indklmn(4,klmn)==2*opt_l)) then
+ + ]
1494 [ - + + + ]: 16115 : do ll=lmin,lmax,2
1495 [ + + ]: 77844 : do ilm=ll**2+1,(ll+1)**2
1496 [ + + ]: 70824 : if (lmselectin(ilm)) then
1497 : 55273 : isel=pawang%gntselect(ilm,klm)
1498 [ + + ]: 55273 : if (isel>0) then
1499 [ + + ]: 24214 : ro_rg(1:cplex)=ro(1:cplex)*pawang%realgnt(isel)
1500 : : ! == rho1(r>0)
1501 [ + + ]: 13102262 : do ir=2,mesh_size
1502 : : rho1(cplex*ir-dplex:ir*cplex,ilm,ispden) =rho1(cplex*ir-dplex:ir*cplex,ilm,ispden)&
1503 [ + + ]: 26192417 : & +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 : 2956087 : jrhoij=jrhoij+pawrhoij%cplex_rhoij
1515 : : end do
1516 : :
1517 : : ! Compute rho1(r=0) and trho1(r=0)
1518 [ + + ]: 85559 : if (cplex==2) then
1519 [ - + - + ]: 1520 : ABI_MALLOC(aa,(5))
1520 [ - + - + ]: 1520 : ABI_MALLOC(bb,(5))
1521 : : end if
1522 [ + + ]: 85559 : if (opt_dens==0.or.opt_dens==1) then
1523 [ + + ]: 360414 : do ll=0,pawtab%lcut_size-1
1524 [ + + ]: 1349762 : do ilm=ll**2+1,(ll+1)**2
1525 [ + + ]: 1264296 : if (lmselectin(ilm)) then
1526 [ + + ]: 753179 : if (cplex==1) then
1527 : 738875 : call pawrad_deducer0(rho1 (:,ilm,ispden),mesh_size,pawrad)
1528 : 738875 : 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 [ + + ]: 85559 : if (cplex==2) then
1564 [ - + ]: 1520 : ABI_FREE(aa)
1565 [ - + ]: 1520 : ABI_FREE(bb)
1566 : : end if
1567 : :
1568 : : ! Add extfpmd contribution
1569 [ + + ]: 85559 : if(present(extfpmd_rho)) then
1570 [ + + + + ]: 102531012 : rho1(:,1,:)=rho1(:,1,:)+sqrt(four_pi)*extfpmd_rho
1571 [ + - ]: 63574 : if(opt_dens<=1) then
1572 [ + + + + ]: 102531012 : 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 [ + + ]: 85559 : if (nzlmopt==-1) then
1578 [ + + ]: 66750 : do ll=0,pawtab%lcut_size-1
1579 [ + + ]: 230960 : do ilm=ll**2+1,(ll+1)**2
1580 : 164210 : m1=zero;mt1=zero
1581 [ + + ]: 164210 : if (cplex==1) then
1582 [ + - - + : 126053112 : m1=maxval(abs(rho1 (1:mesh_size,ilm,ispden)))
- - + + ]
1583 [ + + + - : 125945625 : 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 [ + + ]: 214580 : if (ispden==1) then
1597 [ + + + + : 144630 : if ((ilm>1).and.(m1<tol16).and.(mt1<tol16)) then
+ - ]
1598 : 65682 : lmselectout(ilm)=.false.
1599 : : end if
1600 [ + + ]: 19580 : else if (.not.(lmselectout(ilm))) then
1601 [ + + - + ]: 7576 : 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 [ + + + + : 85559 : if (opt_compch==1.and.ispden==1.and.opt_dens<2) then
+ - ]
1609 [ - + + - : 114492 : ABI_MALLOC(aa,(mesh_size))
+ - + - -
+ - + -
+ ]
1610 : : aa(1:mesh_size)=(rho1(1:mesh_size,1,1)-trho1(1:mesh_size,1,1)) &
1611 [ + + ]: 50437295 : & *pawrad%rad(1:mesh_size)**2
1612 : 38164 : call simp_gen(compchspha(1),aa,pawrad)
1613 : 38164 : compch_sph=compch_sph+compchspha(1)*sqrt(four_pi)
1614 [ - + ]: 38164 : ABI_FREE(aa)
1615 : : end if
1616 : :
1617 : : ! -- Print out moments of densities (if requested)
1618 [ + + + - : 156783 : if (abs(pawprtvol)>=2.and.opt_print==1.and.opt_dens<2) then
+ - ]
1619 [ - + + - : 15810 : ABI_MALLOC(aa,(cplex*mesh_size))
+ - + - -
+ - + -
+ ]
1620 [ - + + - : 10540 : ABI_MALLOC(bb,(cplex*mesh_size))
- + - + -
+ ]
1621 [ + + ]: 5270 : if (opt_dens==0) then
1622 : 4344 : write(msg,'(2a,i3,a,i1,3a)') ch10, &
1623 : 4344 : & ' Atom ',iatom,' (ispden=',ispden,'):',ch10,&
1624 : 8688 : & ' ******* 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 : 5270 : call wrtout(std_out,msg,'PERS')
1631 [ + + ]: 31960 : do ll=0,pawtab%lcut_size-1
1632 [ + + ]: 169470 : do ilm=ll**2+1,(ll+1)**2
1633 [ + + ]: 164200 : if (lmselectin(ilm)) then
1634 [ + + ]: 156564 : do iplex=1,cplex
1635 [ + + ]: 156564 : if (opt_dens==0) then
1636 [ + + ]: 49803364 : do ir=1,mesh_size
1637 : 49738072 : ii=cplex*(ir-1)+iplex
1638 : 49738072 : ro(1)=pawrad%rad(ir)**(2+ll)
1639 : 49738072 : aa(ir)=ro(1)*(rho1(ii,ilm,ispden)-trho1(ii,ilm,ispden))
1640 : 49803364 : bb(ir)=ro(1)*nhat1(ii,ilm,ispden)
1641 : : end do
1642 : 65292 : call simp_gen(compchspha(iplex),aa,pawrad)
1643 : 65292 : call simp_gen(compchsphb(iplex),bb,pawrad)
1644 : : else
1645 [ + + ]: 12453070 : do ir=1,mesh_size
1646 : 12440080 : ii=cplex*(ir-1)+iplex
1647 : 12440080 : ro(1)=pawrad%rad(ir)**(2+ll)
1648 : 12453070 : aa(ir)=ro(1)*(rho1(ii,ilm,ispden)-trho1(ii,ilm,ispden))
1649 : : end do
1650 : 12990 : call simp_gen(compchspha(iplex),aa,pawrad)
1651 : : end if
1652 : : end do
1653 [ + + ]: 78282 : if (opt_dens==0) then
1654 [ + - ]: 65292 : if (cplex==1) then
1655 : : write(msg,'(3x,a,2i2,2(a,es14.7))') &
1656 : 65292 : & 'l,m=',ll,ilm-(ll**2+ll+1),': M=',compchspha(1),&
1657 : 130584 : & ' ** 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 [ + - ]: 12990 : if (cplex==1) then
1665 : : write(msg,'(3x,a,2i2,a,es14.7)') &
1666 : 12990 : & '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 : 78282 : call wrtout(std_out,msg,'PERS')
1673 : : end if
1674 : : end do
1675 : : end do
1676 [ - + ]: 5270 : ABI_FREE(aa)
1677 [ - + ]: 5270 : ABI_FREE(bb)
1678 : : end if
1679 : :
1680 : : ! ----- End loop over spin components
1681 : : end do
1682 : :
1683 [ + + ]: 71224 : if (.not.present(one_over_rad2)) then
1684 : 14042 : ABI_FREE(one_over_rad2_)
1685 : : end if
1686 : :
1687 : : DBG_EXIT("COLL")
1688 : :
1689 [ + + - + : 142448 : 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 : 145259 : 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 : 145259 : real(dp),pointer :: rhoij(:,:)
2065 : :
2066 : : ! *************************************************************************
2067 : :
2068 : : DBG_ENTER("COLL")
2069 : :
2070 : : !Compatibility tests
2071 : 145259 : if (pawrhoij%qphase/=qphase_dij) then
2072 : 0 : msg='pawaccenergy: pawrhoij%qphase/=qphase_dij!'
2073 : 0 : ABI_BUG(msg)
2074 : : end if
2075 [ + + - + ]: 145259 : 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 : 145259 : nspden_rhoij=pawrhoij%nspden
2082 : 145259 : lmn2_size=pawtab%lmn2_size
2083 : :
2084 : : !Special treatment for nspden
2085 : 145259 : nsploop=nspden_rhoij
2086 [ + + + + ]: 145259 : if (nspden_dij==1.and.nspden_rhoij==4) nsploop=1
2087 : :
2088 : : !Non-collinear case: need a temporary rhoij
2089 [ + + + + ]: 11717 : if (nspden_rhoij==4.and.nspden_dij==4) then
2090 : 2379 : cplex_rhoij=2
2091 [ + - - + : 9516 : ABI_MALLOC(rhoij,(2*lmn2_size,4))
+ - + - +
- + - - +
- + ]
2092 : : else
2093 : 142880 : cplex_rhoij=pawrhoij%cplex_rhoij
2094 : 142880 : rhoij => pawrhoij%rhoijp
2095 : : end if
2096 : :
2097 [ + + - + ]: 145259 : add_imaginary=(cplex_dij==2.and.cplex_rhoij==2)
2098 : :
2099 : : !Loop over qphase components
2100 [ + + ]: 291946 : do iq=1,qphase_dij
2101 : 146687 : iq0_rhoij=(iq-1)*lmn2_size*cplex_rhoij
2102 : 146687 : iq0_dij =(iq-1)*lmn2_size*cplex_dij
2103 : :
2104 : : ! Non-collinear case
2105 [ + + + + ]: 146687 : if (nspden_rhoij==4.and.nspden_dij==4) then
2106 [ + + + + ]: 2076983 : rhoij(:,:)=zero
2107 : 2379 : jrhoij=(iq-1)*lmn2_size*pawrhoij%cplex_rhoij+1 ; krhoij=1
2108 [ + + ]: 174724 : do irhoij=1,pawrhoij%nrhoijsel
2109 : 172345 : klmn=pawrhoij%rhoijselect(irhoij)
2110 : 172345 : rhoij(krhoij ,1)= half*(pawrhoij%rhoijp(jrhoij,1)+pawrhoij%rhoijp(jrhoij,4))
2111 : 172345 : 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 : 172345 : rhoij(krhoij ,3)= half*pawrhoij%rhoijp(jrhoij,2)
2115 : 172345 : rhoij(krhoij+1,3)= half*pawrhoij%rhoijp(jrhoij,3)
2116 : 172345 : rhoij(krhoij ,4)= half*pawrhoij%rhoijp(jrhoij,2)
2117 : 172345 : rhoij(krhoij+1,4)=-half*pawrhoij%rhoijp(jrhoij,3)
2118 [ + - ]: 172345 : if (pawrhoij%cplex_rhoij==2) then
2119 : 172345 : rhoij(krhoij+1,1)= half*(pawrhoij%rhoijp(jrhoij+1,1)+pawrhoij%rhoijp(jrhoij+1,4))
2120 : 172345 : 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 : 172345 : rhoij(krhoij ,3)= rhoij(krhoij ,3)-half*pawrhoij%rhoijp(jrhoij+1,3)
2124 : 172345 : rhoij(krhoij+1,3)= rhoij(krhoij+1,3)+half*pawrhoij%rhoijp(jrhoij+1,2)
2125 : 172345 : rhoij(krhoij ,4)= rhoij(krhoij ,4)+half*pawrhoij%rhoijp(jrhoij+1,3)
2126 : 172345 : rhoij(krhoij+1,4)= rhoij(krhoij+1,4)+half*pawrhoij%rhoijp(jrhoij+1,2)
2127 : : end if
2128 : 174724 : 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 [ + + ]: 318333 : do isp_rhoij=1,nsploop
2135 : 171646 : isp_dij=min(isp_rhoij,nspden_dij)
2136 : 171646 : jrhoij=iq0_rhoij+1
2137 [ + + ]: 6049190 : do irhoij=1,pawrhoij%nrhoijsel
2138 : 5730857 : klmn=pawrhoij%rhoijselect(irhoij)
2139 : 5730857 : kklmn=iq0_dij+cplex_dij*(klmn-1)+1
2140 : 5730857 : etmp=rhoij(jrhoij,isp_rhoij)*dij(kklmn,isp_dij)
2141 [ + + ]: 5730857 : if (add_imaginary) etmp=etmp-rhoij(jrhoij+1,isp_rhoij)*dij(kklmn+1,isp_dij)
2142 : 5730857 : epaw=epaw+etmp*pawtab%dltij(klmn)
2143 : 5902503 : jrhoij=jrhoij+cplex_rhoij
2144 : : end do
2145 : : end do ! nsploop
2146 : :
2147 : : ! Contribution to on-site energy (imaginary part)
2148 [ + + + + ]: 291946 : 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 [ + + + + ]: 145259 : if (nspden_rhoij==4.and.nspden_dij==4) then
2173 [ - + ]: 2379 : ABI_FREE(rhoij)
2174 : : end if
2175 : :
2176 : : DBG_EXIT("COLL")
2177 : :
2178 [ - + ]: 145259 : 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 : 141666 : 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 : 141666 : real(dp), contiguous, pointer :: dij_2D(:,:)
2232 : :
2233 : : ! *************************************************************************
2234 : :
2235 : 141666 : size_dij=size(dij)
2236 : :
2237 : : #ifdef HAVE_FC_ISO_C_BINDING
2238 : 141666 : cptr=c_loc(dij(1))
2239 : 141666 : 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 : 141666 : if (present(epaw_im)) then
2246 : 9204 : call pawaccenergy(epaw,pawrhoij,dij_2D,cplex_dij,qphase_dij,1,pawtab,epaw_im=epaw_im)
2247 : : else
2248 : 132462 : 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 [ + + ]: 141666 : 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 : : !!***
|
