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