TABLE OF CONTENTS
ABINIT/psp10nl [ Functions ]
NAME
psp10nl
FUNCTION
Hartwigsen-Goedecker-Hutter nonlocal pseudopotential (from preprint of 1998). Uses Gaussians for fully nonlocal form, analytic expressions.
COPYRIGHT
Copyright (C) 1998-2018 ABINIT group (DCA, FRD, XG, GMR, PT, SC) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt . For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
INPUTS
hij(0:lmax,3,3)=factor defining strength of (max 3) projectors for each angular momentum channel l among 0, 1, ..., lmax lmax=maximum angular momentum mproj=maximum number of projectors in any channel mpsang= 1+maximum angular momentum for nonlocal pseudopotentials mqgrid=number of grid points for qgrid nproj(1:lmax+1)=number of projectors in any channel qgrid(mqgrid)=array of |G| values rr(0:lmax)=core radius for each 0<l<lmax channel (bohr)
OUTPUT
ekb(mpsang,mproj)=Kleinman-Bylander energies ffspl(mqgrid,2,mpssang,mproj)=Kleinman-Bylander form factor f_l(q) and second derivative from spline fit for each angular momentum and each projectors
PARENTS
psp10in
CHILDREN
spline,zhpev
SOURCE
42 #if defined HAVE_CONFIG_H 43 #include "config.h" 44 #endif 45 46 #include "abi_common.h" 47 48 subroutine psp10nl(ekb,ffspl,hij,lmax,mproj,mpsang,mqgrid,nproj,qgrid,rr) 49 50 use defs_basis 51 use m_errors 52 use m_profiling_abi 53 use m_splines 54 55 !This section has been created automatically by the script Abilint (TD). 56 !Do not modify the following lines by hand. 57 #undef ABI_FUNC 58 #define ABI_FUNC 'psp10nl' 59 !End of the abilint section 60 61 implicit none 62 63 !Arguments ------------------------------------ 64 !scalars 65 integer,intent(in) :: lmax,mproj,mpsang,mqgrid 66 !arrays 67 integer,intent(in) :: nproj(mpsang) 68 real(dp),intent(in) :: hij(0:lmax,3,3),qgrid(mqgrid),rr(0:lmax) 69 real(dp),intent(out) :: ekb(mpsang,mproj),ffspl(mqgrid,2,mpsang,mproj) 70 71 !Local variables------------------------------- 72 !scalars 73 integer :: info,ipack,iproj,iqgrid,jproj,ll,numproj 74 real(dp) :: qmax,rrl 75 character(len=500) :: message 76 character :: jobz,uplo 77 !arrays 78 real(dp) :: ap(2,9),rwork1(9),work1(2,9),ww(3),yp1j(3),ypnj(3) 79 real(dp),allocatable :: ppspl(:,:,:,:),uu(:,:),work(:),zz(:,:,:) 80 81 ! ************************************************************************* 82 83 ABI_ALLOCATE(ppspl,(mqgrid,2,mpsang,mproj)) 84 ABI_ALLOCATE(work,(mqgrid)) 85 86 qmax=qgrid(mqgrid) 87 jobz='v' 88 uplo='u' 89 ekb(:,:)=zero 90 91 lloop: do ll=0,lmax 92 ap(:,:)=zero 93 numproj=nproj(ll+1) 94 95 ! Fill up the matrix in packed storage 96 prjloop: do jproj=1,numproj 97 priloop: do iproj=1,jproj 98 ipack=iproj+(jproj-1)*jproj/2 99 if(mod((jproj-1)*jproj,2)/=0) then 100 MSG_ERROR("odd") 101 end if 102 ap(1,ipack)=hij(ll,iproj,jproj) 103 end do priloop 104 end do prjloop 105 106 if(numproj/=0)then 107 108 ABI_ALLOCATE(uu,(numproj,numproj)) 109 ABI_ALLOCATE(zz,(2,numproj,numproj)) 110 111 if (numproj > 1) then 112 call ZHPEV(jobz,uplo,numproj,ap,ww,zz,numproj,work1,rwork1,info) 113 uu(:,:)=zz(1,:,:) 114 else 115 ww(1)=hij(ll,1,1) 116 uu(1,1)=one 117 end if 118 119 ! Initialization of ekb, and spline fitting 120 121 if (ll==0) then ! s channel 122 123 rrl=rr(0) 124 do iproj=1,numproj 125 ekb(1,iproj)=ww(iproj)*32.d0*(rrl**3)*(pi**2.5d0)/(4.d0*pi)**2 126 if(iproj==1)then 127 do iqgrid=1,mqgrid 128 ppspl(iqgrid,1,1,1)=exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrl)**2) 129 end do 130 yp1j(1)=zero 131 ypnj(1)=-(two_pi*rrl)**2*qmax*exp(-0.5d0*(two_pi*qmax*rrl)**2) 132 else if(iproj==2)then 133 do iqgrid=1,mqgrid 134 ppspl(iqgrid,1,1,2)=2.0d0/sqrt(15.0d0) & 135 & *exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrl)**2) & 136 & *( 3.d0-(two_pi*qgrid(iqgrid)*rrl)**2 ) 137 end do 138 yp1j(2)=zero 139 ypnj(2)=2.0d0/sqrt(15.0d0)*(two_pi*rrl)**2*qmax & 140 & *exp(-0.5d0*(two_pi*qmax*rrl)**2) * (-5.d0+(two_pi*qmax*rrl)**2) 141 else if(iproj==3)then 142 do iqgrid=1,mqgrid 143 ppspl(iqgrid,1,1,3)=(4.0d0/3.0d0)/sqrt(105.0d0)*& 144 & exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrl)**2) * & 145 & (15.0d0-10.0d0*(two_pi*qgrid(iqgrid)*rrl)**2 + & 146 & (two_pi*qgrid(iqgrid)*rrl)**4) 147 end do 148 yp1j(3)=zero 149 ypnj(3)=(4.0d0/3.0d0)/sqrt(105.0d0)*exp(-0.5d0*(two_pi*qmax*rrl)**2) * & 150 & (two_pi*rrl)**2*qmax*(-35.0d0+14d0*(two_pi*qmax*rrl)**2-(two_pi*qmax*rrl)**4) 151 end if 152 call spline(qgrid,ppspl(:,1,1,iproj),mqgrid,& 153 & yp1j(iproj),ypnj(iproj),ppspl(:,2,1,iproj)) 154 end do 155 156 else if (ll==1) then ! p channel 157 158 rrl=rr(1) 159 do iproj=1,numproj 160 ekb(2,iproj)=ww(iproj)*64.d0*(rrl**5)*(pi**2.5d0)/(4.d0*pi)**2 161 if(iproj==1)then 162 do iqgrid=1,mqgrid 163 ppspl(iqgrid,1,2,1)=(1.0d0/sqrt(3.0d0))* & 164 & exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrl)**2) * (two_pi*qgrid(iqgrid)) 165 end do 166 yp1j(1)=two_pi*(1.0d0/sqrt(3.0d0)) 167 ypnj(1)=-two_pi*((two_pi*qmax*rrl)**2-1.d0)*exp(-0.5d0*(two_pi*qmax*rrl)**2)*& 168 & (1.0d0/sqrt(3.0d0)) 169 else if(iproj==2)then 170 do iqgrid=1,mqgrid 171 ppspl(iqgrid,1,2,2)=(2.0d0/sqrt(105.0d0))* & 172 & exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrl)**2) * & 173 & (two_pi*qgrid(iqgrid))*(5.0d0-(two_pi*qgrid(iqgrid)*rrl)**2) 174 end do 175 yp1j(2)=(5.0d0*two_pi)*(2.0d0/sqrt(105.0d0)) 176 ypnj(2)=(2.0d0/sqrt(105.0d0))*two_pi*exp(-0.5d0*(two_pi*qmax*rrl)**2)* & 177 & (-8*(two_pi*qmax*rrl)**2 + (two_pi*qmax*rrl)**4 + 5.0d0) 178 else if(iproj==3)then 179 do iqgrid=1,mqgrid 180 ppspl(iqgrid,1,2,3)=(4.0d0/3.0d0)/sqrt(1155d0)*& 181 & exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrl)**2) * & 182 & (two_pi*qgrid(iqgrid))*& 183 & (35.0d0-14.0d0*(two_pi*qgrid(iqgrid)*rrl)**2+(two_pi*qgrid(iqgrid)*rrl)**4) 184 end do 185 yp1j(3)=(35.0d0*two_pi)*(4.0d0/3.0d0)/sqrt(1155.0d0) 186 ypnj(3)=(4.0d0/3.0d0)/sqrt(1155.0d0)*two_pi*exp(-0.5d0*(two_pi*qmax*rrl)**2)* & 187 & (35.0d0-77.0d0*(two_pi*qmax*rrl)**2+19.0d0*(two_pi*qmax*rrl)**4 - & 188 & (two_pi*qmax*rrl)**6) 189 end if 190 call spline(qgrid,ppspl(:,1,2,iproj),mqgrid,& 191 & yp1j(iproj),ypnj(iproj),ppspl(:,2,2,iproj)) 192 end do 193 194 else if (ll==2) then ! d channel 195 196 ! If there is a third projector. Warning : only two projectors are allowed. 197 if ( numproj>2 ) then 198 write(message, '(3a)' )& 199 & ' only two d-projectors are allowed ',ch10,& 200 & ' Action : check your pseudopotential file.' 201 MSG_ERROR(message) 202 end if 203 204 rrl=rr(2) 205 do iproj=1,numproj 206 ekb(3,iproj)=ww(iproj)*128.d0*(rrl**7)*(pi**2.5d0)/(4.d0*pi)**2 207 if(iproj==1)then 208 do iqgrid=1,mqgrid 209 ppspl(iqgrid,1,3,1)=(1.0d0/sqrt(15.0d0))* & 210 & exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrl)**2) * (two_pi*qgrid(iqgrid))**2 211 end do 212 yp1j(1)=zero 213 ypnj(1)=(1.0d0/sqrt(15.0d0))*(two_pi**2)*& 214 & exp(-0.5d0*(two_pi*qmax*rrl)**2)*qmax*(2d0-(two_pi*qmax*rrl)**2) 215 else if(iproj==2)then 216 do iqgrid=1,mqgrid 217 ppspl(iqgrid,1,3,2)=(2.0d0/3.0d0)/sqrt(105.0d0)* & 218 & exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrl)**2) * & 219 & ((two_pi*qgrid(iqgrid))**2)*(7.0d0-(two_pi*qgrid(iqgrid)*rrl)**2) 220 end do 221 yp1j(2)=zero 222 ypnj(2)=(2.0d0/3.0d0)/sqrt(105.0d0)*exp(-0.5d0*(two_pi*qmax*rrl)**2)* & 223 & qmax*(two_pi**2)*( (two_pi*qmax*rrl)**4 - 11.0d0*(two_pi*qmax*rrl)**2 + 14.0d0) 224 end if 225 call spline(qgrid,ppspl(:,1,3,iproj),mqgrid,& 226 & yp1j(iproj),ypnj(iproj),ppspl(:,2,3,iproj)) 227 end do 228 229 else if (ll==3) then ! f channel 230 231 ! If there is a second projector. Warning : only one projector is allowed. 232 if ( numproj>1 ) then 233 write(message, '(a,a,a)' )& 234 & ' only one f-projector is allowed ',ch10,& 235 & ' Action : check your pseudopotential file.' 236 MSG_ERROR(message) 237 end if 238 239 rrl=rr(3) 240 ekb(4,1)=ww(1)*(256.0d0/105.0d0)*(rrl**9)*(pi**2.5d0)/(4.d0*pi)**2 241 do iqgrid=1,mqgrid 242 ppspl(iqgrid,1,4,1)=(two_pi*qgrid(iqgrid))**3* & 243 & exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrl)**2) 244 end do 245 ! Compute yp1,ypn=derivatives of f(q) at q=0, q=qgrid(mqgrid) 246 yp1j(1)=zero 247 ypnj(1)=(two_pi**3)*qmax**2*exp(-0.5d0*(two_pi*qmax*rrl)**2)*& 248 & (3.0d0-(two_pi*qmax*rrl)**2) 249 ! Fit spline to get second derivatives by spline fit 250 call spline(qgrid,ppspl(:,1,4,1),mqgrid,& 251 & yp1j(1),ypnj(1),ppspl(:,2,4,1)) 252 253 else 254 MSG_ERROR("lmax>3?") 255 end if 256 257 ! Linear combination using the eigenvectors 258 ffspl(:,:,ll+1,:)=zero 259 do jproj=1,numproj 260 do iproj=1,numproj 261 do iqgrid=1,mqgrid 262 ffspl(iqgrid,1:2,ll+1,jproj)=ffspl(iqgrid,1:2,ll+1,jproj) & 263 & +uu(iproj,jproj)*ppspl(iqgrid,1:2,ll+1,iproj) 264 end do 265 end do 266 end do 267 268 ABI_DEALLOCATE(uu) 269 ABI_DEALLOCATE(zz) 270 271 ! End condition on numproj(/=0) 272 end if 273 274 end do lloop 275 276 ABI_DEALLOCATE(ppspl) 277 ABI_DEALLOCATE(work) 278 279 end subroutine psp10nl