TABLE OF CONTENTS
ABINIT/psp3nl [ Functions ]
NAME
psp3nl
FUNCTION
Hartwigsen-Goedecker-Hutter nonlocal pseudopotential (from preprint of 1998). Uses Gaussians for fully nonlocal form, analytic expressions.
COPYRIGHT
Copyright (C) 1998-2017 ABINIT group (DCA, FRD, XG, GMR, PT) 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
h11s=factor defining strength of 1st projector for l=0 channel h22s=factor defining strength of 2nd projector for l=0 channel h33s=factor defining strength of 3rd projector for l=0 channel h11p=factor defining strength of 1st projector for l=1 channel h22p=factor defining strength of 2nd projector for l=1 channel h33p=factor defining strength of 2nd projector for l=1 channel h11d=factor defining strength of 1st projector for l=2 channel h22d=factor defining strength of 2nd projector for l=2 channel h33d=factor defining strength of 2nd projector for l=2 channel h11f=factor defining strength of 1st projector for l=3 channel mproj=maximum number of projectors in any channel mpsang= 1+maximum angular momentum for nonlocal pseudopotentials mqgrid=number of grid points for qgrid qgrid(mqgrid)=array of |G| values rrd=core radius for d channel (bohr) rrf=core radius for f channel (bohr) rrp=core radius for p channel (bohr) rrs=core radius for s 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
psp3in
CHILDREN
spline,zhpev
SOURCE
51 #if defined HAVE_CONFIG_H 52 #include "config.h" 53 #endif 54 55 #include "abi_common.h" 56 57 subroutine psp3nl(ekb,ffspl,h11s,h22s,h33s,h11p,h22p,h33p,h11d,h22d,& 58 & h33d,h11f,mproj,mpsang,mqgrid,qgrid,rrd,rrf,rrp,rrs) 59 60 use defs_basis 61 use m_splines 62 use m_errors 63 use m_profiling_abi 64 65 !This section has been created automatically by the script Abilint (TD). 66 !Do not modify the following lines by hand. 67 #undef ABI_FUNC 68 #define ABI_FUNC 'psp3nl' 69 !End of the abilint section 70 71 implicit none 72 73 !Arguments ------------------------------------ 74 !scalars 75 integer,intent(in) :: mproj,mpsang,mqgrid 76 real(dp),intent(in) :: h11d,h11f,h11p,h11s,h22d,h22p,h22s,h33d,h33p,h33s,rrd 77 real(dp),intent(in) :: rrf,rrp,rrs 78 !arrays 79 real(dp),intent(in) :: qgrid(mqgrid) 80 real(dp),intent(out) :: ekb(mpsang,mproj),ffspl(mqgrid,2,mpsang,mproj) 81 82 !Local variables------------------------------- 83 !scalars 84 integer :: info,iproj,iqgrid,ldz,mu,nproj,nu 85 real(dp) :: qmax 86 character(len=500) :: message 87 character :: jobz,uplo 88 !arrays 89 real(dp) :: ap(2,9),rwork1(9),work1(2,9),ww(3),yp1j(3),ypnj(3) 90 real(dp),allocatable :: ppspl(:,:,:,:),uu(:,:),work(:),zz(:,:,:) 91 92 ! ************************************************************************* 93 94 !DEBUG 95 !write(std_out,*)' psp3nl : enter ' 96 !stop 97 !ENDDEBUG 98 99 ABI_ALLOCATE(ppspl,(mqgrid,2,mpsang,mproj)) 100 ABI_ALLOCATE(work,(mqgrid)) 101 102 qmax=qgrid(mqgrid) 103 jobz='v' 104 uplo='u' 105 ekb(:,:)=0.0d0 106 107 !--------------------------------------------------------------- 108 !Treat s channel 109 110 nproj=0 111 ap(:,:)=0.0d0 112 !If there is at least one s-projector 113 if ( abs(h11s) >= 1.0d-8 ) then 114 nproj=1 ; ldz=1 ; ap(1,1)=h11s 115 end if 116 nproj=1 117 !If there is a second projector 118 if ( abs(h22s) >= 1.0d-8 ) then 119 nproj=2 ; ldz=2 ; ap(1,3)=h22s 120 ap(1,2)=-0.5d0*sqrt(0.6d0)*h22s 121 end if 122 !If there is a third projector 123 if ( abs(h33s) >= 1.0d-8 ) then 124 nproj=3 ; ldz=3 ; ap(1,6)=h33s 125 ap(1,4)=0.5d0*sqrt(5.d0/21.d0)*h33s 126 ap(1,5)=-0.5d0*sqrt(100.d0/63.d0)*h33s 127 end if 128 129 if(nproj/=0)then 130 131 ABI_ALLOCATE(uu,(nproj,nproj)) 132 ABI_ALLOCATE(zz,(2,nproj,nproj)) 133 134 if (nproj > 1) then 135 call ZHPEV(jobz,uplo,nproj,ap,ww,zz,ldz,work1,rwork1,info) 136 uu(:,:)=zz(1,:,:) 137 else 138 ww(1)=h11s 139 uu(1,1)=1.0d0 140 end if 141 142 ! Initialization of ekb, and spline fitting 143 do iproj=1,nproj 144 ekb(1,iproj)=ww(iproj)*32.d0*(rrs**3)*(pi**2.5d0)/(4.d0*pi)**2 145 if(iproj==1)then 146 do iqgrid=1,mqgrid 147 ppspl(iqgrid,1,1,1)=exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrs)**2) 148 end do 149 yp1j(1)=0.d0 150 ypnj(1)=-(two_pi*rrs)**2*qmax*exp(-0.5d0*(two_pi*qmax*rrs)**2) 151 else if(iproj==2)then 152 do iqgrid=1,mqgrid 153 ppspl(iqgrid,1,1,2)=2.0d0/sqrt(15.0d0) & 154 & *exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrs)**2) & 155 & *( 3.d0-(two_pi*qgrid(iqgrid)*rrs)**2 ) 156 end do 157 yp1j(2)=0.0d0 158 ypnj(2)=2.0d0/sqrt(15.0d0)*(two_pi*rrs)**2*qmax & 159 & *exp(-0.5d0*(two_pi*qmax*rrs)**2) * (-5.d0+(two_pi*qmax*rrs)**2) 160 else if(iproj==3)then 161 do iqgrid=1,mqgrid 162 ppspl(iqgrid,1,1,3)=(4.0d0/3.0d0)/sqrt(105.0d0)*& 163 & exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrs)**2) * & 164 & (15.0d0-10.0d0*(two_pi*qgrid(iqgrid)*rrs)**2 + & 165 & (two_pi*qgrid(iqgrid)*rrs)**4) 166 end do 167 yp1j(3)=0.0d0 168 ypnj(3)=(4.0d0/3.0d0)/sqrt(105.0d0)*exp(-0.5d0*(two_pi*qmax*rrs)**2) * & 169 & (two_pi*rrs)**2*qmax*(-35.0d0+14d0*(two_pi*qmax*rrs)**2-(two_pi*qmax*rrs)**4) 170 end if 171 call spline(qgrid,ppspl(:,1,1,iproj),mqgrid,& 172 & yp1j(iproj),ypnj(iproj),ppspl(:,2,1,iproj)) 173 end do 174 175 ! Linear combination using the eigenvectors 176 ffspl(:,:,1,:)=0.0d0 177 do mu=1,nproj 178 do nu=1,nproj 179 do iqgrid=1,mqgrid 180 ffspl(iqgrid,1:2,1,mu)=ffspl(iqgrid,1:2,1,mu) & 181 & +uu(nu,mu)*ppspl(iqgrid,1:2,1,nu) 182 end do 183 end do 184 end do 185 186 ABI_DEALLOCATE(uu) 187 ABI_DEALLOCATE(zz) 188 189 ! End condition on nproj(/=0) 190 end if 191 192 !DEBUG 193 !write(std_out,*)' psp3nl : after s channel ' 194 !stop 195 !ENDDEBUG 196 197 !-------------------------------------------------------------------- 198 !Now treat p channel 199 200 nproj=0 201 ap(:,:)=0.0d0 202 !If there is at least one projector 203 if ( abs(h11p) >= 1.0d-8 ) then 204 nproj=1 ; ldz=1 ; ap(1,1)=h11p 205 end if 206 !If there is a second projector 207 if ( abs(h22p) >= 1.0d-8 ) then 208 nproj=2 ; ldz=2 ; ap(1,3)=h22p 209 ap(1,2)=-0.5d0*sqrt(5.d0/7.d0)*h22p 210 end if 211 !If there is a third projector 212 if ( abs(h33p) >= 1.0d-8 ) then 213 nproj=3 ; ldz=3 ; ap(1,6)=h33p 214 ap(1,4)= (1.d0/6.d0)*sqrt(35.d0/11.d0)*h33p 215 ap(1,5)=-(1.d0/6.d0)*(14.d0/sqrt(11.d0))*h33p 216 end if 217 218 if(nproj/=0)then 219 220 ABI_ALLOCATE(uu,(nproj,nproj)) 221 ABI_ALLOCATE(zz,(2,nproj,nproj)) 222 223 if (nproj > 1) then 224 call ZHPEV(jobz,uplo,nproj,ap,ww,zz,ldz,work1,rwork1,info) 225 uu(:,:)=zz(1,:,:) 226 else 227 ww(1)=h11p 228 uu(1,1)=1.0d0 229 end if 230 231 ! Initialization of ekb, and spline fitting 232 do iproj=1,nproj 233 ekb(2,iproj)=ww(iproj)*64.d0*(rrp**5)*(pi**2.5d0)/(4.d0*pi)**2 234 if(iproj==1)then 235 do iqgrid=1,mqgrid 236 ppspl(iqgrid,1,2,1)=(1.0d0/sqrt(3.0d0))* & 237 & exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrp)**2) * (two_pi*qgrid(iqgrid)) 238 end do 239 yp1j(1)=two_pi*(1.0d0/sqrt(3.0d0)) 240 ypnj(1)=-two_pi*((two_pi*qmax*rrp)**2-1.d0)*exp(-0.5d0*(two_pi*qmax*rrp)**2)*& 241 & (1.0d0/sqrt(3.0d0)) 242 else if(iproj==2)then 243 do iqgrid=1,mqgrid 244 ppspl(iqgrid,1,2,2)=(2.0d0/sqrt(105.0d0))* & 245 & exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrp)**2) * & 246 & (two_pi*qgrid(iqgrid))*(5.0d0-(two_pi*qgrid(iqgrid)*rrp)**2) 247 end do 248 yp1j(2)=(5.0d0*two_pi)*(2.0d0/sqrt(105.0d0)) 249 ypnj(2)=(2.0d0/sqrt(105.0d0))*two_pi*exp(-0.5d0*(two_pi*qmax*rrp)**2)* & 250 & (-8*(two_pi*qmax*rrp)**2 + (two_pi*qmax*rrp)**4 + 5.0d0) 251 else if(iproj==3)then 252 do iqgrid=1,mqgrid 253 ppspl(iqgrid,1,2,3)=(4.0d0/3.0d0)/sqrt(1155d0)*& 254 & exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrp)**2) * & 255 & (two_pi*qgrid(iqgrid))*& 256 & (35.0d0-14.0d0*(two_pi*qgrid(iqgrid)*rrp)**2+(two_pi*qgrid(iqgrid)*rrp)**4) 257 end do 258 yp1j(3)=(35.0d0*two_pi)*(4.0d0/3.0d0)/sqrt(1155.0d0) 259 ypnj(3)=(4.0d0/3.0d0)/sqrt(1155.0d0)*two_pi*exp(-0.5d0*(two_pi*qmax*rrp)**2)* & 260 & (35.0d0-77.0d0*(two_pi*qmax*rrp)**2+19.0d0*(two_pi*qmax*rrp)**4 - & 261 & (two_pi*qmax*rrp)**6) 262 end if 263 call spline(qgrid,ppspl(:,1,2,iproj),mqgrid,& 264 & yp1j(iproj),ypnj(iproj),ppspl(:,2,2,iproj)) 265 end do 266 267 ! Linear combination using the eigenvectors 268 ffspl(:,:,2,:)=0.0d0 269 do mu=1,nproj 270 do nu=1,nproj 271 do iqgrid=1,mqgrid 272 ffspl(iqgrid,1:2,2,mu)=ffspl(iqgrid,1:2,2,mu) & 273 & +uu(nu,mu)*ppspl(iqgrid,1:2,2,nu) 274 end do 275 end do 276 end do 277 278 ABI_DEALLOCATE(uu) 279 ABI_DEALLOCATE(zz) 280 281 ! End condition on nproj(/=0) 282 end if 283 284 !DEBUG 285 !write(std_out,*)' psp3nl : after p channel ' 286 !stop 287 !ENDDEBUG 288 289 !----------------------------------------------------------------------- 290 !Now treat d channel. 291 292 nproj=0 293 ap(:,:)=0.0d0 294 !If there is at least one projector 295 if ( abs(h11d) >= 1.0d-8 ) then 296 nproj=1 ; ldz=1 ; ap(1,1)=h11d 297 end if 298 !If there is a second projector 299 if ( abs(h22d) >= 1.0d-8 ) then 300 nproj=2 ; ldz=2 ; ap(1,3)=h22d 301 ap(1,2)=-0.5d0*sqrt(7.d0/9.d0)*h22d 302 end if 303 !If there is a third projector. Warning : only two projectors are allowed. 304 if ( abs(h33d) >= 1.0d-8 ) then 305 write(message, '(a,a,a)' )& 306 & ' only two d-projectors are allowed ',ch10,& 307 & ' Action : check your pseudopotential file.' 308 MSG_ERROR(message) 309 ! nproj=3 ; ldz=3 ; ap(1,6)=h33d 310 ! ap(1,4)= 0.5d0*sqrt(63.d0/143.d0)*h33d 311 ! ap(1,5)= -0.5d0*(18.d0/sqrt(143.d0))*h33d 312 end if 313 314 if(nproj/=0)then 315 316 ABI_ALLOCATE(uu,(nproj,nproj)) 317 ABI_ALLOCATE(zz,(2,nproj,nproj)) 318 319 if (nproj > 1) then 320 call ZHPEV(jobz,uplo,nproj,ap,ww,zz,ldz,work1,rwork1,info) 321 uu(:,:)=zz(1,:,:) 322 else 323 ww(1)=h11d 324 uu(1,1)=1.0d0 325 end if 326 327 ! Initialization of ekb, and spline fitting 328 do iproj=1,nproj 329 ekb(3,iproj)=ww(iproj)*128.d0*(rrd**7)*(pi**2.5d0)/(4.d0*pi)**2 330 if(iproj==1)then 331 do iqgrid=1,mqgrid 332 ppspl(iqgrid,1,3,1)=(1.0d0/sqrt(15.0d0))* & 333 & exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrd)**2) * (two_pi*qgrid(iqgrid))**2 334 end do 335 yp1j(1)=0.0d0 336 ypnj(1)=(1.0d0/sqrt(15.0d0))*(two_pi**2)*& 337 & exp(-0.5d0*(two_pi*qmax*rrd)**2)*qmax*(2d0-(two_pi*qmax*rrd)**2) 338 else if(iproj==2)then 339 do iqgrid=1,mqgrid 340 ppspl(iqgrid,1,3,2)=(2.0d0/3.0d0)/sqrt(105.0d0)* & 341 & exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrd)**2) * & 342 & ((two_pi*qgrid(iqgrid))**2)*(7.0d0-(two_pi*qgrid(iqgrid)*rrd)**2) 343 end do 344 yp1j(2)=0.0d0 345 ypnj(2)=(2.0d0/3.0d0)/sqrt(105.0d0)*exp(-0.5d0*(two_pi*qmax*rrd)**2)* & 346 & qmax*(two_pi**2)*( (two_pi*qmax*rrd)**4 - 11.0d0*(two_pi*qmax*rrd)**2 + 14.0d0) 347 end if 348 call spline(qgrid,ppspl(:,1,3,iproj),mqgrid,& 349 & yp1j(iproj),ypnj(iproj),ppspl(:,2,3,iproj)) 350 end do 351 352 ! Linear combination using the eigenvectors 353 ffspl(:,:,3,:)=0.0d0 354 do mu=1,nproj 355 do nu=1,nproj 356 do iqgrid=1,mqgrid 357 ffspl(iqgrid,1:2,3,mu)=ffspl(iqgrid,1:2,3,mu) & 358 & +uu(nu,mu)*ppspl(iqgrid,1:2,3,nu) 359 end do 360 end do 361 end do 362 363 ABI_DEALLOCATE(uu) 364 ABI_DEALLOCATE(zz) 365 366 ! End condition on nproj(/=0) 367 end if 368 369 !DEBUG 370 !write(std_out,*)' psp3nl : after d channel ' 371 !stop 372 !ENDDEBUG 373 374 !----------------------------------------------------------------------- 375 !Treat now f channel (max one projector ! - so do not use ppspl) 376 377 !l=3 first projector 378 if (abs(h11f)>1.d-12) then 379 ekb(4,1)=h11f*(256.0d0/105.0d0)*(rrf**9)*(pi**2.5d0)/(4.d0*pi)**2 380 do iqgrid=1,mqgrid 381 ffspl(iqgrid,1,4,1)=(two_pi*qgrid(iqgrid))**3* & 382 & exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrf)**2) 383 end do 384 ! Compute yp1,ypn=derivatives of f(q) at q=0, q=qgrid(mqgrid) 385 yp1j(1)=0d0 386 ypnj(1)=(two_pi**3)*qmax**2*exp(-0.5d0*(two_pi*qmax*rrf)**2)*& 387 & (3.0d0-(two_pi*qmax*rrf)**2) 388 ! Fit spline to get second derivatives by spline fit 389 call spline(qgrid,ffspl(:,1,4,1),mqgrid,& 390 & yp1j(1),ypnj(1),ffspl(:,2,4,1)) 391 end if 392 393 !----------------------------------------------------------------------- 394 395 ABI_DEALLOCATE(ppspl) 396 ABI_DEALLOCATE(work) 397 398 !DEBUG 399 !write(std_out,*)' psp3nl : exit ' 400 !stop 401 !ENDDEBUG 402 403 end subroutine psp3nl