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-2018 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

22 #include "libpaw.h"
23 
24 MODULE m_paw_finegrid
25 
26  USE_DEFS
27  USE_MSG_HANDLING
28  USE_MEMORY_PROFILING
29 
30  use m_pawtab,      only : pawtab_type
31  use m_paw_sphharm, only : initylmr
32  use m_paw_numeric, only : paw_jbessel,paw_splint,paw_uniform_splfit,paw_sort_dp
33 
34  implicit none
35 
36  private
37 
38 !public procedures.
39  public :: pawgylm      ! g_l(r-R)*Y_lm(r-R) (and derivatives)
40  public :: pawgylmg     ! Fourier transform of g_l(r-R)*Y_lm(r-R), plane-waves case
41  public :: pawrfgd_fft  ! r-R, plane-waves case
42  public :: pawrfgd_wvl  ! r-R, wavelets case
43  public :: pawexpiqr    ! exp(i.q.(r-R))
44 
45 !declarations for the whole module (were needed to replace the statement functions)
46 !MG: Why this? Global variables are powerful but extremely DANGEROUS
47 integer,private,save :: lambda
48 real(dp),private,save :: pi_over_rshp,sigma
49 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 !

PARENTS

      m_pawdij,pawgrnl,pawmknhat,pawmknhat_psipsi,respfn

CHILDREN

SOURCE

1591 subroutine pawexpiqr(expiqr,gprimd,nfgd,qphon,rfgd,xred)
1592 
1593 
1594 !This section has been created automatically by the script Abilint (TD).
1595 !Do not modify the following lines by hand.
1596 #undef ABI_FUNC
1597 #define ABI_FUNC 'pawexpiqr'
1598 !End of the abilint section
1599 
1600  implicit none
1601 
1602 !Arguments ---------------------------------------------
1603 !scalars
1604  integer,intent(in) :: nfgd
1605 !arrays
1606  real(dp),intent(in) :: gprimd(3,3),qphon(3),xred(3)
1607  real(dp),intent(in) :: rfgd(:,:)
1608  real(dp),intent(out) :: expiqr(2,nfgd)
1609 
1610 !Local variables ------------------------------
1611 !scalars
1612  integer :: ic
1613  logical :: qne0
1614  real(dp) :: phase,phase_xred,qx,qy,qz
1615  character(len=500) :: msg
1616 !arrays
1617 
1618 ! *************************************************************************
1619 
1620  if (size(rfgd)/=3*nfgd) then
1621    msg='rfgd array must be allocated!'
1622    MSG_BUG(msg)
1623  end if
1624 
1625  qne0=(qphon(1)**2+qphon(2)**2+qphon(3)**2>=1.d-15)
1626 
1627 !Compute q in cartesian coordinates
1628  if (qne0) then
1629    qx=gprimd(1,1)*qphon(1)+gprimd(1,2)*qphon(2)+gprimd(1,3)*qphon(3)
1630    qy=gprimd(2,1)*qphon(1)+gprimd(2,2)*qphon(2)+gprimd(2,3)*qphon(3)
1631    qz=gprimd(3,1)*qphon(1)+gprimd(3,2)*qphon(2)+gprimd(3,3)*qphon(3)
1632    phase_xred=two_pi*(qphon(1)*xred(1)+qphon(2)*xred(2)+qphon(3)*xred(3))
1633  end if
1634 
1635 !Compute exp(i.q.r)
1636  if (qne0) then
1637    do ic=1,nfgd
1638      phase=two_pi*(qx*rfgd(1,ic)+qy*rfgd(2,ic)+qz*rfgd(3,ic)) + phase_xred
1639      expiqr(1,ic)=cos(phase)
1640      expiqr(2,ic)=sin(phase)
1641    end do
1642  end if
1643 
1644 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

PARENTS

      m_pawdij,nhatgrid,paw_mknewh0,pawgrnl,pawmknhat,pawmknhat_psipsi
      pawnhatfr,wvl_nhatgrid

CHILDREN

SOURCE

 93 subroutine pawgylm(gylm,gylmgr,gylmgr2,lm_size,nfgd,optgr0,optgr1,optgr2,pawtab,rfgd)
 94 
 95 
 96 !This section has been created automatically by the script Abilint (TD).
 97 !Do not modify the following lines by hand.
 98 #undef ABI_FUNC
 99 #define ABI_FUNC 'pawgylm'
