TABLE OF CONTENTS
- ABINIT/m_paw_atom
- m_paw_atom/atompaw_atompaw_shapebes
- m_paw_atom/atompaw_dij0
- m_paw_atom/atompaw_kij
- m_paw_atom/atompaw_shpfun
- m_paw_atom/atompaw_vhnzc
ABINIT/m_paw_atom [ Modules ]
NAME
m_paw_atom
FUNCTION
atompaw related operations
COPYRIGHT
Copyright (C) 2012-2022 ABINIT group (T. Rangel, MT, JWZ, GJ) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
NOTES
FOR DEVELOPPERS: in order to preserve the portability of libPAW library, please consult ~abinit/src/??_libpaw/libpaw-coding-rules.txt
SOURCE
20 #include "libpaw.h" 21 22 module m_paw_atom 23 24 USE_DEFS 25 USE_MSG_HANDLING 26 USE_MEMORY_PROFILING 27 28 use m_paw_numeric, only : paw_jbessel, paw_solvbes, paw_spline, paw_splint 29 use m_pawrad, only : pawrad_type, simp_gen, poisson, pawrad_deducer0, bound_deriv, pawrad_ifromr 30 use m_pawtab, only : pawtab_type 31 32 implicit none 33 34 private 35 36 public:: atompaw_shpfun 37 public:: atompaw_shapebes 38 public:: atompaw_vhnzc 39 public:: atompaw_dij0 40 public:: atompaw_kij
m_paw_atom/atompaw_atompaw_shapebes [ Functions ]
[ Top ] [ m_paw_atom ] [ Functions ]
NAME
atompaw_shapebes
FUNCTION
Find al and ql parameters for a "Bessel" shape function: Shape(r)=al1.jl(ql1.r)+al2.jl(ql2.r) such as Shape(r) and 2 derivatives are zero at r=rc Intg_0_rc[Shape(r).r^(l+2).dr]=1
INPUTS
ll= l quantum number rc= cut-off radius
OUTPUT
al(2)= al coefficients ql(2)= ql factors
SOURCE
189 subroutine atompaw_shapebes(al,ql,ll,rc) 190 191 !Arguments ------------------------------------ 192 !scalars 193 integer :: ll 194 real(dp) :: rc 195 !arrays 196 real(dp) :: al(2),ql(2) 197 198 !Local variables------------------------------- 199 !scalars 200 integer :: ii 201 real(dp) :: alpha,beta,det,jbes,jbesp,jbespp,qr 202 !arrays 203 real(dp) :: amat(2,2),bb(2) 204 205 ! ************************************************************************* 206 207 alpha=1._dp;beta=0._dp 208 call paw_solvbes(ql,alpha,beta,ll,2) 209 ql(1:2)=ql(1:2)/rc 210 211 do ii=1,2 212 qr=ql(ii)*rc 213 call paw_jbessel(jbes,jbesp,jbespp,ll,1,qr) 214 amat(1,ii)=jbesp*ql(ii) 215 call paw_jbessel(jbes,jbesp,jbespp,ll+1,0,qr) 216 amat(2,ii)=jbes*rc**(ll+2)/ql(ii) ! Intg_0_rc[jl(qr).r^(l+2).dr] 217 end do 218 219 bb(1)=zero;bb(2)=one 220 221 det=amat(1,1)*amat(2,2)-amat(1,2)*amat(2,1) 222 al(1)=(amat(2,2)*bb(1)-amat(1,2)*bb(2))/det 223 al(2)=(amat(1,1)*bb(2)-amat(2,1)*bb(1))/det 224 225 end subroutine atompaw_shapebes
m_paw_atom/atompaw_dij0 [ Functions ]
[ Top ] [ m_paw_atom ] [ Functions ]
NAME
atompaw_dij0
FUNCTION
PAW: Compute "frozen" values of pseudopotential strengths Dij = Dij0
INPUTS
indlmn(6,lmnmax)= array giving l,m,n,lm,ln,s for i=lmn kij(pawtab%lmn2_size)= kinetic part of Dij lmnmax=max number of (l,m,n) components over all type of psps ncore(:)=atomic core density opt_init=flag defining the storage of PAW atomic data 0: PAW atomic data have not been initialized (in pawtab) 1: PAW atomic data have been initialized (in pawtab) pawtab <type(pawtab_type)>=paw tabulated starting data radmesh <type(pawrad_type)>=paw radial mesh (and related data) radmesh_core <type(pawrad_type)>=radial mesh (and related data) for the core densities radmesh_vloc <type(pawrad_type)>=radial mesh (and related data) for the local potential (VH(tnZc)) vhtnzc(:)= local potential VH(tnZc) znucl= valence and total charge of the atomic species
OUTPUT
pawtab%dij0(pawtab%lmn2_size)= Frozen part of the Dij term
SOURCE
310 subroutine atompaw_dij0(indlmn,kij,lmnmax,ncore,opt_init,pawtab,& 311 & radmesh,radmesh_core,radmesh_vloc,vhtnzc,znucl) 312 313 !Arguments --------------------------------------------- 314 !scalars 315 integer,intent(in) :: lmnmax,opt_init 316 real(dp),intent(in) :: znucl 317 type(pawrad_type),intent(in) :: radmesh,radmesh_core,radmesh_vloc 318 type(pawtab_type),intent(inout) :: pawtab 319 !arrays 320 integer,intent(in) :: indlmn(6,lmnmax) 321 real(dp),intent(in) :: kij(pawtab%lmn2_size) 322 real(dp),intent(in) :: ncore(:),vhtnzc(:) 323 !real(dp),optional,intent(in) :: vminushalf(:) 324 325 !Local variables --------------------------------------- 326 integer :: il,ilm,iln,ilmn,j0lmn,jl,jlm,jln,jlmn,klmn,lmn2_size,meshsz,meshsz_core 327 integer :: meshsz_vhtnzc,meshsz_vmh 328 real(dp) :: intg,intvh,yp1,ypn 329 real(dp),allocatable :: ff(:),r2k(:),shpf(:),vhnzc(:),vhtnzc_sph(:),work1(:),work2(:) 330 331 ! ********************************************************************* 332 333 lmn2_size=pawtab%lmn2_size 334 meshsz_vhtnzc=size(vhtnzc) 335 meshsz=min(radmesh%mesh_size,radmesh_core%mesh_size,radmesh_vloc%mesh_size,meshsz_vhtnzc) 336 LIBPAW_ALLOCATE(ff,(meshsz)) 337 338 !Retrieve VH(tnZc) on the correct radial mesh 339 LIBPAW_ALLOCATE(vhtnzc_sph,(meshsz)) 340 if ((radmesh%mesh_type/=radmesh_vloc%mesh_type).or.& 341 & (radmesh%rstep /=radmesh_vloc%rstep) .or.& 342 & (radmesh%lstep /=radmesh_vloc%lstep)) then 343 call bound_deriv(vhtnzc,radmesh_vloc,meshsz_vhtnzc,yp1,ypn) 344 LIBPAW_ALLOCATE(work1,(meshsz_vhtnzc)) 345 LIBPAW_ALLOCATE(work2,(meshsz_vhtnzc)) 346 call paw_spline(radmesh_vloc%rad,vhtnzc,meshsz_vhtnzc,yp1,ypn,work1) 347 call paw_splint(meshsz_vhtnzc,radmesh_vloc%rad,vhtnzc,work1,meshsz,radmesh%rad(1:meshsz),vhtnzc_sph) 348 LIBPAW_DEALLOCATE(work1) 349 LIBPAW_DEALLOCATE(work2) 350 else 351 vhtnzc_sph(1:meshsz)=vhtnzc(1:meshsz) 352 end if 353 354 !Kinetic part of Dij0 355 !==================== 356 pawtab%dij0(1:lmn2_size)=kij(1:lmn2_size) 357 358 !Computation of <phi_i|vh(nZc)|phi_j> on the PAW sphere 359 !====================================================== 360 meshsz_core=size(ncore) 361 LIBPAW_ALLOCATE(vhnzc,(meshsz_core)) 362 call atompaw_vhnzc(ncore,radmesh_core,vhnzc,znucl) 363 do jlmn=1,pawtab%lmn_size 364 j0lmn=jlmn*(jlmn-1)/2 365 jlm=indlmn(4,jlmn);jln=indlmn(5,jlmn) 366 do ilmn=1,jlmn 367 klmn=j0lmn+ilmn 368 ilm=indlmn(4,ilmn);iln=indlmn(5,ilmn) 369 if (jlm==ilm) then 370 ff(1:meshsz)=pawtab%phi(1:meshsz,iln)*pawtab%phi(1:meshsz,jln)*vhnzc(1:meshsz) 371 call simp_gen(intg,ff,radmesh) 372 pawtab%dij0(klmn)=pawtab%dij0(klmn)+intg 373 end if 374 end do 375 end do 376 LIBPAW_DEALLOCATE(vhnzc) 377 378 !Computation of -<tphi_i|vh(tnZc)|tphi_j> on the PAW sphere 379 !========================================================== 380 do jlmn=1,pawtab%lmn_size 381 j0lmn=jlmn*(jlmn-1)/2 382 jlm=indlmn(4,jlmn);jln=indlmn(5,jlmn) 383 do ilmn=1,jlmn 384 klmn=j0lmn+ilmn 385 ilm=indlmn(4,ilmn);iln=indlmn(5,ilmn) 386 if (jlm==ilm) then 387 ff(1:meshsz)=pawtab%tphi(1:meshsz,iln)*pawtab%tphi(1:meshsz,jln)*vhtnzc_sph(1:meshsz) 388 call simp_gen(intg,ff,radmesh) 389 pawtab%dij0(klmn)=pawtab%dij0(klmn)-intg 390 end if 391 end do 392 end do 393 394 !Computation of <phi_i|vminushalf|phi_j> (if any) 395 !================================================= 396 if(pawtab%has_vminushalf==1) then 397 if(size(pawtab%vminushalf)>=1) then 398 meshsz_vmh=min(meshsz,size(pawtab%vminushalf)) 399 do jlmn=1,pawtab%lmn_size 400 j0lmn=jlmn*(jlmn-1)/2 401 jlm=indlmn(4,jlmn);jln=indlmn(5,jlmn) 402 do ilmn=1,jlmn 403 klmn=j0lmn+ilmn 404 ilm=indlmn(4,ilmn);iln=indlmn(5,ilmn) 405 if (jlm==ilm) then 406 ff(1:meshsz_vmh)=pawtab%phi(1:meshsz_vmh,iln)*pawtab%phi(1:meshsz_vmh,jln)*pawtab%vminushalf(1:meshsz_vmh) 407 call simp_gen(intg,ff(1:meshsz_vmh),radmesh) 408 pawtab%dij0(klmn)=pawtab%dij0(klmn)+intg 409 end if 410 end do 411 end do 412 end if 413 end if 414 415 !Computation of -int[vh(tnzc)*Qijhat(r)dr] 416 !========================================= 417 if (opt_init==0) then 418 LIBPAW_ALLOCATE(shpf,(radmesh%mesh_size)) 419 call atompaw_shpfun(0,radmesh,intg,pawtab,shpf) 420 if (pawtab%shape_type==3) then 421 LIBPAW_ALLOCATE(r2k,(radmesh%int_meshsz)) 422 r2k=zero 423 r2k(2:radmesh%int_meshsz)=shpf(2:radmesh%int_meshsz)*radmesh%rad(2:radmesh%int_meshsz)**2 424 if(radmesh%mesh_type==5) then 425 call simp_gen(intg,r2k,radmesh) 426 else 427 call simp_gen(intg,r2k,radmesh,r_for_intg=pawtab%rshp) 428 end if 429 shpf(1:meshsz)=shpf(1:meshsz)/intg 430 LIBPAW_DEALLOCATE(r2k) 431 end if 432 ff(1:meshsz)=vhtnzc_sph(1:meshsz)*shpf(1:meshsz)*radmesh%rad(1:meshsz)**2 433 LIBPAW_DEALLOCATE(shpf) 434 call simp_gen(intvh,ff,radmesh) 435 do jlmn=1,pawtab%lmn_size 436 j0lmn=jlmn*(jlmn-1)/2 437 jl=indlmn(1,jlmn);jln=indlmn(5,jlmn);jlm=indlmn(4,jlmn) 438 do ilmn=1,jlmn 439 klmn=j0lmn+ilmn 440 il=indlmn(1,ilmn);iln=indlmn(5,ilmn);ilm=indlmn(4,ilmn) 441 if (ilm==jlm) then 442 ff(1:meshsz)=(pawtab%phi (1:meshsz,iln)*pawtab%phi (1:meshsz,jln)& 443 & -pawtab%tphi(1:meshsz,iln)*pawtab%tphi(1:meshsz,jln)) 444 call simp_gen(intg,ff,radmesh) 445 pawtab%dij0(klmn)=pawtab%dij0(klmn)-intvh*intg 446 end if 447 end do 448 end do 449 else 450 ff(1:meshsz)=vhtnzc_sph(1:meshsz)*pawtab%shapefunc(1:meshsz,1)*radmesh%rad(1:meshsz)**2 451 call simp_gen(intvh,ff,radmesh) 452 do jlmn=1,pawtab%lmn_size 453 j0lmn=jlmn*(jlmn-1)/2 454 jl=indlmn(1,jlmn);jln=indlmn(5,jlmn);jlm=indlmn(4,jlmn) 455 do ilmn=1,jlmn 456 klmn=j0lmn+ilmn 457 il=indlmn(1,ilmn);iln=indlmn(5,ilmn);ilm=indlmn(4,ilmn) 458 if (ilm==jlm) then 459 intg=pawtab%qijl(1,klmn)*sqrt(four_pi) 460 pawtab%dij0(klmn)=pawtab%dij0(klmn)-intvh*intg 461 end if 462 end do 463 end do 464 end if 465 466 LIBPAW_DEALLOCATE(ff) 467 LIBPAW_DEALLOCATE(vhtnzc_sph) 468 469 end subroutine atompaw_dij0
m_paw_atom/atompaw_kij [ Functions ]
[ Top ] [ m_paw_atom ] [ Functions ]
NAME
atompaw_kij
FUNCTION
PAW: deduce kinetic part of psp strength (Dij) from the knowledge of frozen Dij (Dij0)
INPUTS
indlmn(6,lmnmax)= array giving l,m,n,lm,ln,s for i=lmn lmnmax=max number of (l,m,n) components over all type of psps ncore(:)=atomic core density opt_init=flag defining the storage of PAW atomic data 0: PAW atomic data have not been initialized (in pawtab) 1: PAW atomic data have been initialized (in pawtab) opt_vhnzc=flag defining the inclusion of VH(nZc) in computation 0: VH(nZc) is not taken into account 1: VH(nZc) is taken into account pawtab <type(pawtab_type)>=paw tabulated starting data radmesh <type(pawrad_type)>=paw radial mesh (and related data) radmesh_core <type(pawrad_type)>=radial mesh (and related data) for the core densities radmesh_vloc <type(pawrad_type)>=radial mesh (and related data) for the local potential (VH(tnZc)) vhtnzc(:)= local potential VH(tnZc) znucl= valence and total charge of the atomic species
OUTPUT
kij(pawtab%lmn2_size)= kinetic part of Dij
SOURCE
503 subroutine atompaw_kij(indlmn,kij,lmnmax,ncore,opt_init,opt_vhnzc,pawtab, & 504 & radmesh,radmesh_core,radmesh_vloc,vhtnzc,znucl) 505 506 !Arguments --------------------------------------------- 507 !scalars 508 integer,intent(in) :: lmnmax,opt_init,opt_vhnzc 509 real(dp),intent(in) :: znucl 510 type(pawrad_type),intent(in) :: radmesh,radmesh_core,radmesh_vloc 511 type(pawtab_type),intent(in) :: pawtab 512 !arrays 513 integer,intent(in) :: indlmn(6,lmnmax) 514 real(dp),intent(out) :: kij(pawtab%lmn2_size) 515 real(dp),intent(in) :: ncore(:) 516 real(dp),intent(in) :: vhtnzc(:) 517 518 !Local variables --------------------------------------- 519 integer :: il,ilm,iln,ilmn,j0lmn,jl,jlm,jln,jlmn,klmn,lmn2_size 520 integer :: meshsz,meshsz_core,meshsz_vhtnzc,meshsz_vmh 521 real(dp) :: intg,intvh,yp1,ypn 522 real(dp),allocatable :: ff(:),r2k(:),shpf(:),vhnzc(:),vhtnzc_sph(:),work1(:),work2(:) 523 524 ! ********************************************************************* 525 526 lmn2_size=pawtab%lmn2_size 527 meshsz_vhtnzc=size(vhtnzc) 528 meshsz=min(radmesh%mesh_size,radmesh_core%mesh_size,radmesh_vloc%mesh_size,meshsz_vhtnzc) 529 LIBPAW_ALLOCATE(ff,(meshsz)) 530 531 !Retrieve VH(tnZc) on the correct radial mesh 532 LIBPAW_ALLOCATE(vhtnzc_sph,(meshsz)) 533 if ((radmesh%mesh_type/=radmesh_vloc%mesh_type).or.& 534 & (radmesh%rstep /=radmesh_vloc%rstep) .or.& 535 & (radmesh%lstep /=radmesh_vloc%lstep)) then 536 call bound_deriv(vhtnzc,radmesh_vloc,meshsz_vhtnzc,yp1,ypn) 537 LIBPAW_ALLOCATE(work1,(meshsz_vhtnzc)) 538 LIBPAW_ALLOCATE(work2,(meshsz_vhtnzc)) 539 call paw_spline(radmesh_vloc%rad,vhtnzc,meshsz_vhtnzc,yp1,ypn,work1) 540 call paw_splint(meshsz_vhtnzc,radmesh_vloc%rad,vhtnzc,work1,meshsz,radmesh%rad(1:meshsz),vhtnzc_sph) 541 LIBPAW_DEALLOCATE(work1) 542 LIBPAW_DEALLOCATE(work2) 543 else 544 vhtnzc_sph(1:meshsz)=vhtnzc(1:meshsz) 545 end if 546 547 !Initialize Kij with Dij0 548 !========================= 549 kij(1:lmn2_size)=pawtab%dij0(1:lmn2_size) 550 551 !Substraction of -<phi_i|vh(nZc)|phi_j> on the PAW sphere 552 !======================================================== 553 if (opt_vhnzc/=0) then 554 meshsz_core=size(ncore) 555 LIBPAW_ALLOCATE(vhnzc,(meshsz_core)) 556 call atompaw_vhnzc(ncore,radmesh_core,vhnzc,znucl) 557 do jlmn=1,pawtab%lmn_size 558 j0lmn=jlmn*(jlmn-1)/2 559 jlm=indlmn(4,jlmn);jln=indlmn(5,jlmn) 560 do ilmn=1,jlmn 561 klmn=j0lmn+ilmn 562 ilm=indlmn(4,ilmn);iln=indlmn(5,ilmn) 563 if (jlm==ilm) then 564 ff(1:meshsz)=pawtab%phi(1:meshsz,iln)*pawtab%phi(1:meshsz,jln)*vhnzc(1:meshsz) 565 call simp_gen(intg,ff,radmesh) 566 kij(klmn)=kij(klmn)-intg 567 end if 568 end do 569 end do 570 LIBPAW_DEALLOCATE(vhnzc) 571 end if 572 573 !Substraction of <tphi_i|vh(tnZc)|tphi_j> on the PAW sphere 574 !========================================================== 575 do jlmn=1,pawtab%lmn_size 576 j0lmn=jlmn*(jlmn-1)/2 577 jlm=indlmn(4,jlmn);jln=indlmn(5,jlmn) 578 do ilmn=1,jlmn 579 klmn=j0lmn+ilmn 580 ilm=indlmn(4,ilmn);iln=indlmn(5,ilmn) 581 if (jlm==ilm) then 582 ff(1:meshsz)=pawtab%tphi(1:meshsz,iln)*pawtab%tphi(1:meshsz,jln)*vhtnzc_sph(1:meshsz) 583 call simp_gen(intg,ff,radmesh) 584 kij(klmn)=kij(klmn)+intg 585 end if 586 end do 587 end do 588 589 !Computation of <phi_i|vminushalf|phi_j> (if any) 590 !================================================= 591 if(pawtab%has_vminushalf==1) then 592 if(size(pawtab%vminushalf)>=1) then 593 meshsz_vmh=min(meshsz,size(pawtab%vminushalf)) 594 do jlmn=1,pawtab%lmn_size 595 j0lmn=jlmn*(jlmn-1)/2 596 jlm=indlmn(4,jlmn);jln=indlmn(5,jlmn) 597 do ilmn=1,jlmn 598 klmn=j0lmn+ilmn 599 ilm=indlmn(4,ilmn);iln=indlmn(5,ilmn) 600 if (jlm==ilm) then 601 ff(1:meshsz_vmh)=pawtab%phi(1:meshsz_vmh,iln)*pawtab%phi(1:meshsz_vmh,jln)*pawtab%vminushalf(1:meshsz_vmh) 602 call simp_gen(intg,ff(1:meshsz_vmh),radmesh) 603 kij(klmn)=kij(klmn)-intg 604 end if 605 end do 606 end do 607 end if 608 end if 609 610 !Computation of int[vh(tnzc)*Qijhat(r)dr] 611 !========================================== 612 if (opt_init==0) then 613 LIBPAW_ALLOCATE(shpf,(radmesh%mesh_size)) 614 call atompaw_shpfun(0,radmesh,intg,pawtab,shpf) 615 if (pawtab%shape_type==3) then 616 LIBPAW_ALLOCATE(r2k,(radmesh%int_meshsz)) 617 r2k=zero 618 r2k(2:radmesh%int_meshsz)=shpf(2:radmesh%int_meshsz)*radmesh%rad(2:radmesh%int_meshsz)**2 619 if(radmesh%mesh_type==5) then 620 call simp_gen(intg,r2k,radmesh) 621 else 622 call simp_gen(intg,r2k,radmesh,r_for_intg=pawtab%rshp) 623 end if 624 shpf(1:meshsz)=shpf(1:meshsz)/intg 625 LIBPAW_DEALLOCATE(r2k) 626 end if 627 ff(1:meshsz)=vhtnzc_sph(1:meshsz)*shpf(1:meshsz)*radmesh%rad(1:meshsz)**2 628 LIBPAW_DEALLOCATE(shpf) 629 call simp_gen(intvh,ff,radmesh) 630 do jlmn=1,pawtab%lmn_size 631 j0lmn=jlmn*(jlmn-1)/2 632 jl=indlmn(1,jlmn);jln=indlmn(5,jlmn);jlm=indlmn(4,jlmn) 633 do ilmn=1,jlmn 634 klmn=j0lmn+ilmn 635 il=indlmn(1,ilmn);iln=indlmn(5,ilmn);ilm=indlmn(4,ilmn) 636 if (ilm==jlm) then 637 ff(1:meshsz)=(pawtab%phi (1:meshsz,iln)*pawtab%phi (1:meshsz,jln)& 638 & -pawtab%tphi(1:meshsz,iln)*pawtab%tphi(1:meshsz,jln)) 639 call simp_gen(intg,ff,radmesh) 640 kij(klmn)=kij(klmn)+intvh*intg 641 end if 642 end do 643 end do 644 else 645 ff(1:meshsz)=vhtnzc_sph(1:meshsz)*pawtab%shapefunc(1:meshsz,1)*radmesh%rad(1:meshsz)**2 646 call simp_gen(intvh,ff,radmesh) 647 do jlmn=1,pawtab%lmn_size 648 j0lmn=jlmn*(jlmn-1)/2 649 jl=indlmn(1,jlmn);jln=indlmn(5,jlmn);jlm=indlmn(4,jlmn) 650 do ilmn=1,jlmn 651 klmn=j0lmn+ilmn 652 il=indlmn(1,ilmn);iln=indlmn(5,ilmn);ilm=indlmn(4,ilmn) 653 if (ilm==jlm) then 654 intg=pawtab%qijl(1,klmn)*sqrt(four_pi) 655 kij(klmn)=kij(klmn)+intvh*intg 656 end if 657 end do 658 end do 659 end if 660 661 LIBPAW_DEALLOCATE(ff) 662 LIBPAW_DEALLOCATE(vhtnzc_sph) 663 664 end subroutine atompaw_kij
m_paw_atom/atompaw_shpfun [ Functions ]
[ Top ] [ m_paw_atom ] [ Functions ]
NAME
atompaw_shpfun
FUNCTION
Compute shape function used in the definition of compensation density (PAW)
INPUTS
ll= l quantum number mesh <type(pawrad_type)>=data containing radial grid information pawtab <type(pawtab_type)>=paw tabulated starting data
OUTPUT
norm= factor for shape function normalization SIDE effects shapefunc(:)=shape function g(r) In case of numerical shape function (shape_type=-1), shapefunc array contains the shape function read in psp file at input.
NOTES
Types of shape functions: type -1: numerical shape function, given in psp file type 1: g(r)=k(r).r^l; k(r)=exp(-(r/sigma)^lambda) type 2: g(r)=k(r).r^l; k(r)=[sin(Pi.r/rshp)/(Pi.r/rshp)]^2 type 3: g(r)=alpha1.jl(q1.r)+alpha2.jl(q2.r)
SOURCE
76 subroutine atompaw_shpfun(ll,mesh,norm,pawtab,shapefunc) 77 78 !Arguments --------------------------------------------- 79 !scalars 80 integer,intent(in) :: ll 81 real(dp),intent(out) :: norm 82 type(pawrad_type),intent(in) :: mesh 83 type(pawtab_type),intent(in) :: pawtab 84 !arrays 85 real(dp),intent(inout) :: shapefunc(:) 86 87 !Local variables ------------------------------ 88 !scalars 89 integer :: ir,ishp,mesh_size 90 real(dp) :: arg,besp,bespp,jbes1,jbes2 91 !arrays 92 real(dp) :: alpha(2),qq(2) 93 real(dp),allocatable :: r2k(:) 94 !no_abirules 95 96 !*************************************************************************** 97 98 mesh_size=size(shapefunc) 99 if (mesh_size>mesh%mesh_size) then 100 LIBPAW_BUG('wrong size!') 101 end if 102 103 !Index for shape function cut-off radius 104 ishp=pawrad_ifromr(mesh,pawtab%rshp)-1 105 106 !Computation of non-normalized shape function 107 if (pawtab%shape_type==-1) then 108 shapefunc(1:ishp)=pawtab%shapefunc(1:ishp,1+ll) 109 else if (pawtab%shape_type==1) then 110 if (ll==0) then 111 shapefunc(1)=one 112 do ir=2,ishp 113 shapefunc(ir)=exp(-(mesh%rad(ir)/pawtab%shape_sigma)**pawtab%shape_lambda) 114 end do 115 else 116 shapefunc(1)=zero 117 do ir=2,ishp 118 shapefunc(ir)=exp(-(mesh%rad(ir)/pawtab%shape_sigma)**pawtab%shape_lambda)*mesh%rad(ir)**ll 119 end do 120 end if 121 else if (pawtab%shape_type==2) then 122 if (ll==0) then 123 shapefunc(1)=one 124 do ir=2,ishp 125 arg=pi*mesh%rad(ir)/pawtab%rshp 126 shapefunc(ir)=(sin(arg)/arg)**2 127 end do 128 else 129 shapefunc(1)=zero 130 do ir=2,ishp 131 arg=pi*mesh%rad(ir)/pawtab%rshp 132 shapefunc(ir)=(sin(arg)/(arg))**2 *mesh%rad(ir)**ll 133 end do 134 end if 135 else if (pawtab%shape_type==3) then 136 alpha(1:2)=pawtab%shape_alpha(1:2,1+ll) 137 qq(1:2)=pawtab%shape_q(1:2,1+ll) 138 do ir=1,ishp 139 call paw_jbessel(jbes1,besp,bespp,ll,0,qq(1)*mesh%rad(ir)) 140 call paw_jbessel(jbes2,besp,bespp,ll,0,qq(2)*mesh%rad(ir)) 141 shapefunc(ir)=alpha(1)*jbes1+alpha(2)*jbes2 142 end do 143 end if 144 145 if (ishp<mesh_size) shapefunc(ishp+1:mesh_size)=zero 146 147 !Shape function normalization 148 if (pawtab%shape_type==-1.or.pawtab%shape_type==1.or.pawtab%shape_type==2) then 149 LIBPAW_ALLOCATE(r2k,(mesh_size)) 150 r2k=zero 151 r2k(2:ishp)=shapefunc(2:ishp)*mesh%rad(2:ishp)**(2+ll) 152 if (mesh%mesh_type==5) then 153 call simp_gen(norm,r2k,mesh);norm=one/norm 154 else 155 call simp_gen(norm,r2k,mesh,r_for_intg=pawtab%rshp);norm=one/norm 156 end if 157 shapefunc(1:ishp)=shapefunc(1:ishp)*norm 158 if (pawtab%shape_type==-1) norm=one 159 LIBPAW_DEALLOCATE(r2k) 160 else if (pawtab%shape_type==3) then 161 norm=one 162 end if 163 164 end subroutine atompaw_shpfun
m_paw_atom/atompaw_vhnzc [ Functions ]
[ Top ] [ m_paw_atom ] [ Functions ]
NAME
atompaw_vhnzc
FUNCTION
PAW: compute Hartree potential for n_{Zc}
INPUTS
ncore(:)=atomic core density radmesh_core <type(pawrad_type)>=radial mesh (and related data) for the core densities znucl= valence and total charge of the atomic species
OUTPUT
vhnzc(:)=Hartree potential due to Z_nc
SOURCE
247 subroutine atompaw_vhnzc(ncore,radmesh_core,vhnzc,znucl) 248 249 !Arguments --------------------------------------------- 250 !scalars 251 real(dp),intent(in) :: znucl 252 type(pawrad_type),intent(in) :: radmesh_core 253 !arrays 254 real(dp),intent(in) :: ncore(:) 255 real(dp), intent(out) :: vhnzc(:) 256 257 !Local variables --------------------------------------- 258 integer :: mesh_size 259 real(dp),allocatable :: nwk(:) 260 261 ! ********************************************************************* 262 263 mesh_size=size(ncore) 264 if (mesh_size/=size(vhnzc).or.mesh_size>radmesh_core%mesh_size) then 265 LIBPAW_BUG('wrong sizes!') 266 end if 267 268 LIBPAW_ALLOCATE(nwk,(mesh_size)) 269 270 nwk(:)=ncore(:)*four_pi*radmesh_core%rad(:)**2 271 call poisson(nwk,0,radmesh_core,vhnzc) 272 vhnzc(2:mesh_size)=(vhnzc(2:mesh_size)-znucl)/radmesh_core%rad(2:mesh_size) 273 call pawrad_deducer0(vhnzc,mesh_size,radmesh_core) 274 275 LIBPAW_DEALLOCATE(nwk) 276 277 end subroutine atompaw_vhnzc