TABLE OF CONTENTS

ABINIT/m_paw_finegrid [ Modules ]

[ Top ] [ Modules ]

NAME

  m_paw_finegrid

FUNCTION

  This module contains a set of routines to compute various quantities
  on the fine grid around a given atom.

COPYRIGHT

 Copyright (C) 2013-2024 ABINIT group (MT,FJ)
 This file is distributed under the terms of the
 GNU General Public License, see ~abinit/COPYING
 or http://www.gnu.org/copyleft/gpl.txt .

NOTES

  FOR DEVELOPPERS: in order to preserve the portability of libPAW library,
  please consult ~abinit/src/??_libpaw/libpaw-coding-rules.txt

SOURCE

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(:,:)

m_paw_finegrid/d2shpfunc1 [ Functions ]

[ Top ] [ m_paw_finegrid ] [ Functions ]

m_paw_finegrid/d2shpfunc1_ovr2_0 [ Functions ]

[ Top ] [ m_paw_finegrid ] [ Functions ]

m_paw_finegrid/d2shpfunc1_ovr2_0_2 [ Functions ]

[ Top ] [ m_paw_finegrid ] [ Functions ]

m_paw_finegrid/d2shpfunc1_ovr2_0_3 [ Functions ]

[ Top ] [ m_paw_finegrid ] [ Functions ]

m_paw_finegrid/d2shpfunc1_ovr2_0_4 [ Functions ]

[ Top ] [ m_paw_finegrid ] [ Functions ]

m_paw_finegrid/d2shpfunc2 [ Functions ]

[ Top ] [ m_paw_finegrid ] [ Functions ]

m_paw_finegrid/d2shpfunc2_ovr2_0 [ Functions ]

[ Top ] [ m_paw_finegrid ] [ Functions ]

m_paw_finegrid/d2shpfunc3 [ Functions ]

[ Top ] [ m_paw_finegrid ] [ Functions ]

m_paw_finegrid/dshpfunc1 [ Functions ]

[ Top ] [ m_paw_finegrid ] [ Functions ]

m_paw_finegrid/dshpfunc1_ovr_0 [ Functions ]

[ Top ] [ m_paw_finegrid ] [ Functions ]

m_paw_finegrid/dshpfunc1_ovr_0_2 [ Functions ]

[ Top ] [ m_paw_finegrid ] [ Functions ]

m_paw_finegrid/dshpfunc2 [ Functions ]

[ Top ] [ m_paw_finegrid ] [ Functions ]

m_paw_finegrid/dshpfunc2_ovr_0 [ Functions ]

[ Top ] [ m_paw_finegrid ] [ Functions ]

m_paw_finegrid/dshpfunc3 [ Functions ]

[ Top ] [ m_paw_finegrid ] [ Functions ]

m_paw_finegrid/my_ext_buffers [ Functions ]

[ Top ] [ m_paw_finegrid ] [ Functions ]

m_paw_finegrid/my_ind_positions [ Functions ]

[ Top ] [ m_paw_finegrid ] [ Functions ]

m_paw_finegrid/pawexpiqr [ Functions ]

[ Top ] [ m_paw_finegrid ] [ Functions ]

NAME

 pawexpiqr

FUNCTION

 Compute exp(i.q.r) for each point of the (fine) rectangular grid
 around a given atomic site. R is the position of the atom.
 Used for the determination of phonons at non-zero q wavevector.

INPUTS

  gprimd(3,3)= dimensional primitive translations for reciprocal space
  nfgd= number of (fine grid) FFT points in the paw sphere around current atom
  qphon(3)= wavevector of the phonon
  rfgd(3,nfgd)= coordinates of r-R on the fine grid around current atom
  xred(3)= reduced atomic coordinates

OUTPUT

  expiqr(2,nfgd)= exp(i.q.r) around the current atom
                                 Not allocated if q=0 !

SOURCE

1377 subroutine pawexpiqr(expiqr,gprimd,nfgd,qphon,rfgd,xred)
1378 
1379 !Arguments ---------------------------------------------
1380 !scalars
1381  integer,intent(in) :: nfgd
1382 !arrays
1383  real(dp),intent(in) :: gprimd(3,3),qphon(3),xred(3)
1384  real(dp),intent(in) :: rfgd(:,:)
1385  real(dp),intent(out) :: expiqr(2,nfgd)
1386 
1387 !Local variables ------------------------------
1388 !scalars
1389  integer :: ic
1390  logical :: qne0
1391  real(dp) :: phase,phase_xred,qx,qy,qz
1392  character(len=500) :: msg
1393 !arrays
1394 
1395 ! *************************************************************************
1396 
1397  if (size(rfgd)/=3*nfgd) then
1398    msg='rfgd array must be allocated!'
1399    LIBPAW_BUG(msg)
1400  end if
1401 
1402  qne0=(qphon(1)**2+qphon(2)**2+qphon(3)**2>=1.d-15)
1403 
1404 !Compute q in cartesian coordinates
1405  if (qne0) then
1406    qx=gprimd(1,1)*qphon(1)+gprimd(1,2)*qphon(2)+gprimd(1,3)*qphon(3)
1407    qy=gprimd(2,1)*qphon(1)+gprimd(2,2)*qphon(2)+gprimd(2,3)*qphon(3)
1408    qz=gprimd(3,1)*qphon(1)+gprimd(3,2)*qphon(2)+gprimd(3,3)*qphon(3)
1409    phase_xred=two_pi*(qphon(1)*xred(1)+qphon(2)*xred(2)+qphon(3)*xred(3))
1410  end if
1411 
1412 !Compute exp(i.q.r)
1413  if (qne0) then
1414    do ic=1,nfgd
1415      phase=two_pi*(qx*rfgd(1,ic)+qy*rfgd(2,ic)+qz*rfgd(3,ic)) + phase_xred
1416      expiqr(1,ic)=cos(phase)
1417      expiqr(2,ic)=sin(phase)
1418    end do
1419  end if
1420 
1421 end subroutine pawexpiqr

m_paw_finegrid/pawgylm [ Functions ]

