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-2018 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
21 #include "libpaw.h" 22 23 module m_paw_atom 24 25 USE_DEFS 26 USE_MSG_HANDLING 27 USE_MEMORY_PROFILING 28 29 use m_paw_numeric, only : paw_jbessel, paw_solvbes, paw_spline, paw_splint 30 use m_pawrad, only : pawrad_type, simp_gen, poisson, pawrad_deducer0, bound_deriv, pawrad_ifromr 31 use m_pawtab, only : pawtab_type 32 33 implicit none 34 35 private 36 37 public:: atompaw_shpfun 38 public:: atompaw_shapebes 39 public:: atompaw_vhnzc 40 public:: atompaw_dij0 41 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
PARENTS
m_pawpsp
CHILDREN
atompaw_shpfun,atompaw_vhnzc,bound_deriv,paw_spline,paw_splint,simp_gen
SOURCE
211 subroutine atompaw_shapebes(al,ql,ll,rc) 212 213 214 !This section has been created automatically by the script Abilint (TD). 215 !Do not modify the following lines by hand. 216 #undef ABI_FUNC 217 #define ABI_FUNC 'atompaw_shapebes' 218 !End of the abilint section 219 220 implicit none 221 222 !Arguments ------------------------------------ 223 !scalars 224 integer :: ll 225 real(dp) :: rc 226 !arrays 227 real(dp) :: al(2),ql(2) 228 229 !Local variables------------------------------- 230 !scalars 231 integer :: ii 232 real(dp) :: alpha,beta,det,jbes,jbesp,jbespp,qr 233 !arrays 234 real(dp) :: amat(2,2),bb(2) 235 236 ! ************************************************************************* 237 238 alpha=1._dp;beta=0._dp 239 call paw_solvbes(ql,alpha,beta,ll,2) 240 ql(1:2)=ql(1:2)/rc 241 242 do ii=1,2 243 qr=ql(ii)*rc 244 call paw_jbessel(jbes,jbesp,jbespp,ll,1,qr) 245 amat(1,ii)=jbesp*ql(ii) 246 call paw_jbessel(jbes,jbesp,jbespp,ll+1,0,qr) 247 amat(2,ii)=jbes*rc**(ll+2)/ql(ii) ! Intg_0_rc[jl(qr).r^(l+2).dr] 248 end do 249 250 bb(1)=zero;bb(2)=one 251 252 det=amat(1,1)*amat(2,2)-amat(1,2)*amat(2,1) 253 al(1)=(amat(2,2)*bb(1)-amat(1,2)*bb(2))/det 254 al(2)=(amat(1,1)*bb(2)-amat(2,1)*bb(1))/det 255 256 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
PARENTS
m_pawpsp
CHILDREN
atompaw_shpfun,atompaw_vhnzc,bound_deriv,paw_spline,paw_splint,simp_gen
SOURCE
362 subroutine atompaw_dij0(indlmn,kij,lmnmax,ncore,opt_init,pawtab,& 363 & radmesh,radmesh_core,radmesh_vloc,vhtnzc,znucl) 364 365 366 !This section has been created automatically by the script Abilint (TD). 367 !Do not modify the following lines by hand. 368 #undef ABI_FUNC 369 #define ABI_FUNC 'atompaw_dij0' 370 !End of the abilint section 371 372 implicit none 373 374 !Arguments --------------------------------------------- 375 !scalars 376 integer,intent(in) :: lmnmax,opt_init 377 real(dp),intent(in) :: znucl 378 type(pawrad_type),intent(in) :: radmesh,radmesh_core,radmesh_vloc 379 type(pawtab_type),intent(inout) :: pawtab 380 !arrays 381 integer,intent(in) :: indlmn(6,lmnmax) 382 real(dp),intent(in) :: kij(pawtab%lmn2_size) 383 real(dp),intent(in) :: ncore(:),vhtnzc(:) 384 !real(dp),optional,intent(in) :: vminushalf(:) 385 386 !Local variables --------------------------------------- 387 integer :: il,ilm,iln,ilmn,j0lmn,jl,jlm,jln,jlmn,klmn,lmn2_size,meshsz,meshsz_core 388 integer :: meshsz_vhtnzc,meshsz_vmh 389 real(dp) :: intg,intvh,yp1,ypn 390 real(dp),allocatable :: ff(:),ff1(:),r2k(:),shpf(:),vhnzc(:),vhtnzc_sph(:),work1(:),work2(:) 391 392 ! ********************************************************************* 393 394 lmn2_size=pawtab%lmn2_size 395 meshsz_vhtnzc=size(vhtnzc) 396 meshsz=min(radmesh%mesh_size,radmesh_core%mesh_size,radmesh_vloc%mesh_size,meshsz_vhtnzc) 397 LIBPAW_ALLOCATE(ff,(meshsz)) 398 399 !Retrieve VH(tnZc) on the correct radial mesh 400 LIBPAW_ALLOCATE(vhtnzc_sph,(meshsz)) 401 if ((radmesh%mesh_type/=radmesh_vloc%mesh_type).or.& 402 & (radmesh%rstep /=radmesh_vloc%rstep) .or.& 403 & (radmesh%lstep /=radmesh_vloc%lstep)) then 404 call bound_deriv(vhtnzc,radmesh_vloc,meshsz_vhtnzc,yp1,ypn) 405 LIBPAW_ALLOCATE(work1,(meshsz_vhtnzc)) 406 LIBPAW_ALLOCATE(work2,(meshsz_vhtnzc)) 407 call paw_spline(radmesh_vloc%rad,vhtnzc,meshsz_vhtnzc,yp1,ypn,work1) 408 call paw_splint(meshsz_vhtnzc,radmesh_vloc%rad,vhtnzc,work1,meshsz,radmesh%rad(1:meshsz),vhtnzc_sph) 409 LIBPAW_DEALLOCATE(work1) 410 LIBPAW_DEALLOCATE(work2) 411 else 412 vhtnzc_sph(1:meshsz)=vhtnzc(1:meshsz) 413 end if 414 415 !Kinetic part of Dij0 416 !==================== 417 pawtab%dij0(1:lmn2_size)=kij(1:lmn2_size) 418 419 !Computation of <phi_i|vh(nZc)|phi_j> on the PAW sphere 420 !====================================================== 421 meshsz_core=size(ncore) 422 LIBPAW_ALLOCATE(vhnzc,(meshsz_core)) 423 call atompaw_vhnzc(ncore,radmesh_core,vhnzc,znucl) 424 do jlmn=1,pawtab%lmn_size 425 j0lmn=jlmn*(jlmn-1)/2 426 jlm=indlmn(4,jlmn);jln=indlmn(5,jlmn) 427 do ilmn=1,jlmn 428 klmn=j0lmn+ilmn 429 ilm=indlmn(4,ilmn);iln=indlmn(5,ilmn) 430 if (jlm==ilm) then 431 ff(1:meshsz)=pawtab%phi(1:meshsz,iln)*pawtab%phi(1:meshsz,jln)*vhnzc(1:meshsz) 432 call simp_gen(intg,ff,radmesh) 433 pawtab%dij0(klmn)=pawtab%dij0(klmn)+intg 434 end if 435 end do 436 end do 437 LIBPAW_DEALLOCATE(vhnzc) 438 439 !Computation of -<tphi_i|vh(tnZc)|tphi_j> on the PAW sphere 440 !========================================================== 441 do jlmn=1,pawtab%lmn_size 442 j0lmn=jlmn*(jlmn-1)/2 443 jlm=indlmn(4,jlmn);jln=indlmn(5,jlmn) 444 do ilmn=1,jlmn 445 klmn=j0lmn+ilmn 446 ilm=indlmn(4,ilmn);iln=indlmn(5,ilmn) 447 if (jlm==ilm) then 448 ff(1:meshsz)=pawtab%tphi(1:meshsz,iln)*pawtab%tphi(1:meshsz,jln)*vhtnzc_sph(1:meshsz) 449 call simp_gen(intg,ff,radmesh) 450 pawtab%dij0(klmn)=pawtab%dij0(klmn)-intg 451 end if 452 end do 453 end do 454 455 !Computation of <phi_i|vminushalf|phi_j> (if any) 456 !================================================= 457 if(pawtab%has_vminushalf==1) then 458 if(size(pawtab%vminushalf)>=1) then 459 meshsz_vmh=min(meshsz,size(pawtab%vminushalf)) 460 do jlmn=1,pawtab%lmn_size 461 j0lmn=jlmn*(jlmn-1)/2 462 jlm=indlmn(4,jlmn);jln=indlmn(5,jlmn) 463 do ilmn=1,jlmn 464 klmn=j0lmn+ilmn 465 ilm=indlmn(4,ilmn);iln=indlmn(5,ilmn) 466 if (jlm==ilm) then 467 ff(1:meshsz_vmh)=pawtab%phi(1:meshsz_vmh,iln)*pawtab%phi(1:meshsz_vmh,jln)*pawtab%vminushalf(1:meshsz_vmh) 468 call simp_gen(intg,ff(1:meshsz_vmh),radmesh) 469 pawtab%dij0(klmn)=pawtab%dij0(klmn)+intg 470 end if 471 end do 472 end do 473 end if 474 end if 475 476 !Computation of -int[vh(tnzc)*Qijhat(r)dr] 477 !========================================= 478 if (opt_init==0) then 479 LIBPAW_ALLOCATE(shpf,(radmesh%mesh_size)) 480 call atompaw_shpfun(0,radmesh,intg,pawtab,shpf) 481 if (pawtab%shape_type==3) then 482 LIBPAW_ALLOCATE(r2k,(radmesh%int_meshsz)) 483 r2k=zero 484 r2k(2:radmesh%int_meshsz)=shpf(2:radmesh%int_meshsz)*radmesh%rad(2:radmesh%int_meshsz)**2 485 if(radmesh%mesh_type==5) then 486 call simp_gen(intg,r2k,radmesh) 487 else 488 call simp_gen(intg,r2k,radmesh,r_for_intg=pawtab%rshp) 489 end if 490 shpf(1:meshsz)=shpf(1:meshsz)/intg 491 LIBPAW_DEALLOCATE(r2k) 492 end if 493 ff(1:meshsz)=vhtnzc_sph(1:meshsz)*shpf(1:meshsz)*radmesh%rad(1:meshsz)**2 494 LIBPAW_DEALLOCATE(shpf) 495 call simp_gen(intvh,ff,radmesh) 496 do jlmn=1,pawtab%lmn_size 497 j0lmn=jlmn*(jlmn-1)/2 498 jl=indlmn(1,jlmn);jln=indlmn(5,jlmn);jlm=indlmn(4,jlmn) 499 do ilmn=1,jlmn 500 klmn=j0lmn+ilmn 501 il=indlmn(1,ilmn);iln=indlmn(5,ilmn);ilm=indlmn(4,ilmn) 502 if (ilm==jlm) then 503 ff(1:meshsz)=(pawtab%phi (1:meshsz,iln)*pawtab%phi (1:meshsz,jln)& 504 & -pawtab%tphi(1:meshsz,iln)*pawtab%tphi(1:meshsz,jln)) 505 call simp_gen(intg,ff,radmesh) 506 pawtab%dij0(klmn)=pawtab%dij0(klmn)-intvh*intg 507 end if 508 end do 509 end do 510 else 511 ff(1:meshsz)=vhtnzc_sph(1:meshsz)*pawtab%shapefunc(1:meshsz,1)*radmesh%rad(1:meshsz)**2 512 call simp_gen(intvh,ff,radmesh) 513 do jlmn=1,pawtab%lmn_size 514 j0lmn=jlmn*(jlmn-1)/2 515 jl=indlmn(1,jlmn);jln=indlmn(5,jlmn);jlm=indlmn(4,jlmn) 516 do ilmn=1,jlmn 517 klmn=j0lmn+ilmn 518 il=indlmn(1,ilmn);iln=indlmn(5,ilmn);ilm=indlmn(4,ilmn) 519 if (ilm==jlm) then 520 intg=pawtab%qijl(1,klmn)*sqrt(four_pi) 521 pawtab%dij0(klmn)=pawtab%dij0(klmn)-intvh*intg 522 end if 523 end do 524 end do 525 end if 526 527 LIBPAW_DEALLOCATE(ff) 528 LIBPAW_DEALLOCATE(vhtnzc_sph) 529 530 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
PARENTS
m_pawpsp
CHILDREN
atompaw_shpfun,atompaw_vhnzc,bound_deriv,paw_spline,paw_splint,simp_gen
SOURCE
570 subroutine atompaw_kij(indlmn,kij,lmnmax,ncore,opt_init,opt_vhnzc,pawtab, & 571 & radmesh,radmesh_core,radmesh_vloc,vhtnzc,znucl) 572 573 574 !This section has been created automatically by the script Abilint (TD). 575 !Do not modify the following lines by hand. 576 #undef ABI_FUNC 577 #define ABI_FUNC 'atompaw_kij' 578 !End of the abilint section 579 580 implicit none 581 582 !Arguments --------------------------------------------- 583 !scalars 584 integer,intent(in) :: lmnmax,opt_init,opt_vhnzc 585 real(dp),intent(in) :: znucl 586 type(pawrad_type),intent(in) :: radmesh,radmesh_core,radmesh_vloc 587 type(pawtab_type),intent(in) :: pawtab 588 !arrays 589 integer,intent(in) :: indlmn(6,lmnmax) 590 real(dp),intent(out) :: kij(pawtab%lmn2_size) 591 real(dp),intent(in) :: ncore(:) 592 real(dp),intent(in) :: vhtnzc(:) 593 594 !Local variables --------------------------------------- 595 integer :: il,ilm,iln,ilmn,j0lmn,jl,jlm,jln,jlmn,klmn,lmn2_size 596 integer :: meshsz,meshsz_core,meshsz_vhtnzc,meshsz_vmh 597 real(dp) :: intg,intvh,yp1,ypn 598 real(dp),allocatable :: ff(:),r2k(:),shpf(:),vhnzc(:),vhtnzc_sph(:),work1(:),work2(:) 599 600 ! ********************************************************************* 601 602 lmn2_size=pawtab%lmn2_size 603 meshsz_vhtnzc=size(vhtnzc) 604 meshsz=min(radmesh%mesh_size,radmesh_core%mesh_size,radmesh_vloc%mesh_size,meshsz_vhtnzc) 605 LIBPAW_ALLOCATE(ff,(meshsz)) 606 607 !Retrieve VH(tnZc) on the correct radial mesh 608 LIBPAW_ALLOCATE(vhtnzc_sph,(meshsz)) 609 if ((radmesh%mesh_type/=radmesh_vloc%mesh_type).or.& 610 & (radmesh%rstep /=radmesh_vloc%rstep) .or.& 611 & (radmesh%lstep /=radmesh_vloc%lstep)) then 612 call bound_deriv(vhtnzc,radmesh_vloc,meshsz_vhtnzc,yp1,ypn) 613 LIBPAW_ALLOCATE(work1,(meshsz_vhtnzc)) 614 LIBPAW_ALLOCATE(work2,(meshsz_vhtnzc)) 615 call paw_spline(radmesh_vloc%rad,vhtnzc,meshsz_vhtnzc,yp1,ypn,work1) 616 call paw_splint(meshsz_vhtnzc,radmesh_vloc%rad,vhtnzc,work1,meshsz,radmesh%rad(1:meshsz),vhtnzc_sph) 617 LIBPAW_DEALLOCATE(work1) 618 LIBPAW_DEALLOCATE(work2) 619 else 620 vhtnzc_sph(1:meshsz)=vhtnzc(1:meshsz) 621 end if 622 623 !Initiialize Kij with Dij0 624 !========================= 625 kij(1:lmn2_size)=pawtab%dij0(1:lmn2_size) 626 627 !Substraction of -<phi_i|vh(nZc)|phi_j> on the PAW sphere 628 !======================================================== 629 if (opt_vhnzc/=0) then 630 meshsz_core=size(ncore) 631 LIBPAW_ALLOCATE(vhnzc,(meshsz_core)) 632 call atompaw_vhnzc(ncore,radmesh_core,vhnzc,znucl) 633 do jlmn=1,pawtab%lmn_size 634 j0lmn=jlmn*(jlmn-1)/2 635 jlm=indlmn(4,jlmn);jln=indlmn(5,jlmn) 636 do ilmn=1,jlmn 637 klmn=j0lmn+ilmn 638 ilm=indlmn(4,ilmn);iln=indlmn(5,ilmn) 639 if (jlm==ilm) then 640 ff(1:meshsz)=pawtab%phi(1:meshsz,iln)*pawtab%phi(1:meshsz,jln)*vhnzc(1:meshsz) 641 call simp_gen(intg,ff,radmesh) 642 kij(klmn)=kij(klmn)-intg 643 end if 644 end do 645 end do 646 LIBPAW_DEALLOCATE(vhnzc) 647 end if 648 649 !Substraction of <tphi_i|vh(tnZc)|tphi_j> on the PAW sphere 650 !========================================================== 651 do jlmn=1,pawtab%lmn_size 652 j0lmn=jlmn*(jlmn-1)/2 653 jlm=indlmn(4,jlmn);jln=indlmn(5,jlmn) 654 do ilmn=1,jlmn 655 klmn=j0lmn+ilmn 656 ilm=indlmn(4,ilmn);iln=indlmn(5,ilmn) 657 if (jlm==ilm) then 658 ff(1:meshsz)=pawtab%tphi(1:meshsz,iln)*pawtab%tphi(1:meshsz,jln)*vhtnzc_sph(1:meshsz) 659 call simp_gen(intg,ff,radmesh) 660 kij(klmn)=kij(klmn)+intg 661 end if 662 end do 663 end do 664 665 !Computation of <phi_i|vminushalf|phi_j> (if any) 666 !================================================= 667 if(pawtab%has_vminushalf==1) then 668 if(size(pawtab%vminushalf)>=1) then 669 meshsz_vmh=min(meshsz,size(pawtab%vminushalf)) 670 do jlmn=1,pawtab%lmn_size 671 j0lmn=jlmn*(jlmn-1)/2 672 jlm=indlmn(4,jlmn);jln=indlmn(5,jlmn) 673 do ilmn=1,jlmn 674 klmn=j0lmn+ilmn 675 ilm=indlmn(4,ilmn);iln=indlmn(5,ilmn) 676 if (jlm==ilm) then 677 ff(1:meshsz_vmh)=pawtab%phi(1:meshsz_vmh,iln)*pawtab%phi(1:meshsz_vmh,jln)*pawtab%vminushalf(1:meshsz_vmh) 678 call simp_gen(intg,ff(1:meshsz_vmh),radmesh) 679 kij(klmn)=kij(klmn)-intg 680 end if 681 end do 682 end do 683 end if 684 end if 685 686 !Computation of int[vh(tnzc)*Qijhat(r)dr] 687 !========================================== 688 if (opt_init==0) then 689 LIBPAW_ALLOCATE(shpf,(radmesh%mesh_size)) 690 call atompaw_shpfun(0,radmesh,intg,pawtab,shpf) 691 if (pawtab%shape_type==3) then 692 LIBPAW_ALLOCATE(r2k,(radmesh%int_meshsz)) 693 r2k=zero 694 r2k(2:radmesh%int_meshsz)=shpf(2:radmesh%int_meshsz)*radmesh%rad(2:radmesh%int_meshsz)**2 695 if(radmesh%mesh_type==5) then 696 call simp_gen(intg,r2k,radmesh) 697 else 698 call simp_gen(intg,r2k,radmesh,r_for_intg=pawtab%rshp) 699 end if 700 shpf(1:meshsz)=shpf(1:meshsz)/intg 701 LIBPAW_DEALLOCATE(r2k) 702 end if 703 ff(1:meshsz)=vhtnzc_sph(1:meshsz)*shpf(1:meshsz)*radmesh%rad(1:meshsz)**2 704 LIBPAW_DEALLOCATE(shpf) 705 call simp_gen(intvh,ff,radmesh) 706 do jlmn=1,pawtab%lmn_size 707 j0lmn=jlmn*(jlmn-1)/2 708 jl=indlmn(1,jlmn);jln=indlmn(5,jlmn);jlm=indlmn(4,jlmn) 709 do ilmn=1,jlmn 710 klmn=j0lmn+ilmn 711 il=indlmn(1,ilmn);iln=indlmn(5,ilmn);ilm=indlmn(4,ilmn) 712 if (ilm==jlm) then 713 ff(1:meshsz)=(pawtab%phi (1:meshsz,iln)*pawtab%phi (1:meshsz,jln)& 714 & -pawtab%tphi(1:meshsz,iln)*pawtab%tphi(1:meshsz,jln)) 715 call simp_gen(intg,ff,radmesh) 716 kij(klmn)=kij(klmn)+intvh*intg 717 end if 718 end do 719 end do 720 else 721 ff(1:meshsz)=vhtnzc_sph(1:meshsz)*pawtab%shapefunc(1:meshsz,1)*radmesh%rad(1:meshsz)**2 722 call simp_gen(intvh,ff,radmesh) 723 do jlmn=1,pawtab%lmn_size 724 j0lmn=jlmn*(jlmn-1)/2 725 jl=indlmn(1,jlmn);jln=indlmn(5,jlmn);jlm=indlmn(4,jlmn) 726 do ilmn=1,jlmn 727 klmn=j0lmn+ilmn 728 il=indlmn(1,ilmn);iln=indlmn(5,ilmn);ilm=indlmn(4,ilmn) 729 if (ilm==jlm) then 730 intg=pawtab%qijl(1,klmn)*sqrt(four_pi) 731 kij(klmn)=kij(klmn)+intvh*intg 732 end if 733 end do 734 end do 735 end if 736 737 LIBPAW_DEALLOCATE(ff) 738 LIBPAW_DEALLOCATE(vhtnzc_sph) 739 740 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)
PARENTS
m_paw_atom,m_pawpsp,pawinit
CHILDREN
atompaw_shpfun,atompaw_vhnzc,bound_deriv,paw_spline,paw_splint,simp_gen
SOURCE
83 subroutine atompaw_shpfun(ll,mesh,norm,pawtab,shapefunc) 84 85 86 !This section has been created automatically by the script Abilint (TD). 87 !Do not modify the following lines by hand. 88 #undef ABI_FUNC 89 #define ABI_FUNC 'atompaw_shpfun' 90 !End of the abilint section 91 92 implicit none 93 94 !Arguments --------------------------------------------- 95 !scalars 96 integer,intent(in) :: ll 97 real(dp),intent(out) :: norm 98 type(pawrad_type),intent(in) :: mesh 99 type(pawtab_type),intent(in) :: pawtab 100 !arrays 101 real(dp),intent(inout) :: shapefunc(:) 102 103 !Local variables ------------------------------ 104 !scalars 105 integer :: ir,ishp,mesh_size 106 real(dp) :: arg,besp,bespp,jbes1,jbes2 107 !arrays 108 real(dp) :: alpha(2),qq(2) 109 real(dp),allocatable :: r2k(:) 110 !no_abirules 111 112 !*************************************************************************** 113 114 mesh_size=size(shapefunc) 115 if (mesh_size>mesh%mesh_size) then 116 MSG_BUG('wrong size!') 117 end if 118 119 !Index for shape function cut-off radius 120 ishp=pawrad_ifromr(mesh,pawtab%rshp)-1 121 122 !Computation of non-normalized shape function 123 if (pawtab%shape_type==-1) then 124 shapefunc(1:ishp)=pawtab%shapefunc(1:ishp,1+ll) 125 else if (pawtab%shape_type==1) then 126 if (ll==0) then 127 shapefunc(1)=one 128 do ir=2,ishp 129 shapefunc(ir)=exp(-(mesh%rad(ir)/pawtab%shape_sigma)**pawtab%shape_lambda) 130 end do 131 else 132 shapefunc(1)=zero 133 do ir=2,ishp 134 shapefunc(ir)=exp(-(mesh%rad(ir)/pawtab%shape_sigma)**pawtab%shape_lambda)*mesh%rad(ir)**ll 135 end do 136 end if 137 else if (pawtab%shape_type==2) then 138 if (ll==0) then 139 shapefunc(1)=one 140 do ir=2,ishp 141 arg=pi*mesh%rad(ir)/pawtab%rshp 142 shapefunc(ir)=(sin(arg)/arg)**2 143 end do 144 else 145 shapefunc(1)=zero 146 do ir=2,ishp 147 arg=pi*mesh%rad(ir)/pawtab%rshp 148 shapefunc(ir)=(sin(arg)/(arg))**2 *mesh%rad(ir)**ll 149 end do 150 end if 151 else if (pawtab%shape_type==3) then 152 alpha(1:2)=pawtab%shape_alpha(1:2,1+ll) 153 qq(1:2)=pawtab%shape_q(1:2,1+ll) 154 do ir=1,ishp 155 call paw_jbessel(jbes1,besp,bespp,ll,0,qq(1)*mesh%rad(ir)) 156 call paw_jbessel(jbes2,besp,bespp,ll,0,qq(2)*mesh%rad(ir)) 157 shapefunc(ir)=alpha(1)*jbes1+alpha(2)*jbes2 158 end do 159 end if 160 161 if (ishp<mesh_size) shapefunc(ishp+1:mesh_size)=zero 162 163 !Shape function normalization 164 if (pawtab%shape_type==-1.or.pawtab%shape_type==1.or.pawtab%shape_type==2) then 165 LIBPAW_ALLOCATE(r2k,(mesh_size)) 166 r2k=zero 167 r2k(2:ishp)=shapefunc(2:ishp)*mesh%rad(2:ishp)**(2+ll) 168 if (mesh%mesh_type==5) then 169 call simp_gen(norm,r2k,mesh);norm=one/norm 170 else 171 call simp_gen(norm,r2k,mesh,r_for_intg=pawtab%rshp);norm=one/norm 172 end if 173 shapefunc(1:ishp)=shapefunc(1:ishp)*norm 174 if (pawtab%shape_type==-1) norm=one 175 LIBPAW_DEALLOCATE(r2k) 176 else if (pawtab%shape_type==3) then 177 norm=one 178 end if 179 180 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
PARENTS
m_paw_atom,m_pawpsp
CHILDREN
atompaw_shpfun,atompaw_vhnzc,bound_deriv,paw_spline,paw_splint,simp_gen
SOURCE
284 subroutine atompaw_vhnzc(ncore,radmesh_core,vhnzc,znucl) 285 286 287 !This section has been created automatically by the script Abilint (TD). 288 !Do not modify the following lines by hand. 289 #undef ABI_FUNC 290 #define ABI_FUNC 'atompaw_vhnzc' 291 !End of the abilint section 292 293 implicit none 294 295 !Arguments --------------------------------------------- 296 !scalars 297 real(dp),intent(in) :: znucl 298 type(pawrad_type),intent(in) :: radmesh_core 299 !arrays 300 real(dp),intent(in) :: ncore(:) 301 real(dp), intent(out) :: vhnzc(:) 302 303 !Local variables --------------------------------------- 304 integer :: mesh_size 305 real(dp),allocatable :: nwk(:) 306 307 ! ********************************************************************* 308 309 mesh_size=size(ncore) 310 if (mesh_size/=size(vhnzc).or.mesh_size>radmesh_core%mesh_size) then 311 MSG_BUG('wrong sizes!') 312 end if 313 314 LIBPAW_ALLOCATE(nwk,(mesh_size)) 315 316 nwk(:)=ncore(:)*four_pi*radmesh_core%rad(:)**2 317 call poisson(nwk,0,radmesh_core,vhnzc) 318 vhnzc(2:mesh_size)=(vhnzc(2:mesh_size)-znucl)/radmesh_core%rad(2:mesh_size) 319 call pawrad_deducer0(vhnzc,mesh_size,radmesh_core) 320 321 LIBPAW_DEALLOCATE(nwk) 322 323 end subroutine atompaw_vhnzc