TABLE OF CONTENTS
ABINIT/pspnl_operat_rec [ Functions ]
NAME
pspnl_operat_rec
FUNCTION
It calculates the non-local projectors used in recursion for any psp non-local: The nl interaction in recursion is $$exp{-V_{NL}/beta}=\sum_A\sum_{lm} \sum{ij}Y_{lm}(\hat{r-R_A}')f^l_i(r-R_A)D^l_{i,j}Y_{lm}(\hat{r-R_A})f^l_j{r-R_A}$$ where $D^_{i,j}$ is a matrix previously (see pspnl_operat_rec). In this routine the projectors $Y_{lm}(\hat{r-R_A}')f^l_i(r-R_A)$ are calculated. So an array of dimensions rset%nl%projec(nfftrec,lmnmax,nlrec%npsp)
COPYRIGHT
Copyright (C) 2009-2018 ABINIT group (the_author) 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
metrec<metricrec_type>=contains information concerning metric in recursion: grid_step, metric, infinitesimal volume ngfftrec(18)=is the ngfft grid (truncated if different from ngfft) of recursion debug=debug variable
OUTPUT
nlrec<nlpsprec_type>%projec= array containig the projectors on the real grid nlrec<nlpsprec_type>%intlen= integer linear size of the non-local grid
SIDE EFFECTS
nlrec<nlpsprec_type> data set of non-local pseudo for recursion The better Interaction length (Intlen) is also calculated and printed but recursion use intlen=ngfftrec/2
NOTES
PARENTS
m_rec
CHILDREN
gamma_function,initylmr,wrtout
SOURCE
48 #if defined HAVE_CONFIG_H 49 #include "config.h" 50 #endif 51 52 #include "abi_common.h" 53 54 subroutine pspnl_operat_rec(nlrec,metrec,ngfftrec,debug) 55 56 use defs_basis 57 use defs_rectypes 58 use m_profiling_abi 59 60 use m_paw_sphharm, only : initylmr 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 'pspnl_operat_rec' 66 use interfaces_14_hidewrite 67 use interfaces_28_numeric_noabirule 68 !End of the abilint section 69 70 implicit none 71 72 !Arguments ------------------------------------ 73 !scalars 74 logical,intent(in) :: debug 75 type(metricrec_type),intent(in) ::metrec 76 type(nlpsprec_type),intent(inout) :: nlrec 77 !arrays 78 integer,intent(in) :: ngfftrec(18) 79 !Local variables------------------------------- 80 integer :: ii,intlen 81 integer :: iangol,ipsp,iproj 82 integer :: mpsang,jj,kk,rr 83 integer :: nfftrec 84 integer :: ilmn,il,ilm,in,lmnmax 85 real(dp) :: raggio,rloc,denom,step 86 real(dp) :: delta_out,partial,err 87 character(len=500) :: msg 88 real(dp) :: part_sum(3) 89 real(dp) :: ylmr_gr_dm(0,0,0) 90 real(dp),allocatable :: ylmr(:,:),proj_arr(:,:,:) 91 real(dp),allocatable :: radloc(:,:),nrm(:) 92 93 ! ************************************************************************* 94 95 if(debug)then 96 write(msg,'(80a,a,a)') ('=',ii=1,80),ch10,' pspnl_operat_rec : enter ' 97 call wrtout(std_out,msg,'PERS') 98 end if 99 100 !##################################################################### 101 !--CALCULATE THE (SEMI-)MAXIMUM INTERVAL WHERE ALL THE PROJECTORS ARE 102 !DIFFERENT TO ZERO. 103 !--For any pseudo potential: 104 delta_out = zero 105 step = metrec%tr(1)*half !--Scanning step= grid step/2 106 107 108 do ipsp = 1, nlrec%npsp !--Loop on the pseudos 109 110 ! --For any angular moment: 111 do iangol = 0,maxval(nlrec%indlmn(1,:,ipsp)) !--Loop on the angular moment 112 rloc = nlrec%radii(iangol+1,ipsp) !--Local radius 113 114 ! --For any projector 115 do iproj = 1,nlrec%pspinfo(iangol+1,ipsp) 116 ! --Starting point to searching when the projector goes to zero. 117 ! this correspond to twice the radius wher the projector has its maximum 118 raggio = two*sqrt(real(-2+2*iproj+iangol,dp))*rloc 119 ! --Caculate the gamma func at the denominator 120 call gamma_function(real(iangol+2*iproj,dp)-half,denom) 121 ! --Find the zero 122 ! --The following while cycle should be replaced by a bisection 123 ! --method. Bucause this is calculated only 1 time it is not very 124 ! important. 125 err = one 126 ii=0 127 ! tolloc = 1.d0*abs(minval(nlrec%mat_exp_psp_nl(:nlrec%pspinfo(iangol+1,ipsp),:nlrec%pspinfo(iangol+1,ipsp),iangol+1,ipsp))) 128 do while(abs(err)>1.d-2) 129 raggio = raggio + step 130 err = project_prec(raggio,iproj,iangol,rloc)/sqrt(denom) 131 ii = ii+1 132 end do 133 write(std_out,*)'local delta',raggio,ii 134 delta_out=maxval((/ delta_out,raggio /)) 135 end do !end loop on projectors 136 137 end do !enddo on angular moment 138 end do !enddo on pseudos 139 140 !--CALCULATE how many grid steps correspond to delta_out 141 intlen = int(delta_out/metrec%tr(1)) 142 !--I want that intlen is odd 143 if(mod(intlen,2)==0) intlen = intlen+1 144 145 write(msg,'(2a,i3,a)') ch10,' Interac. length of non-local psp(grid steps)=',intlen,ch10 146 call wrtout(std_out,msg,'COLL') 147 !##################################################################### 148 149 !--Initialisation 150 nfftrec = product(ngfftrec(1:3)) 151 lmnmax = nlrec%lmnmax 152 intlen = ngfftrec(1)/2 153 nlrec%intlen = intlen !--Setted in recursion variables 154 155 !##################################################################### 156 !--CALCULATE E(q,q') 157 !--Cration of the exponential*projectors*ylm matrix 158 159 !--Initialisation 160 ABI_ALLOCATE(nlrec%projec,(nfftrec,lmnmax,nlrec%npsp)) 161 nlrec%projec = zero 162 ABI_ALLOCATE(radloc,(3,nfftrec)) 163 radloc = zero 164 ABI_ALLOCATE(nrm,(nfftrec)) 165 nrm = zero 166 167 !--Loop on pseudo types 168 pseudodo: do ipsp = 1, nlrec%npsp 169 ! --Control if the psp is non-local, else continue 170 if(all(nlrec%pspinfo(:,ipsp)==0)) cycle 171 ! --Vector which stores localy the upper part of symmetrical 172 ! matrix of the exponential of the non-local operator 173 mpsang = maxval(nlrec%indlmn(1,:,ipsp))+1 174 ABI_ALLOCATE(proj_arr,(nfftrec,maxval(nlrec%pspinfo(:,ipsp)),mpsang)) 175 ABI_ALLOCATE(ylmr,(mpsang*mpsang,nfftrec)) 176 proj_arr = zero 177 ylmr = zero 178 179 ! !debug 180 ! write(std_out,*)'mpsang,proj num',mpsang,maxval(nlrec%pspinfo(:,ipsp)) 181 ! !enddebug 182 183 ! --Calculate the projctors 184 do iangol = 0,mpsang-1 185 rloc = nlrec%radii(iangol+1,ipsp) 186 do iproj = 1,nlrec%pspinfo(iangol+1,ipsp) 187 call gamma_function(real(iangol+2*iproj,dp)-half,denom) 188 denom = one/sqrt(denom) 189 do ii = 0,ngfftrec(1)-1 !--3-loop on coordinates 190 do jj = 0,ngfftrec(2)-1 191 do kk = 0,ngfftrec(3)-1 192 ! --Calculate the radii 193 part_sum(:) = real((/ ii,jj,kk /)-intlen,dp)*(metrec%tr) 194 rr = 1+ii+(jj+kk*ngfftrec(2))*ngfftrec(3) 195 radloc(:,rr) = part_sum 196 nrm(rr) = sqrt(sum(part_sum(:)**two)) 197 partial = project_prec(nrm(rr),iproj,iangol,rloc)*denom 198 if(abs(partial)>tol12 ) proj_arr(rr,iproj,iangol+1) = partial 199 end do 200 end do 201 end do !--End 3-loop on coordinates 202 end do 203 end do 204 205 206 ! ------------------------------------------------------------- 207 ! --Calculate the spherical harmonics (Verified: it works well) 208 call initylmr(mpsang,1,nfftrec,nrm(:),1,radloc(:,:),ylmr(:,:),ylmr_gr_dm) 209 ! ------------------------------------------------------------- 210 211 212 do ilmn = 1,lmnmax 213 ilm = nlrec%indlmn(4,ilmn,ipsp) 214 il = nlrec%indlmn(1,ilmn,ipsp)+1 215 in = nlrec%indlmn(3,ilmn,ipsp) 216 write(msg,'(2a,i3,2i2)')ch10,'lm,l,n',ilm,il,in 217 call wrtout(std_out,msg,'COLL') 218 219 nlrec%projec(:,ilmn,ipsp) = ylmr(ilm,:)*proj_arr(:,in,il) 220 end do 221 222 ABI_DEALLOCATE(ylmr) 223 ABI_DEALLOCATE(proj_arr) 224 end do pseudodo !--end loop on pseudo types 225 226 227 ABI_DEALLOCATE(radloc) 228 ABI_DEALLOCATE(nrm) 229 230 if(debug)then 231 write(msg,'(80a,a,a)') ('=',ii=1,80),ch10,' pspnl_operat_rec : exit ' 232 call wrtout(std_out,msg,'PERS') 233 end if 234 235 contains 236 237 function project_prec(raggio,iproj,iangol,rloc) 238 !--Analytical expression of the projectors in hgh-pspeudopotential 239 !--The gamma function at denominator is missing 240 241 !This section has been created automatically by the script Abilint (TD). 242 !Do not modify the following lines by hand. 243 #undef ABI_FUNC 244 #define ABI_FUNC 'project_prec' 245 !End of the abilint section 246 247 real(dp) :: project_prec 248 integer,intent(in) :: iproj,iangol 249 real(dp),intent(in) :: raggio,rloc 250 251 project_prec=sqrt2*(raggio/rloc)**real((iangol+2*(iproj-1)),dp)*& 252 & exp(-((raggio/rloc)**two)*half)/rloc**onehalf 253 end function project_prec 254 255 end subroutine pspnl_operat_rec