Branch data 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 : 1155 : 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 : 1155 : integer, LIBPAW_CONTIGUOUS pointer :: indlmn(:,:)
102 : 2310 : real(dp) :: ang_phipphj(mpsang**2,mpsang**2,8)
103 : 1155 : real(dp),allocatable :: dphi(:),dtphi(:),ff(:),int1(:,:),int2(:,:),rad(:)
104 : :
105 : : ! *************************************************************************
106 : :
107 : 1155 : 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 : 1155 : call setnabla_ylm(ang_phipphj,mpsang)
117 : :
118 [ + + ]: 2316 : do itypat=1,ntypat
119 : :
120 : : ! COMPUTE nabla_ij := <phi_i|nabla|phi_j>-<tphi_i|nabla|tphi_j> for this type
121 : 1161 : mesh_size=pawtab(itypat)%mesh_size
122 : 1161 : lmn_size=pawtab(itypat)%lmn_size
123 : 1161 : nln=pawtab(itypat)%basis_size
124 : :
125 [ + + ]: 1161 : if (allocated(pawtab(itypat)%nabla_ij)) then
126 : 6 : LIBPAW_DEALLOCATE(pawtab(itypat)%nabla_ij)
127 : : end if
128 [ - + + - : 4644 : LIBPAW_ALLOCATE(pawtab(itypat)%nabla_ij,(3,lmn_size,lmn_size))
- + + - +
- + - - +
- + - + ]
129 : 1161 : pawtab(itypat)%has_nabla=1
130 : :
131 [ - + + - : 3483 : LIBPAW_ALLOCATE(ff,(mesh_size))
+ - + - -
+ - + -
+ ]
132 [ - + + - : 2322 : LIBPAW_ALLOCATE(rad,(mesh_size))
- + - + ]
133 [ - + + - : 4644 : LIBPAW_ALLOCATE(int1,(lmn_size,lmn_size))
- + + - +
- + - - +
- + ]
134 [ - + + - : 3483 : LIBPAW_ALLOCATE(int2,(lmn_size,lmn_size))
- + + - -
+ - + ]
135 [ - + + - : 2322 : LIBPAW_ALLOCATE(dphi,(mesh_size))
- + - + -
+ ]
136 [ - + + - : 2322 : LIBPAW_ALLOCATE(dtphi,(mesh_size))
- + - + -
+ ]
137 : :
138 : 1161 : indlmn => pawtab(itypat)%indlmn
139 [ + + ]: 1488679 : 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 [ + + ]: 5808 : do jln=1,nln
145 [ + + ]: 5961099 : ff(1:mesh_size)=pawtab(itypat)%phi(1:mesh_size,jln)
146 : 4647 : call nderiv_gen(dphi,ff,pawrad(itypat))
147 [ + + ]: 5961099 : ff(1:mesh_size)=pawtab(itypat)%tphi(1:mesh_size,jln)
148 : 4647 : call nderiv_gen(dtphi,ff,pawrad(itypat))
149 [ + + ]: 24417 : 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 [ + + ]: 23862498 : & -pawtab(itypat)%tphi(2:mesh_size,iln)*pawtab(itypat)%tphi(2:mesh_size,jln)/rad(2:mesh_size) )
155 : 18609 : call pawrad_deducer0(ff,mesh_size,pawrad(itypat))
156 : 18609 : call simp_gen(intg,ff,pawrad(itypat))
157 : 23256 : 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 [ + + ]: 5808 : do jln=1,nln
165 [ + + ]: 24417 : 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 [ + + ]: 23862498 : & -pawtab(itypat)%tphi(2:mesh_size,iln)*pawtab(itypat)%tphi(2:mesh_size,jln) ) /rad(2:mesh_size)
169 : 18609 : call pawrad_deducer0(ff,mesh_size,pawrad(itypat))
170 : 18609 : call simp_gen(intg,ff,pawrad(itypat))
171 : 23256 : int2(iln,jln)=intg
172 : : end do
173 : : end do
174 : :
175 : : ! Integration of the radial part, Note unpacked loop
176 [ + + ]: 10468 : do jlmn=1,lmn_size
177 : 9307 : jlm=indlmn(4,jlmn)
178 : 9307 : jl =indlmn(5,jlmn)
179 [ + + ]: 85269 : do ilmn=1,lmn_size
180 : 74801 : ilm=indlmn(4,ilmn)
181 : 74801 : il =indlmn(5,ilmn)
182 : :
183 : : pawtab(itypat)%nabla_ij(1,ilmn,jlmn)= &
184 : : & int1(il,jl)* ang_phipphj(ilm,jlm,1) &
185 : 74801 : & +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 : 74801 : & +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 : 84108 : & +int2(il,jl)* ang_phipphj(ilm,jlm,8)
194 : :
195 : : end do !ilmn
196 : : end do !jlmn
197 : :
198 : : ! Symetrization
199 [ + - ]: 1161 : if (lmn_size>1) then
200 [ + + ]: 9307 : do jlmn=2,lmn_size
201 [ + + ]: 42054 : do ilmn=1,jlmn-1
202 [ + + ]: 139134 : do ii=1,3
203 : 98241 : avg=half*(pawtab(itypat)%nabla_ij(ii,ilmn,jlmn)-pawtab(itypat)%nabla_ij(ii,jlmn,ilmn))
204 : 98241 : pawtab(itypat)%nabla_ij(ii,ilmn,jlmn)= avg
205 : 130988 : pawtab(itypat)%nabla_ij(ii,jlmn,ilmn)=-avg
206 : : end do
207 : : end do
208 : : end do
209 : : end if
210 : :
211 : : ! End
212 : 1161 : pawtab(itypat)%has_nabla=2
213 [ - + ]: 1161 : LIBPAW_DEALLOCATE(ff)
214 [ - + ]: 1161 : LIBPAW_DEALLOCATE(rad)
215 [ - + ]: 1161 : LIBPAW_DEALLOCATE(int2)
216 [ - + ]: 1161 : LIBPAW_DEALLOCATE(int1)
217 [ - + ]: 1161 : LIBPAW_DEALLOCATE(dphi)
218 [ - + ]: 2316 : LIBPAW_DEALLOCATE(dtphi)
219 : :
220 : : end do !itypat
221 : :
222 [ - + - + : 2310 : 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 : :
|