100 !End of the abilint section
101 
102  implicit none
103 
104 !Arguments ---------------------------------------------
105 !scalars
106  integer,intent(in) :: lm_size,nfgd,optgr0,optgr1,optgr2
107  type(pawtab_type),intent(in) :: pawtab
108 !arrays
109  real(dp),intent(in) :: rfgd(:,:)
110  real(dp),intent(out) :: gylm(nfgd,optgr0*lm_size)
111  real(dp),intent(out) :: gylmgr(3,nfgd,optgr1*lm_size)
112  real(dp),intent(out) :: gylmgr2(6,nfgd,optgr2*lm_size)
113 
114 !Local variables ------------------------------
115 !scalars
116  integer :: ic,ilm,izero,l_size,ll,normchoice,option,shape_type
117  real(dp) :: arg
118  real(dp) :: jbes1,jbes2,jbesp1,jbesp2,jbespp1,jbespp2,rcut
119  real(dp) :: splfact
120  logical :: compute_gr0,compute_gr1,compute_gr2
121  character(len=500) :: msg
122 !arrays
123  integer,allocatable :: isort(:)
124  real(dp),parameter :: ffact(1:9)=(/1._dp,3._dp,15._dp,105._dp,945._dp,10395._dp,&
125 &                                   135135._dp,2027025._dp,34459425._dp/)
126  real(dp),parameter :: toldev=tol3
127  real(dp) :: ss(3)
128  real(dp),allocatable :: cc(:,:),d2gfact(:,:),d2shpfuncnum(:,:),dgfact(:,:)
129  real(dp),allocatable :: dshpfuncnum(:,:),gfact(:,:)
130  real(dp),allocatable :: rnrm(:),rnrm_inv(:),rnrm_sort(:)
131  real(dp),allocatable :: shpfuncnum(:,:),work(:),ylmr(:,:),ylmrgr(:,:,:)
132 
133 ! *************************************************************************
134 
135  if (optgr0==0.and.optgr1==0.and.optgr2==0) return
136  if (nfgd==0) return
137 
138 !Compatibility test
139 !==========================================================
140  if (size(rfgd)/=3*nfgd) then
141    msg='rfgd array must be allocated at rfgd(3,nfgd)!'
142    MSG_BUG(msg)
143  end if
144  if (pawtab%lcut_size>9) then
145    msg='l_size>10 forbidden!'
146    MSG_BUG(msg)
147  end if
148  if (pawtab%shape_type==1.and.pawtab%shape_lambda<2) then
149    msg='Exponent lambda of gaussian shape function must be > 1!'
150    MSG_ERROR(msg)
151  end if
152 
153 !Initializations
154 !==========================================================
155 !Options for computation
156  compute_gr0=(optgr0==1.or.optgr1==1.or.optgr2==1)
157  compute_gr1=(optgr1==1.or.optgr2==1)
158  compute_gr2=(optgr2==1)
159  l_size=pawtab%lcut_size
160 
161 !Norms of vectors around the atom
162  LIBPAW_ALLOCATE(rnrm,(nfgd))
163  izero=-1
164  do ic=1,nfgd
165    rnrm(ic)=sqrt(rfgd(1,ic)**2+rfgd(2,ic)**2+rfgd(3,ic)**2)
166    if (rnrm(ic)<=tol10) izero=ic  ! Has to be consistent with initylmr !!
167  end do
168 
169 !Initializations
170  if (optgr0==1) gylm=zero
171  if (optgr1==1) gylmgr=zero
172  if (optgr2==1) gylmgr2=zero
173 
174 !Some definitions concerning shape function g_l(r)
175  shape_type=pawtab%shape_type
176  sigma=pawtab%shape_sigma;lambda=pawtab%shape_lambda
177  pi_over_rshp=pi/pawtab%rshp
178  rcut=tol12+pawtab%rshp
179  if (shape_type==3) then
180    LIBPAW_ALLOCATE(alpha,(2,l_size))
181    LIBPAW_ALLOCATE(qq,(2,l_size))
182    do ll=1,l_size
183      alpha(1:2,ll)=pawtab%shape_alpha(1:2,ll)
184      qq(1:2,ll)=pawtab%shape_q(1:2,ll)
185    end do
186  end if
187 
188 !If needed, sort selected radii by increasing norm
189  if (shape_type==-1) then
190    LIBPAW_ALLOCATE(isort,(nfgd))
191    LIBPAW_ALLOCATE(rnrm_sort,(nfgd))
192    do ic=1,nfgd
193      isort(ic)=ic
194    end do
195    rnrm_sort(1:nfgd)=rnrm(1:nfgd)
196    call paw_sort_dp(nfgd,rnrm_sort,isort,tol16)
197  end if
198 
199 !If shape function is "numeric", spline it onto selected radii
200  if (shape_type==-1) then
201    LIBPAW_ALLOCATE(work,(nfgd))
202    if (compute_gr0) then
203      LIBPAW_ALLOCATE(shpfuncnum,(nfgd,l_size))
204      do ll=1,l_size
205        call paw_splint(pawtab%mesh_size,pawtab%rad_for_spline,pawtab%shapefunc(:,ll),&
206 &       pawtab%dshpfunc(:,ll,2),nfgd,rnrm_sort,work)
207        do ic=1,nfgd
208          shpfuncnum(isort(ic),ll)=work(ic)
209        end do
210      end do
211    end if
212    if(compute_gr1) then
213      LIBPAW_ALLOCATE(dshpfuncnum,(nfgd,l_size))
214      do ll=1,l_size
215        call paw_splint(pawtab%mesh_size,pawtab%rad_for_spline,pawtab%dshpfunc(:,ll,1),&
216 &       pawtab%dshpfunc(:,ll,3),nfgd,rnrm_sort,work)
217        do ic=1,nfgd
218          dshpfuncnum(isort(ic),ll)=work(ic)
219        end do
220      end do
221    end if
222    if(compute_gr2) then
223      LIBPAW_ALLOCATE(d2shpfuncnum,(nfgd,l_size))
224      do ll=1,l_size
225        call paw_splint(pawtab%mesh_size,pawtab%rad_for_spline,pawtab%dshpfunc(:,ll,2),&
226 &       pawtab%dshpfunc(:,ll,4),nfgd,rnrm_sort,work)
227        do ic=1,nfgd
228          d2shpfuncnum(isort(ic),ll)=work(ic)
229        end do
230      end do
231    end if
232    LIBPAW_DEALLOCATE(work)
233  end if
234 
235  if (shape_type==-1)  then
236    LIBPAW_DEALLOCATE(isort)
237    LIBPAW_DEALLOCATE(rnrm_sort)
238  end if
239 
240 !If needed, compute limits at r=0 of shape function and derivatives
241  if (izero>0) then
242    LIBPAW_ALLOCATE(cc,(3,min(l_size,3)))
243    cc=zero
244    if (shape_type==-1) then
245      splfact=(pawtab%rad_for_spline(4)-pawtab%rad_for_spline(1))&
246 &     /(pawtab%rad_for_spline(3)-pawtab%rad_for_spline(2))
247    end if
248    do ll=1,min(l_size,3)
249 !    cc(2,l) is g_prime(0)
250      if (optgr0==1.or.optgr1==1.or.optgr2==1) then
251        if (shape_type==-1) then
252          ss(1:3)=pawtab%shapefunc(2:4,ll)/pawtab%rad_for_spline(2:4)**(ll-1)
253          cc(1,ll)=ss(3)+(ss(1)-ss(2))*splfact
254        else if (shape_type==1.or.shape_type==2) then
255          cc(1,ll)=one
256        else if (shape_type==3) then
257          cc(1,ll)=(alpha(1,ll)*qq(1,ll)**(ll-1) &
258 &         +alpha(2,ll)*qq(2,ll)**(ll-1))/ffact(ll)
259        end if
260        cc(1,ll)=cc(1,ll)*pawtab%gnorm(ll)
261      end if
262 !    cc(2,l) is g_prime(0)
263      if (optgr1==1.or.optgr2==1) then
264        if (shape_type==-1) then
265          ss(1:3)=(ss(1:3)-cc(1,ll))/pawtab%rad_for_spline(2:4)
266          cc(2,ll)=ss(3)+(ss(1)-ss(2))*splfact
267        else if (shape_type==1.and.lambda==1) then
268          cc(2,ll)=-one/sigma
269        else
270          cc(2,ll)=zero
271        end if
272        cc(2,ll)=cc(2,ll)*pawtab%gnorm(ll)
273      end if
274 !    cc(3,l) is g_prime_prime(0)
275      if (optgr2==1) then
276        if (shape_type==-1) then
277          ss(1:3)=(ss(1:3)-cc(2,ll))/pawtab%rad_for_spline(2:4)
278          cc(3,ll)=two*(ss(3)+(ss(1)-ss(2))*splfact)
279        else if (shape_type==1) then
280          if (lambda==1) cc(3,ll)=one/sigma**2
281          if (lambda==2) cc(3,ll)=-two/sigma**2
282          if (lambda >2) cc(3,ll)=zero
283        else if (shape_type==2) then
284          cc(3,ll)=-(two/three)*pi_over_rshp**2
285        else if (shape_type==3) then
286          cc(3,ll)=-(alpha(1,ll)*qq(1,ll)**(ll+1) &
287 &         +alpha(2,ll)*qq(2,ll)**(ll+1))/ffact(ll+1)
288        end if
289        cc(3,ll)=cc(3,ll)*pawtab%gnorm(ll)
290      end if
291    end do
292  end if
293 
294 !Y_lm(r-R) calculation
295 !==========================================================
296  normchoice=1 ; option=max(optgr0,2*optgr1,3*optgr2)
297  if(compute_gr0)  then
298    LIBPAW_ALLOCATE(ylmr,(l_size**2,nfgd))
299  end if
300  if(compute_gr1.and.(.not.compute_gr2))  then
301    LIBPAW_ALLOCATE(ylmrgr,(3,l_size**2,nfgd))
302  end if
303  if(compute_gr2)  then
304    LIBPAW_ALLOCATE(ylmrgr,(9,l_size**2,nfgd))
305  end if
306  if (compute_gr0.and.(.not.compute_gr1).and.(.not.compute_gr2)) then
307    call initylmr(l_size,normchoice,nfgd,rnrm,option,rfgd,ylmr)
308  else
309    call initylmr(l_size,normchoice,nfgd,rnrm,option,rfgd,ylmr,ylmrgr)
310  end if
311 
312 !gl(r) and derivatives calculation for l>=0
313 !==========================================================
314 !Compute gl(r), gl_prime(r)/r and (gl_prime_prime(r)-gl_prime(r)/r)/r**2
315  if (compute_gr0)  then
316    LIBPAW_BOUND2_ALLOCATE(gfact,BOUNDS(1,nfgd),BOUNDS(0,l_size-1))
317    gfact(:,:)=zero
318  end if
319  if (compute_gr1)  then
320    LIBPAW_BOUND2_ALLOCATE(dgfact,BOUNDS(1,nfgd),BOUNDS(0,l_size-1))
321    dgfact(:,:)=zero
322  end if
323  if (compute_gr2)  then
324    LIBPAW_BOUND2_ALLOCATE(d2gfact,BOUNDS(1,nfgd),BOUNDS(0,l_size-1))
325    d2gfact(:,:)=zero
326  end if
327  if(compute_gr1) then
328    LIBPAW_ALLOCATE(rnrm_inv,(nfgd))
329    do ic=1,nfgd
330      if (ic/=izero) rnrm_inv(ic)=one/rnrm(ic)
331    end do
332    if (izero>0) rnrm_inv(izero)=zero
333  end if
334 
335 !----- type -1 -----
336  if (shape_type==-1) then
337    if (compute_gr0) then
338      do ll=0,l_size-1
339        do ic=1,nfgd
340          if (rnrm(ic)<=rcut) then
341            gfact(ic,ll)=shpfuncnum(ic,ll+1)
342          end if
343        end do
344      end do
345    end if
346    if (compute_gr1) then
347      do ll=0,l_size-1
348        do ic=1,nfgd
349          if (rnrm(ic)<=rcut) then
350            dgfact(ic,ll)=dshpfuncnum(ic,ll+1)*rnrm_inv(ic)
351          end if
352        end do
353      end do
354    end if
355    if(compute_gr2) then
356      do ll=0,l_size-1
357        do ic=1,nfgd
358          if (rnrm(ic)<=rcut) then
359            d2gfact(ic,ll)=(d2shpfuncnum(ic,ll+1)-dgfact(ic,ll))*rnrm_inv(ic)**2
360          end if
361        end do
362     end do
363    end if
364 
365 !  ----- type 1 or 2 -----
366  else if (shape_type==1.or.shape_type==2) then
367 !  FIRST COMPUTE FACTORS FOR l=0
368    if (optgr0==1.and.optgr1==0.and.optgr2==0) then
369      if (shape_type==1) then
370        do ic=1,nfgd
371          arg=rnrm(ic)
372          if (arg<toldev) then
373            gfact(ic,0)=shapefunc1_0(arg)
374          else if (arg<=rcut) then
375            gfact(ic,0)=shapefunc1(arg)
376          end if
377        end do
378      else ! shape_type==2
379        do ic=1,nfgd
380          arg=rnrm(ic)
381          if (arg<toldev) then
382            gfact(ic,0)=shapefunc2_0(arg)
383          else if (arg<=rcut) then
384            gfact(ic,0)=shapefunc2(arg)
385          end if
386        end do
387      end if
388    else if (optgr1==1.and.optgr2==0) then
389      if (shape_type==1) then
390        do ic=1,nfgd
391          arg=rnrm(ic)
392          if (arg<toldev) then
393            gfact(ic,0)=shapefunc1_0(arg)
394            if (lambda==2) then
395              dgfact(ic,0)=dshpfunc1_ovr_0_2(arg)
396            else ! lambda>2
397              dgfact(ic,0)=dshpfunc1_ovr_0(arg)
398            end if
399          else if (arg<=rcut) then
400            gfact(ic,0)=shapefunc1(arg)
401            dgfact(ic,0)=dshpfunc1(arg)*rnrm_inv(ic)
402          end if
403        end do
404      else ! shape_type==2
405        do ic=1,nfgd
406          arg=rnrm(ic)
407          if (arg<toldev) then
408            gfact(ic,0)=shapefunc2_0(arg)
409            dgfact(ic,0)=dshpfunc2_ovr_0(arg)
410          else if (arg<=rcut) then
411            gfact(ic,0)=shapefunc2(arg)
412            dgfact(ic,0)=dshpfunc2(arg)*rnrm_inv(ic)
413          end if
414        end do
415      end if
416    else if (optgr2==1) then
417      if (shape_type==1) then
418        do ic=1,nfgd
419          arg=rnrm(ic)
420          if (arg<toldev) then
421            gfact(ic,0)=shapefunc1_0(arg)
422            if (lambda==2) then
423              dgfact(ic,0)=dshpfunc1_ovr_0_2(arg)
424              d2gfact(ic,0)=d2shpfunc1_ovr2_0_2(arg)
425            else if (lambda==3) then
426              dgfact(ic,0)=dshpfunc1_ovr_0(arg)
427              if (ic/=izero) then
428                d2gfact(ic,0)=d2shpfunc1_ovr2_0_3(arg)
429              else
430                d2gfact(ic,0)=zero ! Diverging case
431              end if
432            else if (lambda==4) then
433              dgfact(ic,0)=dshpfunc1_ovr_0(arg)
434              d2gfact(ic,0)=d2shpfunc1_ovr2_0_4(arg)
435            else ! lambda>4
436              dgfact(ic,0)=dshpfunc1_ovr_0(arg)
437              d2gfact(ic,0)=d2shpfunc1_ovr2_0(arg)
438            end if
439          else if (arg<=rcut) then
440            gfact(ic,0)=shapefunc1(arg)
441            dgfact(ic,0)=dshpfunc1(arg)*rnrm_inv(ic)
442            d2gfact(ic,0)=(d2shpfunc1(arg)-dgfact(ic,0))*rnrm_inv(ic)**2
443          end if
444        end do
445      else ! shape_type==2
446        do ic=1,nfgd
447          arg=rnrm(ic)
448          if (arg<toldev) then
449            gfact(ic,0)=shapefunc2_0(arg)
450            dgfact(ic,0)=dshpfunc2_ovr_0(arg)
451            d2gfact(ic,0)=d2shpfunc2_ovr2_0(arg)
452          else if (arg<=rcut) then
453            gfact(ic,0)=shapefunc2(arg)
454            dgfact(ic,0)=dshpfunc2(arg)*rnrm_inv(ic)
455            d2gfact(ic,0)=(d2shpfunc2(arg)-dgfact(ic,0))*rnrm_inv(ic)**2
456          end if
457        end do
458      end if
459    end if
460 
461 !  THEN COMPUTE FACTORS FOR l>0 (from l=0)
462    if (compute_gr0) then
463      if (l_size>1) then
464        do ll=1,l_size-1
465          do ic=1,nfgd
466            gfact(ic,ll)=pawtab%gnorm(ll+1)*gfact(ic,0)*rnrm(ic)**ll
467          end do
468        end do
469      end if
470    end if
471    if (compute_gr1) then
472      if (l_size>1) then
473        do ic=1,nfgd
474          dgfact(ic,1)=pawtab%gnorm(2)*(gfact(ic,0)*rnrm_inv(ic)+dgfact(ic,0)*rnrm(ic))
475        end do
476      end if
477      if (l_size>2) then
478        do ic=1,nfgd
479          dgfact(ic,2)=pawtab%gnorm(3)*(two*gfact(ic,0)+dgfact(ic,0)*rnrm(ic)**2)
480        end do
481      end if
482      if (l_size>3) then
483        do ll=3,l_size-1
484          do ic=1,nfgd
485            dgfact(ic,ll)=pawtab%gnorm(ll+1) &
486 &           *(dble(ll)*gfact(ic,0)*rnrm(ic)**(ll-2)+dgfact(ic,0)*rnrm(ic)**ll)
487          end do
488        end do
489      end if
490    end if
491    if (compute_gr2) then
492      if (l_size>1) then
493        do ic=1,nfgd
494          d2gfact(ic,1)=pawtab%gnorm(2) &
495 &         *(-gfact(ic,0)*rnrm_inv(ic)**3+two*dgfact(ic,0)*rnrm_inv(ic)+d2gfact(ic,0)*rnrm(ic))
496        end do
497      end if
498      if (l_size>2) then
499        do ic=1,nfgd
500          d2gfact(ic,2)=pawtab%gnorm(3)*(four*dgfact(ic,0)+d2gfact(ic,0)*rnrm(ic)**2)
501        end do
502      end if
503      if (l_size>3) then
504        do ic=1,nfgd
505          d2gfact(ic,3)=pawtab%gnorm(4) &
506 &         *(three*gfact(ic,0)*rnrm_inv(ic)+6._dp*dgfact(ic,0)*rnrm(ic)+ d2gfact(ic,0)*rnrm(ic)**3)
507        end do
508      end if
509      if (l_size>4) then
510        do ic=1,nfgd
511          d2gfact(ic,4)=pawtab%gnorm(5) &
512 &         *(8._dp*gfact(ic,0)+8._dp*dgfact(ic,0)*rnrm(ic)**2+d2gfact(ic,0)*rnrm(ic)**4)
513        end do
514      end if
515      if (l_size>5) then
516        do ll=5,l_size-1
517          do ic=1,nfgd
518            d2gfact(ic,ll)=pawtab%gnorm(ll+1) &
519 &           *(dble(ll*(ll-2))*gfact(ic,0)*rnrm(ic)**(ll-4) &
520 &            +dble(2*ll)*dgfact(ic,0)*rnrm(ic)**(ll-2) &
521 &            +d2gfact(ic,0)*rnrm(ic)**ll)
522          end do
523        end do
524      end if
525    end if
526    if (compute_gr0) gfact(:,0)=gfact(:,0)*pawtab%gnorm(1)
527    if (compute_gr1) dgfact(:,0)=dgfact(:,0)*pawtab%gnorm(1)
528    if (compute_gr2) d2gfact(:,0)=d2gfact(:,0)*pawtab%gnorm(1)
529 
530 !  ----- type 3 -----
531  else if (shape_type==3) then
532    if (optgr0==1.and.optgr1==0.and.optgr2==0) then
533      do ll=0,l_size-1
534        do ic=1,nfgd
535          arg=rnrm(ic)
536          if (arg<=rcut) then
537            call paw_jbessel(jbes1,jbesp1,jbespp1,ll,0,qq(1,1+ll)*arg)
538            call paw_jbessel(jbes2,jbesp2,jbespp2,ll,0,qq(2,1+ll)*arg)
539            gfact(ic,ll)=shapefunc3(jbes1,jbes2,ll)
540          end if
541        end do
542      end do
543    else if (optgr1==1.and.optgr2==0) then
544      do ll=0,l_size-1
545        do ic=1,nfgd
546          arg=rnrm(ic)
547          if (arg<=rcut) then
548            call paw_jbessel(jbes1,jbesp1,jbespp1,ll,1,qq(1,1+ll)*arg)
549            call paw_jbessel(jbes2,jbesp2,jbespp2,ll,1,qq(2,1+ll)*arg)
550            gfact(ic,ll)=shapefunc3(jbes1,jbes2,ll)
551            dgfact(ic,ll)=dshpfunc3(jbesp1,jbesp2,ll)*rnrm_inv(ic)
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    else if (optgr2==1) then
559      do ll=0,l_size-1
560        do ic=1,nfgd
561          arg=rnrm(ic)
562          if (arg<=rcut) then
563            call paw_jbessel(jbes1,jbesp1,jbespp1,ll,2,qq(1,1+ll)*arg)
564            call paw_jbessel(jbes2,jbesp2,jbespp2,ll,2,qq(2,1+ll)*arg)
565            gfact(ic,ll)=shapefunc3(jbes1,jbes2,ll)
566            dgfact(ic,ll)=dshpfunc3(jbesp1,jbesp2,ll)*rnrm_inv(ic)
567            d2gfact(ic,ll)=(d2shpfunc3(jbespp1,jbespp2,ll)-dgfact(ic,ll))*rnrm_inv(ic)**2
568          end if
569        end do
570      end do
571      if (izero>0.and.l_size>=1)  dgfact(izero,0)=-(alpha(1,1)*qq(1,1)+alpha(2,1)*qq(2,1))/three
572      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))
573 !    Note: for l=1, dgfact is diverging - d2gfact is diverging for l<4
574    end if
575  end if
576 
577 !g_l(r-R)*Y_lm(r-R) calculation
578 !==========================================================
579  if (optgr0==1) then
580 
581    do ll=0,l_size-1
582      do ilm=ll**2+1,min((ll+1)**2,lm_size)
583        do ic=1,nfgd
584          gylm(ic,ilm)=gfact(ic,ll)*ylmr(ilm,ic)
585        end do
586      end do
587    end do
588 
589 !  Special value at r-R=0  (supposing shapefunc(r)->C.r**l when r->0)
590    if (izero>0) then
591      gylm(izero,1:lm_size)=zero
592      if (lm_size>=1) gylm(izero,1)=ylmr(1,izero)*cc(1,1)
593    end if
594 
595  end if
596 
597 !d/dr{g_l(r-R)*Y_lm(r-R)} calculation
598 !==========================================================
599  if(optgr1==1) then
600 
601    do ll=0,l_size-1
602      do ilm=ll**2+1,min((ll+1)**2,lm_size)
603        do ic=1,nfgd
604          gylmgr(1:3,ic,ilm)=gfact(ic,ll)*ylmrgr(1:3,ilm,ic)&
605 &         +dgfact(ic,ll)*rfgd(1:3,ic)*ylmr(ilm,ic)
606        end do
607      end do
608    end do
609 
610 !  Special values at r-R=0  (supposing shapefunc(r)->C.r**l when r->0)
611    if (izero>0) then
612      gylmgr(1:3,izero,1:lm_size)=zero
613      if (lm_size>=1) then
614        arg=cc(2,1)/sqrt(four_pi)
615        gylmgr(1:3,izero,1)=arg
616      end if
617      if (lm_size>=2) then
618        arg=cc(1,2)*sqrt(three/four_pi)
619        gylmgr(2,izero,2)=arg
620        if (lm_size>=3) gylmgr(3,izero,3)=arg
621        if (lm_size>=4) gylmgr(1,izero,4)=arg
622      end if
623    end if
624 
625  end if
626 
627 !d2/dridrj{g_l(r-R)*Y_lm(r-R)} calculation
628 !==========================================================
629  if(optgr2==1) then
630 
631    do ll=0,l_size-1
632      do ilm=ll**2+1,min((ll+1)**2,lm_size)
633        do ic=1,nfgd
634          gylmgr2(1,ic,ilm)=gfact(ic,ll)*ylmrgr(4,ilm,ic) &
635 &         +dgfact(ic,ll)*(ylmr(ilm,ic)+two*rfgd(1,ic)*ylmrgr(1,ilm,ic)) &
636 &         +d2gfact(ic,ll)*ylmr(ilm,ic)*rfgd(1,ic)*rfgd(1,ic)
637          gylmgr2(2,ic,ilm)=gfact(ic,ll)*ylmrgr(5,ilm,ic) &
638 &         +dgfact(ic,ll)*(ylmr(ilm,ic)+two*rfgd(2,ic)*ylmrgr(2,ilm,ic)) &
639 &         +d2gfact(ic,ll)*ylmr(ilm,ic)*rfgd(2,ic)*rfgd(2,ic)
640          gylmgr2(3,ic,ilm)=gfact(ic,ll)*ylmrgr(6,ilm,ic) &
641 &         +dgfact(ic,ll)*(ylmr(ilm,ic)+two*rfgd(3,ic)*ylmrgr(3,ilm,ic)) &
642 &         +d2gfact(ic,ll)*ylmr(ilm,ic)*rfgd(3,ic)*rfgd(3,ic)
643          gylmgr2(4,ic,ilm)=gfact(ic,ll)*ylmrgr(7,ilm,ic) &
644 &         +dgfact(ic,ll)*(rfgd(3,ic)*ylmrgr(2,ilm,ic)+rfgd(2,ic)*ylmrgr(3,ilm,ic)) &
645 &         +d2gfact(ic,ll)*ylmr(ilm,ic)*rfgd(3,ic)*rfgd(2,ic)
646          gylmgr2(5,ic,ilm)=gfact(ic,ll)*ylmrgr(8,ilm,ic) &
647 &         +dgfact(ic,ll)*(rfgd(3,ic)*ylmrgr(1,ilm,ic)+rfgd(1,ic)*ylmrgr(3,ilm,ic)) &
648 &         +d2gfact(ic,ll)*ylmr(ilm,ic)*rfgd(3,ic)*rfgd(1,ic)
649          gylmgr2(6,ic,ilm)=gfact(ic,ll)*ylmrgr(9,ilm,ic) &
650 &         +dgfact(ic,ll)*(rfgd(1,ic)*ylmrgr(2,ilm,ic)+rfgd(2,ic)*ylmrgr(1,ilm,ic)) &
651 &         +d2gfact(ic,ll)*ylmr(ilm,ic)*rfgd(1,ic)*rfgd(2,ic)
652        end do
653      end do
654    end do
655 
656 !  Special values at r-R=0  (supposing shapefunc(r)->C.r**l when r->0)
657    if (izero>0) then
658      gylmgr2(1:6,izero,1:lm_size)=zero
659      if (lm_size>=1) then
660        arg=cc(3,1)/sqrt(four_pi)
661        gylmgr2(1:3,izero,1)=arg
662      end if
663      if (lm_size>=2) then
664        arg=cc(2,2)*sqrt(three/four_pi)
665        gylmgr2(2,izero,2)=two*arg
666        gylmgr2(4,izero,2)=    arg
667        if (lm_size>=3) then
668          gylmgr2(1,izero,3)=two*arg
669          gylmgr2(3,izero,3)=two*arg
670        end if
671        if (lm_size>=4) then
672          gylmgr2(5,izero,4)=arg
673          gylmgr2(6,izero,4)=arg
674        end if
675      end if
676      if (lm_size>=5) then
677        arg=cc(1,3)*sqrt(15._dp/four_pi)
678        gylmgr2(6,izero,5)=arg
679        if (lm_size>=6) gylmgr2(4,izero,6)=arg
680        if (lm_size>=7) then
681          gylmgr2(1,izero,7)=   -arg/sqrt3
682          gylmgr2(2,izero,7)=   -arg/sqrt3
683          gylmgr2(3,izero,7)=two*arg/sqrt3
684        end if
685        if (lm_size>=8) gylmgr2(5,izero,8)=arg
686        if (lm_size>=9) then
687          gylmgr2(1,izero,9)= arg
688          gylmgr2(2,izero,9)=-arg
689        end if
690      end if
691    end if
692 
693  end if
694 
695 !Memory deallocation
696 !==========================================================
697  LIBPAW_DEALLOCATE(rnrm)
698  if (allocated(cc)) then
699    LIBPAW_DEALLOCATE(cc)
700  end if
701  if (compute_gr0)  then
702    LIBPAW_DEALLOCATE(gfact)
703  end if
704  if (compute_gr1)  then
705    LIBPAW_DEALLOCATE(dgfact)
706  end if
707  if (compute_gr2)  then
708    LIBPAW_DEALLOCATE(d2gfact)
709  end if
710  if (compute_gr1)  then
711    LIBPAW_DEALLOCATE(rnrm_inv)
712  end if
713  if (shape_type==3)  then
714    LIBPAW_DEALLOCATE(alpha)
715    LIBPAW_DEALLOCATE(qq)
716  end if
717  if (compute_gr0)  then
718    LIBPAW_DEALLOCATE(ylmr)
719  end if
720  if (compute_gr1)  then
721    LIBPAW_DEALLOCATE(ylmrgr)
722  end if
723  if (shape_type==-1) then
724    if (compute_gr0)  then
725      LIBPAW_DEALLOCATE(shpfuncnum)
726    end if
727    if (compute_gr1)  then
728      LIBPAW_DEALLOCATE(dshpfuncnum)
729    end if
730    if (compute_gr2)  then
731      LIBPAW_DEALLOCATE(d2shpfuncnum)
732    end if
733  end if
734 
735 ! -----------------------------------------------------------------
736 !Small functions related to analytical expression of shape function
737  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

