Line data Source code
1 : !!****m* ABINIT/m_paw_onsite
2 : !! NAME
3 : !! m_paw_onsite
4 : !!
5 : !! FUNCTION
6 : !! This module contains a set of routines to compute various PAW on-site quantities
7 : !! i.e. quantities expressed with <Phi_i|...|Phi_j> and/or <tild_Phi_i|...|tild_Phi_j>.
8 : !!
9 : !! COPYRIGHT
10 : !! Copyright (C) 2013-2026 ABINIT group (MT,FJ)
11 : !! This file is distributed under the terms of the
12 : !! GNU General Public License, see ~abinit/COPYING
13 : !! or http://www.gnu.org/copyleft/gpl.txt .
14 : !!
15 : !! NOTES
16 : !! FOR DEVELOPPERS: in order to preserve the portability of libPAW library,
17 : !! please consult ~abinit/src/??_libpaw/libpaw-coding-rules.txt
18 : !!
19 : !! SOURCE
20 :
21 : #include "libpaw.h"
22 :
23 : MODULE m_paw_onsite
24 :
25 : USE_DEFS
26 : USE_MSG_HANDLING
27 : USE_MEMORY_PROFILING
28 :
29 : use m_paw_atomorb, only : atomorb_type
30 : use m_pawrad, only : pawrad_type, pawrad_deducer0, simp_gen, nderiv_gen
31 : use m_pawtab, only : pawtab_type
32 : use m_paw_sphharm, only : setnabla_ylm
33 :
34 :
35 : implicit none
36 :
37 : private
38 :
39 : !public procedures.
40 : public :: pawnabla_init ! Evaluate valence-valence on-site contribs of the nabla operator in cart. coord.
41 : public :: pawnabla_core_init ! Evaluate core-valence on-site contribs of the nabla operator in cart. coord.
42 :
43 : !!***
44 :
45 : CONTAINS
46 :
47 : !===========================================================
48 : !!***
49 :
50 : !----------------------------------------------------------------------
51 :
52 : !!****f* m_paw_onsite/pawnabla_init
53 : !! NAME
54 : !! pawnabla_init
55 : !!
56 : !! FUNCTION
57 : !! Evaluate all valence-valence onsite contributions of the nabla operator in cartesian coordinates,
58 : !! i.e. <Phi_i|Nabla|Phi_j>-<tPhi_i|Nabla|tPhi_j>.
59 : !!
60 : !! INPUTS
61 : !! mpsang=1+maximum angular momentum
62 : !! ntypat=Number of types of atoms in cell
63 : !! Pawrad(ntypat)<Pawrad_type>=PAW radial mesh and related data:
64 : !! %rad(mesh_size)=The coordinates of all the points of the radial mesh
65 : !! Pawtab(ntypat) <type(pawtab_type>=PAW tabulated starting data:
66 : !! %mesh_size=Dimension of radial mesh
67 : !! %lmn_size=Number of (l,m,n) elements for the PAW basis
68 : !!
69 : !! OUTPUT
70 : !! See side effects
71 : !!
72 : !! SIDE EFFECTS
73 : !! Pawtab(ntypat) <type(pawtab_type>=PAW tabulated starting data:
74 : !! %has_nabla=set to 1 in matrix elements are calculated and stored
75 : !! %nabla_ij(3,lmn_size,lmn_size)= <phi_i|nabla|phi_j>-<tphi_i|nabla|tphi_j>
76 : !!
77 : !! NOTES
78 : !! MG extracted this piece of code from optics_paw.F90 in order to have something more
79 : !! reusable! Note however the storage mode of nabla_ij differs from optics_paw
80 : !! (here Cartesian coordinates run faster). Besides nabla_ij contains the matrix
81 : !! elements of \nabla instead of the elements of the momentum operator p.
82 : !!
83 : !! SOURCE
84 :
85 1154 : subroutine pawnabla_init(mpsang,ntypat,pawrad,pawtab)
86 :
87 : !Arguments ------------------------------------
88 : !scalars
89 : integer,intent(in) :: mpsang,ntypat
90 : !arrays
91 : type(pawtab_type),target,intent(inout) :: pawtab(ntypat)
92 : type(pawrad_type),intent(in) :: pawrad(ntypat)
93 :
94 : !Local variables-------------------------------
95 : !scalars
96 : integer :: ii,nln,il,ilm,ilmn,iln,itypat
97 : integer :: jl,jlm,jlmn,jln,lmn_size,mesh_size
98 : real(dp) :: avg,intg
99 : character(len=500) :: msg
100 : !arrays
101 1154 : integer, LIBPAW_CONTIGUOUS pointer :: indlmn(:,:)
102 2308 : real(dp) :: ang_phipphj(mpsang**2,mpsang**2,8)
103 1154 : real(dp),allocatable :: dphi(:),dtphi(:),ff(:),int1(:,:),int2(:,:),rad(:)
104 :
105 : ! *************************************************************************
106 :
107 1154 : if (mpsang>4)then
108 : write(msg,'(3a)')&
109 0 : & 'Not designed for angular momentum greater than 3 ',ch10,&
110 0 : & 'Modification in the table defined in routine setnabla_ylm is required.'
111 0 : LIBPAW_BUG(msg)
112 : end if
113 :
114 : !Integration of the angular part: all angular integrals have been computed
115 : !outside Abinit and tabulated for each (l,m) value up to l=3
116 1154 : call setnabla_ylm(ang_phipphj,mpsang)
117 :
118 2313 : do itypat=1,ntypat
119 :
120 : ! COMPUTE nabla_ij := <phi_i|nabla|phi_j>-<tphi_i|nabla|tphi_j> for this type
121 1159 : mesh_size=pawtab(itypat)%mesh_size
122 1159 : lmn_size=pawtab(itypat)%lmn_size
123 1159 : nln=pawtab(itypat)%basis_size
124 :
125 1159 : if (allocated(pawtab(itypat)%nabla_ij)) then
126 6 : LIBPAW_DEALLOCATE(pawtab(itypat)%nabla_ij)
127 : end if
128 4636 : LIBPAW_ALLOCATE(pawtab(itypat)%nabla_ij,(3,lmn_size,lmn_size))
129 1159 : pawtab(itypat)%has_nabla=1
130 :
131 3477 : LIBPAW_ALLOCATE(ff,(mesh_size))
132 2318 : LIBPAW_ALLOCATE(rad,(mesh_size))
133 4636 : LIBPAW_ALLOCATE(int1,(lmn_size,lmn_size))
134 3477 : LIBPAW_ALLOCATE(int2,(lmn_size,lmn_size))
135 2318 : LIBPAW_ALLOCATE(dphi,(mesh_size))
136 2318 : LIBPAW_ALLOCATE(dtphi,(mesh_size))
137 :
138 1159 : indlmn => pawtab(itypat)%indlmn
139 1485839 : rad(1:mesh_size)=pawrad(itypat)%rad(1:mesh_size)
140 :
141 : ! int1= \int [ Phi d/dr(Phj) - tPhi d/dr(tPhj) ] r^2 dr
142 : ! = \int [ (phi d/dr(phj) - phi phj/r) - (tphi d/dr(tphj) - tphi tphj/r) ] dr
143 : ! with Phi=phi/r and tPhi=phi/r
144 5796 : do jln=1,nln
145 5946683 : ff(1:mesh_size)=pawtab(itypat)%phi(1:mesh_size,jln)
146 4637 : call nderiv_gen(dphi,ff,pawrad(itypat))
147 5946683 : ff(1:mesh_size)=pawtab(itypat)%tphi(1:mesh_size,jln)
148 4637 : call nderiv_gen(dtphi,ff,pawrad(itypat))
149 24353 : do iln=1,nln
150 : ff(2:mesh_size)= &
151 : & pawtab(itypat)%phi (2:mesh_size,iln)*dphi (2:mesh_size) &
152 : & -pawtab(itypat)%phi (2:mesh_size,iln)*pawtab(itypat)%phi (2:mesh_size,jln)/rad(2:mesh_size) &
153 : & -( pawtab(itypat)%tphi(2:mesh_size,iln)*dtphi(2:mesh_size) &
154 23786550 : & -pawtab(itypat)%tphi(2:mesh_size,iln)*pawtab(itypat)%tphi(2:mesh_size,jln)/rad(2:mesh_size) )
155 18557 : call pawrad_deducer0(ff,mesh_size,pawrad(itypat))
156 18557 : call simp_gen(intg,ff,pawrad(itypat))
157 23194 : int1(iln,jln)=intg
158 : end do
159 : end do
160 :
161 : ! int2= \int [ Phi Phj /r - \int tPhi tPhj /r ] r^2 dr
162 : ! = \int [ phi phj /r - \int tphi tphj /r ] dr
163 : ! with Phi=phi/r and tPhi=phi/r
164 5796 : do jln=1,nln
165 24353 : do iln=1,nln
166 : ff(2:mesh_size)= ( &
167 : & pawtab(itypat)%phi (2:mesh_size,iln)*pawtab(itypat)%phi (2:mesh_size,jln) &
168 23786550 : & -pawtab(itypat)%tphi(2:mesh_size,iln)*pawtab(itypat)%tphi(2:mesh_size,jln) ) /rad(2:mesh_size)
169 18557 : call pawrad_deducer0(ff,mesh_size,pawrad(itypat))
170 18557 : call simp_gen(intg,ff,pawrad(itypat))
171 23194 : int2(iln,jln)=intg
172 : end do
173 : end do
174 :
175 : ! Integration of the radial part, Note unpacked loop
176 10440 : do jlmn=1,lmn_size
177 9281 : jlm=indlmn(4,jlmn)
178 9281 : jl =indlmn(5,jlmn)
179 84853 : do ilmn=1,lmn_size
180 74413 : ilm=indlmn(4,ilmn)
181 74413 : il =indlmn(5,ilmn)
182 :
183 : pawtab(itypat)%nabla_ij(1,ilmn,jlmn)= &
184 : & int1(il,jl)* ang_phipphj(ilm,jlm,1) &
185 74413 : & +int2(il,jl)*(ang_phipphj(ilm,jlm,2)+ang_phipphj(ilm,jlm,3))
186 :
187 : pawtab(itypat)%nabla_ij(2,ilmn,jlmn)= &
188 : & int1(il,jl)* ang_phipphj(ilm,jlm,4) &
189 74413 : & +int2(il,jl)*(ang_phipphj(ilm,jlm,5)+ang_phipphj(ilm,jlm,6))
190 :
191 : pawtab(itypat)%nabla_ij(3,ilmn,jlmn)= &
192 : & int1(il,jl)* ang_phipphj(ilm,jlm,7) &
193 83694 : & +int2(il,jl)* ang_phipphj(ilm,jlm,8)
194 :
195 : end do !ilmn
196 : end do !jlmn
197 :
198 : ! Symetrization
199 1159 : if (lmn_size>1) then
200 9281 : do jlmn=2,lmn_size
201 41847 : do ilmn=1,jlmn-1
202 138386 : do ii=1,3
203 97698 : avg=half*(pawtab(itypat)%nabla_ij(ii,ilmn,jlmn)-pawtab(itypat)%nabla_ij(ii,jlmn,ilmn))
204 97698 : pawtab(itypat)%nabla_ij(ii,ilmn,jlmn)= avg
205 130264 : pawtab(itypat)%nabla_ij(ii,jlmn,ilmn)=-avg
206 : end do
207 : end do
208 : end do
209 : end if
210 :
211 : ! End
212 1159 : pawtab(itypat)%has_nabla=2
213 1159 : LIBPAW_DEALLOCATE(ff)
214 1159 : LIBPAW_DEALLOCATE(rad)
215 1159 : LIBPAW_DEALLOCATE(int2)
216 1159 : LIBPAW_DEALLOCATE(int1)
217 1159 : LIBPAW_DEALLOCATE(dphi)
218 2313 : LIBPAW_DEALLOCATE(dtphi)
219 :
220 : end do !itypat
221 :
222 2308 : end subroutine pawnabla_init
223 : !!***
224 :
225 : !----------------------------------------------------------------------
226 :
227 : !!****f* m_paw_onsite/pawnabla_core_init
228 : !! NAME
229 : !! pawnabla_core_init
230 : !!
231 : !! FUNCTION
232 : !! Evaluate core-valence onsite contributions of the nabla operator in cartesian coordinates,
233 : !! i.e. <Phi_i|Nabla|Phi_core_j>-<tPhi_i|Nabla|tPhi_core_j>.
234 : !! Core wave-functions are only given for one atom type.
235 : !!
236 : !! INPUTS
237 : !! mpsang=1+maximum angular momentum
238 : !! ntypat=Number of types of atoms in cell
239 : !! Pawrad(ntypat)<Pawrad_type>=PAW radial mesh and related data:
240 : !! %rad(mesh_size)=The coordinates of all the points of the radial mesh
241 : !! Pawtab(ntypat) <type(pawtab_type>=PAW tabulated starting data:
242 : !! %mesh_size=Dimension of radial mesh
243 : !! %lmn_size=Number of (l,m,n) elements for the PAW basis
244 : !! atm <type(paw_atomorb_type>= core tabulated data
245 : !!
246 : !! OUTPUT
247 : !! See side effects
248 : !!
249 : !! SIDE EFFECTS
250 : !! Pawtab(ntypat) <type(pawtab_type>=PAW tabulated starting data:
251 : !! %has_nabla=set to 1 in matrix elements are calculated and stored
252 : !! %nabla_ij(3,lmn_size,lmn_size)= <phi_i|nabla|phi_core_j>-<tphi_i|nabla|tphi_core_j>
253 : !!
254 : !! NOTES
255 : !! MG extracted this piece of code from optics_paw.F90 in order to have something more
256 : !! reusable! Note however the storage mode of nabla_ij differs from optics_paw
257 : !! (here Cartesian coordinates run faster). Besides nabla_ij contains the matrix
258 : !! elements of \nabla instead of the elements of the momentum operator p.
259 : !!
260 : !! SOURCE
261 :
262 3 : subroutine pawnabla_core_init(mpsang,ntypat,pawrad,pawtab,atm)
263 :
264 : !Arguments ------------------------------------
265 : !scalars
266 : integer,intent(in) :: mpsang,ntypat
267 : !arrays
268 : type(atomorb_type),intent(in) :: atm(ntypat)
269 : type(pawtab_type),target,intent(inout) :: pawtab(ntypat)
270 : type(pawrad_type),intent(in) :: pawrad(ntypat)
271 :
272 : !Local variables-------------------------------
273 : !scalars
274 : integer :: nln,nln_cor,ilm,ilmn,iln,itypat,sgnkappa
275 : integer :: jl,jm,jlm,jlmn,jln,js,jlm_re,jlm_im,jm_re,jm_im
276 : integer :: lmn_size,lmncmax,lcmax,ltmax,mesh_size,mesh_size_cor
277 : real(dp) :: intg,jmj,cgc
278 : logical :: dirac
279 : character(len=500) :: msg
280 : !arrays
281 3 : integer, LIBPAW_CONTIGUOUS pointer :: indlmn(:,:)
282 3 : real(dp) , allocatable:: ang_phipphj(:,:,:)
283 3 : real(dp),allocatable :: dphi(:),ff(:),int1(:,:),int2(:,:),rad(:)
284 :
285 : ! *************************************************************************
286 :
287 6 : do itypat=1,ntypat
288 3 : if(atm(itypat)%nsppol==2) LIBPAW_ERROR('Work in progress')
289 3 : lcmax=atm(itypat)%l_max
290 6 : ltmax=max(lcmax,mpsang)
291 : enddo
292 :
293 15 : LIBPAW_ALLOCATE(ang_phipphj,(ltmax**2,ltmax**2,8))
294 :
295 3 : if (ltmax>4)then
296 : write(msg,'(3a)')&
297 0 : & 'Not designed for angular momentum greater than 3!',ch10,&
298 0 : & 'Modification in the table defined in routine setnabla_ylm is required.'
299 0 : LIBPAW_BUG(msg)
300 : end if
301 :
302 : !if (mesh_size_cor/=pawrad(1)%mesh_size) then
303 : ! write(msg,'(a)') 'Wrong mesh_size_cor value (1)!'
304 : ! LIBPAW_BUG(msg)
305 : !end if
306 : !if (any(mesh_size_cor/=pawtab(:)%mesh_size)) then
307 : ! write(msg,'(3a)') 'Wrong mesh_size_cor value (2)!',ch10,&
308 : !& 'Should have only one type of atom.'
309 : ! LIBPAW_ERROR(msg)
310 : !end if
311 :
312 : !Integration of the angular part: all angular integrals have been computed
313 : !outside Abinit and tabulated for each (l,m) value up to l=3
314 3 : call setnabla_ylm(ang_phipphj,ltmax)
315 :
316 6 : do itypat=1,ntypat
317 3 : dirac=atm(itypat)%dirac
318 3 : mesh_size_cor=atm(itypat)%mesh_size
319 3 : nln_cor=atm(itypat)%ln_size
320 3 : lmncmax=atm(itypat)%lmn_size
321 :
322 : ! COMPUTE nabla_ij := <phi_i|nabla|phi_cor> for this type
323 3 : mesh_size=min(pawtab(itypat)%partialwave_mesh_size,pawrad(itypat)%mesh_size)
324 3 : mesh_size=min(mesh_size_cor,mesh_size)
325 3 : lmn_size=pawtab(itypat)%lmn_size
326 3 : nln=pawtab(itypat)%basis_size
327 :
328 : if (mesh_size_cor<mesh_size) then
329 : msg='mesh_size and mesh_sier_cor not compatible!'
330 : LIBPAW_BUG(msg)
331 : endif
332 :
333 3 : if (allocated(pawtab(itypat)%nabla_ij)) then
334 0 : LIBPAW_DEALLOCATE(pawtab(itypat)%nabla_ij)
335 : end if
336 12 : LIBPAW_ALLOCATE(pawtab(itypat)%nabla_ij,(3,lmn_size,lmncmax))
337 :
338 3 : if (dirac) then
339 2 : if (allocated(pawtab(itypat)%nabla_im_ij)) then
340 0 : LIBPAW_DEALLOCATE(pawtab(itypat)%nabla_im_ij)
341 : end if
342 6 : LIBPAW_ALLOCATE(pawtab(itypat)%nabla_im_ij,(3,lmn_size,lmncmax))
343 : end if
344 :
345 3 : pawtab(itypat)%has_nabla=1
346 :
347 9 : LIBPAW_ALLOCATE(ff,(mesh_size))
348 6 : LIBPAW_ALLOCATE(rad,(mesh_size))
349 12 : LIBPAW_ALLOCATE(int1,(lmn_size,lmncmax))
350 9 : LIBPAW_ALLOCATE(int2,(lmn_size,lmncmax))
351 6 : LIBPAW_ALLOCATE(dphi,(mesh_size))
352 :
353 3 : indlmn => pawtab(itypat)%indlmn
354 4609 : rad(1:mesh_size)=pawrad(itypat)%rad(1:mesh_size)
355 :
356 : ! int1= \int Phi d/dr(Phi_core) r^2 dr
357 : ! = \int (phi d/dr(phi_core) - phi phj_core/r) dr
358 : ! with Phi=phi/r and Phi_core=phi_core/r
359 14 : do jln=1,nln_cor
360 17831 : ff(1:mesh_size)=atm(itypat)%phi(1:mesh_size,jln,1)
361 11 : call nderiv_gen(dphi,ff,pawrad(itypat))
362 66 : do iln=1,nln
363 : ff(2:mesh_size)=pawtab(itypat)%phi(2:mesh_size,iln)*dphi(2:mesh_size) &
364 87288 : & -pawtab(itypat)%phi(2:mesh_size,iln)*atm(itypat)%phi(2:mesh_size,jln,1)/rad(2:mesh_size)
365 52 : call pawrad_deducer0(ff,mesh_size,pawrad(itypat))
366 52 : call simp_gen(intg,ff,pawrad(itypat))
367 63 : int1(iln,jln)=intg
368 : end do
369 : end do
370 :
371 : ! int2= \int Phi Phi_core /r r^2 dr = \int phi phi_core /r dr
372 : ! with Phi=phi/r and Phi_core=phi_core/r
373 14 : do jln=1,nln_cor
374 66 : do iln=1,nln
375 87288 : ff(2:mesh_size)=(pawtab(itypat)%phi(2:mesh_size,iln)*atm(itypat)%phi(2:mesh_size,jln,1))/rad(2:mesh_size)
376 52 : call pawrad_deducer0(ff,mesh_size,pawrad(itypat))
377 52 : call simp_gen(intg,ff,pawrad(itypat))
378 63 : int2(iln,jln)=intg
379 : end do
380 : end do
381 :
382 : ! ===== FULLY-RELATIVISTIC =====
383 3 : if(dirac) then
384 :
385 : ! Integration of the radial part, Note unpacked loop
386 42 : do jlmn=1,lmncmax
387 40 : jl=atm(itypat)%indlmn(1,jlmn)
388 40 : jm=atm(itypat)%indlmn(2,jlmn)
389 :
390 40 : sgnkappa=atm(itypat)%indlmn(3,jlmn)
391 40 : jmj=half*atm(itypat)%indlmn(8,jlmn) ! 2mj is stored in indlmn_cor
392 40 : js=atm(itypat)%indlmn(6,jlmn) ! 1 is up, 2 is down
393 :
394 : ! Calculate spinor dependend coeffs
395 : ! (Clebsch-Gordan, I guess)
396 40 : cgc=one ! so nothing changes without core spinors
397 40 : if (sgnkappa==1) then
398 8 : if(js==1) then
399 4 : cgc= sqrt((dble(jl)-jmj+half)/dble(2*jl+1))
400 : else
401 4 : cgc=-sqrt((dble(jl)+jmj+half)/dble(2*jl+1))
402 : endif
403 : else
404 32 : if(js==1) then
405 10 : cgc= sqrt((dble(jl)+jmj+half)/dble(2*jl+1))
406 : else
407 22 : cgc= sqrt((dble(jl)-jmj+half)/dble(2*jl+1))
408 : endif
409 : endif
410 :
411 40 : jlm=atm(itypat)%indlmn(4,jlmn)
412 40 : jln=atm(itypat)%indlmn(5,jlmn)
413 :
414 562 : do ilmn=1,lmn_size
415 520 : ilm=indlmn(4,ilmn)
416 520 : iln=indlmn(5,ilmn)
417 :
418 : ! jl was set as a flag for invalid combinations
419 : ! i.e. m=-(l+1) or m=(l+1)
420 : ! In these cases, cgc=0 ; so nabla_ij=0
421 560 : if(jl==-1) then
422 624 : pawtab(itypat)%nabla_ij(1:3,ilmn,jlmn)= zero
423 624 : pawtab(itypat)%nabla_im_ij(1:3,ilmn,jlmn) = zero
424 :
425 : else
426 :
427 : ! if jm<>0, need to convert from complex
428 : ! to real spherical harmonics
429 364 : if(jm<0) then
430 78 : jm_re=abs(jm)
431 78 : jm_im=-abs(jm)
432 78 : jlm_re=jl*(jl+1)+jm_re+1
433 78 : jlm_im=jl*(jl+1)+jm_im+1
434 : pawtab(itypat)%nabla_ij(1,ilmn,jlmn)=half_sqrt2*cgc*( &
435 : & int1(iln,jln)* ang_phipphj(ilm,jlm_re,1) &
436 78 : & +int2(iln,jln)*(ang_phipphj(ilm,jlm_re,2)+ang_phipphj(ilm,jlm_re,3)))
437 : pawtab(itypat)%nabla_ij(2,ilmn,jlmn)=half_sqrt2*cgc*( &
438 : & int1(iln,jln)* ang_phipphj(ilm,jlm_re,4) &
439 78 : & +int2(iln,jln)*(ang_phipphj(ilm,jlm_re,5)+ang_phipphj(ilm,jlm_re,6)))
440 : pawtab(itypat)%nabla_ij(3,ilmn,jlmn)=half_sqrt2*cgc*( &
441 : & int1(iln,jln)* ang_phipphj(ilm,jlm_re,7) &
442 78 : & +int2(iln,jln)* ang_phipphj(ilm,jlm_re,8))
443 : pawtab(itypat)%nabla_im_ij(1,ilmn,jlmn)=-half_sqrt2*cgc*( &
444 : & int1(iln,jln)* ang_phipphj(ilm,jlm_im,1) &
445 78 : & +int2(iln,jln)*(ang_phipphj(ilm,jlm_im,2)+ang_phipphj(ilm,jlm_im,3)))
446 : pawtab(itypat)%nabla_im_ij(2,ilmn,jlmn)=-half_sqrt2*cgc*( &
447 : & int1(iln,jln)* ang_phipphj(ilm,jlm_im,4) &
448 78 : & +int2(iln,jln)*(ang_phipphj(ilm,jlm_im,5)+ang_phipphj(ilm,jlm_im,6)))
449 : pawtab(itypat)%nabla_im_ij(3,ilmn,jlmn)=-half_sqrt2*cgc*( &
450 : & int1(iln,jln)* ang_phipphj(ilm,jlm_im,7) &
451 78 : & +int2(iln,jln)* ang_phipphj(ilm,jlm_im,8))
452 :
453 286 : else if (jm>0) then
454 78 : jm_re=abs(jm)
455 78 : jm_im=-abs(jm)
456 78 : jlm_re=jl*(jl+1)+jm_re+1
457 78 : jlm_im=jl*(jl+1)+jm_im+1
458 : pawtab(itypat)%nabla_ij(1,ilmn,jlmn)=((-1)**jm)*half_sqrt2*cgc*( &
459 : & int1(iln,jln)* ang_phipphj(ilm,jlm_re,1) &
460 78 : & +int2(iln,jln)*(ang_phipphj(ilm,jlm_re,2)+ang_phipphj(ilm,jlm_re,3)))
461 : pawtab(itypat)%nabla_ij(2,ilmn,jlmn)=((-1)**jm)*half_sqrt2*cgc*( &
462 : & int1(iln,jln)* ang_phipphj(ilm,jlm_re,4) &
463 78 : & +int2(iln,jln)*(ang_phipphj(ilm,jlm_re,5)+ang_phipphj(ilm,jlm_re,6)))
464 : pawtab(itypat)%nabla_ij(3,ilmn,jlmn)=((-1)**jm)*half_sqrt2*cgc*( &
465 : & int1(iln,jln)* ang_phipphj(ilm,jlm_re,7) &
466 78 : & +int2(iln,jln)* ang_phipphj(ilm,jlm_re,8))
467 : pawtab(itypat)%nabla_im_ij(1,ilmn,jlmn)=((-1)**jm)*half_sqrt2*cgc*( &
468 : & int1(iln,jln)* ang_phipphj(ilm,jlm_im,1) &
469 78 : & +int2(iln,jln)*(ang_phipphj(ilm,jlm_im,2)+ang_phipphj(ilm,jlm_im,3)))
470 : pawtab(itypat)%nabla_im_ij(2,ilmn,jlmn)=((-1)**jm)*half_sqrt2*cgc*( &
471 : & int1(iln,jln)* ang_phipphj(ilm,jlm_im,4) &
472 78 : & +int2(iln,jln)*(ang_phipphj(ilm,jlm_im,5)+ang_phipphj(ilm,jlm_im,6)))
473 : pawtab(itypat)%nabla_im_ij(3,ilmn,jlmn)=((-1)**jm)*half_sqrt2*cgc*( &
474 : & int1(iln,jln)* ang_phipphj(ilm,jlm_im,7) &
475 78 : & +int2(iln,jln)* ang_phipphj(ilm,jlm_im,8))
476 :
477 : else ! jm=0 : no conversion necessary if m=0
478 : pawtab(itypat)%nabla_ij(1,ilmn,jlmn)=cgc*( &
479 : & int1(iln,jln)* ang_phipphj(ilm,jlm,1) &
480 208 : & +int2(iln,jln)*(ang_phipphj(ilm,jlm,2)+ang_phipphj(ilm,jlm,3)))
481 : pawtab(itypat)%nabla_ij(2,ilmn,jlmn)=cgc*( &
482 : & int1(iln,jln)* ang_phipphj(ilm,jlm,4) &
483 208 : & +int2(iln,jln)*(ang_phipphj(ilm,jlm,5)+ang_phipphj(ilm,jlm,6)))
484 : pawtab(itypat)%nabla_ij(3,ilmn,jlmn)=cgc*( &
485 : & int1(iln,jln)* ang_phipphj(ilm,jlm,7) &
486 208 : & +int2(iln,jln)* ang_phipphj(ilm,jlm,8))
487 208 : pawtab(itypat)%nabla_im_ij(1,ilmn,jlmn)= zero
488 208 : pawtab(itypat)%nabla_im_ij(2,ilmn,jlmn)= zero
489 208 : pawtab(itypat)%nabla_im_ij(3,ilmn,jlmn)= zero
490 : end if
491 :
492 : endif ! jl==-1?
493 :
494 : end do !ilmn
495 : end do !jlmn
496 :
497 2 : pawtab(itypat)%has_nabla=4
498 :
499 : ! ===== NON-RELATIVISTIC OR SCALAR-RELATICISTIC =====
500 : else
501 :
502 : ! Integration of the radial part, Note unpacked loop
503 6 : do jlmn=1,lmncmax
504 5 : jl=atm(itypat)%indlmn(1,jlmn)
505 5 : jlm=atm(itypat)%indlmn(4,jlmn)
506 5 : jln =atm(itypat)%indlmn(5,jlmn)
507 46 : do ilmn=1,lmn_size
508 40 : ilm=indlmn(4,ilmn)
509 40 : iln =indlmn(5,ilmn)
510 : pawtab(itypat)%nabla_ij(1,ilmn,jlmn)=( &
511 : & int1(iln,jln)* ang_phipphj(ilm,jlm,1) &
512 40 : & +int2(iln,jln)*(ang_phipphj(ilm,jlm,2)+ang_phipphj(ilm,jlm,3)))
513 :
514 : pawtab(itypat)%nabla_ij(2,ilmn,jlmn)=( &
515 : & int1(iln,jln)* ang_phipphj(ilm,jlm,4) &
516 40 : & +int2(iln,jln)*(ang_phipphj(ilm,jlm,5)+ang_phipphj(ilm,jlm,6)))
517 :
518 : pawtab(itypat)%nabla_ij(3,ilmn,jlmn)=( &
519 : & int1(iln,jln)* ang_phipphj(ilm,jlm,7) &
520 45 : & +int2(iln,jln)* ang_phipphj(ilm,jlm,8))
521 : end do !ilmn
522 : end do !jlmn
523 :
524 1 : pawtab(itypat)%has_nabla=3
525 :
526 : end if ! Relativistic?
527 3 : LIBPAW_DEALLOCATE(ff)
528 3 : LIBPAW_DEALLOCATE(rad)
529 3 : LIBPAW_DEALLOCATE(int1)
530 3 : LIBPAW_DEALLOCATE(int2)
531 6 : LIBPAW_DEALLOCATE(dphi)
532 :
533 : end do !itypat
534 :
535 3 : LIBPAW_DEALLOCATE(ang_phipphj)
536 :
537 6 : end subroutine pawnabla_core_init
538 : !!***
539 :
540 : !----------------------------------------------------------------------
541 :
542 :
543 : end module m_paw_onsite
544 : !!***
545 :
|
