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 3826 : 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 3826 : real(dp),allocatable :: r2k(:)
95 : !no_abirules
96 :
97 : !***************************************************************************
98 :
99 3826 : mesh_size=size(shapefunc)
100 3826 : 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 3826 : ishp=pawrad_ifromr(mesh,pawtab%rshp)-1
106 :
107 : !Computation of non-normalized shape function
108 3826 : if (pawtab%shape_type==-1) then
109 9424 : shapefunc(1:ishp)=pawtab%shapefunc(1:ishp,1+ll)
110 3810 : 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 3794 : else if (pawtab%shape_type==2) then
123 3078 : if (ll==0) then
124 1316 : shapefunc(1)=one
125 1307437 : do ir=2,ishp
126 1306121 : arg=pi*mesh%rad(ir)/pawtab%rshp
127 1307437 : shapefunc(ir)=(sin(arg)/arg)**2
128 : end do
129 : else
130 1762 : shapefunc(1)=zero
131 1604046 : do ir=2,ishp
132 1602284 : arg=pi*mesh%rad(ir)/pawtab%rshp
133 1604046 : 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 443248 : if (ishp<mesh_size) shapefunc(ishp+1:mesh_size)=zero
147 :
148 : !Shape function normalization
149 3826 : if (pawtab%shape_type==-1.or.pawtab%shape_type==1.or.pawtab%shape_type==2) then
150 9330 : LIBPAW_ALLOCATE(r2k,(mesh_size))
151 3303158 : r2k=zero
152 2926907 : r2k(2:ishp)=shapefunc(2:ishp)*mesh%rad(2:ishp)**(2+ll)
153 3110 : if (mesh%mesh_type==5) then
154 16 : call simp_gen(norm,r2k,mesh);norm=one/norm
155 : else
156 3094 : call simp_gen(norm,r2k,mesh,r_for_intg=pawtab%rshp);norm=one/norm
157 : end if
158 2930017 : shapefunc(1:ishp)=shapefunc(1:ishp)*norm
159 3110 : if (pawtab%shape_type==-1) norm=one
160 3110 : LIBPAW_DEALLOCATE(r2k)
161 716 : else if (pawtab%shape_type==3) then
162 716 : norm=one
163 : end if
164 :
165 3826 : 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 1421 : 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 1421 : mesh_size=size(ncore)
265 1421 : if (mesh_size/=size(vhnzc).or.mesh_size>radmesh_core%mesh_size) then
266 0 : LIBPAW_BUG('wrong sizes!')
267 : end if
268 :
269 4263 : LIBPAW_ALLOCATE(nwk,(mesh_size))
270 :
271 1941392 : nwk(:)=ncore(:)*four_pi*radmesh_core%rad(:)**2
272 1421 : call poisson(nwk,0,radmesh_core,vhnzc)
273 1939971 : vhnzc(2:mesh_size)=(vhnzc(2:mesh_size)-znucl)/radmesh_core%rad(2:mesh_size)
274 1421 : call pawrad_deducer0(vhnzc,mesh_size,radmesh_core)
275 :
276 1421 : LIBPAW_DEALLOCATE(nwk)
277 :
278 1421 : 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 972 : subroutine atompaw_dij0(indlmn,kij,lmnmax,ncore,opt_init,pawtab,&
312 324 : & 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 324 : real(dp),allocatable :: ff(:),r2k(:),shpf(:),vhnzc(:),vhtnzc_sph(:),work1(:),work2(:)
331 :
332 : ! *********************************************************************
333 :
334 324 : lmn2_size=pawtab%lmn2_size
335 324 : meshsz_vhtnzc=size(vhtnzc)
336 324 : meshsz=min(radmesh%mesh_size,radmesh_core%mesh_size,radmesh_vloc%mesh_size,meshsz_vhtnzc)
337 972 : LIBPAW_ALLOCATE(ff,(meshsz))
338 :
339 : !Retrieve VH(tnZc) on the correct radial mesh
340 648 : LIBPAW_ALLOCATE(vhtnzc_sph,(meshsz))
341 : if ((radmesh%mesh_type/=radmesh_vloc%mesh_type).or.&
342 324 : & (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 548207 : vhtnzc_sph(1:meshsz)=vhtnzc(1:meshsz)
353 : end if
354 :
355 : !Kinetic part of Dij0
356 : !====================
357 37232 : 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 324 : meshsz_core=size(ncore)
362 972 : LIBPAW_ALLOCATE(vhnzc,(meshsz_core))
363 324 : call atompaw_vhnzc(ncore,radmesh_core,vhnzc,znucl)
364 4091 : do jlmn=1,pawtab%lmn_size
365 3767 : j0lmn=jlmn*(jlmn-1)/2
366 3767 : jlm=indlmn(4,jlmn);jln=indlmn(5,jlmn)
367 40999 : do ilmn=1,jlmn
368 36908 : klmn=j0lmn+ilmn
369 36908 : ilm=indlmn(4,ilmn);iln=indlmn(5,ilmn)
370 40675 : if (jlm==ilm) then
371 10454577 : ff(1:meshsz)=pawtab%phi(1:meshsz,iln)*pawtab%phi(1:meshsz,jln)*vhnzc(1:meshsz)
372 6009 : call simp_gen(intg,ff,radmesh)
373 6009 : pawtab%dij0(klmn)=pawtab%dij0(klmn)+intg
374 : end if
375 : end do
376 : end do
377 324 : LIBPAW_DEALLOCATE(vhnzc)
378 :
379 : !Computation of -<tphi_i|vh(tnZc)|tphi_j> on the PAW sphere
380 : !==========================================================
381 4091 : do jlmn=1,pawtab%lmn_size
382 3767 : j0lmn=jlmn*(jlmn-1)/2
383 3767 : jlm=indlmn(4,jlmn);jln=indlmn(5,jlmn)
384 40999 : do ilmn=1,jlmn
385 36908 : klmn=j0lmn+ilmn
386 36908 : ilm=indlmn(4,ilmn);iln=indlmn(5,ilmn)
387 40675 : if (jlm==ilm) then
388 10454577 : ff(1:meshsz)=pawtab%tphi(1:meshsz,iln)*pawtab%tphi(1:meshsz,jln)*vhtnzc_sph(1:meshsz)
389 6009 : call simp_gen(intg,ff,radmesh)
390 6009 : 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 324 : 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 324 : if (opt_init==0) then
419 972 : LIBPAW_ALLOCATE(shpf,(radmesh%mesh_size))
420 324 : call atompaw_shpfun(0,radmesh,intg,pawtab,shpf)
421 324 : 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 548207 : ff(1:meshsz)=vhtnzc_sph(1:meshsz)*shpf(1:meshsz)*radmesh%rad(1:meshsz)**2
434 324 : LIBPAW_DEALLOCATE(shpf)
435 324 : call simp_gen(intvh,ff,radmesh)
436 4091 : do jlmn=1,pawtab%lmn_size
437 3767 : j0lmn=jlmn*(jlmn-1)/2
438 3767 : jl=indlmn(1,jlmn);jln=indlmn(5,jlmn);jlm=indlmn(4,jlmn)
439 40999 : do ilmn=1,jlmn
440 36908 : klmn=j0lmn+ilmn
441 36908 : il=indlmn(1,ilmn);iln=indlmn(5,ilmn);ilm=indlmn(4,ilmn)
442 40675 : if (ilm==jlm) then
443 : ff(1:meshsz)=(pawtab%phi (1:meshsz,iln)*pawtab%phi (1:meshsz,jln)&
444 10454577 : & -pawtab%tphi(1:meshsz,iln)*pawtab%tphi(1:meshsz,jln))
445 6009 : call simp_gen(intg,ff,radmesh)
446 6009 : 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 324 : LIBPAW_DEALLOCATE(ff)
468 324 : LIBPAW_DEALLOCATE(vhtnzc_sph)
469 :
470 324 : 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 559 : 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 559 : mesh_size=size(ncore)
707 559 : if (mesh_size/=radmesh_core%mesh_size) then
708 0 : LIBPAW_BUG('wrong sizes!')
709 : end if
710 1677 : LIBPAW_ALLOCATE(vhnzc,(mesh_size))
711 1118 : LIBPAW_ALLOCATE(ff,(mesh_size))
712 559 : ehnzc=zero
713 720017 : vhnzc=zero
714 720017 : ff(:)=ncore(:)*four_pi*radmesh_core%rad(:)**2
715 559 : call poisson(ff,0,radmesh_core,vhnzc)
716 719458 : vhnzc(2:mesh_size)=(vhnzc(2:mesh_size)*half-znucl)
717 559 : call pawrad_deducer0(vhnzc,mesh_size,radmesh_core)
718 720017 : ff(:)=vhnzc(:)*ncore(:)*(four_pi*radmesh_core%rad(:))
719 559 : LIBPAW_DEALLOCATE(vhnzc)
720 559 : call simp_gen(ehnzc,ff,radmesh_core)
721 559 : LIBPAW_DEALLOCATE(ff)
722 :
723 559 : end subroutine atompaw_ehnzc
724 : !!***
725 :
726 : end module m_paw_atom
727 : !!***
|