PARENTS

CHILDREN

      paw_uniform_splfit

SOURCE

1109 subroutine pawgylmg(gprimd,gylmg,kg,kpg,kpt,lmax,nkpg,npw,ntypat,pawtab,ylm)
1110 
1111 
1112 !This section has been created automatically by the script Abilint (TD).
1113 !Do not modify the following lines by hand.
1114 #undef ABI_FUNC
1115 #define ABI_FUNC 'pawgylmg'
1116 !End of the abilint section
1117 
1118  implicit none
1119 
1120 !Arguments ------------------------------------
1121 !scalars
1122  integer,intent(in) :: lmax,nkpg,npw,ntypat
1123 !arrays
1124  integer,intent(in) :: kg(3,npw)
1125  real(dp),intent(in) :: gprimd(3,3),kpg(npw,nkpg),kpt(3)
1126  real(dp),intent(in) :: ylm(npw,lmax**2)
1127  type(pawtab_type),intent(in) :: pawtab(ntypat)
1128 
1129  real(dp),intent(out) :: gylmg(npw,lmax**2,ntypat)
1130 
1131 !Local variables-------------------------------
1132 !scalars
1133  integer :: ig,ilm,itypat,ll,l0,mm,mqgrid
1134  real(dp) :: kpg1,kpg2,kpg3,kpgc1,kpgc2,kpgc3
1135 !arrays
1136  real(dp),allocatable :: glg(:),qgrid(:),kpgnorm(:),shpf(:,:),work(:)
1137 
1138 ! *************************************************************************
1139 
1140 !Get |k+G|:
1141  LIBPAW_ALLOCATE(kpgnorm,(npw))
1142  if (nkpg<3) then
1143    do ig=1,npw
1144      kpg1=kpt(1)+dble(kg(1,ig));kpg2=kpt(2)+dble(kg(2,ig));kpg3=kpt(3)+dble(kg(3,ig))
1145      kpgc1=kpg1*gprimd(1,1)+kpg2*gprimd(1,2)+kpg3*gprimd(1,3)
1146      kpgc2=kpg1*gprimd(2,1)+kpg2*gprimd(2,2)+kpg3*gprimd(2,3)
1147      kpgc3=kpg1*gprimd(3,1)+kpg2*gprimd(3,2)+kpg3*gprimd(3,3)
1148      kpgnorm(ig)=sqrt(kpgc1*kpgc1+kpgc2*kpgc2+kpgc3*kpgc3)
1149    end do
1150  else
1151    do ig=1,npw
1152      kpgc1=kpg(ig,1)*gprimd(1,1)+kpg(ig,2)*gprimd(1,2)+kpg(ig,3)*gprimd(1,3)
1153      kpgc2=kpg(ig,1)*gprimd(2,1)+kpg(ig,2)*gprimd(2,2)+kpg(ig,3)*gprimd(2,3)
1154      kpgc3=kpg(ig,1)*gprimd(3,1)+kpg(ig,2)*gprimd(3,2)+kpg(ig,3)*gprimd(3,3)
1155      kpgnorm(ig)=sqrt(kpgc1*kpgc1+kpgc2*kpgc2+kpgc3*kpgc3)
1156    end do
1157  end if
1158 
1159  LIBPAW_ALLOCATE(glg,(npw))
1160  LIBPAW_ALLOCATE(work,(npw))
1161 
1162 !Loop over types of atoms
1163  do itypat=1,ntypat
1164 
1165    mqgrid=pawtab(itypat)%mqgrid_shp
1166    LIBPAW_ALLOCATE(qgrid,(mqgrid))
1167    LIBPAW_ALLOCATE(shpf,(mqgrid,2))
1168    qgrid(1:mqgrid)=pawtab(itypat)%qgrid_shp(1:mqgrid)
1169 
1170 !  Loops over (l,m) values
1171    do ll=0,pawtab(itypat)%lcut_size-1
1172      l0=ll**2+ll+1
1173 
1174      shpf(1:mqgrid,1:2)=pawtab(itypat)%shapefncg(1:mqgrid,1:2,1+ll)
1175      call paw_uniform_splfit(qgrid,work,shpf,0,kpgnorm,glg,mqgrid,npw)
1176 
1177      do mm=-ll,ll
1178        ilm=l0+mm
1179 
1180        gylmg(1:npw,ilm,itypat)=ylm(1:npw,ilm)*glg(1:npw)
1181 
1182 !      End loops over (l,m) values
1183      end do
1184    end do
1185 
1186 !  End loop over atom types
1187    LIBPAW_DEALLOCATE(qgrid)
1188    LIBPAW_DEALLOCATE(shpf)
1189  end do
1190 
1191  LIBPAW_DEALLOCATE(kpgnorm)
1192  LIBPAW_DEALLOCATE(glg)
1193  LIBPAW_DEALLOCATE(work)
1194 
1195 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.

