TABLE OF CONTENTS
ABINIT/pawnabla_init [ 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