TABLE OF CONTENTS
ABINIT/psp3in [ Functions ]
NAME
psp3in
FUNCTION
Initialize pspcod=3 pseudopotentials (HGH psps PRB58,3641(1998)): continue to read the file, then compute the corresponding local and non-local potentials.
COPYRIGHT
Copyright (C) 1998-2018 ABINIT group (DCA, XG, GMR, FD, 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
dtset <type(dataset_type)>=all input variables in this dataset pspso=spin-orbit characteristics, govern the content of ffspl and ekb if =0 : this input requires NO spin-orbit characteristics of the psp if =2 : this input requires HGH characteristics of the psp if =3 : this input requires HFN characteristics of the psp ipsp=id in the array of the pseudo-potential. zion=nominal valence of atom as specified in psp file
OUTPUT
ekb(lnmax)=Kleinman-Bylander energy, {{\ \begin{equation} \frac{\int_0^\infty [Rl(r)^2 (Vl(r)-Vloc(r))^2 dr]} {\int_0^\infty [Rl(r)^2 (Vl(r)-Vloc(r)) dr]} \end{equation} }} for each (l,n) if any, spin-orbit components begin at l=mpsang+1 epsatm=$ (4\pi)\int_0^\infty [r^2 (V(r)+\frac{Zv}{r}) dr]$ (hartree) ffspl(mqgrid_ff,2,lnmax)=Kleinman-Bylander form factor f_l(q) and second derivative from spline fit for each angular momentum and each projector; if any, spin-orbit components begin at l=mpsang+1 indlmn(6,i)= array giving l,m,n,lm,ln,s for i=ln (if useylm=0) or i=lmn (if useylm=1) nproj(mpssoang)=number of projection functions for each angular momentum vlspl(mqgrid_ff,2)=q^2 Vloc(q) and second derivatives from spline fit
SIDE EFFECTS
Input/output
lmax : at input =value of lmax mentioned at the second line of the psp file
at output= 1
psps <type(pseudopotential_type)>=at output, values depending on the read
pseudo are set.
| lmnmax(IN)=if useylm=1, max number of (l,m,n) comp. over all type of psps
| =if useylm=0, max number of (l,n) comp. over all type of psps
| lnmax(IN)=max. number of (l,n) components over all type of psps
| angular momentum of nonlocal pseudopotential
| mpsang(IN)= 1+maximum angular momentum for nonlocal pseudopotentials
| mpssoang(IN)= 1+maximum (spin*angular momentum) for nonlocal pseudopotentials
| mqgrid_ff(IN)=dimension of q (or G) grid for arrays.
| qgrid_ff(mqgrid_ff)(IN)=values of q on grid from 0 to qmax (bohr^-1) for nl form factors
| useylm(IN)=governs the way the nonlocal operator is to be applied:
| 1=using Ylm, 0=using Legendre polynomials
PARENTS
pspatm
CHILDREN
psp2lo,psp3nl,spline,wrtout,wvl_descr_psp_fill
SOURCE
69 #if defined HAVE_CONFIG_H 70 #include "config.h" 71 #endif 72 73 #include "abi_common.h" 74 75 subroutine psp3in(dtset, ekb, epsatm, ffspl, indlmn, ipsp, lmax, nproj, psps, pspso, vlspl, zion) 76 77 use defs_basis 78 use defs_datatypes 79 use defs_abitypes 80 use m_profiling_abi 81 use m_splines 82 use m_errors 83 #if defined HAVE_BIGDFT 84 use BigDFT_API, only: atomic_info 85 #endif 86 87 !This section has been created automatically by the script Abilint (TD). 88 !Do not modify the following lines by hand. 89 #undef ABI_FUNC 90 #define ABI_FUNC 'psp3in' 91 use interfaces_14_hidewrite 92 use interfaces_43_wvl_wrappers 93 use interfaces_64_psp, except_this_one => psp3in 94 !End of the abilint section 95 96 implicit none 97 98 !Arguments ------------------------------------ 99 !scalars 100 integer,intent(in) :: ipsp,pspso 101 integer,intent(inout) :: lmax 102 real(dp),intent(in) :: zion 103 real(dp),intent(out) :: epsatm 104 type(dataset_type),intent(in) :: dtset 105 type(pseudopotential_type),intent(inout) :: psps 106 !arrays 107 integer,intent(out) :: indlmn(6,psps%lmnmax),nproj(psps%mpssoang) 108 real(dp),intent(inout) :: ekb(psps%lnmax),ffspl(psps%mqgrid_ff,2,psps%lnmax)!vz_i 109 real(dp),intent(out) :: vlspl(psps%mqgrid_ff,2) 110 111 !Local variables------------------------------- 112 !scalars 113 integer :: iln,iln0,index,ipsang,jj,kk,ll,mm,mproj,nn,nso 114 real(dp) :: cc1,cc2,cc3,cc4,h11d,h11f,h11p,h11s,h22d,h22p,h22s,h33d,h33p,h33s 115 real(dp) :: k11d,k11f,k11p,k22d,k22p,k33d,k33p,rloc 116 real(dp) :: rrd,rrf,rrp,rrs,yp1,ypn 117 character(len=500) :: message,errmsg 118 !arrays 119 real(dp),allocatable :: dvlspl(:),ekb_so(:,:),ekb_sr(:,:),ffspl_so(:,:,:,:) 120 real(dp),allocatable :: ffspl_sr(:,:,:,:),work_space(:),work_spl(:) 121 122 ! *************************************************************************** 123 124 !Set various terms to 0 in case not defined below 125 !HGH values 126 rloc=zero ; rrs=zero ; h11p=zero ; k33p=zero ; k11d=zero; 127 cc1=zero ; h11s=zero ; h22p=zero ; rrd=zero ; k22d=zero; 128 cc2=zero ; h22s=zero ; h33p=zero ; h11d=zero ; k33d=zero; 129 cc3=zero ; h33s=zero ; k11p=zero ; h22d=zero ; h11f=zero; 130 cc4=zero ; rrp=zero ; k22p=zero ; h33d=zero ; k11f=zero; 131 rrf=zero 132 nproj(1:psps%mpssoang)=0 133 134 !DEBUG 135 !write(std_out,*) ' psp3in : enter ' 136 !ENDDEBUG 137 138 !Read and write different lines of the pseudopotential file 139 140 read (tmp_unit,*,err=10,iomsg=errmsg) rloc,cc1,cc2,cc3,cc4 141 write(message, '(a,f12.7)' ) ' rloc=',rloc 142 call wrtout(ab_out,message,'COLL') 143 call wrtout(std_out, message,'COLL') 144 write(message, '(a,f12.7,a,f12.7,a,f12.7,a,f12.7)' )& 145 & ' cc1 =',cc1,'; cc2 =',cc2,'; cc3 =',cc3,'; cc4 =',cc4 146 call wrtout(ab_out,message,'COLL') 147 call wrtout(std_out, message,'COLL') 148 149 !For the time being, the s state line must be present and is read, 150 !even for local pseudopotentials (zero must appear) 151 read (tmp_unit,*,err=10,iomsg=errmsg) rrs,h11s,h22s,h33s 152 write(message, '(a,f12.7,a,f12.7,a,f12.7,a,f12.7)' )& 153 & ' rrs =',rrs,'; h11s=',h11s,'; h22s=',h22s,'; h33s=',h33s 154 call wrtout(ab_out,message,'COLL') 155 call wrtout(std_out, message,'COLL') 156 157 if (lmax > 0) then 158 159 read (tmp_unit,*,err=10,iomsg=errmsg) rrp,h11p,h22p,h33p 160 write(message, '(a,f12.7,a,f12.7,a,f12.7,a,f12.7)' )& 161 & ' rrp =',rrp,'; h11p=',h11p,'; h22p=',h22p,'; h33p=',h33p 162 call wrtout(ab_out,message,'COLL') 163 call wrtout(std_out, message,'COLL') 164 165 read (tmp_unit,*,err=10,iomsg=errmsg) k11p,k22p,k33p 166 write(message, '(a,f12.7,a,f12.7,a,f12.7)' )& 167 & ' k11p=',k11p,'; k22p=',k22p,'; k33p=',k33p 168 call wrtout(ab_out,message,'COLL') 169 call wrtout(std_out, message,'COLL') 170 171 end if 172 173 if (lmax > 1) then 174 175 read (tmp_unit,*,err=10,iomsg=errmsg) rrd,h11d,h22d,h33d 176 write(message, '(a,f12.7,a,f12.7,a,f12.7,a,f12.7)' )& 177 & ' rrd =',rrd,'; h11d=',h11d,'; h22d=',h22d,'; h33d=',h33d 178 call wrtout(ab_out,message,'COLL') 179 call wrtout(std_out, message,'COLL') 180 181 read (tmp_unit,*,err=10,iomsg=errmsg) k11d,k22d,k33d 182 write(message, '(a,f12.7,a,f12.7,a,f12.7)' )& 183 & ' k11d=',k11d,'; k22d=',k22d,'; k33d=',k33d 184 call wrtout(ab_out,message,'COLL') 185 call wrtout(std_out, message,'COLL') 186 187 end if 188 189 if (lmax > 2) then 190 191 read (tmp_unit,*,err=10,iomsg=errmsg) rrf,h11f 192 write(message, '(a,f12.7,a,f12.7)' )' rrf =',rrf,'; h11f=',h11f 193 call wrtout(ab_out,message,'COLL') 194 call wrtout(std_out, message,'COLL') 195 196 read (tmp_unit,*,err=10,iomsg=errmsg) k11f 197 write(message, '(a,f12.7)' )' k11f=',k11f 198 call wrtout(ab_out,message,'COLL') 199 call wrtout(std_out, message,'COLL') 200 201 end if 202 203 if (abs(h11s)>1.d-08) nproj(1)=1 204 if (abs(h22s)>1.d-08) nproj(1)=2 205 if (abs(h33s)>1.d-08) nproj(1)=3 206 207 if (abs(h11p)>1.d-08) then 208 nproj(2)=1 209 if (lmax<1) then 210 write(message, '(a,a,a,a,i5,a,e12.4,a,a)' ) ch10,& 211 & ' psp3in : COMMENT -',ch10,& 212 & ' input lmax=',lmax,' does not agree with input h11p=',h11p,ch10,& 213 & ' setting lmax to 1' 214 call wrtout(ab_out,message,'COLL') 215 call wrtout(std_out, message,'COLL') 216 lmax=1 217 end if 218 end if 219 220 if (abs(h22p)>1.d-08) then 221 nproj(2)=2 222 if (lmax<1) then 223 write(message, '(a,a,a,a,i5,a,e12.4,a,a)' ) ch10,& 224 & ' psp3in : COMMENT -',ch10,& 225 & ' input lmax=',lmax,' does not agree with input h22p=',h22p,ch10,& 226 & ' setting lmax to 1' 227 call wrtout(ab_out,message,'COLL') 228 call wrtout(std_out, message,'COLL') 229 lmax=1 230 end if 231 end if 232 233 234 if (abs(h33p)>1.d-08) then 235 nproj(2)=3 236 if (lmax<1) then 237 write(message, '(a,a,a,a,i5,a,e12.4,a,a)' ) ch10,& 238 & ' psp3in : COMMENT -',ch10,& 239 & ' input lmax=',lmax,' does not agree with input h33p=',h33p,ch10,& 240 & ' setting lmax to 1' 241 call wrtout(ab_out,message,'COLL') 242 call wrtout(std_out, message,'COLL') 243 lmax=1 244 end if 245 end if 246 247 if (abs(h11d)>1.d-08) then 248 nproj(3)=1 249 if (lmax<2) then 250 write(message, '(a,a,a,a,i5,a,e12.4,a,a)' ) ch10,& 251 & ' psp3in : COMMENT -',ch10,& 252 & ' input lmax=',lmax,' does not agree with input h11d=',h11d,ch10,& 253 & ' setting lmax to 2' 254 call wrtout(ab_out,message,'COLL') 255 call wrtout(std_out, message,'COLL') 256 lmax=2 257 end if 258 end if 259 260 if (abs(h22d)>1.d-08) then 261 nproj(3)=2 262 if (lmax<2) then 263 write(message, '(a,a,a,a,i5,a,e12.4,a,a)' ) ch10,& 264 & ' psp3in : COMMENT -',ch10,& 265 & ' input lmax=',lmax,' does not agree with input h22d=',h22d,ch10,& 266 & ' setting lmax to 2' 267 call wrtout(ab_out,message,'COLL') 268 call wrtout(std_out, message,'COLL') 269 lmax=2 270 end if 271 end if 272 273 if (abs(h33d)>1.d-08) then 274 nproj(3)=3 275 if (lmax<2) then 276 write(message, '(a,a,a,a,i5,a,e12.4,a,a)' ) ch10,& 277 & ' psp3in : COMMENT -',ch10,& 278 & ' input lmax=',lmax,' does not agree with input h33d=',h33d,ch10,& 279 & ' setting lmax to 2' 280 call wrtout(ab_out,message,'COLL') 281 call wrtout(std_out, message,'COLL') 282 lmax=2 283 end if 284 end if 285 286 if (abs(h11f)>1.d-08) then 287 nproj(4)=1 288 if (lmax<3) then 289 write(message, '(a,a,a,a,i5,a,e12.4,a,a)' ) ch10,& 290 & ' psp3in : COMMENT -',ch10,& 291 & ' input lmax=',lmax,' does not agree with input h11f=',h11f,ch10,& 292 & ' setting lmax to 3' 293 call wrtout(ab_out,message,'COLL') 294 call wrtout(std_out, message,'COLL') 295 lmax=3 296 end if 297 end if 298 299 if(pspso/=0) then 300 301 if (abs(k11p)>1.d-08) then 302 nproj(psps%mpsang+1)=1 303 if (lmax<1) then 304 write(message, '(a,a,a,a,i5,a,e12.4,a,a)' ) ch10,& 305 & ' psp3in : COMMENT -',ch10,& 306 & ' input lmax=',lmax,' does not agree with input k11p=',k11p,ch10,& 307 & ' setting lmax to 1' 308 call wrtout(ab_out,message,'COLL') 309 call wrtout(std_out, message,'COLL') 310 lmax=1 311 end if 312 end if 313 314 if (abs(k22p)>1.d-08) then 315 nproj(psps%mpsang+1)=2 316 if (lmax<1) then 317 write(message, '(a,a,a,a,i5,a,e12.4,a,a)' ) ch10,& 318 & ' psp3in : COMMENT -',ch10,& 319 & ' input lmax=',lmax,' does not agree with input k22p=',k22p,ch10,& 320 & ' setting lmax to 1' 321 call wrtout(ab_out,message,'COLL') 322 call wrtout(std_out, message,'COLL') 323 lmax=1 324 end if 325 end if 326 327 328 if (abs(k33p)>1.d-08) then 329 nproj(psps%mpsang+1)=3 330 if (lmax<1) then 331 write(message, '(a,a,a,a,i5,a,e12.4,a,a)' ) ch10,& 332 & ' psp3in : COMMENT -',ch10,& 333 & ' input lmax=',lmax,' does not agree with input k33p=',k33p,ch10,& 334 & ' setting lmax to 1' 335 call wrtout(ab_out,message,'COLL') 336 call wrtout(std_out, message,'COLL') 337 lmax=1 338 end if 339 end if 340 341 if (abs(k11d)>1.d-08) then 342 nproj(psps%mpsang+2)=1 343 if (lmax<2) then 344 write(message, '(a,a,a,a,i5,a,e12.4,a,a)' ) ch10,& 345 & ' psp3in : COMMENT -',ch10,& 346 & ' input lmax=',lmax,' does not agree with input k11d=',k11d,ch10,& 347 & ' setting lmax to 2' 348 call wrtout(ab_out,message,'COLL') 349 call wrtout(std_out, message,'COLL') 350 lmax=2 351 end if 352 end if 353 354 if (abs(k22d)>1.d-08) then 355 nproj(psps%mpsang+2)=2 356 if (lmax<2) then 357 write(message, '(a,a,a,a,i5,a,e12.4,a,a)' ) ch10,& 358 & ' psp3in : COMMENT -',ch10,& 359 & ' input lmax=',lmax,' does not agree with input k22d=',k22d,ch10,& 360 & ' setting lmax to 2' 361 call wrtout(ab_out,message,'COLL') 362 call wrtout(std_out, message,'COLL') 363 lmax=2 364 end if 365 end if 366 367 if (abs(k33d)>1.d-08) then 368 nproj(psps%mpsang+2)=3 369 if (lmax<2) then 370 write(message, '(a,a,a,a,i5,a,e12.4,a,a)' ) ch10,& 371 & ' psp3in : COMMENT -',ch10,& 372 & ' input lmax=',lmax,' does not agree with input k33d=',k33d,ch10,& 373 & ' setting lmax to 2' 374 call wrtout(ab_out,message,'COLL') 375 call wrtout(std_out, message,'COLL') 376 lmax=2 377 end if 378 end if 379 380 if (abs(k11f)>1.d-08) then 381 nproj(psps%mpsang+3)=1 382 if (lmax<3) then 383 write(message, '(a,a,a,a,i5,a,e12.4,a,a)' ) ch10,& 384 & ' psp3in : COMMENT -',ch10,& 385 & ' input lmax=',lmax,' does not agree with input k11f=',k11f,ch10,& 386 & ' setting lmax to 3' 387 call wrtout(ab_out,message,'COLL') 388 call wrtout(std_out, message,'COLL') 389 lmax=3 390 end if 391 end if 392 393 end if 394 395 !Store the coefficients. 396 psps%gth_params%set(ipsp) = .true. 397 psps%gth_params%psppar(0, :, ipsp) = (/ rloc, cc1, cc2, cc3, cc4, zero, zero /) 398 psps%gth_params%psppar(1, :, ipsp) = (/ rrs, h11s, h22s, h33s, zero, zero, zero /) 399 psps%gth_params%psppar(2, :, ipsp) = (/ rrp, h11p, h22p, h33p, zero, zero, zero /) 400 psps%gth_params%psppar(3, :, ipsp) = (/ rrd, h11d, h22d, h33d, zero, zero, zero /) 401 psps%gth_params%psppar(4, :, ipsp) = (/ rrf, h11f, zero, zero, zero, zero, zero /) 402 403 !Store the k coefficients 404 psps%gth_params%psp_k_par(1, :, ipsp) = (/ zero, zero, zero /) 405 psps%gth_params%psp_k_par(2, :, ipsp) = (/ k11p, k22p, k33p /) 406 psps%gth_params%psp_k_par(3, :, ipsp) = (/ k11d, k22d, k33d /) 407 psps%gth_params%psp_k_par(4, :, ipsp) = (/ k11f, zero, zero /) 408 409 !Additionnal wavelet parameters 410 if (dtset%usewvl == 1) then 411 call wvl_descr_psp_fill(psps%gth_params, ipsp, 0, int(psps%zionpsp(ipsp)), int(psps%znuclpsp(ipsp)), tmp_unit) 412 end if 413 414 !Initialize array indlmn array giving l,m,n,ln,lm,s for i=lmn 415 nso=1;if(pspso/=0) nso=2 416 index=0;iln=0;indlmn(:,:)=0 417 do nn=1,nso 418 do ipsang=1+(nn-1)*(lmax+1),nn*lmax+1 419 if (nproj(ipsang)>0) then 420 ll=ipsang-(nn-1)*lmax-1 421 do kk=1,nproj(ipsang) 422 iln=iln+1 423 do mm=1,2*ll*psps%useylm+1 424 index=index+1 425 indlmn(1,index)=ll 426 indlmn(2,index)=mm-ll*psps%useylm-1 427 indlmn(3,index)=kk 428 indlmn(4,index)=ll*ll+(1-psps%useylm)*ll+mm 429 indlmn(5,index)=iln 430 indlmn(6,index)=nn 431 end do 432 end do 433 end if 434 end do 435 end do 436 437 ABI_ALLOCATE(dvlspl,(psps%mqgrid_ff)) 438 !First, the local potential -- compute on q grid and fit spline 439 call psp2lo(cc1,cc2,cc3,cc4,dvlspl,epsatm,psps%mqgrid_ff,psps%qgrid_ff,& 440 & vlspl(:,1),rloc,.true.,yp1,ypn,zion) 441 ABI_DEALLOCATE(dvlspl) 442 443 !DEBUG 444 !write(std_out,*)' psp3in : after psp2lo ' 445 !ENDDEBUG 446 447 !Fit spline to q^2 V(q) (Numerical Recipes subroutine) 448 ABI_ALLOCATE(work_space,(psps%mqgrid_ff)) 449 ABI_ALLOCATE(work_spl,(psps%mqgrid_ff)) 450 call spline (psps%qgrid_ff,vlspl(:,1),psps%mqgrid_ff,yp1,ypn,work_spl) 451 vlspl(:,2)=work_spl(:) 452 ABI_DEALLOCATE(work_space) 453 ABI_DEALLOCATE(work_spl) 454 455 !Second, compute KB energies and form factors and fit splines 456 ekb(:)=zero 457 458 !Check if any nonlocal projectors are being used 459 mproj=maxval(nproj) 460 if (mproj>0) then 461 462 ABI_ALLOCATE(ekb_sr,(psps%mpsang,mproj)) 463 ABI_ALLOCATE(ffspl_sr,(psps%mqgrid_ff,2,psps%mpsang,mproj)) 464 ABI_ALLOCATE(ekb_so,(psps%mpsang,mproj)) 465 ABI_ALLOCATE(ffspl_so,(psps%mqgrid_ff,2,psps%mpsang,mproj)) 466 467 call psp3nl(ekb_sr,ffspl_sr,h11s,h22s,h33s,h11p,h22p,h33p,h11d,h22d,& 468 & h33d,h11f,mproj,psps%mpsang,psps%mqgrid_ff,psps%qgrid_ff,rrd,rrf,rrp,rrs) 469 if(pspso/=0) then 470 call psp3nl(ekb_so,ffspl_so,zero,zero,zero,k11p,k22p,k33p,k11d,& 471 & k22d,k33d,k11f,mproj,psps%mpsang,psps%mqgrid_ff,psps%qgrid_ff,rrd,rrf,rrp,rrs) 472 end if 473 474 475 ! Convert ekb and ffspl 476 iln0=0 477 do jj=1,psps%lmnmax 478 iln=indlmn(5,jj) 479 if (iln>iln0) then 480 iln0=iln 481 if (indlmn(6,jj)<=1) then 482 ekb(iln)=ekb_sr(1+indlmn(1,jj),indlmn(3,jj)) 483 ffspl(:,:,iln)=ffspl_sr(:,:,1+indlmn(1,jj),indlmn(3,jj)) 484 else 485 ekb(iln)=ekb_so(1+indlmn(1,jj),indlmn(3,jj)) 486 ffspl(:,:,iln)=ffspl_so(:,:,1+indlmn(1,jj),indlmn(3,jj)) 487 end if 488 end if 489 end do 490 491 ABI_DEALLOCATE(ekb_sr) 492 ABI_DEALLOCATE(ffspl_sr) 493 ABI_DEALLOCATE(ekb_so) 494 ABI_DEALLOCATE(ffspl_so) 495 496 end if 497 498 return 499 500 ! Handle IO error 501 10 continue 502 MSG_ERROR(errmsg) 503 504 end subroutine psp3in
ABINIT/psp3nl [ Functions ]
NAME
psp3nl
FUNCTION
Hartwigsen-Goedecker-Hutter nonlocal pseudopotential (from preprint of 1998). Uses Gaussians for fully nonlocal form, analytic expressions.
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
549 subroutine psp3nl(ekb,ffspl,h11s,h22s,h33s,h11p,h22p,h33p,h11d,h22d,& 550 & h33d,h11f,mproj,mpsang,mqgrid,qgrid,rrd,rrf,rrp,rrs) 551 552 use defs_basis 553 use m_splines 554 use m_errors 555 use m_profiling_abi 556 557 !This section has been created automatically by the script Abilint (TD). 558 !Do not modify the following lines by hand. 559 #undef ABI_FUNC 560 #define ABI_FUNC 'psp3nl' 561 !End of the abilint section 562 563 implicit none 564 565 !Arguments ------------------------------------ 566 !scalars 567 integer,intent(in) :: mproj,mpsang,mqgrid 568 real(dp),intent(in) :: h11d,h11f,h11p,h11s,h22d,h22p,h22s,h33d,h33p,h33s,rrd 569 real(dp),intent(in) :: rrf,rrp,rrs 570 !arrays 571 real(dp),intent(in) :: qgrid(mqgrid) 572 real(dp),intent(out) :: ekb(mpsang,mproj),ffspl(mqgrid,2,mpsang,mproj) 573 574 !Local variables------------------------------- 575 !scalars 576 integer :: info,iproj,iqgrid,ldz,mu,nproj,nu 577 real(dp) :: qmax 578 character(len=500) :: message 579 character :: jobz,uplo 580 !arrays 581 real(dp) :: ap(2,9),rwork1(9),work1(2,9),ww(3),yp1j(3),ypnj(3) 582 real(dp),allocatable :: ppspl(:,:,:,:),uu(:,:),work(:),zz(:,:,:) 583 584 ! ************************************************************************* 585 586 !DEBUG 587 !write(std_out,*)' psp3nl : enter ' 588 !stop 589 !ENDDEBUG 590 591 ABI_ALLOCATE(ppspl,(mqgrid,2,mpsang,mproj)) 592 ABI_ALLOCATE(work,(mqgrid)) 593 594 qmax=qgrid(mqgrid) 595 jobz='v' 596 uplo='u' 597 ekb(:,:)=0.0d0 598 599 !--------------------------------------------------------------- 600 !Treat s channel 601 602 nproj=0 603 ap(:,:)=0.0d0 604 !If there is at least one s-projector 605 if ( abs(h11s) >= 1.0d-8 ) then 606 nproj=1 ; ldz=1 ; ap(1,1)=h11s 607 end if 608 nproj=1 609 !If there is a second projector 610 if ( abs(h22s) >= 1.0d-8 ) then 611 nproj=2 ; ldz=2 ; ap(1,3)=h22s 612 ap(1,2)=-0.5d0*sqrt(0.6d0)*h22s 613 end if 614 !If there is a third projector 615 if ( abs(h33s) >= 1.0d-8 ) then 616 nproj=3 ; ldz=3 ; ap(1,6)=h33s 617 ap(1,4)=0.5d0*sqrt(5.d0/21.d0)*h33s 618 ap(1,5)=-0.5d0*sqrt(100.d0/63.d0)*h33s 619 end if 620 621 if(nproj/=0)then 622 623 ABI_ALLOCATE(uu,(nproj,nproj)) 624 ABI_ALLOCATE(zz,(2,nproj,nproj)) 625 626 if (nproj > 1) then 627 call ZHPEV(jobz,uplo,nproj,ap,ww,zz,ldz,work1,rwork1,info) 628 uu(:,:)=zz(1,:,:) 629 else 630 ww(1)=h11s 631 uu(1,1)=1.0d0 632 end if 633 634 ! Initialization of ekb, and spline fitting 635 do iproj=1,nproj 636 ekb(1,iproj)=ww(iproj)*32.d0*(rrs**3)*(pi**2.5d0)/(4.d0*pi)**2 637 if(iproj==1)then 638 do iqgrid=1,mqgrid 639 ppspl(iqgrid,1,1,1)=exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrs)**2) 640 end do 641 yp1j(1)=0.d0 642 ypnj(1)=-(two_pi*rrs)**2*qmax*exp(-0.5d0*(two_pi*qmax*rrs)**2) 643 else if(iproj==2)then 644 do iqgrid=1,mqgrid 645 ppspl(iqgrid,1,1,2)=2.0d0/sqrt(15.0d0) & 646 & *exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrs)**2) & 647 & *( 3.d0-(two_pi*qgrid(iqgrid)*rrs)**2 ) 648 end do 649 yp1j(2)=0.0d0 650 ypnj(2)=2.0d0/sqrt(15.0d0)*(two_pi*rrs)**2*qmax & 651 & *exp(-0.5d0*(two_pi*qmax*rrs)**2) * (-5.d0+(two_pi*qmax*rrs)**2) 652 else if(iproj==3)then 653 do iqgrid=1,mqgrid 654 ppspl(iqgrid,1,1,3)=(4.0d0/3.0d0)/sqrt(105.0d0)*& 655 & exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrs)**2) * & 656 & (15.0d0-10.0d0*(two_pi*qgrid(iqgrid)*rrs)**2 + & 657 & (two_pi*qgrid(iqgrid)*rrs)**4) 658 end do 659 yp1j(3)=0.0d0 660 ypnj(3)=(4.0d0/3.0d0)/sqrt(105.0d0)*exp(-0.5d0*(two_pi*qmax*rrs)**2) * & 661 & (two_pi*rrs)**2*qmax*(-35.0d0+14d0*(two_pi*qmax*rrs)**2-(two_pi*qmax*rrs)**4) 662 end if 663 call spline(qgrid,ppspl(:,1,1,iproj),mqgrid,& 664 & yp1j(iproj),ypnj(iproj),ppspl(:,2,1,iproj)) 665 end do 666 667 ! Linear combination using the eigenvectors 668 ffspl(:,:,1,:)=0.0d0 669 do mu=1,nproj 670 do nu=1,nproj 671 do iqgrid=1,mqgrid 672 ffspl(iqgrid,1:2,1,mu)=ffspl(iqgrid,1:2,1,mu) & 673 & +uu(nu,mu)*ppspl(iqgrid,1:2,1,nu) 674 end do 675 end do 676 end do 677 678 ABI_DEALLOCATE(uu) 679 ABI_DEALLOCATE(zz) 680 681 ! End condition on nproj(/=0) 682 end if 683 684 !DEBUG 685 !write(std_out,*)' psp3nl : after s channel ' 686 !stop 687 !ENDDEBUG 688 689 !-------------------------------------------------------------------- 690 !Now treat p channel 691 692 nproj=0 693 ap(:,:)=0.0d0 694 !If there is at least one projector 695 if ( abs(h11p) >= 1.0d-8 ) then 696 nproj=1 ; ldz=1 ; ap(1,1)=h11p 697 end if 698 !If there is a second projector 699 if ( abs(h22p) >= 1.0d-8 ) then 700 nproj=2 ; ldz=2 ; ap(1,3)=h22p 701 ap(1,2)=-0.5d0*sqrt(5.d0/7.d0)*h22p 702 end if 703 !If there is a third projector 704 if ( abs(h33p) >= 1.0d-8 ) then 705 nproj=3 ; ldz=3 ; ap(1,6)=h33p 706 ap(1,4)= (1.d0/6.d0)*sqrt(35.d0/11.d0)*h33p 707 ap(1,5)=-(1.d0/6.d0)*(14.d0/sqrt(11.d0))*h33p 708 end if 709 710 if(nproj/=0)then 711 712 ABI_ALLOCATE(uu,(nproj,nproj)) 713 ABI_ALLOCATE(zz,(2,nproj,nproj)) 714 715 if (nproj > 1) then 716 call ZHPEV(jobz,uplo,nproj,ap,ww,zz,ldz,work1,rwork1,info) 717 uu(:,:)=zz(1,:,:) 718 else 719 ww(1)=h11p 720 uu(1,1)=1.0d0 721 end if 722 723 ! Initialization of ekb, and spline fitting 724 do iproj=1,nproj 725 ekb(2,iproj)=ww(iproj)*64.d0*(rrp**5)*(pi**2.5d0)/(4.d0*pi)**2 726 if(iproj==1)then 727 do iqgrid=1,mqgrid 728 ppspl(iqgrid,1,2,1)=(1.0d0/sqrt(3.0d0))* & 729 & exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrp)**2) * (two_pi*qgrid(iqgrid)) 730 end do 731 yp1j(1)=two_pi*(1.0d0/sqrt(3.0d0)) 732 ypnj(1)=-two_pi*((two_pi*qmax*rrp)**2-1.d0)*exp(-0.5d0*(two_pi*qmax*rrp)**2)*& 733 & (1.0d0/sqrt(3.0d0)) 734 else if(iproj==2)then 735 do iqgrid=1,mqgrid 736 ppspl(iqgrid,1,2,2)=(2.0d0/sqrt(105.0d0))* & 737 & exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrp)**2) * & 738 & (two_pi*qgrid(iqgrid))*(5.0d0-(two_pi*qgrid(iqgrid)*rrp)**2) 739 end do 740 yp1j(2)=(5.0d0*two_pi)*(2.0d0/sqrt(105.0d0)) 741 ypnj(2)=(2.0d0/sqrt(105.0d0))*two_pi*exp(-0.5d0*(two_pi*qmax*rrp)**2)* & 742 & (-8*(two_pi*qmax*rrp)**2 + (two_pi*qmax*rrp)**4 + 5.0d0) 743 else if(iproj==3)then 744 do iqgrid=1,mqgrid 745 ppspl(iqgrid,1,2,3)=(4.0d0/3.0d0)/sqrt(1155d0)*& 746 & exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrp)**2) * & 747 & (two_pi*qgrid(iqgrid))*& 748 & (35.0d0-14.0d0*(two_pi*qgrid(iqgrid)*rrp)**2+(two_pi*qgrid(iqgrid)*rrp)**4) 749 end do 750 yp1j(3)=(35.0d0*two_pi)*(4.0d0/3.0d0)/sqrt(1155.0d0) 751 ypnj(3)=(4.0d0/3.0d0)/sqrt(1155.0d0)*two_pi*exp(-0.5d0*(two_pi*qmax*rrp)**2)* & 752 & (35.0d0-77.0d0*(two_pi*qmax*rrp)**2+19.0d0*(two_pi*qmax*rrp)**4 - & 753 & (two_pi*qmax*rrp)**6) 754 end if 755 call spline(qgrid,ppspl(:,1,2,iproj),mqgrid,& 756 & yp1j(iproj),ypnj(iproj),ppspl(:,2,2,iproj)) 757 end do 758 759 ! Linear combination using the eigenvectors 760 ffspl(:,:,2,:)=0.0d0 761 do mu=1,nproj 762 do nu=1,nproj 763 do iqgrid=1,mqgrid 764 ffspl(iqgrid,1:2,2,mu)=ffspl(iqgrid,1:2,2,mu) & 765 & +uu(nu,mu)*ppspl(iqgrid,1:2,2,nu) 766 end do 767 end do 768 end do 769 770 ABI_DEALLOCATE(uu) 771 ABI_DEALLOCATE(zz) 772 773 ! End condition on nproj(/=0) 774 end if 775 776 !DEBUG 777 !write(std_out,*)' psp3nl : after p channel ' 778 !stop 779 !ENDDEBUG 780 781 !----------------------------------------------------------------------- 782 !Now treat d channel. 783 784 nproj=0 785 ap(:,:)=0.0d0 786 !If there is at least one projector 787 if ( abs(h11d) >= 1.0d-8 ) then 788 nproj=1 ; ldz=1 ; ap(1,1)=h11d 789 end if 790 !If there is a second projector 791 if ( abs(h22d) >= 1.0d-8 ) then 792 nproj=2 ; ldz=2 ; ap(1,3)=h22d 793 ap(1,2)=-0.5d0*sqrt(7.d0/9.d0)*h22d 794 end if 795 !If there is a third projector. Warning : only two projectors are allowed. 796 if ( abs(h33d) >= 1.0d-8 ) then 797 write(message, '(a,a,a)' )& 798 & ' only two d-projectors are allowed ',ch10,& 799 & ' Action : check your pseudopotential file.' 800 MSG_ERROR(message) 801 ! nproj=3 ; ldz=3 ; ap(1,6)=h33d 802 ! ap(1,4)= 0.5d0*sqrt(63.d0/143.d0)*h33d 803 ! ap(1,5)= -0.5d0*(18.d0/sqrt(143.d0))*h33d 804 end if 805 806 if(nproj/=0)then 807 808 ABI_ALLOCATE(uu,(nproj,nproj)) 809 ABI_ALLOCATE(zz,(2,nproj,nproj)) 810 811 if (nproj > 1) then 812 call ZHPEV(jobz,uplo,nproj,ap,ww,zz,ldz,work1,rwork1,info) 813 uu(:,:)=zz(1,:,:) 814 else 815 ww(1)=h11d 816 uu(1,1)=1.0d0 817 end if 818 819 ! Initialization of ekb, and spline fitting 820 do iproj=1,nproj 821 ekb(3,iproj)=ww(iproj)*128.d0*(rrd**7)*(pi**2.5d0)/(4.d0*pi)**2 822 if(iproj==1)then 823 do iqgrid=1,mqgrid 824 ppspl(iqgrid,1,3,1)=(1.0d0/sqrt(15.0d0))* & 825 & exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrd)**2) * (two_pi*qgrid(iqgrid))**2 826 end do 827 yp1j(1)=0.0d0 828 ypnj(1)=(1.0d0/sqrt(15.0d0))*(two_pi**2)*& 829 & exp(-0.5d0*(two_pi*qmax*rrd)**2)*qmax*(2d0-(two_pi*qmax*rrd)**2) 830 else if(iproj==2)then 831 do iqgrid=1,mqgrid 832 ppspl(iqgrid,1,3,2)=(2.0d0/3.0d0)/sqrt(105.0d0)* & 833 & exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrd)**2) * & 834 & ((two_pi*qgrid(iqgrid))**2)*(7.0d0-(two_pi*qgrid(iqgrid)*rrd)**2) 835 end do 836 yp1j(2)=0.0d0 837 ypnj(2)=(2.0d0/3.0d0)/sqrt(105.0d0)*exp(-0.5d0*(two_pi*qmax*rrd)**2)* & 838 & qmax*(two_pi**2)*( (two_pi*qmax*rrd)**4 - 11.0d0*(two_pi*qmax*rrd)**2 + 14.0d0) 839 end if 840 call spline(qgrid,ppspl(:,1,3,iproj),mqgrid,& 841 & yp1j(iproj),ypnj(iproj),ppspl(:,2,3,iproj)) 842 end do 843 844 ! Linear combination using the eigenvectors 845 ffspl(:,:,3,:)=0.0d0 846 do mu=1,nproj 847 do nu=1,nproj 848 do iqgrid=1,mqgrid 849 ffspl(iqgrid,1:2,3,mu)=ffspl(iqgrid,1:2,3,mu) & 850 & +uu(nu,mu)*ppspl(iqgrid,1:2,3,nu) 851 end do 852 end do 853 end do 854 855 ABI_DEALLOCATE(uu) 856 ABI_DEALLOCATE(zz) 857 858 ! End condition on nproj(/=0) 859 end if 860 861 !DEBUG 862 !write(std_out,*)' psp3nl : after d channel ' 863 !stop 864 !ENDDEBUG 865 866 !----------------------------------------------------------------------- 867 !Treat now f channel (max one projector ! - so do not use ppspl) 868 869 !l=3 first projector 870 if (abs(h11f)>1.d-12) then 871 ekb(4,1)=h11f*(256.0d0/105.0d0)*(rrf**9)*(pi**2.5d0)/(4.d0*pi)**2 872 do iqgrid=1,mqgrid 873 ffspl(iqgrid,1,4,1)=(two_pi*qgrid(iqgrid))**3* & 874 & exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrf)**2) 875 end do 876 ! Compute yp1,ypn=derivatives of f(q) at q=0, q=qgrid(mqgrid) 877 yp1j(1)=0d0 878 ypnj(1)=(two_pi**3)*qmax**2*exp(-0.5d0*(two_pi*qmax*rrf)**2)*& 879 & (3.0d0-(two_pi*qmax*rrf)**2) 880 ! Fit spline to get second derivatives by spline fit 881 call spline(qgrid,ffspl(:,1,4,1),mqgrid,& 882 & yp1j(1),ypnj(1),ffspl(:,2,4,1)) 883 end if 884 885 !----------------------------------------------------------------------- 886 887 ABI_DEALLOCATE(ppspl) 888 ABI_DEALLOCATE(work) 889 890 !DEBUG 891 !write(std_out,*)' psp3nl : exit ' 892 !stop 893 !ENDDEBUG 894 895 end subroutine psp3nl