PARENTS

      nhatgrid,pawgrnl

CHILDREN

SOURCE

1232 subroutine pawrfgd_fft(ifftsph,gmet,n1,n2,n3,nfgd,rcut,rfgd,rprimd,ucvol,xred, &
1233 &                      fft_distrib,fft_index,me_fft) ! optional arguments
1234 
1235 
1236 !This section has been created automatically by the script Abilint (TD).
1237 !Do not modify the following lines by hand.
1238 #undef ABI_FUNC
1239 #define ABI_FUNC 'pawrfgd_fft'
1240 !End of the abilint section
1241 
1242  implicit none
1243 
1244 !Arguments ---------------------------------------------
1245 
1246 !scalars
1247  integer,intent(in) :: n1,n2,n3
1248  integer,intent(out) :: nfgd
1249  integer,optional,intent(in) :: me_fft
1250  real(dp),intent(in) :: rcut,ucvol
1251 !arrays
1252  integer,target,optional,intent(in) :: fft_distrib(n3),fft_index(n3)
1253  integer,allocatable,intent(out) :: ifftsph(:)
1254  real(dp),intent(in) :: gmet(3,3),rprimd(3,3),xred(3)
1255  real(dp),allocatable,intent(out) :: rfgd(:,:)
1256 
1257 !Local variables ------------------------------
1258 !scalars
1259  integer,parameter :: ishift=15
1260  integer :: i1,i2,i3,ifft_local,ix,iy,iz,izloc,me_fft_,n1a,n1b,n2a,n2b,n3a,n3b,ncmax
1261  real(dp),parameter :: delta=0.99_dp
1262  real(dp) :: difx,dify,difz,rr1,rr2,rr3,r2,r2cut,rx,ry,rz
1263  character(len=500) :: msg
1264 !arrays
1265  integer,allocatable :: ifftsph_tmp(:)
1266  integer,pointer :: fft_distrib_(:),fft_index_(:)
1267  real(dp),allocatable :: rfgd_tmp(:,:)
1268 
1269 ! *************************************************************************
1270 
1271 !Define a "box" around the atom
1272  r2cut=1.0000001_dp*rcut**2
1273  rr1=sqrt(r2cut*gmet(1,1))
1274  rr2=sqrt(r2cut*gmet(2,2))
1275  rr3=sqrt(r2cut*gmet(3,3))
1276  n1a=int((xred(1)-rr1+ishift)*n1+delta)-ishift*n1
1277  n1b=int((xred(1)+rr1+ishift)*n1      )-ishift*n1
1278  n2a=int((xred(2)-rr2+ishift)*n2+delta)-ishift*n2
1279  n2b=int((xred(2)+rr2+ishift)*n2      )-ishift*n2
1280  n3a=int((xred(3)-rr3+ishift)*n3+delta)-ishift*n3
1281  n3b=int((xred(3)+rr3+ishift)*n3      )-ishift*n3
1282 
1283 !Get the distrib associated with this fft_grid
1284  if (present(fft_distrib).and.present(fft_index).and.present(me_fft)) then
1285    me_fft_=me_fft ; fft_distrib_ => fft_distrib ; fft_index_ => fft_index
1286  else
1287    me_fft_=0
1288    LIBPAW_POINTER_ALLOCATE(fft_distrib_,(n3))
1289    LIBPAW_POINTER_ALLOCATE(fft_index_,(n3))
1290    fft_distrib_=0;fft_index_=(/(i3,i3=1,n3)/)
1291  end if
1292 
1293 !Temporary allocate "large" arrays
1294  ncmax=1+int(1.5_dp*(n1*n2*n3)*four_pi/(three*ucvol)*rcut**3)
1295  LIBPAW_ALLOCATE(ifftsph_tmp,(ncmax))
1296  LIBPAW_ALLOCATE(rfgd_tmp,(3,ncmax))
1297 
1298 !Set number of points to zero
1299  nfgd=0
1300 
1301 !Loop over FFT points
1302  do i3=n3a,n3b
1303    iz=mod(i3+ishift*n3,n3)
1304    if (fft_distrib_(iz+1)==me_fft_) then
1305      izloc=fft_index_(iz+1) - 1
1306      difz=dble(i3)/dble(n3)-xred(3)
1307      do i2=n2a,n2b
1308        iy=mod(i2+ishift*n2,n2)
1309        dify=dble(i2)/dble(n2)-xred(2)
1310        do i1=n1a,n1b
1311          ix=mod(i1+ishift*n1,n1)
1312          difx=dble(i1)/dble(n1)-xred(1)
1313 
1314 !        Compute r-R
1315          rx=difx*rprimd(1,1)+dify*rprimd(1,2)+difz*rprimd(1,3)
1316          ry=difx*rprimd(2,1)+dify*rprimd(2,2)+difz*rprimd(2,3)
1317          rz=difx*rprimd(3,1)+dify*rprimd(3,2)+difz*rprimd(3,3)
1318          r2=rx**2+ry**2+rz**2
1319 
1320 !        Select matching points
1321          if (r2 <= r2cut) then
1322            ifft_local=1+ix+n1*(iy+n2*izloc)
1323            if (ifft_local>0) then
1324              nfgd=nfgd+1
1325              if (nfgd>ncmax) then
1326                msg='Number of fft points around atom exceeds max. allowed!'
1327                MSG_BUG(msg)
1328              end if
1329              rfgd_tmp(1,nfgd)=rx
1330              rfgd_tmp(2,nfgd)=ry
1331              rfgd_tmp(3,nfgd)=rz
1332              ifftsph_tmp(nfgd)=ifft_local
1333            end if
1334          end if
1335 
1336 !      End of loops
1337        end do
1338      end do
1339    end if
1340  end do
1341 
1342 !Now fill output arrays
1343  if (allocated(ifftsph)) then
1344    LIBPAW_DEALLOCATE(ifftsph)
1345  end if
1346  if (allocated(rfgd)) then
1347    LIBPAW_DEALLOCATE(rfgd)
1348  end if
1349  LIBPAW_ALLOCATE(ifftsph,(nfgd))
1350  LIBPAW_ALLOCATE(rfgd,(3,nfgd))
1351  ifftsph(1:nfgd)=ifftsph_tmp(1:nfgd)
1352  rfgd(1:3,1:nfgd)=rfgd_tmp(1:3,1:nfgd)
1353 
1354 !Release temporary memory
1355  LIBPAW_DEALLOCATE(ifftsph_tmp)
1356  LIBPAW_DEALLOCATE(rfgd_tmp)
1357  if (.not.present(fft_distrib).or..not.present(fft_index)) then
1358    LIBPAW_POINTER_DEALLOCATE(fft_distrib_)
1359    LIBPAW_POINTER_DEALLOCATE(fft_index_)
1360  end if
1361 
1362 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.

