Branch data Line data Source code
1 : : !!****m* ABINIT/m_paw_atom
2 : : !! NAME
3 : : !! m_paw_atom
4 : : !!
5 : : !! FUNCTION
6 : : !! atompaw related operations
7 : : !!
8 : : !! COPYRIGHT
9 : : !! Copyright (C) 2012-2026 ABINIT group (T. Rangel, MT, JWZ, GJ)
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 : : !! NOTES
15 : : !! FOR DEVELOPPERS: in order to preserve the portability of libPAW library,
16 : : !! please consult ~abinit/src/??_libpaw/libpaw-coding-rules.txt
17 : : !!
18 : : !! SOURCE
19 : :
20 : : #include "libpaw.h"
21 : :
22 : : module m_paw_atom
23 : :
24 : : USE_DEFS
25 : : USE_MSG_HANDLING
26 : : USE_MEMORY_PROFILING
27 : :
28 : : use m_paw_numeric, only : paw_jbessel, paw_solvbes, paw_spline, paw_splint
29 : : use m_pawrad, only : pawrad_type, simp_gen, poisson, pawrad_deducer0, bound_deriv, pawrad_ifromr
30 : : use m_pawtab, only : pawtab_type
31 : :
32 : : implicit none
33 : :
34 : : private
35 : :
36 : : public:: atompaw_shpfun
37 : : public:: atompaw_shapebes
38 : : public:: atompaw_vhnzc
39 : : public:: atompaw_dij0
40 : : public:: atompaw_kij
41 : : public:: atompaw_ehnzc
42 : : !!***
43 : :
44 : : CONTAINS !===========================================================
45 : : !!***
46 : :
47 : : !!****f* m_paw_atom/atompaw_shpfun
48 : : !! NAME
49 : : !! atompaw_shpfun
50 : : !!
51 : : !! FUNCTION
52 : : !! Compute shape function used in the definition
53 : : !! of compensation density (PAW)
54 : : !!
55 : : !! INPUTS
56 : : !! ll= l quantum number
57 : : !! mesh <type(pawrad_type)>=data containing radial grid information
58 : : !! pawtab <type(pawtab_type)>=paw tabulated starting data
59 : : !!
60 : : !! OUTPUT
61 : : !! norm= factor for shape function normalization
62 : : !!
63 : : !! SIDE effects
64 : : !! shapefunc(:)=shape function g(r)
65 : : !! In case of numerical shape function (shape_type=-1), shapefunc
66 : : !! array contains the shape function read in psp file at input.
67 : : !!
68 : : !! NOTES
69 : : !! Types of shape functions:
70 : : !! type -1: numerical shape function, given in psp file
71 : : !! type 1: g(r)=k(r).r^l; k(r)=exp(-(r/sigma)^lambda)
72 : : !! type 2: g(r)=k(r).r^l; k(r)=[sin(Pi.r/rshp)/(Pi.r/rshp)]^2
73 : : !! type 3: g(r)=alpha1.jl(q1.r)+alpha2.jl(q2.r)
74 : : !!
75 : : !! SOURCE
76 : :
77 [ - + ]: 3971 : subroutine atompaw_shpfun(ll,mesh,norm,pawtab,shapefunc)
78 : :
79 : : !Arguments ---------------------------------------------
80 : : !scalars
81 : : integer,intent(in) :: ll
82 : : real(dp),intent(out) :: norm
83 : : type(pawrad_type),intent(in) :: mesh
84 : : type(pawtab_type),intent(in) :: pawtab
85 : : !arrays
86 : : real(dp),intent(inout) :: shapefunc(:)
87 : :
88 : : !Local variables ------------------------------
89 : : !scalars
90 : : integer :: ir,ishp,mesh_size
91 : : real(dp) :: arg,besp,bespp,jbes1,jbes2
92 : : !arrays
93 : : real(dp) :: alpha(2),qq(2)
94 : 3971 : real(dp),allocatable :: r2k(:)
95 : : !no_abirules
96 : :
97 : : !***************************************************************************
98 : :
99 : 3971 : mesh_size=size(shapefunc)
100 [ - + ]: 3971 : if (mesh_size>mesh%mesh_size) then
101 : 0 : LIBPAW_BUG('wrong size!')
102 : : end if
103 : :
104 : : !Index for shape function cut-off radius
105 : 3971 : ishp=pawrad_ifromr(mesh,pawtab%rshp)-1
106 : :
107 : : !Computation of non-normalized shape function
108 [ + + ]: 3971 : if (pawtab%shape_type==-1) then
109 [ + + ]: 9424 : shapefunc(1:ishp)=pawtab%shapefunc(1:ishp,1+ll)
110 [ + + ]: 3955 : else if (pawtab%shape_type==1) then
111 [ + + ]: 16 : if (ll==0) then
112 : 8 : shapefunc(1)=one
113 [ + + ]: 3008 : do ir=2,ishp
114 : 3008 : shapefunc(ir)=exp(-(mesh%rad(ir)/pawtab%shape_sigma)**pawtab%shape_lambda)
115 : : end do
116 : : else
117 : 8 : shapefunc(1)=zero
118 [ + + ]: 3008 : do ir=2,ishp
119 : 3008 : shapefunc(ir)=exp(-(mesh%rad(ir)/pawtab%shape_sigma)**pawtab%shape_lambda)*mesh%rad(ir)**ll
120 : : end do
121 : : end if
122 [ + + ]: 3939 : else if (pawtab%shape_type==2) then
123 [ + + ]: 3223 : if (ll==0) then
124 : 1371 : shapefunc(1)=one
125 [ + + ]: 1371100 : do ir=2,ishp
126 : 1369729 : arg=pi*mesh%rad(ir)/pawtab%rshp
127 : 1371100 : shapefunc(ir)=(sin(arg)/arg)**2
128 : : end do
129 : : else
130 : 1852 : shapefunc(1)=zero
131 [ + + ]: 1683448 : do ir=2,ishp
132 : 1681596 : arg=pi*mesh%rad(ir)/pawtab%rshp
133 : 1683448 : shapefunc(ir)=(sin(arg)/(arg))**2 *mesh%rad(ir)**ll
134 : : end do
135 : : end if
136 [ + - ]: 716 : else if (pawtab%shape_type==3) then
137 [ + + ]: 2148 : alpha(1:2)=pawtab%shape_alpha(1:2,1+ll)
138 [ + + ]: 2148 : qq(1:2)=pawtab%shape_q(1:2,1+ll)
139 [ + + ]: 640577 : do ir=1,ishp
140 : 639861 : call paw_jbessel(jbes1,besp,bespp,ll,0,qq(1)*mesh%rad(ir))
141 : 639861 : call paw_jbessel(jbes2,besp,bespp,ll,0,qq(2)*mesh%rad(ir))
142 : 640577 : shapefunc(ir)=alpha(1)*jbes1+alpha(2)*jbes2
143 : : end do
144 : : end if
145 : :
146 [ + - + + ]: 461213 : if (ishp<mesh_size) shapefunc(ishp+1:mesh_size)=zero
147 : :
148 : : !Shape function normalization
149 [ + + + + ]: 3971 : if (pawtab%shape_type==-1.or.pawtab%shape_type==1.or.pawtab%shape_type==2) then
150 [ - + + - : 9765 : LIBPAW_ALLOCATE(r2k,(mesh_size))
+ - + - -
+ - + ]
151 [ + + ]: 3464188 : r2k=zero
152 [ + + ]: 3069972 : r2k(2:ishp)=shapefunc(2:ishp)*mesh%rad(2:ishp)**(2+ll)
153 [ + + ]: 3255 : if (mesh%mesh_type==5) then
154 : 16 : call simp_gen(norm,r2k,mesh);norm=one/norm
155 : : else
156 : 3239 : call simp_gen(norm,r2k,mesh,r_for_intg=pawtab%rshp);norm=one/norm
157 : : end if
158 [ + + ]: 3073227 : shapefunc(1:ishp)=shapefunc(1:ishp)*norm
159 [ + + ]: 3255 : if (pawtab%shape_type==-1) norm=one
160 [ - + ]: 3255 : LIBPAW_DEALLOCATE(r2k)
161 [ + - ]: 716 : else if (pawtab%shape_type==3) then
162 : 716 : norm=one
163 : : end if
164 : :
165 [ - + ]: 3971 : end subroutine atompaw_shpfun
166 : : !!***
167 : :
168 : : !----------------------------------------------------------------------
169 : :
170 : : !!****f* m_paw_atom/atompaw_atompaw_shapebes
171 : : !! NAME
172 : : !! atompaw_shapebes
173 : : !!
174 : : !! FUNCTION
175 : : !! Find al and ql parameters for a "Bessel" shape function:
176 : : !! Shape(r)=al1.jl(ql1.r)+al2.jl(ql2.r)
177 : : !! such as Shape(r) and 2 derivatives are zero at r=rc
178 : : !! Intg_0_rc[Shape(r).r^(l+2).dr]=1
179 : : !!
180 : : !! INPUTS
181 : : !! ll= l quantum number
182 : : !! rc= cut-off radius
183 : : !!
184 : : !! OUTPUT
185 : : !! al(2)= al coefficients
186 : : !! ql(2)= ql factors
187 : : !!
188 : : !! SOURCE
189 : :
190 : 430 : subroutine atompaw_shapebes(al,ql,ll,rc)
191 : :
192 : : !Arguments ------------------------------------
193 : : !scalars
194 : : integer :: ll
195 : : real(dp) :: rc
196 : : !arrays
197 : : real(dp) :: al(2),ql(2)
198 : :
199 : : !Local variables-------------------------------
200 : : !scalars
201 : : integer :: ii
202 : : real(dp) :: alpha,beta,det,jbes,jbesp,jbespp,qr
203 : : !arrays
204 : : real(dp) :: amat(2,2),bb(2)
205 : :
206 : : ! *************************************************************************
207 : :
208 : 430 : alpha=1._dp;beta=0._dp
209 : 430 : call paw_solvbes(ql,alpha,beta,ll,2)
210 [ + + ]: 1290 : ql(1:2)=ql(1:2)/rc
211 : :
212 [ + + ]: 1290 : do ii=1,2
213 : 860 : qr=ql(ii)*rc
214 : 860 : call paw_jbessel(jbes,jbesp,jbespp,ll,1,qr)
215 : 860 : amat(1,ii)=jbesp*ql(ii)
216 : 860 : call paw_jbessel(jbes,jbesp,jbespp,ll+1,0,qr)
217 : 1290 : amat(2,ii)=jbes*rc**(ll+2)/ql(ii) ! Intg_0_rc[jl(qr).r^(l+2).dr]
218 : : end do
219 : :
220 : 430 : bb(1)=zero;bb(2)=one
221 : :
222 : 430 : det=amat(1,1)*amat(2,2)-amat(1,2)*amat(2,1)
223 : 430 : al(1)=(amat(2,2)*bb(1)-amat(1,2)*bb(2))/det
224 : 430 : al(2)=(amat(1,1)*bb(2)-amat(2,1)*bb(1))/det
225 : :
226 : 430 : end subroutine atompaw_shapebes
227 : : !!***
228 : :
229 : : !----------------------------------------------------------------------
230 : :
231 : : !!****f* m_paw_atom/atompaw_vhnzc
232 : : !! NAME
233 : : !! atompaw_vhnzc
234 : : !!
235 : : !! FUNCTION
236 : : !! PAW: compute Hartree potential for n_{Zc}
237 : : !!
238 : : !! INPUTS
239 : : !! ncore(:)=atomic core density
240 : : !! radmesh_core <type(pawrad_type)>=radial mesh (and related data) for the core densities
241 : : !! znucl= valence and total charge of the atomic species
242 : : !!
243 : : !! OUTPUT
244 : : !! vhnzc(:)=Hartree potential due to Z_nc
245 : : !!
246 : : !! SOURCE
247 : :
248 [ - + - + ]: 1465 : subroutine atompaw_vhnzc(ncore,radmesh_core,vhnzc,znucl)
249 : :
250 : : !Arguments ---------------------------------------------
251 : : !scalars
252 : : real(dp),intent(in) :: znucl
253 : : type(pawrad_type),intent(in) :: radmesh_core
254 : : !arrays
255 : : real(dp),intent(in) :: ncore(:)
256 : : real(dp), intent(out) :: vhnzc(:)
257 : :
258 : : !Local variables ---------------------------------------
259 : : integer :: mesh_size
260 : : real(dp),allocatable :: nwk(:)
261 : :
262 : : ! *********************************************************************
263 : :
264 : 1465 : mesh_size=size(ncore)
265 [ + - - + ]: 1465 : if (mesh_size/=size(vhnzc).or.mesh_size>radmesh_core%mesh_size) then
266 : 0 : LIBPAW_BUG('wrong sizes!')
267 : : end if
268 : :
269 [ - + + - : 4395 : LIBPAW_ALLOCATE(nwk,(mesh_size))
+ - + - -
+ - + ]
270 : :
271 [ + + ]: 2017729 : nwk(:)=ncore(:)*four_pi*radmesh_core%rad(:)**2
272 : 1465 : call poisson(nwk,0,radmesh_core,vhnzc)
273 [ + + ]: 2016264 : vhnzc(2:mesh_size)=(vhnzc(2:mesh_size)-znucl)/radmesh_core%rad(2:mesh_size)
274 [ - + - - : 1465 : call pawrad_deducer0(vhnzc,mesh_size,radmesh_core)
- - - + -
- ]
275 : :
276 [ - + ]: 1465 : LIBPAW_DEALLOCATE(nwk)
277 : :
278 : 1465 : end subroutine atompaw_vhnzc
279 : : !!***
280 : :
281 : : !----------------------------------------------------------------------
282 : :
283 : : !!****f* m_paw_atom/atompaw_dij0
284 : : !! NAME
285 : : !! atompaw_dij0
286 : : !!
287 : : !! FUNCTION
288 : : !! PAW: Compute "frozen" values of pseudopotential strengths Dij = Dij0
289 : : !!
290 : : !! INPUTS
291 : : !! indlmn(6,lmnmax)= array giving l,m,n,lm,ln,s for i=lmn
292 : : !! kij(pawtab%lmn2_size)= kinetic part of Dij
293 : : !! lmnmax=max number of (l,m,n) components over all type of psps
294 : : !! ncore(:)=atomic core density
295 : : !! opt_init=flag defining the storage of PAW atomic data
296 : : !! 0: PAW atomic data have not been initialized (in pawtab)
297 : : !! 1: PAW atomic data have been initialized (in pawtab)
298 : : !! pawtab <type(pawtab_type)>=paw tabulated starting data
299 : : !! radmesh <type(pawrad_type)>=paw radial mesh (and related data)
300 : : !! radmesh_core <type(pawrad_type)>=radial mesh (and related data) for the core densities
301 : : !! radmesh_vloc <type(pawrad_type)>=radial mesh (and related data) for the local potential (VH(tnZc))
302 : : !! vhtnzc(:)= local potential VH(tnZc)
303 : : !! znucl= valence and total charge of the atomic species
304 : : !!
305 : : !! OUTPUT
306 : : !! pawtab%dij0(pawtab%lmn2_size)= Frozen part of the Dij term
307 : : !!
308 : : !! SOURCE
309 : :
310 : :
311 : 1008 : subroutine atompaw_dij0(indlmn,kij,lmnmax,ncore,opt_init,pawtab,&
312 [ - + - + ]: 336 : & radmesh,radmesh_core,radmesh_vloc,vhtnzc,znucl)
313 : :
314 : : !Arguments ---------------------------------------------
315 : : !scalars
316 : : integer,intent(in) :: lmnmax,opt_init
317 : : real(dp),intent(in) :: znucl
318 : : type(pawrad_type),intent(in) :: radmesh,radmesh_core,radmesh_vloc
319 : : type(pawtab_type),intent(inout) :: pawtab
320 : : !arrays
321 : : integer,intent(in) :: indlmn(6,lmnmax)
322 : : real(dp),intent(in) :: kij(pawtab%lmn2_size)
323 : : real(dp),intent(in) :: ncore(:),vhtnzc(:)
324 : : !real(dp),optional,intent(in) :: vminushalf(:)
325 : :
326 : : !Local variables ---------------------------------------
327 : : integer :: il,ilm,iln,ilmn,j0lmn,jl,jlm,jln,jlmn,klmn,lmn2_size,meshsz,meshsz_core
328 : : integer :: meshsz_vhtnzc,meshsz_vmh
329 : : real(dp) :: intg,intvh,yp1,ypn
330 : 336 : real(dp),allocatable :: ff(:),r2k(:),shpf(:),vhnzc(:),vhtnzc_sph(:),work1(:),work2(:)
331 : :
332 : : ! *********************************************************************
333 : :
334 : 336 : lmn2_size=pawtab%lmn2_size
335 : 336 : meshsz_vhtnzc=size(vhtnzc)
336 : 336 : meshsz=min(radmesh%mesh_size,radmesh_core%mesh_size,radmesh_vloc%mesh_size,meshsz_vhtnzc)
337 [ - + + - : 1008 : LIBPAW_ALLOCATE(ff,(meshsz))
+ - + - -
+ - + ]
338 : :
339 : : !Retrieve VH(tnZc) on the correct radial mesh
340 [ - + + - : 672 : LIBPAW_ALLOCATE(vhtnzc_sph,(meshsz))
- + - + ]
341 : : if ((radmesh%mesh_type/=radmesh_vloc%mesh_type).or.&
342 [ + - + - : 336 : & (radmesh%rstep /=radmesh_vloc%rstep) .or.&
+ - ]
343 : : & (radmesh%lstep /=radmesh_vloc%lstep)) then
344 [ # # # # : 0 : call bound_deriv(vhtnzc,radmesh_vloc,meshsz_vhtnzc,yp1,ypn)
# # # # ]
345 [ # # # # : 0 : LIBPAW_ALLOCATE(work1,(meshsz_vhtnzc))
# # # # #
# # # ]
346 [ # # # # : 0 : LIBPAW_ALLOCATE(work2,(meshsz_vhtnzc))
# # # # ]
347 [ # # # # : 0 : call paw_spline(radmesh_vloc%rad,vhtnzc,meshsz_vhtnzc,yp1,ypn,work1)
# # # # ]
348 [ # # # # : 0 : call paw_splint(meshsz_vhtnzc,radmesh_vloc%rad,vhtnzc,work1,meshsz,radmesh%rad(1:meshsz),vhtnzc_sph)
# # # # ]
349 [ # # ]: 0 : LIBPAW_DEALLOCATE(work1)
350 [ # # ]: 0 : LIBPAW_DEALLOCATE(work2)
351 : : else
352 [ + + ]: 572330 : vhtnzc_sph(1:meshsz)=vhtnzc(1:meshsz)
353 : : end if
354 : :
355 : : !Kinetic part of Dij0
356 : : !====================
357 [ + + ]: 37811 : pawtab%dij0(1:lmn2_size)=kij(1:lmn2_size)
358 : :
359 : : !Computation of <phi_i|vh(nZc)|phi_j> on the PAW sphere
360 : : !======================================================
361 : 336 : meshsz_core=size(ncore)
362 [ - + + - : 1008 : LIBPAW_ALLOCATE(vhnzc,(meshsz_core))
+ - + - -
+ - + ]
363 : 336 : call atompaw_vhnzc(ncore,radmesh_core,vhnzc,znucl)
364 [ + + ]: 4209 : do jlmn=1,pawtab%lmn_size
365 : 3873 : j0lmn=jlmn*(jlmn-1)/2
366 : 3873 : jlm=indlmn(4,jlmn);jln=indlmn(5,jlmn)
367 [ + + ]: 41684 : do ilmn=1,jlmn
368 : 37475 : klmn=j0lmn+ilmn
369 : 37475 : ilm=indlmn(4,ilmn);iln=indlmn(5,ilmn)
370 [ + + ]: 41348 : if (jlm==ilm) then
371 [ + + ]: 10775568 : ff(1:meshsz)=pawtab%phi(1:meshsz,iln)*pawtab%phi(1:meshsz,jln)*vhnzc(1:meshsz)
372 : 6168 : call simp_gen(intg,ff,radmesh)
373 : 6168 : pawtab%dij0(klmn)=pawtab%dij0(klmn)+intg
374 : : end if
375 : : end do
376 : : end do
377 [ - + ]: 336 : LIBPAW_DEALLOCATE(vhnzc)
378 : :
379 : : !Computation of -<tphi_i|vh(tnZc)|tphi_j> on the PAW sphere
380 : : !==========================================================
381 [ + + ]: 4209 : do jlmn=1,pawtab%lmn_size
382 : 3873 : j0lmn=jlmn*(jlmn-1)/2
383 : 3873 : jlm=indlmn(4,jlmn);jln=indlmn(5,jlmn)
384 [ + + ]: 41684 : do ilmn=1,jlmn
385 : 37475 : klmn=j0lmn+ilmn
386 : 37475 : ilm=indlmn(4,ilmn);iln=indlmn(5,ilmn)
387 [ + + ]: 41348 : if (jlm==ilm) then
388 [ + + ]: 10775568 : ff(1:meshsz)=pawtab%tphi(1:meshsz,iln)*pawtab%tphi(1:meshsz,jln)*vhtnzc_sph(1:meshsz)
389 : 6168 : call simp_gen(intg,ff,radmesh)
390 : 6168 : pawtab%dij0(klmn)=pawtab%dij0(klmn)-intg
391 : : end if
392 : : end do
393 : : end do
394 : :
395 : : !Computation of <phi_i|vminushalf|phi_j> (if any)
396 : : !=================================================
397 [ + + ]: 336 : if(pawtab%has_vminushalf==1) then
398 [ + - ]: 1 : if(size(pawtab%vminushalf)>=1) then
399 : 1 : meshsz_vmh=min(meshsz,size(pawtab%vminushalf))
400 [ + + ]: 9 : do jlmn=1,pawtab%lmn_size
401 : 8 : j0lmn=jlmn*(jlmn-1)/2
402 : 8 : jlm=indlmn(4,jlmn);jln=indlmn(5,jlmn)
403 [ + + ]: 45 : do ilmn=1,jlmn
404 : 36 : klmn=j0lmn+ilmn
405 : 36 : ilm=indlmn(4,ilmn);iln=indlmn(5,ilmn)
406 [ + + ]: 44 : if (jlm==ilm) then
407 [ + + ]: 24024 : ff(1:meshsz_vmh)=pawtab%phi(1:meshsz_vmh,iln)*pawtab%phi(1:meshsz_vmh,jln)*pawtab%vminushalf(1:meshsz_vmh)
408 : 12 : call simp_gen(intg,ff(1:meshsz_vmh),radmesh)
409 : 12 : pawtab%dij0(klmn)=pawtab%dij0(klmn)+intg
410 : : end if
411 : : end do
412 : : end do
413 : : end if
414 : : end if
415 : :
416 : : !Computation of -int[vh(tnzc)*Qijhat(r)dr]
417 : : !=========================================
418 [ + - ]: 336 : if (opt_init==0) then
419 [ - + + - : 1008 : LIBPAW_ALLOCATE(shpf,(radmesh%mesh_size))
+ - + - -
+ - + ]
420 : 336 : call atompaw_shpfun(0,radmesh,intg,pawtab,shpf)
421 [ + + ]: 336 : if (pawtab%shape_type==3) then
422 [ - + + - : 165 : LIBPAW_ALLOCATE(r2k,(radmesh%int_meshsz))
+ - + - -
+ - + ]
423 [ + + ]: 58071 : r2k=zero
424 [ + + ]: 58016 : r2k(2:radmesh%int_meshsz)=shpf(2:radmesh%int_meshsz)*radmesh%rad(2:radmesh%int_meshsz)**2
425 [ - + ]: 55 : if(radmesh%mesh_type==5) then
426 : 0 : call simp_gen(intg,r2k,radmesh)
427 : : else
428 : 55 : call simp_gen(intg,r2k,radmesh,r_for_intg=pawtab%rshp)
429 : : end if
430 [ + + ]: 83459 : shpf(1:meshsz)=shpf(1:meshsz)/intg
431 [ - + ]: 55 : LIBPAW_DEALLOCATE(r2k)
432 : : end if
433 [ + + ]: 572330 : ff(1:meshsz)=vhtnzc_sph(1:meshsz)*shpf(1:meshsz)*radmesh%rad(1:meshsz)**2
434 [ - + ]: 336 : LIBPAW_DEALLOCATE(shpf)
435 : 336 : call simp_gen(intvh,ff,radmesh)
436 [ + + ]: 4209 : do jlmn=1,pawtab%lmn_size
437 : 3873 : j0lmn=jlmn*(jlmn-1)/2
438 : 3873 : jl=indlmn(1,jlmn);jln=indlmn(5,jlmn);jlm=indlmn(4,jlmn)
439 [ + + ]: 41684 : do ilmn=1,jlmn
440 : 37475 : klmn=j0lmn+ilmn
441 : 37475 : il=indlmn(1,ilmn);iln=indlmn(5,ilmn);ilm=indlmn(4,ilmn)
442 [ + + ]: 41348 : if (ilm==jlm) then
443 : : ff(1:meshsz)=(pawtab%phi (1:meshsz,iln)*pawtab%phi (1:meshsz,jln)&
444 [ + + ]: 10775568 : & -pawtab%tphi(1:meshsz,iln)*pawtab%tphi(1:meshsz,jln))
445 : 6168 : call simp_gen(intg,ff,radmesh)
446 : 6168 : pawtab%dij0(klmn)=pawtab%dij0(klmn)-intvh*intg
447 : : end if
448 : : end do
449 : : end do
450 : : else
451 [ # # ]: 0 : ff(1:meshsz)=vhtnzc_sph(1:meshsz)*pawtab%shapefunc(1:meshsz,1)*radmesh%rad(1:meshsz)**2
452 : 0 : call simp_gen(intvh,ff,radmesh)
453 [ # # ]: 0 : do jlmn=1,pawtab%lmn_size
454 : 0 : j0lmn=jlmn*(jlmn-1)/2
455 : 0 : jl=indlmn(1,jlmn);jln=indlmn(5,jlmn);jlm=indlmn(4,jlmn)
456 [ # # ]: 0 : do ilmn=1,jlmn
457 : 0 : klmn=j0lmn+ilmn
458 : 0 : il=indlmn(1,ilmn);iln=indlmn(5,ilmn);ilm=indlmn(4,ilmn)
459 [ # # ]: 0 : if (ilm==jlm) then
460 : 0 : intg=pawtab%qijl(1,klmn)*sqrt(four_pi)
461 : 0 : pawtab%dij0(klmn)=pawtab%dij0(klmn)-intvh*intg
462 : : end if
463 : : end do
464 : : end do
465 : : end if
466 : :
467 [ - + ]: 336 : LIBPAW_DEALLOCATE(ff)
468 [ - + ]: 336 : LIBPAW_DEALLOCATE(vhtnzc_sph)
469 : :
470 [ - + - + ]: 336 : end subroutine atompaw_dij0
471 : : !!***
472 : :
473 : : !----------------------------------------------------------------------
474 : :
475 : : !!****f* m_paw_atom/atompaw_kij
476 : : !! NAME
477 : : !! atompaw_kij
478 : : !!
479 : : !! FUNCTION
480 : : !! PAW: deduce kinetic part of psp strength (Dij) from the knowledge of frozen Dij (Dij0)
481 : : !!
482 : : !! INPUTS
483 : : !! indlmn(6,lmnmax)= array giving l,m,n,lm,ln,s for i=lmn
484 : : !! lmnmax=max number of (l,m,n) components over all type of psps
485 : : !! ncore(:)=atomic core density
486 : : !! opt_init=flag defining the storage of PAW atomic data
487 : : !! 0: PAW atomic data have not been initialized (in pawtab)
488 : : !! 1: PAW atomic data have been initialized (in pawtab)
489 : : !! opt_vhnzc=flag defining the inclusion of VH(nZc) in computation
490 : : !! 0: VH(nZc) is not taken into account
491 : : !! 1: VH(nZc) is taken into account
492 : : !! pawtab <type(pawtab_type)>=paw tabulated starting data
493 : : !! radmesh <type(pawrad_type)>=paw radial mesh (and related data)
494 : : !! radmesh_core <type(pawrad_type)>=radial mesh (and related data) for the core densities
495 : : !! radmesh_vloc <type(pawrad_type)>=radial mesh (and related data) for the local potential (VH(tnZc))
496 : : !! vhtnzc(:)= local potential VH(tnZc)
497 : : !! znucl= valence and total charge of the atomic species
498 : : !!
499 : : !! OUTPUT
500 : : !! kij(pawtab%lmn2_size)= kinetic part of Dij
501 : : !!
502 : : !! SOURCE
503 : :
504 : 66 : subroutine atompaw_kij(indlmn,kij,lmnmax,ncore,opt_init,opt_vhnzc,pawtab, &
505 [ - + - + ]: 22 : & radmesh,radmesh_core,radmesh_vloc,vhtnzc,znucl)
506 : :
507 : : !Arguments ---------------------------------------------
508 : : !scalars
509 : : integer,intent(in) :: lmnmax,opt_init,opt_vhnzc
510 : : real(dp),intent(in) :: znucl
511 : : type(pawrad_type),intent(in) :: radmesh,radmesh_core,radmesh_vloc
512 : : type(pawtab_type),intent(in) :: pawtab
513 : : !arrays
514 : : integer,intent(in) :: indlmn(6,lmnmax)
515 : : real(dp),intent(out) :: kij(pawtab%lmn2_size)
516 : : real(dp),intent(in) :: ncore(:)
517 : : real(dp),intent(in) :: vhtnzc(:)
518 : :
519 : : !Local variables ---------------------------------------
520 : : integer :: il,ilm,iln,ilmn,j0lmn,jl,jlm,jln,jlmn,klmn,lmn2_size
521 : : integer :: meshsz,meshsz_core,meshsz_vhtnzc,meshsz_vmh
522 : : real(dp) :: intg,intvh,yp1,ypn
523 : 22 : real(dp),allocatable :: ff(:),r2k(:),shpf(:),vhnzc(:),vhtnzc_sph(:),work1(:),work2(:)
524 : :
525 : : ! *********************************************************************
526 : :
527 : 22 : lmn2_size=pawtab%lmn2_size
528 : 22 : meshsz_vhtnzc=size(vhtnzc)
529 : 22 : meshsz=min(radmesh%mesh_size,radmesh_core%mesh_size,radmesh_vloc%mesh_size,meshsz_vhtnzc)
530 [ - + + - : 66 : LIBPAW_ALLOCATE(ff,(meshsz))
+ - + - -
+ - + ]
531 : :
532 : : !Retrieve VH(tnZc) on the correct radial mesh
533 [ - + + - : 44 : LIBPAW_ALLOCATE(vhtnzc_sph,(meshsz))
- + - + ]
534 : : if ((radmesh%mesh_type/=radmesh_vloc%mesh_type).or.&
535 [ + - + - : 22 : & (radmesh%rstep /=radmesh_vloc%rstep) .or.&
+ - ]
536 : : & (radmesh%lstep /=radmesh_vloc%lstep)) then
537 [ # # # # : 0 : call bound_deriv(vhtnzc,radmesh_vloc,meshsz_vhtnzc,yp1,ypn)
# # # # ]
538 [ # # # # : 0 : LIBPAW_ALLOCATE(work1,(meshsz_vhtnzc))
# # # # #
# # # ]
539 [ # # # # : 0 : LIBPAW_ALLOCATE(work2,(meshsz_vhtnzc))
# # # # ]
540 [ # # # # : 0 : call paw_spline(radmesh_vloc%rad,vhtnzc,meshsz_vhtnzc,yp1,ypn,work1)
# # # # ]
541 [ # # # # : 0 : call paw_splint(meshsz_vhtnzc,radmesh_vloc%rad,vhtnzc,work1,meshsz,radmesh%rad(1:meshsz),vhtnzc_sph)
# # # # ]
542 [ # # ]: 0 : LIBPAW_DEALLOCATE(work1)
543 [ # # ]: 0 : LIBPAW_DEALLOCATE(work2)
544 : : else
545 [ + + ]: 12823 : vhtnzc_sph(1:meshsz)=vhtnzc(1:meshsz)
546 : : end if
547 : :
548 : : !Initialize Kij with Dij0
549 : : !=========================
550 [ + + ]: 1084 : kij(1:lmn2_size)=pawtab%dij0(1:lmn2_size)
551 : :
552 : : !Substraction of -<phi_i|vh(nZc)|phi_j> on the PAW sphere
553 : : !========================================================
554 [ + + ]: 22 : if (opt_vhnzc/=0) then
555 : 5 : meshsz_core=size(ncore)
556 [ - + + - : 15 : LIBPAW_ALLOCATE(vhnzc,(meshsz_core))
+ - + - -
+ - + ]
557 : 5 : call atompaw_vhnzc(ncore,radmesh_core,vhnzc,znucl)
558 [ + + ]: 45 : do jlmn=1,pawtab%lmn_size
559 : 40 : j0lmn=jlmn*(jlmn-1)/2
560 : 40 : jlm=indlmn(4,jlmn);jln=indlmn(5,jlmn)
561 [ + + ]: 225 : do ilmn=1,jlmn
562 : 180 : klmn=j0lmn+ilmn
563 : 180 : ilm=indlmn(4,ilmn);iln=indlmn(5,ilmn)
564 [ + + ]: 220 : if (jlm==ilm) then
565 [ + + ]: 36300 : ff(1:meshsz)=pawtab%phi(1:meshsz,iln)*pawtab%phi(1:meshsz,jln)*vhnzc(1:meshsz)
566 : 60 : call simp_gen(intg,ff,radmesh)
567 : 60 : kij(klmn)=kij(klmn)-intg
568 : : end if
569 : : end do
570 : : end do
571 [ - + ]: 5 : LIBPAW_DEALLOCATE(vhnzc)
572 : : end if
573 : :
574 : : !Substraction of <tphi_i|vh(tnZc)|tphi_j> on the PAW sphere
575 : : !==========================================================
576 [ + + ]: 218 : do jlmn=1,pawtab%lmn_size
577 : 196 : j0lmn=jlmn*(jlmn-1)/2
578 : 196 : jlm=indlmn(4,jlmn);jln=indlmn(5,jlmn)
579 [ + + ]: 1280 : do ilmn=1,jlmn
580 : 1062 : klmn=j0lmn+ilmn
581 : 1062 : ilm=indlmn(4,ilmn);iln=indlmn(5,ilmn)
582 [ + + ]: 1258 : if (jlm==ilm) then
583 [ + + ]: 193866 : ff(1:meshsz)=pawtab%tphi(1:meshsz,iln)*pawtab%tphi(1:meshsz,jln)*vhtnzc_sph(1:meshsz)
584 : 294 : call simp_gen(intg,ff,radmesh)
585 : 294 : kij(klmn)=kij(klmn)+intg
586 : : end if
587 : : end do
588 : : end do
589 : :
590 : : !Computation of <phi_i|vminushalf|phi_j> (if any)
591 : : !=================================================
592 [ - + ]: 22 : if(pawtab%has_vminushalf==1) then
593 [ # # ]: 0 : if(size(pawtab%vminushalf)>=1) then
594 : 0 : meshsz_vmh=min(meshsz,size(pawtab%vminushalf))
595 [ # # ]: 0 : do jlmn=1,pawtab%lmn_size
596 : 0 : j0lmn=jlmn*(jlmn-1)/2
597 : 0 : jlm=indlmn(4,jlmn);jln=indlmn(5,jlmn)
598 [ # # ]: 0 : do ilmn=1,jlmn
599 : 0 : klmn=j0lmn+ilmn
600 : 0 : ilm=indlmn(4,ilmn);iln=indlmn(5,ilmn)
601 [ # # ]: 0 : if (jlm==ilm) then
602 [ # # ]: 0 : ff(1:meshsz_vmh)=pawtab%phi(1:meshsz_vmh,iln)*pawtab%phi(1:meshsz_vmh,jln)*pawtab%vminushalf(1:meshsz_vmh)
603 : 0 : call simp_gen(intg,ff(1:meshsz_vmh),radmesh)
604 : 0 : kij(klmn)=kij(klmn)-intg
605 : : end if
606 : : end do
607 : : end do
608 : : end if
609 : : end if
610 : :
611 : : !Computation of int[vh(tnzc)*Qijhat(r)dr]
612 : : !==========================================
613 [ + - ]: 22 : if (opt_init==0) then
614 [ - + + - : 66 : LIBPAW_ALLOCATE(shpf,(radmesh%mesh_size))
+ - + - -
+ - + ]
615 : 22 : call atompaw_shpfun(0,radmesh,intg,pawtab,shpf)
616 [ - + ]: 22 : if (pawtab%shape_type==3) then
617 [ # # # # : 0 : LIBPAW_ALLOCATE(r2k,(radmesh%int_meshsz))
# # # # #
# # # ]
618 [ # # ]: 0 : r2k=zero
619 [ # # ]: 0 : r2k(2:radmesh%int_meshsz)=shpf(2:radmesh%int_meshsz)*radmesh%rad(2:radmesh%int_meshsz)**2
620 [ # # ]: 0 : if(radmesh%mesh_type==5) then
621 : 0 : call simp_gen(intg,r2k,radmesh)
622 : : else
623 : 0 : call simp_gen(intg,r2k,radmesh,r_for_intg=pawtab%rshp)
624 : : end if
625 [ # # ]: 0 : shpf(1:meshsz)=shpf(1:meshsz)/intg
626 [ # # ]: 0 : LIBPAW_DEALLOCATE(r2k)
627 : : end if
628 [ + + ]: 12823 : ff(1:meshsz)=vhtnzc_sph(1:meshsz)*shpf(1:meshsz)*radmesh%rad(1:meshsz)**2
629 [ - + ]: 22 : LIBPAW_DEALLOCATE(shpf)
630 : 22 : call simp_gen(intvh,ff,radmesh)
631 [ + + ]: 218 : do jlmn=1,pawtab%lmn_size
632 : 196 : j0lmn=jlmn*(jlmn-1)/2
633 : 196 : jl=indlmn(1,jlmn);jln=indlmn(5,jlmn);jlm=indlmn(4,jlmn)
634 [ + + ]: 1280 : do ilmn=1,jlmn
635 : 1062 : klmn=j0lmn+ilmn
636 : 1062 : il=indlmn(1,ilmn);iln=indlmn(5,ilmn);ilm=indlmn(4,ilmn)
637 [ + + ]: 1258 : if (ilm==jlm) then
638 : : ff(1:meshsz)=(pawtab%phi (1:meshsz,iln)*pawtab%phi (1:meshsz,jln)&
639 [ + + ]: 193866 : & -pawtab%tphi(1:meshsz,iln)*pawtab%tphi(1:meshsz,jln))
640 : 294 : call simp_gen(intg,ff,radmesh)
641 : 294 : kij(klmn)=kij(klmn)+intvh*intg
642 : : end if
643 : : end do
644 : : end do
645 : : else
646 [ # # ]: 0 : ff(1:meshsz)=vhtnzc_sph(1:meshsz)*pawtab%shapefunc(1:meshsz,1)*radmesh%rad(1:meshsz)**2
647 : 0 : call simp_gen(intvh,ff,radmesh)
648 [ # # ]: 0 : do jlmn=1,pawtab%lmn_size
649 : 0 : j0lmn=jlmn*(jlmn-1)/2
650 : 0 : jl=indlmn(1,jlmn);jln=indlmn(5,jlmn);jlm=indlmn(4,jlmn)
651 [ # # ]: 0 : do ilmn=1,jlmn
652 : 0 : klmn=j0lmn+ilmn
653 : 0 : il=indlmn(1,ilmn);iln=indlmn(5,ilmn);ilm=indlmn(4,ilmn)
654 [ # # ]: 0 : if (ilm==jlm) then
655 : 0 : intg=pawtab%qijl(1,klmn)*sqrt(four_pi)
656 : 0 : kij(klmn)=kij(klmn)+intvh*intg
657 : : end if
658 : : end do
659 : : end do
660 : : end if
661 : :
662 [ - + ]: 22 : LIBPAW_DEALLOCATE(ff)
663 [ - + ]: 22 : LIBPAW_DEALLOCATE(vhtnzc_sph)
664 : :
665 [ - + - + : 22 : end subroutine atompaw_kij
- + ]
666 : : !!***
667 : :
668 : :
669 : : !----------------------------------------------------------------------
670 : :
671 : : !!****f* m_paw_atom/atompaw_ehnzc
672 : : !! NAME
673 : : !! atompaw_ehnzc
674 : : !!
675 : : !! FUNCTION
676 : : !! Computes the contribution of nZc to the Hartree energy
677 : : !!
678 : : !! INPUTS
679 : : !! ncore(:) = radial core density
680 : : !! radmesh_core<type(pawrad_type)> = radial mesh for core density
681 : : !! znucl = nucleus Z number
682 : : !!
683 : : !! OUTPUT
684 : : !! ehnzc = nZc Hartree energy
685 : : !!
686 : : !! SOURCE
687 : :
688 [ - + ]: 575 : subroutine atompaw_ehnzc(ncore,radmesh_core,ehnzc,znucl)
689 : :
690 : : !Arguments ---------------------------------------------
691 : : !scalars
692 : : real(dp),intent(out) :: ehnzc
693 : : real(dp),intent(in) :: znucl
694 : : type(pawrad_type),intent(in) :: radmesh_core
695 : : !arrays
696 : : real(dp),intent(in) :: ncore(:)
697 : :
698 : : !Local variables ---------------------------------------
699 : : !scalars
700 : : integer :: mesh_size
701 : : !arrays
702 : : real(dp),allocatable :: vhnzc(:),ff(:)
703 : :
704 : : ! *********************************************************************
705 : :
706 : 575 : mesh_size=size(ncore)
707 [ - + ]: 575 : if (mesh_size/=radmesh_core%mesh_size) then
708 : 0 : LIBPAW_BUG('wrong sizes!')
709 : : end if
710 [ - + + - : 1725 : LIBPAW_ALLOCATE(vhnzc,(mesh_size))
+ - + - -
+ - + ]
711 [ - + + - : 1150 : LIBPAW_ALLOCATE(ff,(mesh_size))
- + - + ]
712 : 575 : ehnzc=zero
713 [ + + ]: 746124 : vhnzc=zero
714 [ + + ]: 746124 : ff(:)=ncore(:)*four_pi*radmesh_core%rad(:)**2
715 : 575 : call poisson(ff,0,radmesh_core,vhnzc)
716 [ + + ]: 745549 : vhnzc(2:mesh_size)=(vhnzc(2:mesh_size)*half-znucl)
717 : 575 : call pawrad_deducer0(vhnzc,mesh_size,radmesh_core)
718 [ + + ]: 746124 : ff(:)=vhnzc(:)*ncore(:)*(four_pi*radmesh_core%rad(:))
719 [ - + ]: 575 : LIBPAW_DEALLOCATE(vhnzc)
720 : 575 : call simp_gen(ehnzc,ff,radmesh_core)
721 [ - + ]: 575 : LIBPAW_DEALLOCATE(ff)
722 : :
723 : 575 : end subroutine atompaw_ehnzc
724 : : !!***
725 : :
726 : : end module m_paw_atom
727 : : !!***
|