[ Top ] [ m_paw_finegrid ] [ Functions ]

NAME

 pawgylm

FUNCTION

 Compute g_l(r-R)*Y_lm(r-R) (and derivatives) on the fine (rectangular) grid
 around one atom (g_l=radial shape function).
 R is the position of the atom

INPUTS

  lm_size=number of lm components to be calculated
  nfgd= number of (fine grid) FFT points in the paw sphere around current atom
  optgr0= 1 if g_l(r-R)*Y_lm(r-R) are computed
  optgr1= 1 if first derivatives of g_l(r-R)*Y_lm(r-R) are computed
  optgr2= 1 if second derivatives of g_l(r-R)*Y_lm(r-R) are computed
  pawtab <type(pawtab_type)>=paw tabulated starting data for current atom
  rfgd(3,nfgd)= coordinates of r-R on the fine rect. grid around current atom

OUTPUT

  if (optgr0==1)
    gylm(nfgd,lm_size)= g_l(r-R)*Y_lm(r-R) around current atom
  if (optgr1==1)
    gylmgr(3,nfgd,lm_size)= derivatives of g_l(r-R)*Y_lm(r-R) wrt cart. coordinates
  if (optgr2==1)
    gylmgr2(6,nfgd,lm_size)= second derivatives of g_l(r-R)*Y_lm(r-R) wrt cart. coordinates

SOURCE

 86 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  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  real(dp),allocatable :: cc(:,:),d2gfact(:,:),d2shpfuncnum(:,:),dgfact(:,:)
