Branch data 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 [ - + ]: 43685 : 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 : 43685 : 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 : 43685 : real(dp),allocatable :: cc(:,:),d2gfact(:,:),d2shpfuncnum(:,:),dgfact(:,:)
113 : 43685 : real(dp),allocatable :: dshpfuncnum(:,:),gfact(:,:)
114 : 43685 : real(dp),allocatable :: rnrm(:),rnrm_inv(:),rnrm_sort(:)
115 : 43685 : real(dp),allocatable :: shpfuncnum(:,:),work(:),ylmr(:,:),ylmrgr(:,:,:)
116 : :
117 : : ! *************************************************************************
118 : :
119 [ + + - + : 43685 : if (optgr0==0.and.optgr1==0.and.optgr2==0) return
- - ]
120 [ + + ]: 43685 : if (nfgd==0) return
121 : :
122 : : !Compatibility test
123 : : !==========================================================
124 [ + + - + ]: 131007 : 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 [ + + - + ]: 43669 : 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 [ + + + - : 43669 : compute_gr0=(optgr0==1.or.optgr1==1.or.optgr2==1)
- - ]
141 [ + + - + ]: 43669 : compute_gr1=(optgr1==1.or.optgr2==1)
142 : 43669 : compute_gr2=(optgr2==1)
143 : 43669 : l_size=pawtab%lcut_size
144 : :
145 : : !Norms of vectors around the atom
146 [ - + + - : 131007 : LIBPAW_ALLOCATE(rnrm,(nfgd))
+ - + - -
+ - + ]
147 : 41171474 : izero=-1
148 [ + + ]: 41171474 : do ic=1,nfgd
149 : 41127805 : rnrm(ic)=sqrt(rfgd(1,ic)**2+rfgd(2,ic)**2+rfgd(3,ic)**2)
150 [ + + ]: 41171474 : if (rnrm(ic)<=tol10) izero=ic ! Has to be consistent with initylmr !!
151 : : end do
152 : :
153 : : !Initializations
154 [ + + + + : 399377235 : if (optgr0==1) gylm=zero
+ + ]
155 [ + + + + : 633253611 : if (optgr1==1) gylmgr=zero
+ + + + ]
156 [ + + + + : 23256252 : if (optgr2==1) gylmgr2=zero
+ + + + ]
157 : :
158 : : !Some definitions concerning shape function g_l(r)
159 : 43669 : shape_type=pawtab%shape_type
160 : 43669 : sigma=pawtab%shape_sigma;lambda=pawtab%shape_lambda
161 : 43669 : pi_over_rshp=pi/pawtab%rshp
162 : 43669 : rcut=tol12+pawtab%rshp
163 [ + + ]: 43669 : 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 [ + + ]: 43669 : 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 [ + + ]: 43669 : 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 [ + + ]: 43669 : if (izero>0) then
226 [ - + + - : 111042 : LIBPAW_ALLOCATE(cc,(3,l_size))
+ - + - -
+ - + ]
227 [ + + + + ]: 488606 : cc=zero
228 [ + + ]: 37014 : 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 [ + + ]: 149912 : do ll=1,l_size
233 : : ! cc(2,l) is g_prime(0)
234 [ + + - + : 112898 : if (optgr0==1.or.optgr1==1.or.optgr2==1) then
- - ]
235 [ + + ]: 112898 : 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 [ + + ]: 112883 : else if (shape_type==1.or.shape_type==2) then
239 : 6694 : 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 : 112898 : cc(1,ll)=cc(1,ll)*pawtab%gnorm(ll)
245 : : end if
246 : : ! cc(2,l) is g_prime(0)
247 [ + + - + ]: 112898 : if (optgr1==1.or.optgr2==1) then
248 [ + + ]: 15992 : 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 [ + + - + ]: 15977 : else if (shape_type==1.and.lambda==1) then
252 : 0 : cc(2,ll)=-one/sigma
253 : : else
254 : 15977 : cc(2,ll)=zero
255 : : end if
256 : 15992 : cc(2,ll)=cc(2,ll)*pawtab%gnorm(ll)
257 : : end if
258 : : ! cc(3,l) is g_prime_prime(0)
259 [ + + ]: 149912 : if (optgr2==1) then
260 [ - + ]: 219 : 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 [ - + ]: 219 : 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 [ + + ]: 219 : else if (shape_type==2) then
268 : 202 : 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 : 219 : 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 : 43669 : normchoice=1 ; option=max(optgr0,2*optgr1,3*optgr2)
281 [ + - ]: 43669 : if(compute_gr0) then
282 [ - + + - : 174676 : LIBPAW_ALLOCATE(ylmr,(l_size**2,nfgd))
- + + - +
- + - - +
- + ]
283 : : end if
284 [ + + ]: 43669 : if(compute_gr1.and.(.not.compute_gr2)) then
285 [ - + + - : 35496 : LIBPAW_ALLOCATE(ylmrgr,(3,l_size**2,nfgd))
- + + - +
- + - - +
- + ]
286 : : end if
287 [ + + ]: 43669 : if(compute_gr2) then
288 [ - + + - : 404 : LIBPAW_ALLOCATE(ylmrgr,(9,l_size**2,nfgd))
- + + - +
- + - - +
- + - + ]
289 : : end if
290 [ + + + - ]: 43669 : if (compute_gr0.and.(.not.compute_gr1).and.(.not.compute_gr2)) then
291 [ - + - - : 34694 : call initylmr(l_size,normchoice,nfgd,rnrm,option,rfgd,ylmr)
- - - - -
+ ]
292 : : else
293 [ - + - - : 8975 : 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 [ + - ]: 43669 : if (compute_gr0) then
300 [ - + + - : 174676 : LIBPAW_BOUND2_ALLOCATE(gfact,BOUNDS(1,nfgd),BOUNDS(0,l_size-1))
- + + - +
- + - - +
- + ]
301 [ + + + + ]: 135389559 : gfact(:,:)=zero
302 : : end if
303 [ + + ]: 43669 : if (compute_gr1) then
304 [ - + + - : 35900 : LIBPAW_BOUND2_ALLOCATE(dgfact,BOUNDS(1,nfgd),BOUNDS(0,l_size-1))
- + + - +
- + - - +
- + ]
305 [ + + + + ]: 39287671 : dgfact(:,:)=zero
306 : : end if
307 [ + + ]: 43669 : if (compute_gr2) then
308 [ - + + - : 404 : LIBPAW_BOUND2_ALLOCATE(d2gfact,BOUNDS(1,nfgd),BOUNDS(0,l_size-1))
- + + - +
- + - - +
- + ]
309 [ + + + + ]: 784050 : d2gfact(:,:)=zero
310 : : end if
311 [ + + ]: 43669 : if(compute_gr1) then
312 [ - + + - : 17950 : LIBPAW_ALLOCATE(rnrm_inv,(nfgd))
- + - + ]
313 [ + + ]: 10973562 : do ic=1,nfgd
314 [ + + ]: 10973562 : if (ic/=izero) rnrm_inv(ic)=one/rnrm(ic)
315 : : end do
316 [ + + ]: 8975 : if (izero>0) rnrm_inv(izero)=zero
317 : : end if
318 : :
319 : : !----- type -1 -----
320 [ + + ]: 43669 : 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 [ + + ]: 43666 : else if (shape_type==1.or.shape_type==2) then
351 : : ! FIRST COMPUTE FACTORS FOR l=0
352 [ + + + + : 7876 : if (optgr0==1.and.optgr1==0.and.optgr2==0) then
+ - ]
353 [ + + ]: 3393 : 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 [ + + ]: 5479031 : do ic=1,nfgd
364 : 5475640 : arg=rnrm(ic)
365 [ + + ]: 5479031 : if (arg<toldev) then
366 : 822 : gfact(ic,0)=shapefunc2_0(arg)
367 [ + - ]: 5474818 : else if (arg<=rcut) then
368 : 5474818 : gfact(ic,0)=shapefunc2(arg)
369 : : end if
370 : : end do
371 : : end if
372 [ + - + + ]: 4483 : else if (optgr1==1.and.optgr2==0) then
373 [ + + ]: 4405 : 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 [ + + ]: 6538649 : do ic=1,nfgd
390 : 6534246 : arg=rnrm(ic)
391 [ + + ]: 6538649 : if (arg<toldev) then
392 : 920 : gfact(ic,0)=shapefunc2_0(arg)
393 : 920 : dgfact(ic,0)=dshpfunc2_ovr_0(arg)
394 [ + - ]: 6533326 : else if (arg<=rcut) then
395 : 6533326 : gfact(ic,0)=shapefunc2(arg)
396 : 6533326 : dgfact(ic,0)=dshpfunc2(arg)*rnrm_inv(ic)
397 : : end if
398 : : end do
399 : : end if
400 [ + - ]: 78 : else if (optgr2==1) then
401 [ - + ]: 78 : 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 [ + + ]: 172904 : do ic=1,nfgd
431 : 172826 : arg=rnrm(ic)
432 [ + + ]: 172904 : if (arg<toldev) then
433 : 60 : gfact(ic,0)=shapefunc2_0(arg)
434 : 60 : dgfact(ic,0)=dshpfunc2_ovr_0(arg)
435 : 60 : d2gfact(ic,0)=d2shpfunc2_ovr2_0(arg)
436 [ + - ]: 172766 : else if (arg<=rcut) then
437 : 172766 : gfact(ic,0)=shapefunc2(arg)
438 : 172766 : dgfact(ic,0)=dshpfunc2(arg)*rnrm_inv(ic)
439 : 172766 : 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 [ + - ]: 7876 : if (compute_gr0) then
447 [ + + ]: 7876 : if (l_size>1) then
448 [ + + ]: 19508 : do ll=1,l_size-1
449 [ + + ]: 33804756 : do ic=1,nfgd
450 : 33798968 : 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 [ + + ]: 7876 : if (compute_gr1) then
456 [ + + ]: 4483 : if (l_size>1) then
457 [ + + ]: 6312649 : do ic=1,nfgd
458 : 6312649 : dgfact(ic,1)=pawtab%gnorm(2)*(gfact(ic,0)*rnrm_inv(ic)+dgfact(ic,0)*rnrm(ic))
459 : : end do
460 : : end if
461 [ + + ]: 4483 : if (l_size>2) then
462 [ + + ]: 6312649 : do ic=1,nfgd
463 : 6312649 : dgfact(ic,2)=pawtab%gnorm(3)*(two*gfact(ic,0)+dgfact(ic,0)*rnrm(ic)**2)
464 : : end do
465 : : end if
466 [ + + ]: 4483 : if (l_size>3) then
467 [ + + ]: 1710 : do ll=3,l_size-1
468 [ + + ]: 5864120 : do ic=1,nfgd
469 : : dgfact(ic,ll)=pawtab%gnorm(ll+1) &
470 : 5863570 : & *(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 [ + + ]: 7876 : if (compute_gr2) then
476 [ + - ]: 78 : if (l_size>1) then
477 [ + + ]: 172904 : do ic=1,nfgd
478 : : d2gfact(ic,1)=pawtab%gnorm(2) &
479 : 172904 : & *(-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 [ + - ]: 78 : if (l_size>2) then
483 [ + + ]: 172904 : do ic=1,nfgd
484 : 172904 : d2gfact(ic,2)=pawtab%gnorm(3)*(four*dgfact(ic,0)+d2gfact(ic,0)*rnrm(ic)**2)
485 : : end do
486 : : end if
487 [ + + ]: 78 : if (l_size>3) then
488 [ + + ]: 92822 : do ic=1,nfgd
489 : : d2gfact(ic,3)=pawtab%gnorm(4) &
490 : 92822 : & *(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 [ + + ]: 78 : if (l_size>4) then
494 [ + + ]: 92822 : do ic=1,nfgd
495 : : d2gfact(ic,4)=pawtab%gnorm(5) &
496 : 92822 : & *(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 [ - + ]: 78 : 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 [ + - + + ]: 12201224 : if (compute_gr0) gfact(:,0)=gfact(:,0)*pawtab%gnorm(1)
511 [ + + + + ]: 6720266 : if (compute_gr1) dgfact(:,0)=dgfact(:,0)*pawtab%gnorm(1)
512 [ + + + + ]: 180702 : 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 [ + + ]: 43669 : if (optgr0==1) then
564 : :
565 [ + + ]: 166779 : do ll=0,l_size-1
566 [ + + ]: 544387 : do ilm=ll**2+1,min((ll+1)**2,lm_size)
567 [ + + ]: 399457960 : do ic=1,nfgd
568 : 399333566 : 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 [ + + ]: 42385 : if (izero>0) then
575 [ + + ]: 372389 : gylm(izero,1:lm_size)=zero
576 [ + - ]: 36283 : 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 [ + + ]: 43669 : if(optgr1==1) then
584 : :
585 [ + + ]: 34351 : do ll=0,l_size-1
586 [ + + ]: 115157 : do ilm=ll**2+1,min((ll+1)**2,lm_size)
587 [ + + ]: 158388466 : do ic=1,nfgd
588 : : gylmgr(1:3,ic,ilm)=gfact(ic,ll)*ylmrgr(1:3,ilm,ic)&
589 [ + + ]: 633209942 : & +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 [ + + ]: 8975 : if (izero>0) then
596 [ + + + + ]: 218085 : gylmgr(1:3,izero,1:lm_size)=zero
597 [ + - ]: 5013 : if (lm_size>=1) then
598 : 5013 : arg=cc(2,1)/sqrt(four_pi)
599 [ + + ]: 20052 : gylmgr(1:3,izero,1)=arg
600 : : end if
601 [ + + ]: 5013 : if (lm_size>=2) then
602 : 4975 : arg=cc(1,2)*sqrt(three/four_pi)
603 : 4975 : gylmgr(2,izero,2)=arg
604 [ + - ]: 4975 : if (lm_size>=3) gylmgr(3,izero,3)=arg
605 [ + - ]: 4975 : 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 [ + + ]: 43669 : if(optgr2==1) then
614 : :
615 [ + + ]: 428 : do ll=0,l_size-1
616 [ + + ]: 1529 : do ilm=ll**2+1,min((ll+1)**2,lm_size)
617 [ + + ]: 3317354 : 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 : 3315926 : & +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 : 3315926 : & +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 : 3315926 : & +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 : 3315926 : & +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 : 3315926 : & +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 : 3317027 : & +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 [ + + ]: 101 : if (izero>0) then
642 [ + + + + ]: 5504 : gylmgr2(1:6,izero,1:lm_size)=zero
643 [ + - ]: 65 : if (lm_size>=1) then
644 : 65 : arg=cc(3,1)/sqrt(four_pi)
645 [ + + ]: 260 : gylmgr2(1:3,izero,1)=arg
646 : : end if
647 [ + - ]: 65 : if (lm_size>=2) then
648 : 65 : arg=cc(2,2)*sqrt(three/four_pi)
649 : 65 : gylmgr2(2,izero,2)=two*arg
650 : 65 : gylmgr2(4,izero,2)= arg
651 [ + - ]: 65 : if (lm_size>=3) then
652 : 65 : gylmgr2(1,izero,3)=two*arg
653 : 65 : gylmgr2(3,izero,3)=two*arg
654 : : end if
655 [ + - ]: 65 : if (lm_size>=4) then
656 : 65 : gylmgr2(5,izero,4)=arg
657 : 65 : gylmgr2(6,izero,4)=arg
658 : : end if
659 : : end if
660 [ + - ]: 65 : if (lm_size>=5) then
661 : 65 : arg=cc(1,3)*sqrt(15._dp/four_pi)
662 : 65 : gylmgr2(6,izero,5)=arg
663 [ + - ]: 65 : if (lm_size>=6) gylmgr2(4,izero,6)=arg
664 [ + - ]: 65 : if (lm_size>=7) then
665 : 65 : gylmgr2(1,izero,7)= -arg/sqrt3
666 : 65 : gylmgr2(2,izero,7)= -arg/sqrt3
667 : 65 : gylmgr2(3,izero,7)=two*arg/sqrt3
668 : : end if
669 [ + - ]: 65 : if (lm_size>=8) gylmgr2(5,izero,8)=arg
670 [ + - ]: 65 : if (lm_size>=9) then
671 : 65 : gylmgr2(1,izero,9)= arg
672 : 65 : gylmgr2(2,izero,9)=-arg
673 : : end if
674 : : end if
675 : : end if
676 : :
677 : : end if
678 : :
679 : : !Memory deallocation
680 : : !==========================================================
681 [ - + ]: 43669 : LIBPAW_DEALLOCATE(rnrm)
682 [ + + ]: 43669 : if (allocated(cc)) then
683 : 37014 : LIBPAW_DEALLOCATE(cc)
684 : : end if
685 [ + - ]: 43669 : if (compute_gr0) then
686 [ - + ]: 43669 : LIBPAW_DEALLOCATE(gfact)
687 : : end if
688 [ + + ]: 43669 : if (compute_gr1) then
689 [ - + ]: 8975 : LIBPAW_DEALLOCATE(dgfact)
690 : : end if
691 [ + + ]: 43669 : if (compute_gr2) then
692 [ - + ]: 101 : LIBPAW_DEALLOCATE(d2gfact)
693 : : end if
694 [ + + ]: 43669 : if (compute_gr1) then
695 [ - + ]: 8975 : LIBPAW_DEALLOCATE(rnrm_inv)
696 : : end if
697 [ + + ]: 43669 : if (shape_type==3) then
698 [ - + ]: 35790 : LIBPAW_DEALLOCATE(alpha)
699 [ - + ]: 35790 : LIBPAW_DEALLOCATE(qq)
700 : : end if
701 [ + - ]: 43669 : if (compute_gr0) then
702 [ - + ]: 43669 : LIBPAW_DEALLOCATE(ylmr)
703 : : end if
704 [ + + ]: 43669 : if (compute_gr1) then
705 [ - + ]: 8975 : LIBPAW_DEALLOCATE(ylmrgr)
706 : : end if
707 [ + + - + : 43688 : 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 : 12180910 : function shapefunc2(arg)
747 : :
748 : : real(dp) :: shapefunc2
749 : : real(dp),intent(in) :: arg
750 : 12180910 : shapefunc2=(sin(pi_over_rshp*arg)/(pi_over_rshp*arg))**2
751 : 12180910 : end function shapefunc2
752 : : !!***
753 : : ! ------------------------------------------------
754 : : !!****f* m_paw_finegrid/shapefunc2_0
755 : : ! shapefunc2_0 is g(x) (sinc2) for small x
756 : 1802 : function shapefunc2_0(arg)
757 : :
758 : : real(dp) :: shapefunc2_0
759 : : real(dp),intent(in) :: arg
760 : 1802 : shapefunc2_0=one-(pi_over_rshp*arg)**2/three+two*(pi_over_rshp*arg)**4/45._dp
761 : 1802 : 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 : 6706092 : 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 : 6706092 : & *(pi_over_rshp*arg*cos(pi_over_rshp*arg)-sin(pi_over_rshp*arg))
813 : 6706092 : 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 : 980 : function dshpfunc2_ovr_0(arg)
819 : :
820 : : real(dp) :: dshpfunc2_ovr_0
821 : : real(dp),intent(in) :: arg
822 : 980 : dshpfunc2_ovr_0=-two*pi_over_rshp**2/3._dp+8._dp*pi_over_rshp**4*arg**2/45._dp
823 : 980 : 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 : 172766 : 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 : 172766 : & -four*pi_over_rshp*arg*cos(pi_over_rshp*arg)*sin(pi_over_rshp*arg))
900 : 172766 : 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 : 60 : 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 : 60 : & +41._dp/6300._dp*pi_over_rshp**8*arg**4
911 : 60 : 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 : 5033 : subroutine pawrfgd_fft(ifftsph,gmet,n1,n2,n3,nfgd,rcut,rfgd,rprimd,ucvol,xred, &
1065 : 5033 : & 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 : 5033 : integer,allocatable :: ifftsph_tmp(:)
1086 : 5033 : integer,pointer :: fft_distrib_(:),fft_index_(:)
1087 : 5033 : real(dp),allocatable :: rfgd_tmp(:,:)
1088 : :
1089 : : ! *************************************************************************
1090 : :
1091 : : !Define a "box" around the atom
1092 : 5033 : r2cut=1.0000001_dp*rcut**2
1093 : 5033 : rr1=sqrt(r2cut*gmet(1,1))
1094 : 5033 : rr2=sqrt(r2cut*gmet(2,2))
1095 : 5033 : 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 : 5033 : dif=xred(1)-rr1
1100 : 5033 : n1a=int(dif*n1)+1
1101 : 3372 : if (dif<0) n1a = n1a-1
1102 : 5033 : dif=xred(2)-rr2
1103 : 5033 : n2a=int(dif*n2)+1
1104 [ + + ]: 5033 : if (dif<0) n2a = n2a-1
1105 : 5033 : dif=xred(3)-rr3
1106 : 5033 : n3a=int(dif*n3)+1
1107 [ + + ]: 5033 : 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 : 5033 : dif=xred(1)+rr1
1112 : 5033 : n1b=int(dif*n1)+1
1113 [ + + ]: 5033 : if (dif>0) n1b = n1b+1
1114 : 5033 : dif=xred(2)+rr2
1115 : 5033 : n2b=int(dif*n2)+1
1116 [ + + ]: 5033 : if (dif>0) n2b = n2b+1
1117 : 5033 : dif=xred(3)+rr3
1118 : 5033 : n3b=int(dif*n3)+1
1119 [ + + ]: 5033 : if (dif>0) n3b = n3b+1
1120 : :
1121 : : !Get the distrib associated with this fft_grid
1122 [ + - + - ]: 5033 : if (present(fft_distrib).and.present(fft_index).and.present(me_fft)) then
1123 : 5033 : 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 : 5033 : ncmax=1+int(1.5_dp*(n1*n2*n3)*four_pi/(three*ucvol)*rcut**3)
1133 [ - + + - : 15099 : LIBPAW_ALLOCATE(ifftsph_tmp,(ncmax))
+ - + - -
+ - + ]
1134 [ - + + - : 15099 : LIBPAW_ALLOCATE(rfgd_tmp,(3,ncmax))
+ - + - -
+ - + ]
1135 : :
1136 : : !Set number of points to zero
1137 : 5033 : nfgd=0
1138 : :
1139 : : !Loop over FFT points
1140 [ + + ]: 87992 : do i3=n3a,n3b
1141 : 82959 : iz=modulo(i3,n3)
1142 [ + - - + ]: 82959 : 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 [ + + ]: 87992 : if (fft_distrib_(iz+1)==me_fft_) then
1147 : 78181 : izloc=fft_index_(iz+1) - 1
1148 [ + - - + ]: 78181 : 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 : 78181 : difz=dble(i3)/dble(n3)-xred(3)
1153 [ + + ]: 1578575 : do i2=n2a,n2b
1154 : 1500394 : iy=modulo(i2,n2)
1155 [ + - - + ]: 1500394 : 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 : 1500394 : dify=dble(i2)/dble(n2)-xred(2)
1160 [ + + ]: 35196212 : do i1=n1a,n1b
1161 : 33617637 : ix=modulo(i1,n1)
1162 [ + - - + ]: 33617637 : 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 : 33617637 : difx=dble(i1)/dble(n1)-xred(1)
1167 : :
1168 : : ! Compute r-R
1169 : 33617637 : rx=difx*rprimd(1,1)+dify*rprimd(1,2)+difz*rprimd(1,3)
1170 : 33617637 : ry=difx*rprimd(2,1)+dify*rprimd(2,2)+difz*rprimd(2,3)
1171 : 33617637 : rz=difx*rprimd(3,1)+dify*rprimd(3,2)+difz*rprimd(3,3)
1172 : 33617637 : r2=rx**2+ry**2+rz**2
1173 : :
1174 : : ! Select matching points
1175 [ + + ]: 35118031 : if (r2 <= r2cut) then
1176 : 10260059 : ifft_local=1+ix+n1*(iy+n2*izloc)
1177 [ + - - + ]: 10260059 : 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 [ + - ]: 10260059 : if (ifft_local>0) then
1182 : 10260059 : nfgd=nfgd+1
1183 [ - + ]: 10260059 : 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 : 10260059 : rfgd_tmp(1,nfgd)=rx
1188 : 10260059 : rfgd_tmp(2,nfgd)=ry
1189 : 10260059 : rfgd_tmp(3,nfgd)=rz
1190 : 10260059 : 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 [ - + ]: 5033 : if (allocated(ifftsph)) then
1202 : 0 : LIBPAW_DEALLOCATE(ifftsph)
1203 : : end if
1204 [ - + ]: 5033 : if (allocated(rfgd)) then
1205 : 0 : LIBPAW_DEALLOCATE(rfgd)
1206 : : end if
1207 [ + + + - : 15099 : LIBPAW_ALLOCATE(ifftsph,(nfgd))
+ - + + -
+ - + -
+ ]
1208 [ + + + - : 15099 : LIBPAW_ALLOCATE(rfgd,(3,nfgd))
+ - + + -
+ - + -
+ ]
1209 [ + + ]: 10265092 : ifftsph(1:nfgd)=ifftsph_tmp(1:nfgd)
1210 [ + + + + ]: 41045269 : rfgd(1:3,1:nfgd)=rfgd_tmp(1:3,1:nfgd)
1211 : :
1212 : : !Release temporary memory
1213 [ - + ]: 5033 : LIBPAW_DEALLOCATE(ifftsph_tmp)
1214 [ - + ]: 5033 : LIBPAW_DEALLOCATE(rfgd_tmp)
1215 [ - + ]: 5033 : 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 [ + + ]: 5033 : 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 : : !!***
|