PARENTS

      wvl_nhatgrid

CHILDREN

SOURCE

1399 subroutine pawrfgd_wvl(geocode,hh,ifftsph,i3s,n1,n1i,n2,n2i,n3,n3pi,&
1400 &                      nfgd,rcut,rloc,rfgd,shift,xcart)
1401 
1402 
1403 !This section has been created automatically by the script Abilint (TD).
1404 !Do not modify the following lines by hand.
1405 #undef ABI_FUNC
1406 #define ABI_FUNC 'pawrfgd_wvl'
1407 !End of the abilint section
1408 
1409  implicit none
1410 
1411 !Arguments ---------------------------------------------
1412 
1413 !scalars
1414  integer,intent(in) :: i3s,n1,n1i,n2,n2i,n3,n3pi,shift
1415  integer,intent(out) :: nfgd
1416  real(dp),intent(in) :: rcut,rloc
1417  character(1),intent(in) :: geocode
1418 !arrays
1419  integer,allocatable,intent(out) :: ifftsph(:)
1420  real(dp),intent(in) :: hh(3),xcart(3)
1421  real(dp),allocatable,intent(out) :: rfgd(:,:)
1422 
1423 !Local variables ------------------------------
1424 !scalars
1425  integer :: i1,i2,i3,iex,iey,iez,ind,isx,isy,isz,j1,j2,j3
1426  integer :: nbl1,nbr1,nbl2,nbr2,nbl3,nbr3,ncmax
1427  logical :: gox,goy,goz,perx,pery,perz
1428  real(dp) :: cutoff,r2,r2cut,rx,ry,rz,xx,yy,zz
1429 !arrays
1430  integer,allocatable :: ifftsph_tmp(:)
1431  real(dp),allocatable :: rfgd_tmp(:,:)
1432 
1433 ! *************************************************************************
1434 
1435 !Data for periodicity in the three directions
1436  perx=(geocode/='F')
1437  pery=(geocode=='P')
1438  perz=(geocode/='F')
1439  call my_ext_buffers(perx,nbl1,nbr1)
1440  call my_ext_buffers(pery,nbl2,nbr2)
1441  call my_ext_buffers(perz,nbl3,nbr3)
1442 
1443 !Define a "box" around the atom
1444  cutoff=10.d0*rloc
1445  r2cut=1.0000001_dp*rcut**2
1446  rx=xcart(1)
1447  ry=xcart(2)
1448  rz=xcart(3)
1449  isx=floor((rx-cutoff)/hh(1))
1450  isy=floor((ry-cutoff)/hh(2))
1451  isz=floor((rz-cutoff)/hh(3))
1452  iex=ceiling((rx+cutoff)/hh(1))
1453  iey=ceiling((ry+cutoff)/hh(2))
1454  iez=ceiling((rz+cutoff)/hh(3))
1455 
1456 !Temporary allocate "large" arrays
1457 !  use factor 1+int(1.1*, for safety reasons
1458  ncmax=1
1459  if (n3pi>0) ncmax=1+int((rcut/hh(1)+1.0)*(rcut/hh(2)+1.0)*(rcut/hh(3)+1.0)*four_pi/three)
1460  LIBPAW_ALLOCATE(ifftsph_tmp,(ncmax))
1461  LIBPAW_ALLOCATE(rfgd_tmp,(3,ncmax))
1462 
1463 !Set number of points to zero
1464  nfgd=0
1465 
1466 !Loop over WVL points
1467  do i3=isz,iez
1468    zz=real(i3,kind=8)*hh(3)-rz
1469    call my_ind_positions(perz,i3,n3,j3,goz)
1470    j3=j3+nbl3+1
1471    do i2=isy,iey
1472      yy=real(i2,kind=8)*hh(2)-ry
1473      call my_ind_positions(pery,i2,n2,j2,goy)
1474      do i1=isx,iex
1475        xx=real(i1,kind=8)*hh(1)-rx
1476        call my_ind_positions(perx,i1,n1,j1,gox)
1477        r2=xx**2+yy**2+zz**2
1478        if (j3>=i3s.and.j3<=i3s+n3pi-1.and.goy.and.gox) then
1479 
1480 !        Select matching points
1481          if (r2<=r2cut) then
1482            ind=j1+1+nbl1+(j2+nbl2)*n1i+(j3-i3s)*n1i*n2i
1483            nfgd=nfgd+1
1484            rfgd_tmp(:,nfgd)=[xx,yy,zz]
1485            ifftsph_tmp(nfgd)=shift+ind
1486          end if
1487 
1488 !      End of loops
1489        end if
1490      end do
1491    end do
1492  end do
1493 
1494 !Now fill output arrays
1495  if (allocated(ifftsph)) then
1496    LIBPAW_DEALLOCATE(ifftsph)
1497  end if
1498  if (allocated(rfgd)) then
1499    LIBPAW_DEALLOCATE(rfgd)
1500  end if
1501  LIBPAW_ALLOCATE(ifftsph,(nfgd))
1502  LIBPAW_ALLOCATE(rfgd,(3,nfgd))
1503  ifftsph(1:nfgd)=ifftsph_tmp(1:nfgd)
1504  rfgd(1:3,1:nfgd)=rfgd_tmp(1:3,1:nfgd)
1505 
1506 !Release temporary memory
1507  LIBPAW_DEALLOCATE(ifftsph_tmp)
1508  LIBPAW_DEALLOCATE(rfgd_tmp)
1509 
1510 !*********************************************************************
1511 !Small functions related to boundary conditions
1512  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 ]