113  real(dp),allocatable :: dshpfuncnum(:,:),gfact(:,:)
114  real(dp),allocatable :: rnrm(:),rnrm_inv(:),rnrm_sort(:)
115  real(dp),allocatable :: shpfuncnum(:,:),work(:),ylmr(:,:),ylmrgr(:,:,:)
116 
117 ! *************************************************************************
118 
119  if (optgr0==0.and.optgr1==0.and.optgr2==0) return
120  if (nfgd==0) return
121 
122 !Compatibility test
123 !==========================================================
124  if (size(rfgd)/=3*nfgd) then
125    msg='rfgd array must be allocated at rfgd(3,nfgd)!'
126    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  if (pawtab%shape_type==1.and.pawtab%shape_lambda<2) then
133    msg='Exponent lambda of gaussian shape function must be > 1!'
134    LIBPAW_ERROR(msg)
135  end if
136 
137 !Initializations
138 !==========================================================
139 !Options for computation
140  compute_gr0=(optgr0==1.or.optgr1==1.or.optgr2==1)
141  compute_gr1=(optgr1==1.or.optgr2==1)
142  compute_gr2=(optgr2==1)
143  l_size=pawtab%lcut_size
144 
145 !Norms of vectors around the atom
146  LIBPAW_ALLOCATE(rnrm,(nfgd))
147  izero=-1
148  do ic=1,nfgd
149    rnrm(ic)=sqrt(rfgd(1,ic)**2+rfgd(2,ic)**2+rfgd(3,ic)**2)
150    if (rnrm(ic)<=tol10) izero=ic  ! Has to be consistent with initylmr !!
151  end do
152 
153 !Initializations
154  if (optgr0==1) gylm=zero
155  if (optgr1==1) gylmgr=zero
156  if (optgr2==1) gylmgr2=zero
157 
158 !Some definitions concerning shape function g_l(r)
159  shape_type=pawtab%shape_type
160  sigma=pawtab%shape_sigma;lambda=pawtab%shape_lambda
161  pi_over_rshp=pi/pawtab%rshp
162  rcut=tol12+pawtab%rshp
163  if (shape_type==3) then
164    LIBPAW_ALLOCATE(alpha,(2,l_size))
165    LIBPAW_ALLOCATE(qq,(2,l_size))
166    do ll=1,l_size
167      alpha(1:2,ll)=pawtab%shape_alpha(1:2,ll)
168      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  if (shape_type==-1) then
174    LIBPAW_ALLOCATE(isort,(nfgd))
175    LIBPAW_ALLOCATE(rnrm_sort,(nfgd))
176    do ic=1,nfgd
177      isort(ic)=ic
178    end do
179    rnrm_sort(1:nfgd)=rnrm(1:nfgd)
180    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    LIBPAW_ALLOCATE(work,(nfgd))
186    if (compute_gr0) then
187      LIBPAW_ALLOCATE(shpfuncnum,(nfgd,l_size))
188      do ll=1,l_size
189        call paw_splint(pawtab%mesh_size,pawtab%rad_for_spline,pawtab%shapefunc(:,ll),&
190 &       pawtab%dshpfunc(:,ll,2),nfgd,rnrm_sort,work)
191        do ic=1,nfgd
192          shpfuncnum(isort(ic),ll)=work(ic)
193        end do
194      end do
195    end if
196    if(compute_gr1) then
197      LIBPAW_ALLOCATE(dshpfuncnum,(nfgd,l_size))
198      do ll=1,l_size
199        call paw_splint(pawtab%mesh_size,pawtab%rad_for_spline,pawtab%dshpfunc(:,ll,1),&
200 &       pawtab%dshpfunc(:,ll,3),nfgd,rnrm_sort,work)
201        do ic=1,nfgd
202          dshpfuncnum(isort(ic),ll)=work(ic)
203        end do
204      end do
205    end if
206    if(compute_gr2) then
207      LIBPAW_ALLOCATE(d2shpfuncnum,(nfgd,l_size))
208      do ll=1,l_size
209        call paw_splint(pawtab%mesh_size,pawtab%rad_for_spline,pawtab%dshpfunc(:,ll,2),&
210 &       pawtab%dshpfunc(:,ll,4),nfgd,rnrm_sort,work)
211        do ic=1,nfgd
212          d2shpfuncnum(isort(ic),ll)=work(ic)
213        end do
214      end do
215    end if
216    LIBPAW_DEALLOCATE(work)
217  end if
218 
219  if (shape_type==-1)  then
220    LIBPAW_DEALLOCATE(isort)
221    LIBPAW_DEALLOCATE(rnrm_sort)
222  end if
223 
224 !If needed, compute limits at r=0 of shape function and derivatives
225  if (izero>0) then
226    LIBPAW_ALLOCATE(cc,(3,min(l_size,3)))
227    cc=zero
228    if (shape_type==-1) then
229      splfact=(pawtab%rad_for_spline(4)-pawtab%rad_for_spline(1))&
230 &     /(pawtab%rad_for_spline(3)-pawtab%rad_for_spline(2))
231    end if
232    do ll=1,min(l_size,3)
233 !    cc(2,l) is g_prime(0)
234      if (optgr0==1.or.optgr1==1.or.optgr2==1) then
235        if (shape_type==-1) then
236          ss(1:3)=pawtab%shapefunc(2:4,ll)/pawtab%rad_for_spline(2:4)**(ll-1)
237          cc(1,ll)=ss(3)+(ss(1)-ss(2))*splfact
238        else if (shape_type==1.or.shape_type==2) then
239          cc(1,ll)=one
240        else if (shape_type==3) then
241          cc(1,ll)=(alpha(1,ll)*qq(1,ll)**(ll-1) &
242 &         +alpha(2,ll)*qq(2,ll)**(ll-1))/ffact(ll)
243        end if
244        cc(1,ll)=cc(1,ll)*pawtab%gnorm(ll)
245      end if
246 !    cc(2,l) is g_prime(0)
247      if (optgr1==1.or.optgr2==1) then
248        if (shape_type==-1) then
249          ss(1:3)=(ss(1:3)-cc(1,ll))/pawtab%rad_for_spline(2:4)
250          cc(2,ll)=ss(3)+(ss(1)-ss(2))*splfact
251        else if (shape_type==1.and.lambda==1) then
252          cc(2,ll)=-one/sigma
253        else
254          cc(2,ll)=zero
255        end if
256        cc(2,ll)=cc(2,ll)*pawtab%gnorm(ll)
257      end if
258 !    cc(3,l) is g_prime_prime(0)
259      if (optgr2==1) then
260        if (shape_type==-1) then
261          ss(1:3)=(ss(1:3)-cc(2,ll))/pawtab%rad_for_spline(2:4)
262          cc(3,ll)=two*(ss(3)+(ss(1)-ss(2))*splfact)
263        else if (shape_type==1) then
264          if (lambda==1) cc(3,ll)=one/sigma**2
265          if (lambda==2) cc(3,ll)=-two/sigma**2
266          if (lambda >2) cc(3,ll)=zero
267        else if (shape_type==2) then
268          cc(3,ll)=-(two/three)*pi_over_rshp**2
269        else if (shape_type==3) then
270          cc(3,ll)=-(alpha(1,ll)*qq(1,ll)**(ll+1) &
271 &         +alpha(2,ll)*qq(2,ll)**(ll+1))/ffact(ll+1)
272        end if
273        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  normchoice=1 ; option=max(optgr0,2*optgr1,3*optgr2)
281  if(compute_gr0)  then
282    LIBPAW_ALLOCATE(ylmr,(l_size**2,nfgd))
283  end if
284  if(compute_gr1.and.(.not.compute_gr2))  then
285    LIBPAW_ALLOCATE(ylmrgr,(3,l_size**2,nfgd))
286  end if
287  if(compute_gr2)  then
288    LIBPAW_ALLOCATE(ylmrgr,(9,l_size**2,nfgd))
289  end if
290  if (compute_gr0.and.(.not.compute_gr1).and.(.not.compute_gr2)) then
291    call initylmr(l_size,normchoice,nfgd,rnrm,option,rfgd,ylmr)
292  else
293    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  if (compute_gr0)  then
300    LIBPAW_BOUND2_ALLOCATE(gfact,BOUNDS(1,nfgd),BOUNDS(0,l_size-1))
301    gfact(:,:)=zero
302  end if
303  if (compute_gr1)  then
304    LIBPAW_BOUND2_ALLOCATE(dgfact,BOUNDS(1,nfgd),BOUNDS(0,l_size-1))
305    dgfact(:,:)=zero
306  end if
307  if (compute_gr2)  then
308    LIBPAW_BOUND2_ALLOCATE(d2gfact,BOUNDS(1,nfgd),BOUNDS(0,l_size-1))
309    d2gfact(:,:)=zero
310  end if
311  if(compute_gr1) then
312    LIBPAW_ALLOCATE(rnrm_inv,(nfgd))
313    do ic=1,nfgd
314      if (ic/=izero) rnrm_inv(ic)=one/rnrm(ic)
315    end do
316    if (izero>0) rnrm_inv(izero)=zero
317  end if
318 
319 !----- type -1 -----
320  if (shape_type==-1) then
321    if (compute_gr0) then
322      do ll=0,l_size-1
323        do ic=1,nfgd
324          if (rnrm(ic)<=rcut) then
325            gfact(ic,ll)=shpfuncnum(ic,ll+1)
326          end if
327        end do
328      end do
329    end if
330    if (compute_gr1) then
331      do ll=0,l_size-1
332        do ic=1,nfgd
333          if (rnrm(ic)<=rcut) then
334            dgfact(ic,ll)=dshpfuncnum(ic,ll+1)*rnrm_inv(ic)
335          end if
336        end do
337      end do
338    end if
339    if(compute_gr2) then
340      do ll=0,l_size-1
341        do ic=1,nfgd
342          if (rnrm(ic)<=rcut) then
343            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  else if (shape_type==1.or.shape_type==2) then
351 !  FIRST COMPUTE FACTORS FOR l=0
352    if (optgr0==1.and.optgr1==0.and.optgr2==0) then
353      if (shape_type==1) then
354        do ic=1,nfgd
355          arg=rnrm(ic)
356          if (arg<toldev) then
357            gfact(ic,0)=shapefunc1_0(arg)
358          else if (arg<=rcut) then
359            gfact(ic,0)=shapefunc1(arg)
360          end if
361        end do
362      else ! shape_type==2
363        do ic=1,nfgd
364          arg=rnrm(ic)
365          if (arg<toldev) then
366            gfact(ic,0)=shapefunc2_0(arg)
367          else if (arg<=rcut) then
368            gfact(ic,0)=shapefunc2(arg)
369          end if
370        end do
371      end if
372    else if (optgr1==1.and.optgr2==0) then
373      if (shape_type==1) then
374        do ic=1,nfgd
375          arg=rnrm(ic)
376          if (arg<toldev) then
377            gfact(ic,0)=shapefunc1_0(arg)
378            if (lambda==2) then
379              dgfact(ic,0)=dshpfunc1_ovr_0_2(arg)
380            else ! lambda>2
381              dgfact(ic,0)=dshpfunc1_ovr_0(arg)
382            end if
383          else if (arg<=rcut) then
384            gfact(ic,0)=shapefunc1(arg)
385            dgfact(ic,0)=dshpfunc1(arg)*rnrm_inv(ic)
386          end if
387        end do
388      else ! shape_type==2
389        do ic=1,nfgd
390          arg=rnrm(ic)
391          if (arg<toldev) then
392            gfact(ic,0)=shapefunc2_0(arg)
393            dgfact(ic,0)=dshpfunc2_ovr_0(arg)
394          else if (arg<=rcut) then
395            gfact(ic,0)=shapefunc2(arg)
396            dgfact(ic,0)=dshpfunc2(arg)*rnrm_inv(ic)
397          end if
398        end do
399      end if
400    else if (optgr2==1) then
401      if (shape_type==1) then
402        do ic=1,nfgd
403          arg=rnrm(ic)
404          if (arg<toldev) then
405            gfact(ic,0)=shapefunc1_0(arg)
406            if (lambda==2) then
407              dgfact(ic,0)=dshpfunc1_ovr_0_2(arg)
408              d2gfact(ic,0)=d2shpfunc1_ovr2_0_2(arg)
409            else if (lambda==3) then
410              dgfact(ic,0)=dshpfunc1_ovr_0(arg)
411              if (ic/=izero) then
412                d2gfact(ic,0)=d2shpfunc1_ovr2_0_3(arg)
413              else
414                d2gfact(ic,0)=zero ! Diverging case
415              end if
416            else if (lambda==4) then
417              dgfact(ic,0)=dshpfunc1_ovr_0(arg)
418              d2gfact(ic,0)=d2shpfunc1_ovr2_0_4(arg)
419            else ! lambda>4
420              dgfact(ic,0)=dshpfunc1_ovr_0(arg)
421              d2gfact(ic,0)=d2shpfunc1_ovr2_0(arg)
422            end if
423          else if (arg<=rcut) then
424            gfact(ic,0)=shapefunc1(arg)
425            dgfact(ic,0)=dshpfunc1(arg)*rnrm_inv(ic)
426            d2gfact(ic,0)=(d2shpfunc1(arg)-dgfact(ic,0))*rnrm_inv(ic)**2
427          end if
428        end do
429      else ! shape_type==2
430        do ic=1,nfgd
431          arg=rnrm(ic)
432          if (arg<toldev) then
433            gfact(ic,0)=shapefunc2_0(arg)
434            dgfact(ic,0)=dshpfunc2_ovr_0(arg)
435            d2gfact(ic,0)=d2shpfunc2_ovr2_0(arg)
436          else if (arg<=rcut) then
437            gfact(ic,0)=shapefunc2(arg)
438            dgfact(ic,0)=dshpfunc2(arg)*rnrm_inv(ic)
439            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    if (compute_gr0) then
447      if (l_size>1) then
448        do ll=1,l_size-1
449          do ic=1,nfgd
450            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    if (compute_gr1) then
456      if (l_size>1) then
457        do ic=1,nfgd
458          dgfact(ic,1)=pawtab%gnorm(2)*(gfact(ic,0)*rnrm_inv(ic)+dgfact(ic,0)*rnrm(ic))
459        end do
460      end if
461      if (l_size>2) then
462        do ic=1,nfgd
463          dgfact(ic,2)=pawtab%gnorm(3)*(two*gfact(ic,0)+dgfact(ic,0)*rnrm(ic)**2)
464        end do
465      end if
466      if (l_size>3) then
467        do ll=3,l_size-1
468          do ic=1,nfgd
469            dgfact(ic,ll)=pawtab%gnorm(ll+1) &
470 &           *(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    if (compute_gr2) then
476      if (l_size>1) then
477        do ic=1,nfgd
478          d2gfact(ic,1)=pawtab%gnorm(2) &
479 &         *(-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      if (l_size>2) then
483        do ic=1,nfgd
484          d2gfact(ic,2)=pawtab%gnorm(3)*(four*dgfact(ic,0)+d2gfact(ic,0)*rnrm(ic)**2)
485        end do
486      end if
487      if (l_size>3) then
488        do ic=1,nfgd
489          d2gfact(ic,3)=pawtab%gnorm(4) &
490 &         *(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      if (l_size>4) then
494        do ic=1,nfgd
495          d2gfact(ic,4)=pawtab%gnorm(5) &
496 &         *(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      if (l_size>5) then
500        do ll=5,l_size-1
501          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 &            +d2gfact(ic,0)*rnrm(ic)**ll)
506          end do
507        end do
508      end if
509    end if
510    if (compute_gr0) gfact(:,0)=gfact(:,0)*pawtab%gnorm(1)
511    if (compute_gr1) dgfact(:,0)=dgfact(:,0)*pawtab%gnorm(1)
512    if (compute_gr2) d2gfact(:,0)=d2gfact(:,0)*pawtab%gnorm(1)
513 
514 !  ----- type 3 -----
515  else if (shape_type==3) then
516    if (optgr0==1.and.optgr1==0.and.optgr2==0) then
517      do ll=0,l_size-1
518        do ic=1,nfgd
519          arg=rnrm(ic)
520          if (arg<=rcut) then
521            call paw_jbessel(jbes1,jbesp1,jbespp1,ll,0,qq(1,1+ll)*arg)
522            call paw_jbessel(jbes2,jbesp2,jbespp2,ll,0,qq(2,1+ll)*arg)
523            gfact(ic,ll)=shapefunc3(jbes1,jbes2,ll)
524          end if
525        end do
526      end do
527    else if (optgr1==1.and.optgr2==0) then
528      do ll=0,l_size-1
529        do ic=1,nfgd
530          arg=rnrm(ic)
531          if (arg<=rcut) then
532            call paw_jbessel(jbes1,jbesp1,jbespp1,ll,1,qq(1,1+ll)*arg)
533            call paw_jbessel(jbes2,jbesp2,jbespp2,ll,1,qq(2,1+ll)*arg)
534            gfact(ic,ll)=shapefunc3(jbes1,jbes2,ll)
535            dgfact(ic,ll)=dshpfunc3(jbesp1,jbesp2,ll)*rnrm_inv(ic)
536          end if
537        end do
538      end do
539      if (izero>0.and.l_size>=1)  dgfact(izero,0)=-(alpha(1,1)*qq(1,1)+alpha(2,1)*qq(2,1))/three
540      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    else if (optgr2==1) then
543      do ll=0,l_size-1
544        do ic=1,nfgd
545          arg=rnrm(ic)
546          if (arg<=rcut) then
547            call paw_jbessel(jbes1,jbesp1,jbespp1,ll,2,qq(1,1+ll)*arg)
548            call paw_jbessel(jbes2,jbesp2,jbespp2,ll,2,qq(2,1+ll)*arg)
549            gfact(ic,ll)=shapefunc3(jbes1,jbes2,ll)
550            dgfact(ic,ll)=dshpfunc3(jbesp1,jbesp2,ll)*rnrm_inv(ic)
551            d2gfact(ic,ll)=(d2shpfunc3(jbespp1,jbespp2,ll)-dgfact(ic,ll))*rnrm_inv(ic)**2
552          end if
553        end do
554      end do
555      if (izero>0.and.l_size>=1)  dgfact(izero,0)=-(alpha(1,1)*qq(1,1)+alpha(2,1)*qq(2,1))/three
556      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  if (optgr0==1) then
564 
565    do ll=0,l_size-1
566      do ilm=ll**2+1,min((ll+1)**2,lm_size)
567        do ic=1,nfgd
568          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    if (izero>0) then
575      gylm(izero,1:lm_size)=zero
576      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  if(optgr1==1) then
584 
585    do ll=0,l_size-1
586      do ilm=ll**2+1,min((ll+1)**2,lm_size)
587        do ic=1,nfgd
588          gylmgr(1:3,ic,ilm)=gfact(ic,ll)*ylmrgr(1:3,ilm,ic)&
589 &         +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    if (izero>0) then
596      gylmgr(1:3,izero,1:lm_size)=zero
597      if (lm_size>=1) then
598        arg=cc(2,1)/sqrt(four_pi)
599        gylmgr(1:3,izero,1)=arg
600      end if
601      if (lm_size>=2) then
602        arg=cc(1,2)*sqrt(three/four_pi)
603        gylmgr(2,izero,2)=arg
604        if (lm_size>=3) gylmgr(3,izero,3)=arg
605        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  if(optgr2==1) then
614 
615    do ll=0,l_size-1
616      do ilm=ll**2+1,min((ll+1)**2,lm_size)
617        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 &         +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 &         +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 &         +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 &         +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 &         +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 &         +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    if (izero>0) then
642      gylmgr2(1:6,izero,1:lm_size)=zero
643      if (lm_size>=1) then
644        arg=cc(3,1)/sqrt(four_pi)
645        gylmgr2(1:3,izero,1)=arg
646      end if
647      if (lm_size>=2) then
648        arg=cc(2,2)*sqrt(three/four_pi)
649        gylmgr2(2,izero,2)=two*arg
650        gylmgr2(4,izero,2)=    arg
651        if (lm_size>=3) then
652          gylmgr2(1,izero,3)=two*arg
653          gylmgr2(3,izero,3)=two*arg
654        end if
655        if (lm_size>=4) then
656          gylmgr2(5,izero,4)=arg
657          gylmgr2(6,izero,4)=arg
658        end if
659      end if
660      if (lm_size>=5) then
661        arg=cc(1,3)*sqrt(15._dp/four_pi)
662        gylmgr2(6,izero,5)=arg
663        if (lm_size>=6) gylmgr2(4,izero,6)=arg
664        if (lm_size>=7) then
665          gylmgr2(1,izero,7)=   -arg/sqrt3
666          gylmgr2(2,izero,7)=   -arg/sqrt3
667          gylmgr2(3,izero,7)=two*arg/sqrt3
668        end if
669        if (lm_size>=8) gylmgr2(5,izero,8)=arg
670        if (lm_size>=9) then
671          gylmgr2(1,izero,9)= arg
672          gylmgr2(2,izero,9)=-arg
673        end if
674      end if
675    end if
676 
677  end if
678 
679 !Memory deallocation
680 !==========================================================
681  LIBPAW_DEALLOCATE(rnrm)
682  if (allocated(cc)) then
683    LIBPAW_DEALLOCATE(cc)
684  end if
685  if (compute_gr0)  then
686    LIBPAW_DEALLOCATE(gfact)
687  end if
688  if (compute_gr1)  then
689    LIBPAW_DEALLOCATE(dgfact)
690  end if
691  if (compute_gr2)  then
692    LIBPAW_DEALLOCATE(d2gfact)
693  end if
694  if (compute_gr1)  then
695    LIBPAW_DEALLOCATE(rnrm_inv)
696  end if
697  if (shape_type==3)  then
698    LIBPAW_DEALLOCATE(alpha)
699    LIBPAW_DEALLOCATE(qq)
700  end if
701  if (compute_gr0)  then
702    LIBPAW_DEALLOCATE(ylmr)
703  end if
704  if (compute_gr1)  then
705    LIBPAW_DEALLOCATE(ylmrgr)
706  end if
707  if (shape_type==-1) then
708    if (compute_gr0)  then
709      LIBPAW_DEALLOCATE(shpfuncnum)
710    end if
711    if (compute_gr1)  then
712      LIBPAW_DEALLOCATE(dshpfuncnum)
713    end if
714    if (compute_gr2)  then
715      LIBPAW_DEALLOCATE(d2shpfuncnum)
716    end if
717  end if
718 
719 ! -----------------------------------------------------------------
720 !Small functions related to analytical expression of shape function
721  CONTAINS

m_paw_finegrid/pawgylmg [ Functions ]

[ Top ] [ m_paw_finegrid ] [ Functions ]

NAME

 pawgylmg

FUNCTION

 PAW: Compute Fourier transform of each g_l(r).Y_lm(r) function

INPUTS

  gprimd(3,3)=dimensional reciprocal space primitive translations
  kg(3,npw)=integer coordinates of planewaves in basis sphere for this k point.
  kpg(npw,nkpg)= (k+G) components (only if useylm=1)
  kpt(3)=reduced coordinates of k point
  lmax=1+max. value of l angular momentum
  nkpg=second dimension of kpg_k (0 if useylm=0)
  npw=number of planewaves in basis sphere
  ntypat=number of types of atoms
  pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data
  ylm(npw,lmax**2)=real spherical harmonics for each G and k point

OUTPUT

  gylmg(npw,lmax**2,ntypat)=Fourier transform of each g_l(r).Y_lm(r) function

SOURCE

 955 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  real(dp),allocatable :: glg(:),qgrid(:),kpgnorm(:),shpf(:,:),work(:)
 974 
 975 ! *************************************************************************
 976 
 977 !Get |k+G|:
 978  LIBPAW_ALLOCATE(kpgnorm,(npw))
 979  if (nkpg<3) then
 980    do ig=1,npw
 981      kpg1=kpt(1)+dble(kg(1,ig));kpg2=kpt(2)+dble(kg(2,ig));kpg3=kpt(3)+dble(kg(3,ig))
 982      kpgc1=kpg1*gprimd(1,1)+kpg2*gprimd(1,2)+kpg3*gprimd(1,3)
 983      kpgc2=kpg1*gprimd(2,1)+kpg2*gprimd(2,2)+kpg3*gprimd(2,3)
 984      kpgc3=kpg1*gprimd(3,1)+kpg2*gprimd(3,2)+kpg3*gprimd(3,3)
 985      kpgnorm(ig)=sqrt(kpgc1*kpgc1+kpgc2*kpgc2+kpgc3*kpgc3)
 986    end do
 987  else
 988    do ig=1,npw
 989      kpgc1=kpg(ig,1)*gprimd(1,1)+kpg(ig,2)*gprimd(1,2)+kpg(ig,3)*gprimd(1,3)
 990      kpgc2=kpg(ig,1)*gprimd(2,1)+kpg(ig,2)*gprimd(2,2)+kpg(ig,3)*gprimd(2,3)
 991      kpgc3=kpg(ig,1)*gprimd(3,1)+kpg(ig,2)*gprimd(3,2)+kpg(ig,3)*gprimd(3,3)
 992      kpgnorm(ig)=sqrt(kpgc1*kpgc1+kpgc2*kpgc2+kpgc3*kpgc3)
 993    end do
 994  end if
 995 
 996  LIBPAW_ALLOCATE(glg,(npw))
 997  LIBPAW_ALLOCATE(work,(npw))
 998 
 999 !Loop over types of atoms
1000  do itypat=1,ntypat
1001 
1002    mqgrid=pawtab(itypat)%mqgrid_shp
1003    LIBPAW_ALLOCATE(qgrid,(mqgrid))
1004    LIBPAW_ALLOCATE(shpf,(mqgrid,2))
1005    qgrid(1:mqgrid)=pawtab(itypat)%qgrid_shp(1:mqgrid)
1006 
1007 !  Loops over (l,m) values
1008    do ll=0,pawtab(itypat)%lcut_size-1
1009      l0=ll**2+ll+1
1010 
1011      shpf(1:mqgrid,1:2)=pawtab(itypat)%shapefncg(1:mqgrid,1:2,1+ll)
1012      call paw_uniform_splfit(qgrid,work,shpf,0,kpgnorm,glg,mqgrid,npw)
1013 
1014      do mm=-ll,ll
1015        ilm=l0+mm
1016 
1017        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    LIBPAW_DEALLOCATE(qgrid)
1025    LIBPAW_DEALLOCATE(shpf)
1026  end do
1027 
1028  LIBPAW_DEALLOCATE(kpgnorm)
1029  LIBPAW_DEALLOCATE(glg)
1030  LIBPAW_DEALLOCATE(work)
1031 
1032 end subroutine pawgylmg

m_paw_finegrid/pawrfgd_fft [ Functions ]

[ Top ] [ m_paw_finegrid ] [ Functions ]

NAME

 pawrfgd_fft

FUNCTION

 Determine each point of the (fine) rectangular grid
 around a given atom and compute r-R vectors.
 R is the position of the atom.

INPUTS

  [fft_distrib(n3)]= (optional) index of processes which own fft planes in 3rd dimension
  [fft_index(n3)]= (optional) local fft indexes for current process
  gmet(3,3)=reciprocal space metric tensor in bohr**-2
  [me_fft]= (optional) my rank in the FFT MPI communicator
  n1,n2,n3= sizes of the FFT grid (entire simulation cell)
  rcut= radius of the sphere around the atom
  rprimd(3,3)=dimensional primitive translations in real space (bohr)
  ucvol= unit cell volume
  xred(3)= reduced coordinates of the atom

OUTPUT

  ifftsph(nfgd)= FFT index (fine grid) of the points in the sphere around current atom
  nfgd= number of points in the sphere around current atom
  rfgd(3,nfgd)= cartesian coordinates of r-R.

SOURCE

1064 subroutine pawrfgd_fft(ifftsph,gmet,n1,n2,n3,nfgd,rcut,rfgd,rprimd,ucvol,xred, &
1065 &                      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,parameter :: ishift=15
1082  integer :: i1,i2,i3,ifft_local,ix,iy,iz,izloc,me_fft_,n1a,n1b,n2a,n2b,n3a,n3b,ncmax
1083  real(dp),parameter :: delta=0.99_dp
1084  real(dp) :: difx,dify,difz,rr1,rr2,rr3,r2,r2cut,rx,ry,rz
1085  character(len=500) :: msg
1086 !arrays
1087  integer,allocatable :: ifftsph_tmp(:)
1088  integer,pointer :: fft_distrib_(:),fft_index_(:)
1089  real(dp),allocatable :: rfgd_tmp(:,:)
1090 
1091 ! *************************************************************************
1092 
1093 !Define a "box" around the atom
1094  r2cut=1.0000001_dp*rcut**2
1095  rr1=sqrt(r2cut*gmet(1,1))
1096  rr2=sqrt(r2cut*gmet(2,2))
1097  rr3=sqrt(r2cut*gmet(3,3))
1098  n1a=int((xred(1)-rr1+ishift)*n1+delta)-ishift*n1
1099  n1b=int((xred(1)+rr1+ishift)*n1      )-ishift*n1
1100  n2a=int((xred(2)-rr2+ishift)*n2+delta)-ishift*n2
1101  n2b=int((xred(2)+rr2+ishift)*n2      )-ishift*n2
1102  n3a=int((xred(3)-rr3+ishift)*n3+delta)-ishift*n3
1103  n3b=int((xred(3)+rr3+ishift)*n3      )-ishift*n3
1104 
1105 !Get the distrib associated with this fft_grid
1106  if (present(fft_distrib).and.present(fft_index).and.present(me_fft)) then
1107    me_fft_=me_fft ; fft_distrib_ => fft_distrib ; fft_index_ => fft_index
1108  else
1109    me_fft_=0
1110    LIBPAW_POINTER_ALLOCATE(fft_distrib_,(n3))
1111    LIBPAW_POINTER_ALLOCATE(fft_index_,(n3))
1112    fft_distrib_=0;fft_index_=(/(i3,i3=1,n3)/)
1113  end if
1114 
1115 !Temporary allocate "large" arrays
1116  ncmax=1+int(1.5_dp*(n1*n2*n3)*four_pi/(three*ucvol)*rcut**3)
1117  LIBPAW_ALLOCATE(ifftsph_tmp,(ncmax))
1118  LIBPAW_ALLOCATE(rfgd_tmp,(3,ncmax))
1119 
1120 !Set number of points to zero
1121  nfgd=0
1122 
1123 !Loop over FFT points
1124  do i3=n3a,n3b
1125    iz=mod(i3+ishift*n3,n3)
1126    if (fft_distrib_(iz+1)==me_fft_) then
1127      izloc=fft_index_(iz+1) - 1
1128      difz=dble(i3)/dble(n3)-xred(3)
1129      do i2=n2a,n2b
1130        iy=mod(i2+ishift*n2,n2)
1131        dify=dble(i2)/dble(n2)-xred(2)
1132        do i1=n1a,n1b
1133          ix=mod(i1+ishift*n1,n1)
1134          difx=dble(i1)/dble(n1)-xred(1)
1135 
1136 !        Compute r-R
1137          rx=difx*rprimd(1,1)+dify*rprimd(1,2)+difz*rprimd(1,3)
1138          ry=difx*rprimd(2,1)+dify*rprimd(2,2)+difz*rprimd(2,3)
1139          rz=difx*rprimd(3,1)+dify*rprimd(3,2)+difz*rprimd(3,3)
1140          r2=rx**2+ry**2+rz**2
1141 
1142 !        Select matching points
1143          if (r2 <= r2cut) then
1144            ifft_local=1+ix+n1*(iy+n2*izloc)
1145            if (ifft_local>0) then
1146              nfgd=nfgd+1
1147              if (nfgd>ncmax) then
1148                msg='Number of fft points around atom exceeds max. allowed!'
1149                LIBPAW_BUG(msg)
1150              end if
1151              rfgd_tmp(1,nfgd)=rx
1152              rfgd_tmp(2,nfgd)=ry
1153              rfgd_tmp(3,nfgd)=rz
1154              ifftsph_tmp(nfgd)=ifft_local
1155            end if
1156          end if
1157 
1158 !      End of loops
1159        end do
1160      end do
1161    end if
1162  end do
1163 
1164 !Now fill output arrays
1165  if (allocated(ifftsph)) then
1166    LIBPAW_DEALLOCATE(ifftsph)
1167  end if
1168  if (allocated(rfgd)) then
1169    LIBPAW_DEALLOCATE(rfgd)
1170  end if
1171  LIBPAW_ALLOCATE(ifftsph,(nfgd))
1172  LIBPAW_ALLOCATE(rfgd,(3,nfgd))
1173  ifftsph(1:nfgd)=ifftsph_tmp(1:nfgd)
1174  rfgd(1:3,1:nfgd)=rfgd_tmp(1:3,1:nfgd)
1175 
1176 !Release temporary memory
1177  LIBPAW_DEALLOCATE(ifftsph_tmp)
1178  LIBPAW_DEALLOCATE(rfgd_tmp)
1179  if (.not.present(fft_distrib).or..not.present(fft_index)) then
1180    LIBPAW_POINTER_DEALLOCATE(fft_distrib_)
1181    LIBPAW_POINTER_DEALLOCATE(fft_index_)
1182  end if
1183 
1184 end subroutine pawrfgd_fft

m_paw_finegrid/pawrfgd_wvl [ Functions ]

[ Top ] [ m_paw_finegrid ] [ Functions ]

NAME

 pawrfgd_wvl

FUNCTION

 Determine each point of the (fine) rectangular grid
 around a given atom and compute r-R vectors.
 R is the position of the atom.

INPUTS

  geocode= code for geometry (boundary conditions)
  hh(3)=fine grid spacing
  i3s= TO BE COMPLETED
  n1,n2,n3= TO BE COMPLETED
  n1i,n2i,n3pi= TO BE COMPLETED
  rcut= radius of the sphere around the atom
  rloc= cut-off radius for local psp?
  shift= TO BE COMPLETED
  xred(3)= cartesian coordinates of the atom

OUTPUT

  ifftsph(nfgd)= FFT index (fine grid) of the points in the sphere around current atom
  nfgd= number of points in the sphere around current atom
  rfgd(3,nfgd)= cartesian coordinates of r-R.

SOURCE

1216 subroutine pawrfgd_wvl(geocode,hh,ifftsph,i3s,n1,n1i,n2,n2i,n3,n3pi,&
1217 &                      nfgd,rcut,rloc,rfgd,shift,xcart)
1218 
1219 !Arguments ---------------------------------------------
1220 !scalars
1221  integer,intent(in) :: i3s,n1,n1i,n2,n2i,n3,n3pi,shift
1222  integer,intent(out) :: nfgd
1223  real(dp),intent(in) :: rcut,rloc
1224  character(1),intent(in) :: geocode
1225 !arrays
1226  integer,allocatable,intent(out) :: ifftsph(:)
1227  real(dp),intent(in) :: hh(3),xcart(3)
1228  real(dp),allocatable,intent(out) :: rfgd(:,:)
1229 
1230 !Local variables ------------------------------
1231 !scalars
1232  integer :: i1,i2,i3,iex,iey,iez,ind,isx,isy,isz,j1,j2,j3
1233  integer :: nbl1,nbr1,nbl2,nbr2,nbl3,nbr3,ncmax
1234  logical :: gox,goy,goz,perx,pery,perz
1235  real(dp) :: cutoff,r2,r2cut,rx,ry,rz,xx,yy,zz
1236 !arrays
1237  integer,allocatable :: ifftsph_tmp(:)
1238  real(dp),allocatable :: rfgd_tmp(:,:)
1239 
1240 ! *************************************************************************
1241 
1242 !Data for periodicity in the three directions
1243  perx=(geocode/='F')
1244  pery=(geocode=='P')
1245  perz=(geocode/='F')
1246  call my_ext_buffers(perx,nbl1,nbr1)
1247  call my_ext_buffers(pery,nbl2,nbr2)
1248  call my_ext_buffers(perz,nbl3,nbr3)
1249 
1250 !Define a "box" around the atom
1251  cutoff=10.d0*rloc
1252  r2cut=1.0000001_dp*rcut**2
1253  rx=xcart(1)
1254  ry=xcart(2)
1255  rz=xcart(3)
1256  isx=floor((rx-cutoff)/hh(1))
1257  isy=floor((ry-cutoff)/hh(2))
1258  isz=floor((rz-cutoff)/hh(3))
1259  iex=ceiling((rx+cutoff)/hh(1))
1260  iey=ceiling((ry+cutoff)/hh(2))
1261  iez=ceiling((rz+cutoff)/hh(3))
1262 
1263 !Temporary allocate "large" arrays
1264 !  use factor 1+int(1.1*, for safety reasons
1265  ncmax=1
1266  if (n3pi>0) ncmax=1+int((rcut/hh(1)+1.0)*(rcut/hh(2)+1.0)*(rcut/hh(3)+1.0)*four_pi/three)
1267  LIBPAW_ALLOCATE(ifftsph_tmp,(ncmax))
1268  LIBPAW_ALLOCATE(rfgd_tmp,(3,ncmax))
1269 
1270 !Set number of points to zero
1271  nfgd=0
1272 
1273 !Loop over WVL points
1274  do i3=isz,iez
1275    zz=real(i3,kind=8)*hh(3)-rz
1276    call my_ind_positions(perz,i3,n3,j3,goz)
1277    j3=j3+nbl3+1
1278    do i2=isy,iey
1279      yy=real(i2,kind=8)*hh(2)-ry
1280      call my_ind_positions(pery,i2,n2,j2,goy)
1281      do i1=isx,iex
1282        xx=real(i1,kind=8)*hh(1)-rx
1283        call my_ind_positions(perx,i1,n1,j1,gox)
1284        r2=xx**2+yy**2+zz**2
1285        if (j3>=i3s.and.j3<=i3s+n3pi-1.and.goy.and.gox) then
1286 
1287 !        Select matching points
1288          if (r2<=r2cut) then
1289            ind=j1+1+nbl1+(j2+nbl2)*n1i+(j3-i3s)*n1i*n2i
1290            nfgd=nfgd+1
1291            rfgd_tmp(:,nfgd)=[xx,yy,zz]
1292            ifftsph_tmp(nfgd)=shift+ind
1293          end if
1294 
1295 !      End of loops
1296        end if
1297      end do
1298    end do
1299  end do
1300 
1301 !Now fill output arrays
1302  if (allocated(ifftsph)) then
1303    LIBPAW_DEALLOCATE(ifftsph)
1304  end if
1305  if (allocated(rfgd)) then
1306    LIBPAW_DEALLOCATE(rfgd)
1307  end if
1308  LIBPAW_ALLOCATE(ifftsph,(nfgd))
1309  LIBPAW_ALLOCATE(rfgd,(3,nfgd))
1310  ifftsph(1:nfgd)=ifftsph_tmp(1:nfgd)
1311  rfgd(1:3,1:nfgd)=rfgd_tmp(1:3,1:nfgd)
1312 
1313 !Release temporary memory
1314  LIBPAW_DEALLOCATE(ifftsph_tmp)
1315  LIBPAW_DEALLOCATE(rfgd_tmp)
1316 
1317 !*********************************************************************
1318 !Small functions related to boundary conditions
1319  contains

m_paw_finegrid/shapefunc1 [ Functions ]

[ Top ] [ m_paw_finegrid ] [ Functions ]

m_paw_finegrid/shapefunc1_0 [ Functions ]

[ Top ] [ m_paw_finegrid ] [ Functions ]

m_paw_finegrid/shapefunc2 [ Functions ]

[ Top ] [ m_paw_finegrid ] [ Functions ]

m_paw_finegrid/shapefunc2_0 [ Functions ]

[ Top ] [ m_paw_finegrid ] [ Functions ]

m_paw_finegrid/shapefunc3 [ Functions ]

[ Top ] [ m_paw_finegrid ] [ Functions ]