TABLE OF CONTENTS


ABINIT/pawnabla_init [ Functions ]

[ Top ] [ Functions ]

NAME

 pawnabla_init

FUNCTION

 Evaluate onsite contributions of the nabla operator in cartesian coordinates.
 Store values in the pawtab% data structure.

COPYRIGHT

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

INPUTS

  mpsang=1+maximum angular momentum
  ntypat=Number of types of atoms in cell
  Pawrad(ntypat)<Pawrad_type>=PAW radial mesh and related data:
    %rad(mesh_size)=The coordinates of all the points of the radial mesh
  Pawtab(ntypat) <type(pawtab_type>=PAW tabulated starting data:
    %mesh_size=Dimension of radial mesh
    %lmn_size=Number of (l,m,n) elements for the PAW basis

OUTPUT

  See side effects

SIDE EFFECTS

  Pawtab(ntypat) <type(pawtab_type>=PAW tabulated starting data:
    %has_nabla=set to 1 in matrix elements are calculated and stored
    %nabla_ij(3,lmn_size,lmn_size)= <phi_i|nabla|phi_j>-<tphi_i|nabla|tphi_j>

NOTES

  MG extracted this piece of code from optics_paw.F90 in order to have something more 
  reusable! Note however the storage mode of nabla_ij differs from optics_paw 
  (here Cartesian coordinates run faster). Besides nabla_ij contains the matrix 
  elements of \nabla instead of the elements of the momentum operator p.

PARENTS

      bethe_salpeter,screening

CHILDREN

      int_ang,nderiv_gen,pawrad_deducer0,simp_gen

SOURCE

 47 #if defined HAVE_CONFIG_H
 48 #include "config.h"
 49 #endif
 50 
 51 #include "abi_common.h"
 52 
 53 subroutine pawnabla_init(mpsang,ntypat,pawrad,pawtab)
 54     
 55  use defs_basis
 56  use m_errors
 57  use m_profiling_abi
 58 
 59  use m_pawrad, only : pawrad_type, pawrad_deducer0, simp_gen, nderiv_gen
 60  use m_pawtab, only : pawtab_type
 61 
 62 !This section has been created automatically by the script Abilint (TD).
 63 !Do not modify the following lines by hand.
 64 #undef ABI_FUNC
 65 #define ABI_FUNC 'pawnabla_init'
 66  use interfaces_65_paw, except_this_one => pawnabla_init
 67 !End of the abilint section
 68 
 69  implicit none
 70 
 71 !Arguments ------------------------------------
 72 !scalars
 73  integer,intent(in) :: mpsang,ntypat
 74 !arrays
 75  type(pawtab_type),target,intent(inout) :: pawtab(ntypat)
 76  type(pawrad_type),intent(in) :: pawrad(ntypat)
 77 
 78 !Local variables-------------------------------
 79 !scalars
 80  integer :: nln,il,ilm,ilmn,iln,itypat
 81  integer :: jl,jlm,jlmn,jln,lmn_size,mesh_size 
 82  real(dp) :: intg
 83  character(len=500) :: msg      
 84 !arrays
 85  integer,ABI_CONTIGUOUS pointer :: indlmn(:,:)
 86  real(dp) :: ang_phipphj(mpsang**2,mpsang**2,8)
 87  real(dp),allocatable :: dphi(:),dtphi(:),ff(:),int1(:,:),int2(:,:),rad(:)
 88  
 89 ! *************************************************************************
 90 
 91  DBG_ENTER("COLL")
 92 
 93  if (mpsang>4)then
 94    write(msg,'(3a)')&
 95 &   'Not designed for angular momentum greater than 3 ',ch10,&
 96 &   'Modification in the table defined in ang_int.F90 is required. '
 97    MSG_ERROR(msg)
 98  end if
 99 !
100 !Integration of the angular part: all angular integrals have been computed 
101 !outside Abinit and tabulated for each (l,m) value up to l=2
102  call int_ang(ang_phipphj,mpsang)
103 
104  do itypat=1,ntypat
105 !  
106 !  COMPUTE nabla_ij := <phi_i|nabla|phi_j>-<tphi_i|nabla|tphi_j> for this type
107    mesh_size=pawtab(itypat)%mesh_size
108    lmn_size=pawtab(itypat)%lmn_size
109    nln=pawtab(itypat)%basis_size
110 
111    if (allocated(pawtab(itypat)%nabla_ij)) then
112      ABI_DEALLOCATE(pawtab(itypat)%nabla_ij)
113    end if
114    ABI_ALLOCATE(pawtab(itypat)%nabla_ij,(3,lmn_size,lmn_size))
115    pawtab(itypat)%has_nabla=1
116 
117    ABI_ALLOCATE(ff,(mesh_size))
118    ABI_ALLOCATE(rad,(mesh_size))
119    ABI_ALLOCATE(int2,(lmn_size,lmn_size))
120    ABI_ALLOCATE(int1,(lmn_size,lmn_size))
121    ABI_ALLOCATE(dphi,(mesh_size))
122    ABI_ALLOCATE(dtphi,(mesh_size))
123 
124    indlmn => pawtab(itypat)%indlmn
125    rad(1:mesh_size)=pawrad(itypat)%rad(1:mesh_size)
126 !
127 !  int1=\int phi phj/r dr - \int tphi tphj /r dr
128    do jln=1,nln
129      do iln=1,nln
130        ff(2:mesh_size)= ( &
131 &       pawtab(itypat)%phi (2:mesh_size,iln)*pawtab(itypat)%phi (2:mesh_size,jln) &
132 &       -pawtab(itypat)%tphi(2:mesh_size,iln)*pawtab(itypat)%tphi(2:mesh_size,jln) ) /rad(2:mesh_size)
133        call pawrad_deducer0(ff,mesh_size,pawrad(itypat))
134        call simp_gen(intg,ff,pawrad(itypat))
135        int1(iln,jln)=intg
136      end do
137    end do
138 !
139 !  int2=\int phi/r d/dr(phj/r) r^2dr - \int tphi/r d/dr(tphj/r)r^2 dr
140    do jln=1,nln
141      ff(1:mesh_size)=pawtab(itypat)%phi(1:mesh_size,jln)
142      call nderiv_gen(dphi,ff,pawrad(itypat))
143      ff(1:mesh_size)=pawtab(itypat)%tphi(1:mesh_size,jln)
144      call nderiv_gen(dtphi,ff,pawrad(itypat))
145 
146      do iln=1,nln
147        ff(2:mesh_size)= &
148 &       pawtab(itypat)%phi (2:mesh_size,iln)*dphi (2:mesh_size) &
149 &       -pawtab(itypat)%phi (2:mesh_size,iln)*pawtab(itypat)%phi (2:mesh_size,jln)/rad(2:mesh_size) &
150 &       -( pawtab(itypat)%tphi(2:mesh_size,iln)*dtphi(2:mesh_size) &
151 &       -pawtab(itypat)%tphi(2:mesh_size,iln)*pawtab(itypat)%tphi(2:mesh_size,jln)/rad(2:mesh_size) )
152        call pawrad_deducer0(ff,mesh_size,pawrad(itypat))
153        call simp_gen(intg,ff,pawrad(itypat))
154        int2(iln,jln)=intg
155      end do
156    end do
157 !  
158 !  1-c Integration of the radial part, Note unpacked loop
159    do jlmn=1,lmn_size
160      jlm=indlmn(4,jlmn)
161      jl =indlmn(5,jlmn)
162      do ilmn=1,lmn_size
163        ilm=indlmn(4,ilmn)
164        il =indlmn(5,ilmn)
165 
166        pawtab(itypat)%nabla_ij(1,ilmn,jlmn)= &
167 &       int2(il,jl)* ang_phipphj(ilm,jlm,1) &
168 &       +int1(il,jl)*(ang_phipphj(ilm,jlm,2)+ang_phipphj(ilm,jlm,3))
169 
170        pawtab(itypat)%nabla_ij(2,ilmn,jlmn)= &
171 &       int2(il,jl)* ang_phipphj(ilm,jlm,4) &
172 &       +int1(il,jl)*(ang_phipphj(ilm,jlm,5)+ang_phipphj(ilm,jlm,6))
173 
174        pawtab(itypat)%nabla_ij(3,ilmn,jlmn)= &
175 &       int2(il,jl)* ang_phipphj(ilm,jlm,7) &
176 &       +int1(il,jl)* ang_phipphj(ilm,jlm,8)
177 
178      end do !ilmn
179    end do !jlmn
180 
181    pawtab(itypat)%has_nabla=2
182    ABI_DEALLOCATE(ff)
183    ABI_DEALLOCATE(rad)
184    ABI_DEALLOCATE(int2)
185    ABI_DEALLOCATE(int1)
186    ABI_DEALLOCATE(dphi)
187    ABI_DEALLOCATE(dtphi)
188 
189  end do !itypat
190 
191  DBG_EXIT("COLL")
192 
193 end subroutine pawnabla_init