Line data Source code
1 : !!****m* ABINIT/m_paw_finegrid
2 : !! NAME
3 : !! m_paw_finegrid
4 : !!
5 : !! FUNCTION
6 : !! This module contains a set of routines to compute various quantities
7 : !! on the fine grid around a given atom.
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_finegrid
24 :
25 : USE_DEFS
26 : USE_MSG_HANDLING
27 : USE_MEMORY_PROFILING
28 :
29 : use m_pawtab, only : pawtab_type
30 : use m_paw_sphharm, only : initylmr
31 : use m_paw_numeric, only : paw_jbessel,paw_splint,paw_uniform_splfit,paw_sort_dp
32 :
33 : implicit none
34 :
35 : private
36 :
37 : !public procedures.
38 : public :: pawgylm ! g_l(r-R)*Y_lm(r-R) (and derivatives)
39 : public :: pawgylmg ! Fourier transform of g_l(r-R)*Y_lm(r-R), plane-waves case
40 : public :: pawrfgd_fft ! r-R, plane-waves case
41 : public :: pawrfgd_wvl ! r-R, wavelets case
42 : public :: pawexpiqr ! exp(i.q.(r-R))
43 :
44 : !declarations for the whole module (were needed to replace the statement functions)
45 : !MG: Why this? Global variables are powerful but extremely DANGEROUS
46 : integer,private,save :: lambda
47 : real(dp),private,save :: pi_over_rshp,sigma
48 : real(dp),private,save,allocatable :: alpha(:,:),qq(:,:)
49 : !!***
50 :
51 : CONTAINS
52 :
53 : !===========================================================
54 : !!***
55 :
56 : !----------------------------------------------------------------------
57 :
58 : !!****f* m_paw_finegrid/pawgylm
59 : !! NAME
60 : !! pawgylm
61 : !!
62 : !! FUNCTION
63 : !! Compute g_l(r-R)*Y_lm(r-R) (and derivatives) on the fine (rectangular) grid
64 : !! around one atom (g_l=radial shape function).
65 : !! R is the position of the atom
66 : !!
67 : !! INPUTS
68 : !! lm_size=number of lm components to be calculated
69 : !! nfgd= number of (fine grid) FFT points in the paw sphere around current atom
70 : !! optgr0= 1 if g_l(r-R)*Y_lm(r-R) are computed
71 : !! optgr1= 1 if first derivatives of g_l(r-R)*Y_lm(r-R) are computed
72 : !! optgr2= 1 if second derivatives of g_l(r-R)*Y_lm(r-R) are computed
73 : !! pawtab <type(pawtab_type)>=paw tabulated starting data for current atom
74 : !! rfgd(3,nfgd)= coordinates of r-R on the fine rect. grid around current atom
75 : !!
76 : !! OUTPUT
77 : !! if (optgr0==1)
78 : !! gylm(nfgd,lm_size)= g_l(r-R)*Y_lm(r-R) around current atom
79 : !! if (optgr1==1)
80 : !! gylmgr(3,nfgd,lm_size)= derivatives of g_l(r-R)*Y_lm(r-R) wrt cart. coordinates
81 : !! if (optgr2==1)
82 : !! gylmgr2(6,nfgd,lm_size)= second derivatives of g_l(r-R)*Y_lm(r-R) wrt cart. coordinates
83 : !!
84 : !! SOURCE
85 :
86 43589 : subroutine pawgylm(gylm,gylmgr,gylmgr2,lm_size,nfgd,optgr0,optgr1,optgr2,pawtab,rfgd)
87 :
88 : !Arguments ---------------------------------------------
89 : !scalars
90 : integer,intent(in) :: lm_size,nfgd,optgr0,optgr1,optgr2
91 : type(pawtab_type),intent(in) :: pawtab
92 : !arrays
93 : real(dp),intent(in) :: rfgd(:,:)
94 : real(dp),intent(out) :: gylm(nfgd,optgr0*lm_size)
95 : real(dp),intent(out) :: gylmgr(3,nfgd,optgr1*lm_size)
96 : real(dp),intent(out) :: gylmgr2(6,nfgd,optgr2*lm_size)
97 :
98 : !Local variables ------------------------------
99 : !scalars
100 : integer :: ic,ilm,izero,l_size,ll,normchoice,option,shape_type
101 : real(dp) :: arg
102 : real(dp) :: jbes1,jbes2,jbesp1,jbesp2,jbespp1,jbespp2,rcut
103 : real(dp) :: splfact
104 : logical :: compute_gr0,compute_gr1,compute_gr2
105 : character(len=500) :: msg
106 : !arrays
107 43589 : integer,allocatable :: isort(:)
108 : real(dp),parameter :: ffact(1:9)=(/1._dp,3._dp,15._dp,105._dp,945._dp,10395._dp,&
109 : & 135135._dp,2027025._dp,34459425._dp/)
110 : real(dp),parameter :: toldev=tol3
111 : real(dp) :: ss(3)
112 43589 : real(dp),allocatable :: cc(:,:),d2gfact(:,:),d2shpfuncnum(:,:),dgfact(:,:)
113 43589 : real(dp),allocatable :: dshpfuncnum(:,:),gfact(:,:)
114 43589 : real(dp),allocatable :: rnrm(:),rnrm_inv(:),rnrm_sort(:)
115 43589 : real(dp),allocatable :: shpfuncnum(:,:),work(:),ylmr(:,:),ylmrgr(:,:,:)
116 :
117 : ! *************************************************************************
118 :
119 43589 : if (optgr0==0.and.optgr1==0.and.optgr2==0) return
120 43589 : if (nfgd==0) return
121 :
122 : !Compatibility test
123 : !==========================================================
124 130719 : if (size(rfgd)/=3*nfgd) then
125 0 : msg='rfgd array must be allocated at rfgd(3,nfgd)!'
126 0 : LIBPAW_BUG(msg)
127 : end if
128 : !if (pawtab%lcut_size>9) then
129 : ! msg='l_size>10 forbidden!'
130 : ! LIBPAW_BUG(msg)
131 : !end if
132 43573 : if (pawtab%shape_type==1.and.pawtab%shape_lambda<2) then
133 0 : msg='Exponent lambda of gaussian shape function must be > 1!'
134 0 : LIBPAW_ERROR(msg)
135 : end if
136 :
137 : !Initializations
138 : !==========================================================
139 : !Options for computation
140 43573 : compute_gr0=(optgr0==1.or.optgr1==1.or.optgr2==1)
141 43573 : compute_gr1=(optgr1==1.or.optgr2==1)
142 43573 : compute_gr2=(optgr2==1)
143 43573 : l_size=pawtab%lcut_size
144 :
145 : !Norms of vectors around the atom
146 130719 : LIBPAW_ALLOCATE(rnrm,(nfgd))
147 40587004 : izero=-1
148 40587004 : do ic=1,nfgd
149 40543431 : rnrm(ic)=sqrt(rfgd(1,ic)**2+rfgd(2,ic)**2+rfgd(3,ic)**2)
150 40587004 : if (rnrm(ic)<=tol10) izero=ic ! Has to be consistent with initylmr !!
151 : end do
152 :
153 : !Initializations
154 388522548 : if (optgr0==1) gylm=zero
155 591658563 : if (optgr1==1) gylmgr=zero
156 14522104 : if (optgr2==1) gylmgr2=zero
157 :
158 : !Some definitions concerning shape function g_l(r)
159 43573 : shape_type=pawtab%shape_type
160 43573 : sigma=pawtab%shape_sigma;lambda=pawtab%shape_lambda
161 43573 : pi_over_rshp=pi/pawtab%rshp
162 43573 : rcut=tol12+pawtab%rshp
163 43573 : if (shape_type==3) then
164 107370 : LIBPAW_ALLOCATE(alpha,(2,l_size))
165 71580 : LIBPAW_ALLOCATE(qq,(2,l_size))
166 143300 : do ll=1,l_size
167 322530 : alpha(1:2,ll)=pawtab%shape_alpha(1:2,ll)
168 358320 : qq(1:2,ll)=pawtab%shape_q(1:2,ll)
169 : end do
170 : end if
171 :
172 : !If needed, sort selected radii by increasing norm
173 43573 : if (shape_type==-1) then
174 9 : LIBPAW_ALLOCATE(isort,(nfgd))
175 6 : LIBPAW_ALLOCATE(rnrm_sort,(nfgd))
176 4380 : do ic=1,nfgd
177 4380 : isort(ic)=ic
178 : end do
179 4380 : rnrm_sort(1:nfgd)=rnrm(1:nfgd)
180 3 : call paw_sort_dp(nfgd,rnrm_sort,isort,tol16)
181 : end if
182 :
183 : !If shape function is "numeric", spline it onto selected radii
184 : if (shape_type==-1) then
185 6 : LIBPAW_ALLOCATE(work,(nfgd))
186 3 : if (compute_gr0) then
187 12 : LIBPAW_ALLOCATE(shpfuncnum,(nfgd,l_size))
188 18 : do ll=1,l_size
189 : call paw_splint(pawtab%mesh_size,pawtab%rad_for_spline,pawtab%shapefunc(:,ll),&
190 15 : & pawtab%dshpfunc(:,ll,2),nfgd,rnrm_sort,work)
191 21903 : do ic=1,nfgd
192 21900 : shpfuncnum(isort(ic),ll)=work(ic)
193 : end do
194 : end do
195 : end if
196 3 : if(compute_gr1) then
197 12 : LIBPAW_ALLOCATE(dshpfuncnum,(nfgd,l_size))
198 18 : do ll=1,l_size
199 : call paw_splint(pawtab%mesh_size,pawtab%rad_for_spline,pawtab%dshpfunc(:,ll,1),&
200 15 : & pawtab%dshpfunc(:,ll,3),nfgd,rnrm_sort,work)
201 21903 : do ic=1,nfgd
202 21900 : dshpfuncnum(isort(ic),ll)=work(ic)
203 : end do
204 : end do
205 : end if
206 3 : if(compute_gr2) then
207 0 : LIBPAW_ALLOCATE(d2shpfuncnum,(nfgd,l_size))
208 0 : do ll=1,l_size
209 : call paw_splint(pawtab%mesh_size,pawtab%rad_for_spline,pawtab%dshpfunc(:,ll,2),&
210 0 : & pawtab%dshpfunc(:,ll,4),nfgd,rnrm_sort,work)
211 0 : do ic=1,nfgd
212 0 : d2shpfuncnum(isort(ic),ll)=work(ic)
213 : end do
214 : end do
215 : end if
216 3 : LIBPAW_DEALLOCATE(work)
217 : end if
218 :
219 43573 : if (shape_type==-1) then
220 3 : LIBPAW_DEALLOCATE(isort)
221 3 : LIBPAW_DEALLOCATE(rnrm_sort)
222 : end if
223 :
224 : !If needed, compute limits at r=0 of shape function and derivatives
225 43573 : if (izero>0) then
226 110838 : LIBPAW_ALLOCATE(cc,(3,l_size))
227 487538 : cc=zero
228 36946 : if (shape_type==-1) then
229 : splfact=(pawtab%rad_for_spline(4)-pawtab%rad_for_spline(1))&
230 3 : & /(pawtab%rad_for_spline(3)-pawtab%rad_for_spline(2))
231 : end if
232 149594 : do ll=1,l_size
233 : ! cc(2,l) is g_prime(0)
234 112648 : if (optgr0==1.or.optgr1==1.or.optgr2==1) then
235 112648 : if (shape_type==-1) then
236 60 : ss(1:3)=pawtab%shapefunc(2:4,ll)/pawtab%rad_for_spline(2:4)**(ll-1)
237 15 : cc(1,ll)=ss(3)+(ss(1)-ss(2))*splfact
238 112633 : else if (shape_type==1.or.shape_type==2) then
239 6444 : cc(1,ll)=one
240 106189 : else if (shape_type==3) then
241 : cc(1,ll)=(alpha(1,ll)*qq(1,ll)**(ll-1) &
242 106189 : & +alpha(2,ll)*qq(2,ll)**(ll-1))/ffact(ll)
243 : end if
244 112648 : cc(1,ll)=cc(1,ll)*pawtab%gnorm(ll)
245 : end if
246 : ! cc(2,l) is g_prime(0)
247 112648 : if (optgr1==1.or.optgr2==1) then
248 15802 : if (shape_type==-1) then
249 60 : ss(1:3)=(ss(1:3)-cc(1,ll))/pawtab%rad_for_spline(2:4)
250 15 : cc(2,ll)=ss(3)+(ss(1)-ss(2))*splfact
251 15787 : else if (shape_type==1.and.lambda==1) then
252 0 : cc(2,ll)=-one/sigma
253 : else
254 15787 : cc(2,ll)=zero
255 : end if
256 15802 : cc(2,ll)=cc(2,ll)*pawtab%gnorm(ll)
257 : end if
258 : ! cc(3,l) is g_prime_prime(0)
259 149594 : if (optgr2==1) then
260 203 : if (shape_type==-1) then
261 0 : ss(1:3)=(ss(1:3)-cc(2,ll))/pawtab%rad_for_spline(2:4)
262 0 : cc(3,ll)=two*(ss(3)+(ss(1)-ss(2))*splfact)
263 203 : else if (shape_type==1) then
264 0 : if (lambda==1) cc(3,ll)=one/sigma**2
265 0 : if (lambda==2) cc(3,ll)=-two/sigma**2
266 0 : if (lambda >2) cc(3,ll)=zero
267 203 : else if (shape_type==2) then
268 186 : cc(3,ll)=-(two/three)*pi_over_rshp**2
269 17 : else if (shape_type==3) then
270 : cc(3,ll)=-(alpha(1,ll)*qq(1,ll)**(ll+1) &
271 17 : & +alpha(2,ll)*qq(2,ll)**(ll+1))/ffact(ll+1)
272 : end if
273 203 : cc(3,ll)=cc(3,ll)*pawtab%gnorm(ll)
274 : end if
275 : end do
276 : end if
277 :
278 : !Y_lm(r-R) calculation
279 : !==========================================================
280 43573 : normchoice=1 ; option=max(optgr0,2*optgr1,3*optgr2)
281 43573 : if(compute_gr0) then
282 174292 : LIBPAW_ALLOCATE(ylmr,(l_size**2,nfgd))
283 : end if
284 43573 : if(compute_gr1.and.(.not.compute_gr2)) then
285 35224 : LIBPAW_ALLOCATE(ylmrgr,(3,l_size**2,nfgd))
286 : end if
287 43573 : if(compute_gr2) then
288 388 : LIBPAW_ALLOCATE(ylmrgr,(9,l_size**2,nfgd))
289 : end if
290 43573 : if (compute_gr0.and.(.not.compute_gr1).and.(.not.compute_gr2)) then
291 34670 : call initylmr(l_size,normchoice,nfgd,rnrm,option,rfgd,ylmr)
292 : else
293 8903 : call initylmr(l_size,normchoice,nfgd,rnrm,option,rfgd,ylmr,ylmrgr)
294 : end if
295 :
296 : !gl(r) and derivatives calculation for l>=0
297 : !==========================================================
298 : !Compute gl(r), gl_prime(r)/r and (gl_prime_prime(r)-gl_prime(r)/r)/r**2
299 43573 : if (compute_gr0) then
300 174292 : LIBPAW_BOUND2_ALLOCATE(gfact,BOUNDS(1,nfgd),BOUNDS(0,l_size-1))
301 132765529 : gfact(:,:)=zero
302 : end if
303 43573 : if (compute_gr1) then
304 35612 : LIBPAW_BOUND2_ALLOCATE(dgfact,BOUNDS(1,nfgd),BOUNDS(0,l_size-1))
305 37063783 : dgfact(:,:)=zero
306 : end if
307 43573 : if (compute_gr2) then
308 388 : LIBPAW_BOUND2_ALLOCATE(d2gfact,BOUNDS(1,nfgd),BOUNDS(0,l_size-1))
309 521746 : d2gfact(:,:)=zero
310 : end if
311 43573 : if(compute_gr1) then
312 17806 : LIBPAW_ALLOCATE(rnrm_inv,(nfgd))
313 10480822 : do ic=1,nfgd
314 10480822 : if (ic/=izero) rnrm_inv(ic)=one/rnrm(ic)
315 : end do
316 8903 : if (izero>0) rnrm_inv(izero)=zero
317 : end if
318 :
319 : !----- type -1 -----
320 43573 : if (shape_type==-1) then
321 3 : if (compute_gr0) then
322 18 : do ll=0,l_size-1
323 21903 : do ic=1,nfgd
324 21900 : if (rnrm(ic)<=rcut) then
325 21885 : gfact(ic,ll)=shpfuncnum(ic,ll+1)
326 : end if
327 : end do
328 : end do
329 : end if
330 3 : if (compute_gr1) then
331 18 : do ll=0,l_size-1
332 21903 : do ic=1,nfgd
333 21900 : if (rnrm(ic)<=rcut) then
334 21885 : dgfact(ic,ll)=dshpfuncnum(ic,ll+1)*rnrm_inv(ic)
335 : end if
336 : end do
337 : end do
338 : end if
339 3 : if(compute_gr2) then
340 0 : do ll=0,l_size-1
341 0 : do ic=1,nfgd
342 0 : if (rnrm(ic)<=rcut) then
343 0 : d2gfact(ic,ll)=(d2shpfuncnum(ic,ll+1)-dgfact(ic,ll))*rnrm_inv(ic)**2
344 : end if
345 : end do
346 : end do
347 : end if
348 :
349 : ! ----- type 1 or 2 -----
350 43570 : else if (shape_type==1.or.shape_type==2) then
351 : ! FIRST COMPUTE FACTORS FOR l=0
352 7780 : if (optgr0==1.and.optgr1==0.and.optgr2==0) then
353 3369 : if (shape_type==1) then
354 5320 : do ic=1,nfgd
355 5318 : arg=rnrm(ic)
356 5320 : if (arg<toldev) then
357 2 : gfact(ic,0)=shapefunc1_0(arg)
358 5316 : else if (arg<=rcut) then
359 5316 : gfact(ic,0)=shapefunc1(arg)
360 : end if
361 : end do
362 : else ! shape_type==2
363 5387301 : do ic=1,nfgd
364 5383934 : arg=rnrm(ic)
365 5387301 : if (arg<toldev) then
366 806 : gfact(ic,0)=shapefunc2_0(arg)
367 5383128 : else if (arg<=rcut) then
368 5383128 : gfact(ic,0)=shapefunc2(arg)
369 : end if
370 : end do
371 : end if
372 4411 : else if (optgr1==1.and.optgr2==0) then
373 4337 : if (shape_type==1) then
374 5320 : do ic=1,nfgd
375 5318 : arg=rnrm(ic)
376 5320 : if (arg<toldev) then
377 2 : gfact(ic,0)=shapefunc1_0(arg)
378 2 : if (lambda==2) then
379 2 : dgfact(ic,0)=dshpfunc1_ovr_0_2(arg)
380 : else ! lambda>2
381 0 : dgfact(ic,0)=dshpfunc1_ovr_0(arg)
382 : end if
383 5316 : else if (arg<=rcut) then
384 5316 : gfact(ic,0)=shapefunc1(arg)
385 5316 : dgfact(ic,0)=dshpfunc1(arg)*rnrm_inv(ic)
386 : end if
387 : end do
388 : else ! shape_type==2
389 6102617 : do ic=1,nfgd
390 6098282 : arg=rnrm(ic)
391 6102617 : if (arg<toldev) then
392 872 : gfact(ic,0)=shapefunc2_0(arg)
393 872 : dgfact(ic,0)=dshpfunc2_ovr_0(arg)
394 6097410 : else if (arg<=rcut) then
395 6097410 : gfact(ic,0)=shapefunc2(arg)
396 6097410 : dgfact(ic,0)=dshpfunc2(arg)*rnrm_inv(ic)
397 : end if
398 : end do
399 : end if
400 74 : else if (optgr2==1) then
401 74 : if (shape_type==1) then
402 0 : do ic=1,nfgd
403 0 : arg=rnrm(ic)
404 0 : if (arg<toldev) then
405 0 : gfact(ic,0)=shapefunc1_0(arg)
406 0 : if (lambda==2) then
407 0 : dgfact(ic,0)=dshpfunc1_ovr_0_2(arg)
408 0 : d2gfact(ic,0)=d2shpfunc1_ovr2_0_2(arg)
409 0 : else if (lambda==3) then
410 0 : dgfact(ic,0)=dshpfunc1_ovr_0(arg)
411 0 : if (ic/=izero) then
412 0 : d2gfact(ic,0)=d2shpfunc1_ovr2_0_3(arg)
413 : else
414 0 : d2gfact(ic,0)=zero ! Diverging case
415 : end if
416 0 : else if (lambda==4) then
417 0 : dgfact(ic,0)=dshpfunc1_ovr_0(arg)
418 0 : d2gfact(ic,0)=d2shpfunc1_ovr2_0_4(arg)
419 : else ! lambda>4
420 0 : dgfact(ic,0)=dshpfunc1_ovr_0(arg)
421 0 : d2gfact(ic,0)=d2shpfunc1_ovr2_0(arg)
422 : end if
423 0 : else if (arg<=rcut) then
424 0 : gfact(ic,0)=shapefunc1(arg)
425 0 : dgfact(ic,0)=dshpfunc1(arg)*rnrm_inv(ic)
426 0 : d2gfact(ic,0)=(d2shpfunc1(arg)-dgfact(ic,0))*rnrm_inv(ic)**2
427 : end if
428 : end do
429 : else ! shape_type==2
430 116196 : do ic=1,nfgd
431 116122 : arg=rnrm(ic)
432 116196 : if (arg<toldev) then
433 56 : gfact(ic,0)=shapefunc2_0(arg)
434 56 : dgfact(ic,0)=dshpfunc2_ovr_0(arg)
435 56 : d2gfact(ic,0)=d2shpfunc2_ovr2_0(arg)
436 116066 : else if (arg<=rcut) then
437 116066 : gfact(ic,0)=shapefunc2(arg)
438 116066 : dgfact(ic,0)=dshpfunc2(arg)*rnrm_inv(ic)
439 116066 : d2gfact(ic,0)=(d2shpfunc2(arg)-dgfact(ic,0))*rnrm_inv(ic)**2
440 : end if
441 : end do
442 : end if
443 : end if
444 :
445 : ! THEN COMPUTE FACTORS FOR l>0 (from l=0)
446 7780 : if (compute_gr0) then
447 7780 : if (l_size>1) then
448 19174 : do ll=1,l_size-1
449 31765196 : do ic=1,nfgd
450 31759504 : gfact(ic,ll)=pawtab%gnorm(ll+1)*gfact(ic,0)*rnrm(ic)**ll
451 : end do
452 : end do
453 : end if
454 : end if
455 7780 : if (compute_gr1) then
456 4411 : if (l_size>1) then
457 5819909 : do ic=1,nfgd
458 5819909 : dgfact(ic,1)=pawtab%gnorm(2)*(gfact(ic,0)*rnrm_inv(ic)+dgfact(ic,0)*rnrm(ic))
459 : end do
460 : end if
461 4411 : if (l_size>2) then
462 5819909 : do ic=1,nfgd
463 5819909 : dgfact(ic,2)=pawtab%gnorm(3)*(two*gfact(ic,0)+dgfact(ic,0)*rnrm(ic)**2)
464 : end do
465 : end if
466 4411 : if (l_size>3) then
467 1659 : do ll=3,l_size-1
468 5118507 : do ic=1,nfgd
469 : dgfact(ic,ll)=pawtab%gnorm(ll+1) &
470 5117974 : & *(dble(ll)*gfact(ic,0)*rnrm(ic)**(ll-2)+dgfact(ic,0)*rnrm(ic)**ll)
471 : end do
472 : end do
473 : end if
474 : end if
475 7780 : if (compute_gr2) then
476 74 : if (l_size>1) then
477 116196 : do ic=1,nfgd
478 : d2gfact(ic,1)=pawtab%gnorm(2) &
479 116196 : & *(-gfact(ic,0)*rnrm_inv(ic)**3+two*dgfact(ic,0)*rnrm_inv(ic)+d2gfact(ic,0)*rnrm(ic))
480 : end do
481 : end if
482 74 : if (l_size>2) then
483 116196 : do ic=1,nfgd
484 116196 : d2gfact(ic,2)=pawtab%gnorm(3)*(four*dgfact(ic,0)+d2gfact(ic,0)*rnrm(ic)**2)
485 : end do
486 : end if
487 74 : if (l_size>3) then
488 46734 : do ic=1,nfgd
489 : d2gfact(ic,3)=pawtab%gnorm(4) &
490 46734 : & *(three*gfact(ic,0)*rnrm_inv(ic)+6._dp*dgfact(ic,0)*rnrm(ic)+ d2gfact(ic,0)*rnrm(ic)**3)
491 : end do
492 : end if
493 74 : if (l_size>4) then
494 46734 : do ic=1,nfgd
495 : d2gfact(ic,4)=pawtab%gnorm(5) &
496 46734 : & *(8._dp*gfact(ic,0)+8._dp*dgfact(ic,0)*rnrm(ic)**2+d2gfact(ic,0)*rnrm(ic)**4)
497 : end do
498 : end if
499 74 : if (l_size>5) then
500 0 : do ll=5,l_size-1
501 0 : do ic=1,nfgd
502 : d2gfact(ic,ll)=pawtab%gnorm(ll+1) &
503 : & *(dble(ll*(ll-2))*gfact(ic,0)*rnrm(ic)**(ll-4) &
504 : & +dble(2*ll)*dgfact(ic,0)*rnrm(ic)**(ll-2) &
505 0 : & +d2gfact(ic,0)*rnrm(ic)**ll)
506 : end do
507 : end do
508 : end if
509 : end if
510 11616754 : if (compute_gr0) gfact(:,0)=gfact(:,0)*pawtab%gnorm(1)
511 6227502 : if (compute_gr1) dgfact(:,0)=dgfact(:,0)*pawtab%gnorm(1)
512 123902 : if (compute_gr2) d2gfact(:,0)=d2gfact(:,0)*pawtab%gnorm(1)
513 :
514 : ! ----- type 3 -----
515 35790 : else if (shape_type==3) then
516 35790 : if (optgr0==1.and.optgr1==0.and.optgr2==0) then
517 125459 : do ll=0,l_size-1
518 75304044 : do ic=1,nfgd
519 75178585 : arg=rnrm(ic)
520 75272743 : if (arg<=rcut) then
521 75178585 : call paw_jbessel(jbes1,jbesp1,jbespp1,ll,0,qq(1,1+ll)*arg)
522 75178585 : call paw_jbessel(jbes2,jbesp2,jbespp2,ll,0,qq(2,1+ll)*arg)
523 75178585 : gfact(ic,ll)=shapefunc3(jbes1,jbes2,ll)
524 : end if
525 : end do
526 : end do
527 4489 : else if (optgr1==1.and.optgr2==0) then
528 17747 : do ll=0,l_size-1
529 13975928 : do ic=1,nfgd
530 13958181 : arg=rnrm(ic)
531 13971462 : if (arg<=rcut) then
532 13958181 : call paw_jbessel(jbes1,jbesp1,jbespp1,ll,1,qq(1,1+ll)*arg)
533 13958181 : call paw_jbessel(jbes2,jbesp2,jbespp2,ll,1,qq(2,1+ll)*arg)
534 13958181 : gfact(ic,ll)=shapefunc3(jbes1,jbes2,ll)
535 13958181 : dgfact(ic,ll)=dshpfunc3(jbesp1,jbesp2,ll)*rnrm_inv(ic)
536 : end if
537 : end do
538 : end do
539 4466 : if (izero>0.and.l_size>=1) dgfact(izero,0)=-(alpha(1,1)*qq(1,1)+alpha(2,1)*qq(2,1))/three
540 4038 : if (izero>0.and.l_size>=3) dgfact(izero,2)=two/15._dp*(alpha(1,1)*qq(1,1)+alpha(2,1)*qq(2,1))
541 : ! Note: for l=1, dgfact is diverging - d2gfact is diverging for l<4
542 23 : else if (optgr2==1) then
543 94 : do ll=0,l_size-1
544 79616 : do ic=1,nfgd
545 79522 : arg=rnrm(ic)
546 79593 : if (arg<=rcut) then
547 79522 : call paw_jbessel(jbes1,jbesp1,jbespp1,ll,2,qq(1,1+ll)*arg)
548 79522 : call paw_jbessel(jbes2,jbesp2,jbespp2,ll,2,qq(2,1+ll)*arg)
549 79522 : gfact(ic,ll)=shapefunc3(jbes1,jbes2,ll)
550 79522 : dgfact(ic,ll)=dshpfunc3(jbesp1,jbesp2,ll)*rnrm_inv(ic)
551 79522 : d2gfact(ic,ll)=(d2shpfunc3(jbespp1,jbespp2,ll)-dgfact(ic,ll))*rnrm_inv(ic)**2
552 : end if
553 : end do
554 : end do
555 23 : if (izero>0.and.l_size>=1) dgfact(izero,0)=-(alpha(1,1)*qq(1,1)+alpha(2,1)*qq(2,1))/three
556 5 : if (izero>0.and.l_size>=3) dgfact(izero,2)=two/15._dp*(alpha(1,1)*qq(1,1)+alpha(2,1)*qq(2,1))
557 : ! Note: for l=1, dgfact is diverging - d2gfact is diverging for l<4
558 : end if
559 : end if
560 :
561 : !g_l(r-R)*Y_lm(r-R) calculation
562 : !==========================================================
563 43573 : if (optgr0==1) then
564 :
565 166395 : do ll=0,l_size-1
566 542909 : do ilm=ll**2+1,min((ll+1)**2,lm_size)
567 388603071 : do ic=1,nfgd
568 388478975 : gylm(ic,ilm)=gfact(ic,ll)*ylmr(ilm,ic)
569 : end do
570 : end do
571 : end do
572 :
573 : ! Special value at r-R=0 (supposing shapefunc(r)->C.r**l when r->0)
574 42299 : if (izero>0) then
575 371459 : gylm(izero,1:lm_size)=zero
576 36222 : if (lm_size>=1) gylm(izero,1)=ylmr(1,izero)*cc(1,1)
577 : end if
578 :
579 : end if
580 :
581 : !d/dr{g_l(r-R)*Y_lm(r-R)} calculation
582 : !==========================================================
583 43573 : if(optgr1==1) then
584 :
585 34029 : do ll=0,l_size-1
586 113915 : do ilm=ll**2+1,min((ll+1)**2,lm_size)
587 147988788 : do ic=1,nfgd
588 : gylmgr(1:3,ic,ilm)=gfact(ic,ll)*ylmrgr(1:3,ilm,ic)&
589 591614990 : & +dgfact(ic,ll)*rfgd(1:3,ic)*ylmr(ilm,ic)
590 : end do
591 : end do
592 : end do
593 :
594 : ! Special values at r-R=0 (supposing shapefunc(r)->C.r**l when r->0)
595 8903 : if (izero>0) then
596 215073 : gylmgr(1:3,izero,1:lm_size)=zero
597 4961 : if (lm_size>=1) then
598 4961 : arg=cc(2,1)/sqrt(four_pi)
599 19844 : gylmgr(1:3,izero,1)=arg
600 : end if
601 4961 : if (lm_size>=2) then
602 4923 : arg=cc(1,2)*sqrt(three/four_pi)
603 4923 : gylmgr(2,izero,2)=arg
604 4923 : if (lm_size>=3) gylmgr(3,izero,3)=arg
605 4923 : if (lm_size>=4) gylmgr(1,izero,4)=arg
606 : end if
607 : end if
608 :
609 : end if
610 :
611 : !d2/dridrj{g_l(r-R)*Y_lm(r-R)} calculation
612 : !==========================================================
613 43573 : if(optgr2==1) then
614 :
615 408 : do ll=0,l_size-1
616 1441 : do ilm=ll**2+1,min((ll+1)**2,lm_size)
617 2069558 : do ic=1,nfgd
618 : gylmgr2(1,ic,ilm)=gfact(ic,ll)*ylmrgr(4,ilm,ic) &
619 : & +dgfact(ic,ll)*(ylmr(ilm,ic)+two*rfgd(1,ic)*ylmrgr(1,ilm,ic)) &
620 2068214 : & +d2gfact(ic,ll)*ylmr(ilm,ic)*rfgd(1,ic)*rfgd(1,ic)
621 : gylmgr2(2,ic,ilm)=gfact(ic,ll)*ylmrgr(5,ilm,ic) &
622 : & +dgfact(ic,ll)*(ylmr(ilm,ic)+two*rfgd(2,ic)*ylmrgr(2,ilm,ic)) &
623 2068214 : & +d2gfact(ic,ll)*ylmr(ilm,ic)*rfgd(2,ic)*rfgd(2,ic)
624 : gylmgr2(3,ic,ilm)=gfact(ic,ll)*ylmrgr(6,ilm,ic) &
625 : & +dgfact(ic,ll)*(ylmr(ilm,ic)+two*rfgd(3,ic)*ylmrgr(3,ilm,ic)) &
626 2068214 : & +d2gfact(ic,ll)*ylmr(ilm,ic)*rfgd(3,ic)*rfgd(3,ic)
627 : gylmgr2(4,ic,ilm)=gfact(ic,ll)*ylmrgr(7,ilm,ic) &
628 : & +dgfact(ic,ll)*(rfgd(3,ic)*ylmrgr(2,ilm,ic)+rfgd(2,ic)*ylmrgr(3,ilm,ic)) &
629 2068214 : & +d2gfact(ic,ll)*ylmr(ilm,ic)*rfgd(3,ic)*rfgd(2,ic)
630 : gylmgr2(5,ic,ilm)=gfact(ic,ll)*ylmrgr(8,ilm,ic) &
631 : & +dgfact(ic,ll)*(rfgd(3,ic)*ylmrgr(1,ilm,ic)+rfgd(1,ic)*ylmrgr(3,ilm,ic)) &
632 2068214 : & +d2gfact(ic,ll)*ylmr(ilm,ic)*rfgd(3,ic)*rfgd(1,ic)
633 : gylmgr2(6,ic,ilm)=gfact(ic,ll)*ylmrgr(9,ilm,ic) &
634 : & +dgfact(ic,ll)*(rfgd(1,ic)*ylmrgr(2,ilm,ic)+rfgd(2,ic)*ylmrgr(1,ilm,ic)) &
635 2069247 : & +d2gfact(ic,ll)*ylmr(ilm,ic)*rfgd(1,ic)*rfgd(2,ic)
636 : end do
637 : end do
638 : end do
639 :
640 : ! Special values at r-R=0 (supposing shapefunc(r)->C.r**l when r->0)
641 97 : if (izero>0) then
642 5024 : gylmgr2(1:6,izero,1:lm_size)=zero
643 61 : if (lm_size>=1) then
644 61 : arg=cc(3,1)/sqrt(four_pi)
645 244 : gylmgr2(1:3,izero,1)=arg
646 : end if
647 61 : if (lm_size>=2) then
648 61 : arg=cc(2,2)*sqrt(three/four_pi)
649 61 : gylmgr2(2,izero,2)=two*arg
650 61 : gylmgr2(4,izero,2)= arg
651 61 : if (lm_size>=3) then
652 61 : gylmgr2(1,izero,3)=two*arg
653 61 : gylmgr2(3,izero,3)=two*arg
654 : end if
655 61 : if (lm_size>=4) then
656 61 : gylmgr2(5,izero,4)=arg
657 61 : gylmgr2(6,izero,4)=arg
658 : end if
659 : end if
660 61 : if (lm_size>=5) then
661 61 : arg=cc(1,3)*sqrt(15._dp/four_pi)
662 61 : gylmgr2(6,izero,5)=arg
663 61 : if (lm_size>=6) gylmgr2(4,izero,6)=arg
664 61 : if (lm_size>=7) then
665 61 : gylmgr2(1,izero,7)= -arg/sqrt3
666 61 : gylmgr2(2,izero,7)= -arg/sqrt3
667 61 : gylmgr2(3,izero,7)=two*arg/sqrt3
668 : end if
669 61 : if (lm_size>=8) gylmgr2(5,izero,8)=arg
670 61 : if (lm_size>=9) then
671 61 : gylmgr2(1,izero,9)= arg
672 61 : gylmgr2(2,izero,9)=-arg
673 : end if
674 : end if
675 : end if
676 :
677 : end if
678 :
679 : !Memory deallocation
680 : !==========================================================
681 43573 : LIBPAW_DEALLOCATE(rnrm)
682 43573 : if (allocated(cc)) then
683 36946 : LIBPAW_DEALLOCATE(cc)
684 : end if
685 43573 : if (compute_gr0) then
686 43573 : LIBPAW_DEALLOCATE(gfact)
687 : end if
688 43573 : if (compute_gr1) then
689 8903 : LIBPAW_DEALLOCATE(dgfact)
690 : end if
691 43573 : if (compute_gr2) then
692 97 : LIBPAW_DEALLOCATE(d2gfact)
693 : end if
694 43573 : if (compute_gr1) then
695 8903 : LIBPAW_DEALLOCATE(rnrm_inv)
696 : end if
697 43573 : if (shape_type==3) then
698 35790 : LIBPAW_DEALLOCATE(alpha)
699 35790 : LIBPAW_DEALLOCATE(qq)
700 : end if
701 43573 : if (compute_gr0) then
702 43573 : LIBPAW_DEALLOCATE(ylmr)
703 : end if
704 43573 : if (compute_gr1) then
705 8903 : LIBPAW_DEALLOCATE(ylmrgr)
706 : end if
707 43592 : if (shape_type==-1) then
708 3 : if (compute_gr0) then
709 3 : LIBPAW_DEALLOCATE(shpfuncnum)
710 : end if
711 3 : if (compute_gr1) then
712 3 : LIBPAW_DEALLOCATE(dshpfuncnum)
713 : end if
714 3 : if (compute_gr2) then
715 0 : LIBPAW_DEALLOCATE(d2shpfuncnum)
716 : end if
717 : end if
718 :
719 : ! -----------------------------------------------------------------
720 : !Small functions related to analytical expression of shape function
721 : CONTAINS
722 : !!***
723 : ! ------------------------------------------------
724 : !!****f* m_paw_finegrid/shapefunc1
725 : ! shapefunc1 is g(x) (gaussian)
726 10632 : function shapefunc1(arg)
727 :
728 : real(dp) :: shapefunc1
729 : real(dp),intent(in) :: arg
730 10632 : shapefunc1=exp(-(arg/sigma)**lambda)
731 10632 : end function shapefunc1
732 : !!***
733 : ! ------------------------------------------------
734 : !!****f* m_paw_finegrid/shapefunc1_0
735 : ! shapefunc1_0 is g(x) (gaussian) for small x
736 4 : function shapefunc1_0(arg)
737 :
738 : real(dp) :: shapefunc1_0
739 : real(dp),intent(in) :: arg
740 4 : shapefunc1_0=one-(arg/sigma)**lambda+half*(arg/sigma)**(2*lambda)-(arg/sigma)**(3*lambda)/6._dp
741 4 : end function shapefunc1_0
742 : !!***
743 : ! ------------------------------------------------
744 : !!****f* m_paw_finegrid/shapefunc2
745 : ! shapefunc2 is g(x) (sinc2)
746 11596604 : function shapefunc2(arg)
747 :
748 : real(dp) :: shapefunc2
749 : real(dp),intent(in) :: arg
750 11596604 : shapefunc2=(sin(pi_over_rshp*arg)/(pi_over_rshp*arg))**2
751 11596604 : end function shapefunc2
752 : !!***
753 : ! ------------------------------------------------
754 : !!****f* m_paw_finegrid/shapefunc2_0
755 : ! shapefunc2_0 is g(x) (sinc2) for small x
756 1734 : function shapefunc2_0(arg)
757 :
758 : real(dp) :: shapefunc2_0
759 : real(dp),intent(in) :: arg
760 1734 : shapefunc2_0=one-(pi_over_rshp*arg)**2/three+two*(pi_over_rshp*arg)**4/45._dp
761 1734 : end function shapefunc2_0
762 : !!***
763 : ! ------------------------------------------------
764 : !!****f* m_paw_finegrid/shapefunc3
765 : ! shapefunc3 is g(x) (Bessel)
766 89216288 : function shapefunc3(jbes1,jbes2,argl)
767 :
768 : integer,intent(in) :: argl
769 : real(dp) :: shapefunc3
770 : real(dp),intent(in) :: jbes1,jbes2
771 89216288 : shapefunc3= alpha(1,1+argl)*jbes1+alpha(2,1+argl)*jbes2
772 89216288 : end function shapefunc3
773 : !!***
774 : ! ------------------------------------------------
775 : !!****f* m_paw_finegrid/dshpfunc1
776 : ! dshpfunc1(x) is g_prime(x) (gaussian)
777 5316 : function dshpfunc1(arg)
778 :
779 : real(dp) :: dshpfunc1
780 : real(dp),intent(in) :: arg
781 5316 : dshpfunc1=-lambda/sigma*(arg/sigma)**(lambda-1)*exp(-(arg/sigma)**lambda)
782 5316 : end function dshpfunc1
783 : !!***
784 : ! ------------------------------------------------
785 : !!****f* m_paw_finegrid/dshpfunc1_ovr_0
786 : ! dshpfunc1_ovr_0(x) is g_prime(x)/x (gaussian) for small x and lambda>2
787 0 : function dshpfunc1_ovr_0(arg)
788 :
789 : real(dp) :: dshpfunc1_ovr_0
790 : real(dp),intent(in) :: arg
791 0 : dshpfunc1_ovr_0=-lambda/sigma**2*((arg/sigma)**(lambda-2)-(arg/sigma)**(2*lambda-2))
792 0 : end function dshpfunc1_ovr_0
793 : !!***
794 : ! ------------------------------------------------
795 : !!****f* m_paw_finegrid/dshpfunc1_ovr_0_2
796 : ! dshpfunc1_ovr_0_2(x) is g_prime(x)/x (gaussian) for small x and lambda=2
797 2 : function dshpfunc1_ovr_0_2(arg)
798 :
799 : real(dp) :: dshpfunc1_ovr_0_2
800 : real(dp),intent(in) :: arg
801 2 : dshpfunc1_ovr_0_2=-two/sigma**2*(one-(arg/sigma)**2+half*(arg/sigma)**4)
802 2 : end function dshpfunc1_ovr_0_2
803 : !!***
804 : ! ------------------------------------------------
805 : !!****f* m_paw_finegrid/dshpfunc2
806 : ! dshpfunc2(x) is g_prime(x) (sinc2)
807 6213476 : function dshpfunc2(arg)
808 :
809 : real(dp) :: dshpfunc2
810 : real(dp),intent(in) :: arg
811 : dshpfunc2=two*pi_over_rshp*sin(pi_over_rshp*arg)/(pi_over_rshp*arg)**3&
812 6213476 : & *(pi_over_rshp*arg*cos(pi_over_rshp*arg)-sin(pi_over_rshp*arg))
813 6213476 : end function dshpfunc2
814 : !!***
815 : ! ------------------------------------------------
816 : !!****f* m_paw_finegrid/dshpfunc2_ovr_0
817 : ! dshpfunc2_ovr_0(x) is g_prime(x)/x (sinc2) for small x
818 928 : function dshpfunc2_ovr_0(arg)
819 :
820 : real(dp) :: dshpfunc2_ovr_0
821 : real(dp),intent(in) :: arg
822 928 : dshpfunc2_ovr_0=-two*pi_over_rshp**2/3._dp+8._dp*pi_over_rshp**4*arg**2/45._dp
823 928 : end function dshpfunc2_ovr_0
824 : !!***
825 : ! ------------------------------------------------
826 : !!****f* m_paw_finegrid/dshpfunc3
827 : ! dshpfunc3(x) is g_prime(x) (Bessel)
828 14037703 : function dshpfunc3(jbesp1,jbesp2,argl)
829 :
830 : integer :: argl
831 : real(dp) :: dshpfunc3
832 : real(dp),intent(in) :: jbesp1,jbesp2
833 : dshpfunc3= alpha(1,1+argl)*qq(1,1+argl)*jbesp1 &
834 14037703 : & +alpha(2,1+argl)*qq(2,1+argl)*jbesp2
835 14037703 : end function dshpfunc3
836 : !!***
837 : ! ------------------------------------------------
838 : !!****f* m_paw_finegrid/d2shpfunc1
839 : ! d2shpfunc1(x) is g_prime_prime(x) (gaussian)
840 0 : function d2shpfunc1(arg)
841 :
842 : real(dp) :: d2shpfunc1
843 : real(dp),intent(in) :: arg
844 : d2shpfunc1=lambda/(sigma**2)*(lambda*(arg/sigma)**(2*lambda-2) &
845 0 : & -(lambda-1)*(arg/sigma)**(lambda-2))*exp(-(arg/sigma)**lambda)
846 0 : end function d2shpfunc1
847 : !!***
848 : ! ------------------------------------------------
849 : !!****f* m_paw_finegrid/d2shpfunc1_ovr2_0
850 : ! d2shpfunc1_ovr2_0(x) is (g_prime_prime(x)-g_prime(x)/x)/x**2 (gaussian) for small x and lambda>4
851 0 : function d2shpfunc1_ovr2_0(arg)
852 :
853 : real(dp) :: d2shpfunc1_ovr2_0
854 : real(dp),intent(in) :: arg
855 : d2shpfunc1_ovr2_0=-lambda/(sigma**4)*((lambda-2)*(arg/sigma)**(lambda-4) &
856 0 : & -(lambda-1)*two*(arg/sigma)**(2*lambda-4))
857 0 : end function d2shpfunc1_ovr2_0
858 : !!***
859 : ! ------------------------------------------------
860 : !!****f* m_paw_finegrid/d2shpfunc1_ovr2_0_2
861 : ! d2shpfunc1_ovr2_0_2(x) is (g_prime_prime(x)-g_prime(x)/x)/x**2 (gaussian) for small x and lambda==2
862 0 : function d2shpfunc1_ovr2_0_2(arg)
863 :
864 : real(dp) :: d2shpfunc1_ovr2_0_2
865 : real(dp),intent(in) :: arg
866 0 : d2shpfunc1_ovr2_0_2=four/(sigma**4)*(one-(arg/sigma)**2)
867 0 : end function d2shpfunc1_ovr2_0_2
868 : !!***
869 : ! ------------------------------------------------
870 : !!****f* m_paw_finegrid/d2shpfunc1_ovr2_0_3
871 : ! d2shpfunc1_ovr2_0_3(x) is (g_prime_prime(x)-g_prime(x)/x)/x**2 (gaussian) for small x and lambda==3
872 0 : function d2shpfunc1_ovr2_0_3(arg)
873 :
874 : real(dp) :: d2shpfunc1_ovr2_0_3
875 : real(dp),intent(in) :: arg
876 0 : d2shpfunc1_ovr2_0_3=-three/arg/sigma**3+12._dp*arg**2/sigma**6-half*21._dp*arg**5/sigma**9
877 0 : end function d2shpfunc1_ovr2_0_3
878 : !!***
879 : ! ------------------------------------------------
880 : !!****f* m_paw_finegrid/d2shpfunc1_ovr2_0_4
881 : ! d2shpfunc1_ovr2_0_4(x) is (g_prime_prime(x)-g_prime(x)/x)/x**2 (gaussian) for small x and lambda==4
882 0 : function d2shpfunc1_ovr2_0_4(arg)
883 :
884 : real(dp) :: d2shpfunc1_ovr2_0_4
885 : real(dp),intent(in) :: arg
886 0 : d2shpfunc1_ovr2_0_4=-8._dp/(sigma**4)*(one-three*(arg/sigma)**4)
887 0 : end function d2shpfunc1_ovr2_0_4
888 : !!***
889 : ! ------------------------------------------------
890 : !!****f* m_paw_finegrid/d2shpfunc2
891 : ! d2shpfunc2(x) is g_prime_prime(x) (sinc2)
892 116066 : function d2shpfunc2(arg)
893 :
894 : real(dp) :: d2shpfunc2
895 : real(dp),intent(in) :: arg
896 : d2shpfunc2=two/(pi_over_rshp**2*arg**4)* &
897 : & (pi_over_rshp**2*arg**2*(cos(pi_over_rshp*arg))**2 &
898 : & +(three-pi_over_rshp**2*arg**2)*(sin(pi_over_rshp*arg))**2 &
899 116066 : & -four*pi_over_rshp*arg*cos(pi_over_rshp*arg)*sin(pi_over_rshp*arg))
900 116066 : end function d2shpfunc2
901 : !!***
902 : ! ------------------------------------------------
903 : !!****f* m_paw_finegrid/d2shpfunc2_ovr2_0
904 : ! d2shpfunc2_ovr2_0(x) is (g_prime_prime(x)-g_prime(x)/x)/x**2 (sinc2) for small x
905 56 : function d2shpfunc2_ovr2_0(arg)
906 :
907 : real(dp) :: d2shpfunc2_ovr2_0
908 : real(dp),intent(in) :: arg
909 : d2shpfunc2_ovr2_0=16._dp/45._dp*pi_over_rshp**4-8._dp/105._dp*pi_over_rshp**6*arg**2 &
910 56 : & +41._dp/6300._dp*pi_over_rshp**8*arg**4
911 56 : end function d2shpfunc2_ovr2_0
912 : !!***
913 : ! ------------------------------------------------
914 : !!****f* m_paw_finegrid/d2shpfunc3
915 : ! d2shpfunc3(x) is g_prime_prime(x) (Bessel)
916 79522 : function d2shpfunc3(jbespp1,jbespp2,argl)
917 :
918 : integer,intent(in) :: argl
919 : real(dp) :: d2shpfunc3
920 : real(dp),intent(in) :: jbespp1,jbespp2
921 : d2shpfunc3= alpha(1,1+argl)*(qq(1,1+argl)**2)*jbespp1 &
922 79522 : & +alpha(2,1+argl)*(qq(2,1+argl)**2)*jbespp2
923 79522 : end function d2shpfunc3
924 : ! ------------------------------------------------
925 :
926 : end subroutine pawgylm
927 : !!***
928 :
929 : !----------------------------------------------------------------------
930 :
931 : !!****f* m_paw_finegrid/pawgylmg
932 : !! NAME
933 : !! pawgylmg
934 : !!
935 : !! FUNCTION
936 : !! PAW: Compute Fourier transform of each g_l(r).Y_lm(r) function
937 : !!
938 : !! INPUTS
939 : !! gprimd(3,3)=dimensional reciprocal space primitive translations
940 : !! kg(3,npw)=integer coordinates of planewaves in basis sphere for this k point.
941 : !! kpg(npw,nkpg)= (k+G) components (only if useylm=1)
942 : !! kpt(3)=reduced coordinates of k point
943 : !! lmax=1+max. value of l angular momentum
944 : !! nkpg=second dimension of kpg_k (0 if useylm=0)
945 : !! npw=number of planewaves in basis sphere
946 : !! ntypat=number of types of atoms
947 : !! pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data
948 : !! ylm(npw,lmax**2)=real spherical harmonics for each G and k point
949 : !!
950 : !! OUTPUT
951 : !! gylmg(npw,lmax**2,ntypat)=Fourier transform of each g_l(r).Y_lm(r) function
952 : !!
953 : !! SOURCE
954 :
955 0 : subroutine pawgylmg(gprimd,gylmg,kg,kpg,kpt,lmax,nkpg,npw,ntypat,pawtab,ylm)
956 :
957 : !Arguments ------------------------------------
958 : !scalars
959 : integer,intent(in) :: lmax,nkpg,npw,ntypat
960 : !arrays
961 : integer,intent(in) :: kg(3,npw)
962 : real(dp),intent(in) :: gprimd(3,3),kpg(npw,nkpg),kpt(3)
963 : real(dp),intent(in) :: ylm(npw,lmax**2)
964 : type(pawtab_type),intent(in) :: pawtab(ntypat)
965 :
966 : real(dp),intent(out) :: gylmg(npw,lmax**2,ntypat)
967 :
968 : !Local variables-------------------------------
969 : !scalars
970 : integer :: ig,ilm,itypat,ll,l0,mm,mqgrid
971 : real(dp) :: kpg1,kpg2,kpg3,kpgc1,kpgc2,kpgc3
972 : !arrays
973 0 : real(dp),allocatable :: glg(:),qgrid(:),kpgnorm(:),shpf(:,:),work(:)
974 :
975 : ! *************************************************************************
976 :
977 : !Get |k+G|:
978 0 : LIBPAW_ALLOCATE(kpgnorm,(npw))
979 0 : if (nkpg<3) then
980 0 : do ig=1,npw
981 0 : kpg1=kpt(1)+dble(kg(1,ig));kpg2=kpt(2)+dble(kg(2,ig));kpg3=kpt(3)+dble(kg(3,ig))
982 0 : kpgc1=kpg1*gprimd(1,1)+kpg2*gprimd(1,2)+kpg3*gprimd(1,3)
983 0 : kpgc2=kpg1*gprimd(2,1)+kpg2*gprimd(2,2)+kpg3*gprimd(2,3)
984 0 : kpgc3=kpg1*gprimd(3,1)+kpg2*gprimd(3,2)+kpg3*gprimd(3,3)
985 0 : kpgnorm(ig)=sqrt(kpgc1*kpgc1+kpgc2*kpgc2+kpgc3*kpgc3)
986 : end do
987 : else
988 0 : do ig=1,npw
989 0 : kpgc1=kpg(ig,1)*gprimd(1,1)+kpg(ig,2)*gprimd(1,2)+kpg(ig,3)*gprimd(1,3)
990 0 : kpgc2=kpg(ig,1)*gprimd(2,1)+kpg(ig,2)*gprimd(2,2)+kpg(ig,3)*gprimd(2,3)
991 0 : kpgc3=kpg(ig,1)*gprimd(3,1)+kpg(ig,2)*gprimd(3,2)+kpg(ig,3)*gprimd(3,3)
992 0 : kpgnorm(ig)=sqrt(kpgc1*kpgc1+kpgc2*kpgc2+kpgc3*kpgc3)
993 : end do
994 : end if
995 :
996 0 : LIBPAW_ALLOCATE(glg,(npw))
997 0 : LIBPAW_ALLOCATE(work,(npw))
998 :
999 : !Loop over types of atoms
1000 0 : do itypat=1,ntypat
1001 :
1002 0 : mqgrid=pawtab(itypat)%mqgrid_shp
1003 0 : LIBPAW_ALLOCATE(qgrid,(mqgrid))
1004 0 : LIBPAW_ALLOCATE(shpf,(mqgrid,2))
1005 0 : qgrid(1:mqgrid)=pawtab(itypat)%qgrid_shp(1:mqgrid)
1006 :
1007 : ! Loops over (l,m) values
1008 0 : do ll=0,pawtab(itypat)%lcut_size-1
1009 0 : l0=ll**2+ll+1
1010 :
1011 0 : shpf(1:mqgrid,1:2)=pawtab(itypat)%shapefncg(1:mqgrid,1:2,1+ll)
1012 0 : call paw_uniform_splfit(qgrid,work,shpf,0,kpgnorm,glg,mqgrid,npw)
1013 :
1014 0 : do mm=-ll,ll
1015 0 : ilm=l0+mm
1016 :
1017 0 : gylmg(1:npw,ilm,itypat)=ylm(1:npw,ilm)*glg(1:npw)
1018 :
1019 : ! End loops over (l,m) values
1020 : end do
1021 : end do
1022 :
1023 : ! End loop over atom types
1024 0 : LIBPAW_DEALLOCATE(qgrid)
1025 0 : LIBPAW_DEALLOCATE(shpf)
1026 : end do
1027 :
1028 0 : LIBPAW_DEALLOCATE(kpgnorm)
1029 0 : LIBPAW_DEALLOCATE(glg)
1030 0 : LIBPAW_DEALLOCATE(work)
1031 :
1032 0 : end subroutine pawgylmg
1033 : !!***
1034 :
1035 : !----------------------------------------------------------------------
1036 :
1037 : !!****f* m_paw_finegrid/pawrfgd_fft
1038 : !! NAME
1039 : !! pawrfgd_fft
1040 : !!
1041 : !! FUNCTION
1042 : !! Determine each point of the (fine) rectangular grid
1043 : !! around a given atom and compute r-R vectors.
1044 : !! R is the position of the atom.
1045 : !!
1046 : !! INPUTS
1047 : !! [fft_distrib(n3)]= (optional) index of processes which own fft planes in 3rd dimension
1048 : !! [fft_index(n3)]= (optional) local fft indexes for current process
1049 : !! gmet(3,3)=reciprocal space metric tensor in bohr**-2
1050 : !! [me_fft]= (optional) my rank in the FFT MPI communicator
1051 : !! n1,n2,n3= sizes of the FFT grid (entire simulation cell)
1052 : !! rcut= radius of the sphere around the atom
1053 : !! rprimd(3,3)=dimensional primitive translations in real space (bohr)
1054 : !! ucvol= unit cell volume
1055 : !! xred(3)= reduced coordinates of the atom
1056 : !!
1057 : !! OUTPUT
1058 : !! ifftsph(nfgd)= FFT index (fine grid) of the points in the sphere around current atom
1059 : !! nfgd= number of points in the sphere around current atom
1060 : !! rfgd(3,nfgd)= cartesian coordinates of r-R.
1061 : !!
1062 : !! SOURCE
1063 :
1064 4937 : subroutine pawrfgd_fft(ifftsph,gmet,n1,n2,n3,nfgd,rcut,rfgd,rprimd,ucvol,xred, &
1065 4937 : & fft_distrib,fft_index,me_fft) ! optional arguments
1066 :
1067 : !Arguments ---------------------------------------------
1068 : !scalars
1069 : integer,intent(in) :: n1,n2,n3
1070 : integer,intent(out) :: nfgd
1071 : integer,optional,intent(in) :: me_fft
1072 : real(dp),intent(in) :: rcut,ucvol
1073 : !arrays
1074 : integer,target,optional,intent(in) :: fft_distrib(n3),fft_index(n3)
1075 : integer,allocatable,intent(out) :: ifftsph(:)
1076 : real(dp),intent(in) :: gmet(3,3),rprimd(3,3),xred(3)
1077 : real(dp),allocatable,intent(out) :: rfgd(:,:)
1078 :
1079 : !Local variables ------------------------------
1080 : !scalars
1081 : integer :: i1,i2,i3,ifft_local,ix,iy,iz,izloc,me_fft_,n1a,n1b,n2a,n2b,n3a,n3b,ncmax
1082 : real(dp) :: dif,difx,dify,difz,rr1,rr2,rr3,r2,r2cut,rx,ry,rz
1083 : character(len=500) :: msg
1084 : !arrays
1085 4937 : integer,allocatable :: ifftsph_tmp(:)
1086 4937 : integer,pointer :: fft_distrib_(:),fft_index_(:)
1087 4937 : real(dp),allocatable :: rfgd_tmp(:,:)
1088 :
1089 : ! *************************************************************************
1090 :
1091 : !Define a "box" around the atom
1092 4937 : r2cut=1.0000001_dp*rcut**2
1093 4937 : rr1=sqrt(r2cut*gmet(1,1))
1094 4937 : rr2=sqrt(r2cut*gmet(2,2))
1095 4937 : rr3=sqrt(r2cut*gmet(3,3))
1096 :
1097 : !Computing lower bounds of the box
1098 : !If dif<0, int(dif*n) will give an integer above what we want, so we subtract one
1099 4937 : dif=xred(1)-rr1
1100 4937 : n1a=int(dif*n1)+1
1101 3307 : if (dif<0) n1a = n1a-1
1102 4937 : dif=xred(2)-rr2
1103 4937 : n2a=int(dif*n2)+1
1104 4937 : if (dif<0) n2a = n2a-1
1105 4937 : dif=xred(3)-rr3
1106 4937 : n3a=int(dif*n3)+1
1107 4937 : if (dif<0) n3a = n3a-1
1108 :
1109 : !Computing upper bounds of the box
1110 : !If dif>0, int(dif*n) will give an integer below what we want, so we add one
1111 4937 : dif=xred(1)+rr1
1112 4937 : n1b=int(dif*n1)+1
1113 4937 : if (dif>0) n1b = n1b+1
1114 4937 : dif=xred(2)+rr2
1115 4937 : n2b=int(dif*n2)+1
1116 4937 : if (dif>0) n2b = n2b+1
1117 4937 : dif=xred(3)+rr3
1118 4937 : n3b=int(dif*n3)+1
1119 4937 : if (dif>0) n3b = n3b+1
1120 :
1121 : !Get the distrib associated with this fft_grid
1122 4937 : if (present(fft_distrib).and.present(fft_index).and.present(me_fft)) then
1123 4937 : me_fft_=me_fft ; fft_distrib_ => fft_distrib ; fft_index_ => fft_index
1124 : else
1125 0 : me_fft_=0
1126 0 : LIBPAW_POINTER_ALLOCATE(fft_distrib_,(n3))
1127 0 : LIBPAW_POINTER_ALLOCATE(fft_index_,(n3))
1128 0 : fft_distrib_=0;fft_index_=(/(i3,i3=1,n3)/)
1129 : end if
1130 :
1131 : !Temporary allocate "large" arrays
1132 4937 : ncmax=1+int(1.5_dp*(n1*n2*n3)*four_pi/(three*ucvol)*rcut**3)
1133 14811 : LIBPAW_ALLOCATE(ifftsph_tmp,(ncmax))
1134 14811 : LIBPAW_ALLOCATE(rfgd_tmp,(3,ncmax))
1135 :
1136 : !Set number of points to zero
1137 4937 : nfgd=0
1138 :
1139 : !Loop over FFT points
1140 85515 : do i3=n3a,n3b
1141 80578 : iz=modulo(i3,n3)
1142 80578 : if (iz<0.or.iz>n3-1) then
1143 0 : msg='iz<0 or iz>n3-1'
1144 0 : LIBPAW_ERROR(msg)
1145 : end if
1146 85515 : if (fft_distrib_(iz+1)==me_fft_) then
1147 75838 : izloc=fft_index_(iz+1) - 1
1148 75838 : if (izloc<0.or.izloc>n3-1) then
1149 0 : msg='izloc<0 or izloc>n3-1'
1150 0 : LIBPAW_ERROR(msg)
1151 : end if
1152 75838 : difz=dble(i3)/dble(n3)-xred(3)
1153 1501051 : do i2=n2a,n2b
1154 1425213 : iy=modulo(i2,n2)
1155 1425213 : if (iy<0.or.iy>n2-1) then
1156 0 : msg='iy<0 or iy>n2-1'
1157 0 : LIBPAW_ERROR(msg)
1158 : end if
1159 1425213 : dify=dble(i2)/dble(n2)-xred(2)
1160 32217825 : do i1=n1a,n1b
1161 30716774 : ix=modulo(i1,n1)
1162 30716774 : if (ix<0.or.ix>n1-1) then
1163 0 : msg='ix<0 or ix>n1-1'
1164 0 : LIBPAW_ERROR(msg)
1165 : end if
1166 30716774 : difx=dble(i1)/dble(n1)-xred(1)
1167 :
1168 : ! Compute r-R
1169 30716774 : rx=difx*rprimd(1,1)+dify*rprimd(1,2)+difz*rprimd(1,3)
1170 30716774 : ry=difx*rprimd(2,1)+dify*rprimd(2,2)+difz*rprimd(2,3)
1171 30716774 : rz=difx*rprimd(3,1)+dify*rprimd(3,2)+difz*rprimd(3,3)
1172 30716774 : r2=rx**2+ry**2+rz**2
1173 :
1174 : ! Select matching points
1175 32141987 : if (r2 <= r2cut) then
1176 9675685 : ifft_local=1+ix+n1*(iy+n2*izloc)
1177 9675685 : if (ifft_local<1.or.ifft_local>n1*n2*n3) then
1178 0 : msg='ifft_local<1 or ifft_local>n1*n2*n3'
1179 0 : LIBPAW_ERROR(msg)
1180 : end if
1181 9675685 : if (ifft_local>0) then
1182 9675685 : nfgd=nfgd+1
1183 9675685 : if (nfgd>ncmax) then
1184 0 : msg='Number of fft points around atom exceeds max. allowed!'
1185 0 : LIBPAW_BUG(msg)
1186 : end if
1187 9675685 : rfgd_tmp(1,nfgd)=rx
1188 9675685 : rfgd_tmp(2,nfgd)=ry
1189 9675685 : rfgd_tmp(3,nfgd)=rz
1190 9675685 : ifftsph_tmp(nfgd)=ifft_local
1191 : end if
1192 : end if
1193 :
1194 : ! End of loops
1195 : end do
1196 : end do
1197 : end if
1198 : end do
1199 :
1200 : !Now fill output arrays
1201 4937 : if (allocated(ifftsph)) then
1202 0 : LIBPAW_DEALLOCATE(ifftsph)
1203 : end if
1204 4937 : if (allocated(rfgd)) then
1205 0 : LIBPAW_DEALLOCATE(rfgd)
1206 : end if
1207 14811 : LIBPAW_ALLOCATE(ifftsph,(nfgd))
1208 14811 : LIBPAW_ALLOCATE(rfgd,(3,nfgd))
1209 9680622 : ifftsph(1:nfgd)=ifftsph_tmp(1:nfgd)
1210 38707677 : rfgd(1:3,1:nfgd)=rfgd_tmp(1:3,1:nfgd)
1211 :
1212 : !Release temporary memory
1213 4937 : LIBPAW_DEALLOCATE(ifftsph_tmp)
1214 4937 : LIBPAW_DEALLOCATE(rfgd_tmp)
1215 4937 : if (.not.present(fft_distrib).or..not.present(fft_index)) then
1216 0 : LIBPAW_POINTER_DEALLOCATE(fft_distrib_)
1217 0 : LIBPAW_POINTER_DEALLOCATE(fft_index_)
1218 : end if
1219 :
1220 4937 : end subroutine pawrfgd_fft
1221 : !!***
1222 :
1223 : !----------------------------------------------------------------------
1224 :
1225 : !!****f* m_paw_finegrid/pawrfgd_wvl
1226 : !! NAME
1227 : !! pawrfgd_wvl
1228 : !!
1229 : !! FUNCTION
1230 : !! Determine each point of the (fine) rectangular grid
1231 : !! around a given atom and compute r-R vectors.
1232 : !! R is the position of the atom.
1233 : !!
1234 : !! INPUTS
1235 : !! geocode= code for geometry (boundary conditions)
1236 : !! hh(3)=fine grid spacing
1237 : !! i3s= TO BE COMPLETED
1238 : !! n1,n2,n3= TO BE COMPLETED
1239 : !! n1i,n2i,n3pi= TO BE COMPLETED
1240 : !! rcut= radius of the sphere around the atom
1241 : !! rloc= cut-off radius for local psp?
1242 : !! shift= TO BE COMPLETED
1243 : !! xred(3)= cartesian coordinates of the atom
1244 : !!
1245 : !! OUTPUT
1246 : !! ifftsph(nfgd)= FFT index (fine grid) of the points in the sphere around current atom
1247 : !! nfgd= number of points in the sphere around current atom
1248 : !! rfgd(3,nfgd)= cartesian coordinates of r-R.
1249 : !!
1250 : !! SOURCE
1251 :
1252 0 : subroutine pawrfgd_wvl(geocode,hh,ifftsph,i3s,n1,n1i,n2,n2i,n3,n3pi,&
1253 : & nfgd,rcut,rloc,rfgd,shift,xcart)
1254 :
1255 : !Arguments ---------------------------------------------
1256 : !scalars
1257 : integer,intent(in) :: i3s,n1,n1i,n2,n2i,n3,n3pi,shift
1258 : integer,intent(out) :: nfgd
1259 : real(dp),intent(in) :: rcut,rloc
1260 : character(1),intent(in) :: geocode
1261 : !arrays
1262 : integer,allocatable,intent(out) :: ifftsph(:)
1263 : real(dp),intent(in) :: hh(3),xcart(3)
1264 : real(dp),allocatable,intent(out) :: rfgd(:,:)
1265 :
1266 : !Local variables ------------------------------
1267 : !scalars
1268 : integer :: i1,i2,i3,iex,iey,iez,ind,isx,isy,isz,j1,j2,j3
1269 : integer :: nbl1,nbr1,nbl2,nbr2,nbl3,nbr3,ncmax
1270 : logical :: gox,goy,goz,perx,pery,perz
1271 : real(dp) :: cutoff,r2,r2cut,rx,ry,rz,xx,yy,zz
1272 : !arrays
1273 0 : integer,allocatable :: ifftsph_tmp(:)
1274 0 : real(dp),allocatable :: rfgd_tmp(:,:)
1275 :
1276 : ! *************************************************************************
1277 :
1278 : !Data for periodicity in the three directions
1279 0 : perx=(geocode/='F')
1280 0 : pery=(geocode=='P')
1281 0 : perz=(geocode/='F')
1282 0 : call my_ext_buffers(perx,nbl1,nbr1)
1283 0 : call my_ext_buffers(pery,nbl2,nbr2)
1284 0 : call my_ext_buffers(perz,nbl3,nbr3)
1285 :
1286 : !Define a "box" around the atom
1287 0 : cutoff=10.d0*rloc
1288 0 : r2cut=1.0000001_dp*rcut**2
1289 0 : rx=xcart(1)
1290 0 : ry=xcart(2)
1291 0 : rz=xcart(3)
1292 0 : isx=floor((rx-cutoff)/hh(1))
1293 0 : isy=floor((ry-cutoff)/hh(2))
1294 0 : isz=floor((rz-cutoff)/hh(3))
1295 0 : iex=ceiling((rx+cutoff)/hh(1))
1296 0 : iey=ceiling((ry+cutoff)/hh(2))
1297 0 : iez=ceiling((rz+cutoff)/hh(3))
1298 :
1299 : !Temporary allocate "large" arrays
1300 : ! use factor 1+int(1.1*, for safety reasons
1301 0 : ncmax=1
1302 0 : if (n3pi>0) ncmax=1+int((rcut/hh(1)+1.0)*(rcut/hh(2)+1.0)*(rcut/hh(3)+1.0)*four_pi/three)
1303 0 : LIBPAW_ALLOCATE(ifftsph_tmp,(ncmax))
1304 0 : LIBPAW_ALLOCATE(rfgd_tmp,(3,ncmax))
1305 :
1306 : !Set number of points to zero
1307 0 : nfgd=0
1308 :
1309 : !Loop over WVL points
1310 0 : do i3=isz,iez
1311 0 : zz=real(i3,kind=8)*hh(3)-rz
1312 0 : call my_ind_positions(perz,i3,n3,j3,goz)
1313 0 : j3=j3+nbl3+1
1314 0 : do i2=isy,iey
1315 0 : yy=real(i2,kind=8)*hh(2)-ry
1316 0 : call my_ind_positions(pery,i2,n2,j2,goy)
1317 0 : do i1=isx,iex
1318 0 : xx=real(i1,kind=8)*hh(1)-rx
1319 0 : call my_ind_positions(perx,i1,n1,j1,gox)
1320 0 : r2=xx**2+yy**2+zz**2
1321 0 : if (j3>=i3s.and.j3<=i3s+n3pi-1.and.goy.and.gox) then
1322 :
1323 : ! Select matching points
1324 0 : if (r2<=r2cut) then
1325 0 : ind=j1+1+nbl1+(j2+nbl2)*n1i+(j3-i3s)*n1i*n2i
1326 0 : nfgd=nfgd+1
1327 0 : rfgd_tmp(:,nfgd)=[xx,yy,zz]
1328 0 : ifftsph_tmp(nfgd)=shift+ind
1329 : end if
1330 :
1331 : ! End of loops
1332 : end if
1333 : end do
1334 : end do
1335 : end do
1336 :
1337 : !Now fill output arrays
1338 0 : if (allocated(ifftsph)) then
1339 0 : LIBPAW_DEALLOCATE(ifftsph)
1340 : end if
1341 0 : if (allocated(rfgd)) then
1342 0 : LIBPAW_DEALLOCATE(rfgd)
1343 : end if
1344 0 : LIBPAW_ALLOCATE(ifftsph,(nfgd))
1345 0 : LIBPAW_ALLOCATE(rfgd,(3,nfgd))
1346 0 : ifftsph(1:nfgd)=ifftsph_tmp(1:nfgd)
1347 0 : rfgd(1:3,1:nfgd)=rfgd_tmp(1:3,1:nfgd)
1348 :
1349 : !Release temporary memory
1350 0 : LIBPAW_DEALLOCATE(ifftsph_tmp)
1351 0 : LIBPAW_DEALLOCATE(rfgd_tmp)
1352 :
1353 : !*********************************************************************
1354 : !Small functions related to boundary conditions
1355 : contains
1356 : !!***
1357 : ! ------------------------------------------------
1358 : !!****f* m_paw_finegrid/my_ind_positions
1359 0 : subroutine my_ind_positions(periodic,i,n,j,go)
1360 :
1361 : integer,intent(in) :: i,n
1362 : logical,intent(in) :: periodic
1363 : integer,intent(out) :: j
1364 : logical,intent(out) :: go
1365 0 : if (periodic) then
1366 0 : j=modulo(i,2*n+2) ; go=.true.
1367 : else
1368 0 : j=i ; go=(i>=-14.and.i<=2*n+16)
1369 : end if
1370 0 : end subroutine my_ind_positions
1371 : !!***
1372 : ! ------------------------------------------------
1373 : !!****f* m_paw_finegrid/my_ext_buffers
1374 0 : subroutine my_ext_buffers(periodic,nl,nr)
1375 :
1376 : logical, intent(in) :: periodic
1377 : integer, intent(out) :: nl,nr
1378 0 : if (periodic) then
1379 0 : nl=0 ; nr=0
1380 : else
1381 0 : nl=14 ; nr=15
1382 : end if
1383 : end subroutine my_ext_buffers
1384 : ! ------------------------------------------------
1385 :
1386 : end subroutine pawrfgd_wvl
1387 : !!***
1388 :
1389 : !----------------------------------------------------------------------
1390 :
1391 : !!****f* m_paw_finegrid/pawexpiqr
1392 : !! NAME
1393 : !! pawexpiqr
1394 : !!
1395 : !! FUNCTION
1396 : !! Compute exp(i.q.r) for each point of the (fine) rectangular grid
1397 : !! around a given atomic site. R is the position of the atom.
1398 : !! Used for the determination of phonons at non-zero q wavevector.
1399 : !!
1400 : !! INPUTS
1401 : !! gprimd(3,3)= dimensional primitive translations for reciprocal space
1402 : !! nfgd= number of (fine grid) FFT points in the paw sphere around current atom
1403 : !! qphon(3)= wavevector of the phonon
1404 : !! rfgd(3,nfgd)= coordinates of r-R on the fine grid around current atom
1405 : !! xred(3)= reduced atomic coordinates
1406 : !!
1407 : !! OUTPUT
1408 : !! expiqr(2,nfgd)= exp(i.q.r) around the current atom
1409 : !! Not allocated if q=0 !
1410 : !!
1411 : !! SOURCE
1412 :
1413 31310 : subroutine pawexpiqr(expiqr,gprimd,nfgd,qphon,rfgd,xred)
1414 :
1415 : !Arguments ---------------------------------------------
1416 : !scalars
1417 : integer,intent(in) :: nfgd
1418 : !arrays
1419 : real(dp),intent(in) :: gprimd(3,3),qphon(3),xred(3)
1420 : real(dp),intent(in) :: rfgd(:,:)
1421 : real(dp),intent(out) :: expiqr(2,nfgd)
1422 :
1423 : !Local variables ------------------------------
1424 : !scalars
1425 : integer :: ic
1426 : logical :: qne0
1427 : real(dp) :: phase,phase_xred,qx,qy,qz
1428 : character(len=500) :: msg
1429 : !arrays
1430 :
1431 : ! *************************************************************************
1432 :
1433 93930 : if (size(rfgd)/=3*nfgd) then
1434 0 : msg='rfgd array must be allocated!'
1435 0 : LIBPAW_BUG(msg)
1436 : end if
1437 :
1438 31310 : qne0=(qphon(1)**2+qphon(2)**2+qphon(3)**2>=1.d-15)
1439 :
1440 : !Compute q in cartesian coordinates
1441 31310 : if (qne0) then
1442 31310 : qx=gprimd(1,1)*qphon(1)+gprimd(1,2)*qphon(2)+gprimd(1,3)*qphon(3)
1443 31310 : qy=gprimd(2,1)*qphon(1)+gprimd(2,2)*qphon(2)+gprimd(2,3)*qphon(3)
1444 31310 : qz=gprimd(3,1)*qphon(1)+gprimd(3,2)*qphon(2)+gprimd(3,3)*qphon(3)
1445 31310 : phase_xred=two_pi*(qphon(1)*xred(1)+qphon(2)*xred(2)+qphon(3)*xred(3))
1446 : end if
1447 :
1448 : !Compute exp(i.q.r)
1449 : if (qne0) then
1450 24069430 : do ic=1,nfgd
1451 24038120 : phase=two_pi*(qx*rfgd(1,ic)+qy*rfgd(2,ic)+qz*rfgd(3,ic)) + phase_xred
1452 24038120 : expiqr(1,ic)=cos(phase)
1453 24069430 : expiqr(2,ic)=sin(phase)
1454 : end do
1455 : end if
1456 :
1457 31310 : end subroutine pawexpiqr
1458 : !!***
1459 :
1460 : !----------------------------------------------------------------------
1461 :
1462 : END MODULE m_paw_finegrid
1463 : !!***
|
