TABLE OF CONTENTS
- ABINIT/m_psp_hgh
- ABINIT/psp2lo
- ABINIT/psp3in
- m_psp_hgh/psp10in
- m_psp_hgh/psp10nl
- m_psp_hgh/psp2in
- m_psp_hgh/psp2nl
- m_psp_hgh/psp3nl
ABINIT/m_psp_hgh [ Modules ]
NAME
m_psp_hgh
FUNCTION
Initialize pspcod=2, 3, 10 pseudopotentials (GTH)
COPYRIGHT
Copyright (C) 1998-2022 ABINIT group (DCA, XG, GMR, MT, 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 .
SOURCE
16 #if defined HAVE_CONFIG_H 17 #include "config.h" 18 #endif 19 20 #include "abi_common.h" 21 22 module m_psp_hgh 23 24 use defs_basis 25 use m_splines 26 use m_abicore 27 use m_errors 28 use m_dtset 29 30 use defs_datatypes, only : pseudopotential_type 31 use m_special_funcs, only : abi_derfc 32 #if defined HAVE_BIGDFT 33 use BigDFT_API, only: atomic_info 34 #endif 35 use m_wvl_descr_psp, only : wvl_descr_psp_fill 36 37 implicit none 38 39 private
ABINIT/psp2lo [ Functions ]
NAME
psp2lo
FUNCTION
Treat local part of Goedecker-Teter-Hutter pseudopotentials (pspcod=2), as well as Hartwigsen-Goedecker-Hutter pseudopotentials (pspcod=3)
INPUTS
cc1,2,3,4=parameters from analytic pseudopotential form mqgrid=number of grid points in q from 0 to qmax. qgrid(mqgrid)=values of q (or G) on grid from 0 to qmax (bohr^-1) if vlspl_recipSpace is .true. else values of r on grid from 0 to 2pi / qmax * mqgrid_ff (bohr). rloc=local pseudopotential core radius (bohr) vlspl_recipSpace= .true. if computation of vlspl is done in reciprocal space zion=valence charge of atom parameters for local potential: rloc,c1,c2,c3,c4
OUTPUT
dvloc(mqgrid)=dVloc(r)/dr (only allocated if vlspl_recipSpace is false). epsatm=$4\pi\int[r^2 (v(r)+\frac{Zv}{r} dr]$ q2vq(mqgrid)&=&q^2 v(q) \nonumber \\ &=&-Zv/\pi +q^2 4\pi\int[(\frac{\sin(2\pi qr)}{2\pi qr})(r^2 v(r)+r Zv)dr]\nonumber\\ &=&\exp(-K^2*rloc^2/2) \nonumber \\ && *(-\frac{zion}{\pi}+(\frac{K^2*rloc^3}{\sqrt{2*\pi}}* (c1+c2*(3-(rloc*K)^2) \nonumber \\ && +c3*(15-10(rloc*K)^2+(rloc*K)^4) \nonumber \\ && +c4*(105-105*(rloc*K)^2+21*(rloc*K)^4-(rloc*K)^6)) \nonumber for GTH vloc with $K=(2\pi q)$. yp1,ypn=derivative of q^2 v(q) wrt q at q=0 and q=qmax (needed for spline fitter).
SOURCE
420 subroutine psp2lo(cc1,cc2,cc3,cc4,dvloc,epsatm,mqgrid,qgrid,q2vq,& 421 & rloc,vlspl_recipSpace,yp1,ypn,zion) 422 423 !Arguments ------------------------------------ 424 !scalars 425 integer,intent(in) :: mqgrid 426 real(dp),intent(in) :: cc1,cc2,cc3,cc4,rloc,zion 427 real(dp),intent(out) :: epsatm,yp1,ypn 428 logical,intent(in) :: vlspl_recipSpace 429 !arrays 430 real(dp),intent(in) :: qgrid(mqgrid) 431 real(dp),intent(out) :: dvloc(mqgrid),q2vq(mqgrid) 432 433 !Local variables------------------------------- 434 !scalars 435 integer :: iqgrid 436 real(dp) :: erfValue,gaussValue,polyValue,qmax,rq,rq2 437 character(len=500) :: message 438 439 ! ************************************************************************* 440 441 !Compute epsatm = lim(q->0) [Vloc(q) + zion/(Pi*q^2)] 442 epsatm=2.d0*pi*rloc**2*zion+(2.d0*pi)**(1.5d0)*rloc**3*& 443 & (cc1+3.d0*cc2+15.d0*cc3+105.d0*cc4) 444 445 !If vlspl_recipSpace is .true., we compute V(q)*q^2 in reciprocal space, 446 !else we compute V(r) in real space. 447 if (vlspl_recipSpace) then 448 write(message, '(a)' ) '- Local part computed in reciprocal space.' 449 call wrtout(ab_out,message,'COLL') 450 call wrtout(std_out, message,'COLL') 451 452 ! d(q^2*V(q))/d(q) at q=0 and q=qmax 453 qmax=qgrid(mqgrid) 454 rq2=(2.d0*pi*qmax*rloc)**2 455 yp1=0.d0 456 ypn= (2.d0*pi*qmax*rloc**2)*exp(-0.5d0*rq2)* & 457 & (2.d0*zion + sqrt(2.d0*pi)*rloc*& 458 & (cc1*(2.d0-rq2) + cc2*(6.d0-7.d0*rq2+rq2**2) +& 459 & cc3*(30.d0-55.d0*rq2+16.d0*rq2**2-rq2**3) +& 460 & cc4*(210.d0-525.d0*rq2+231.d0*rq2**2-29.d0*rq2**3+rq2**4))) 461 ! ypn has been tested against Maple-derived expression. 462 463 ! Compute q^2*vloc(q) on uniform grid 464 do iqgrid=1,mqgrid 465 rq2=(2.d0*pi*qgrid(iqgrid)*rloc)**2 466 q2vq(iqgrid)=exp(-0.5d0*rq2)*(-zion/pi+rq2*(rloc/sqrt(2.d0*pi)) *& 467 & ( cc1 + cc2*(3.d0-rq2) + cc3*(15.d0-10.d0*rq2+rq2**2) +& 468 & cc4*(105.d0-rq2*(105.d0-rq2*(21.d0-rq2))) )) 469 end do 470 else 471 write(message, '(a)' ) '- Local part computed in real space.' 472 call wrtout(ab_out,message,'COLL') 473 call wrtout(std_out, message,'COLL') 474 475 ! Compute derivatives for splines computations 476 yp1 = 0.d0 477 rq2 = (qgrid(mqgrid) / rloc) ** 2 478 erfValue = abi_derfc(sqrt(0.5d0 * rq2)) 479 ypn = - 2.0d0 * zion / sqrt(2.d0 * pi) / qgrid(mqgrid) / rloc 480 ypn = ypn - rq2 * (cc1 + cc2 * rq2 + cc3 * rq2 ** 2 + cc4 * rq2 ** 3) / qgrid(mqgrid) 481 ypn = ypn + (2.d0 * cc2 * rq2 + 4.d0 * cc3 * rq2 ** 2 + 6.d0 * cc4 * rq2 ** 3) / qgrid(mqgrid) 482 ypn = ypn * exp(-0.5d0 * rq2) 483 ypn = ypn + zion / qgrid(mqgrid) ** 2 * erfValue 484 ! Note that ypn has been calculated on a full-proof a4 paper sheet. 485 486 ! Compute local potential and its first derivatives. 487 do iqgrid = 1, mqgrid, 1 488 rq2 = (qgrid(iqgrid) / rloc) ** 2 489 ! Compute erf() part 490 ! Case r = 0 491 gaussValue = exp(-0.5d0 * rq2) 492 if (qgrid(iqgrid) == 0.d0) then 493 q2vq(iqgrid) = -zion / rloc * sqrt(2.d0 / pi) 494 dvloc(iqgrid) = 0.d0 495 else 496 erfValue = abi_derfc(sqrt(0.5d0 * rq2)) 497 q2vq(iqgrid) = -zion / qgrid(iqgrid) * (1.0d0 - erfValue) 498 dvloc(iqgrid) = - sqrt(2.d0 / pi) * zion * gaussValue / (qgrid(iqgrid) * rloc) - & 499 & q2vq(iqgrid) / qgrid(iqgrid) 500 end if 501 ! Add the gaussian part 502 polyValue = cc1 + cc2 * rq2 + cc3 * rq2 ** 2 + cc4 * rq2 ** 3 503 q2vq(iqgrid) = q2vq(iqgrid) + gaussValue * polyValue 504 rq = qgrid(iqgrid) / rloc 505 dvloc(iqgrid) = dvloc(iqgrid) - qgrid(iqgrid) / rloc ** 2 * gaussValue * polyValue + & 506 & gaussValue * (2.0d0 * cc2 * rq / rloc + 3.0d0 * cc3 * rq ** 3 / rloc + & 507 & 6.0d0 * cc4 * rq ** 5 / rloc) 508 end do 509 510 write(message, '(a,f12.7,a,a,f12.7,a,a,a,f12.7)' ) & 511 & ' | dr spline step is : ', qgrid(2), ch10, & 512 & ' | r > ', qgrid(mqgrid) ,' is set to 0.', ch10, & 513 & ' | last non-nul potential value is : ', q2vq(mqgrid) 514 call wrtout(ab_out,message,'COLL') 515 call wrtout(std_out, message,'COLL') 516 end if 517 518 end subroutine psp2lo
ABINIT/psp3in [ Functions ]
NAME
psp3in
FUNCTION
Initialize pspcod=3 pseudopotentials (HGH psps PRB58,3641(1998) [[cite:Hartwigsen1998]]): continue to read the file, then compute the corresponding local and non-local potentials.
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
SOURCE
575 subroutine psp3in(dtset, ekb, epsatm, ffspl, indlmn, ipsp, lmax, nproj, psps, pspso, vlspl, zion) 576 577 !Arguments ------------------------------------ 578 !scalars 579 integer,intent(in) :: ipsp,pspso 580 integer,intent(inout) :: lmax 581 real(dp),intent(in) :: zion 582 real(dp),intent(out) :: epsatm 583 type(dataset_type),intent(in) :: dtset 584 type(pseudopotential_type),intent(inout) :: psps 585 !arrays 586 integer,intent(out) :: indlmn(6,psps%lmnmax),nproj(psps%mpssoang) 587 real(dp),intent(inout) :: ekb(psps%lnmax),ffspl(psps%mqgrid_ff,2,psps%lnmax)!vz_i 588 real(dp),intent(out) :: vlspl(psps%mqgrid_ff,2) 589 590 !Local variables------------------------------- 591 !scalars 592 integer :: iln,iln0,index,ipsang,jj,kk,ll,mm,mproj,nn,nso 593 real(dp) :: cc1,cc2,cc3,cc4,h11d,h11f,h11p,h11s,h22d,h22p,h22s,h33d,h33p,h33s 594 real(dp) :: k11d,k11f,k11p,k22d,k22p,k33d,k33p,rloc 595 real(dp) :: rrd,rrf,rrp,rrs,yp1,ypn 596 character(len=500) :: message,errmsg 597 !arrays 598 real(dp),allocatable :: dvlspl(:),ekb_so(:,:),ekb_sr(:,:),ffspl_so(:,:,:,:) 599 real(dp),allocatable :: ffspl_sr(:,:,:,:),work_space(:),work_spl(:) 600 601 ! *************************************************************************** 602 603 !Set various terms to 0 in case not defined below 604 !HGH values 605 rloc=zero ; rrs=zero ; h11p=zero ; k33p=zero ; k11d=zero; 606 cc1=zero ; h11s=zero ; h22p=zero ; rrd=zero ; k22d=zero; 607 cc2=zero ; h22s=zero ; h33p=zero ; h11d=zero ; k33d=zero; 608 cc3=zero ; h33s=zero ; k11p=zero ; h22d=zero ; h11f=zero; 609 cc4=zero ; rrp=zero ; k22p=zero ; h33d=zero ; k11f=zero; 610 rrf=zero 611 nproj(1:psps%mpssoang)=0 612 613 !Read and write different lines of the pseudopotential file 614 615 read (tmp_unit,*,err=10,iomsg=errmsg) rloc,cc1,cc2,cc3,cc4 616 write(message, '(a,f12.7)' ) ' rloc=',rloc 617 call wrtout(ab_out,message,'COLL') 618 call wrtout(std_out, message,'COLL') 619 write(message, '(a,f12.7,a,f12.7,a,f12.7,a,f12.7)' )' cc1 =',cc1,'; cc2 =',cc2,'; cc3 =',cc3,'; cc4 =',cc4 620 call wrtout(ab_out,message,'COLL') 621 call wrtout(std_out, message,'COLL') 622 623 !For the time being, the s state line must be present and is read, 624 !even for local pseudopotentials (zero must appear) 625 read (tmp_unit,*,err=10,iomsg=errmsg) rrs,h11s,h22s,h33s 626 write(message, '(a,f12.7,a,f12.7,a,f12.7,a,f12.7)' )' rrs =',rrs,'; h11s=',h11s,'; h22s=',h22s,'; h33s=',h33s 627 call wrtout(ab_out,message,'COLL') 628 call wrtout(std_out, message,'COLL') 629 630 if (lmax > 0) then 631 632 read (tmp_unit,*,err=10,iomsg=errmsg) rrp,h11p,h22p,h33p 633 write(message, '(a,f12.7,a,f12.7,a,f12.7,a,f12.7)' )' rrp =',rrp,'; h11p=',h11p,'; h22p=',h22p,'; h33p=',h33p 634 call wrtout(ab_out,message,'COLL') 635 call wrtout(std_out, message,'COLL') 636 637 read (tmp_unit,*,err=10,iomsg=errmsg) k11p,k22p,k33p 638 write(message, '(a,f12.7,a,f12.7,a,f12.7)' )' k11p=',k11p,'; k22p=',k22p,'; k33p=',k33p 639 call wrtout(ab_out,message,'COLL') 640 call wrtout(std_out, message,'COLL') 641 642 end if 643 644 if (lmax > 1) then 645 read (tmp_unit,*,err=10,iomsg=errmsg) rrd,h11d,h22d,h33d 646 write(message, '(a,f12.7,a,f12.7,a,f12.7,a,f12.7)' )' rrd =',rrd,'; h11d=',h11d,'; h22d=',h22d,'; h33d=',h33d 647 call wrtout(ab_out,message,'COLL') 648 call wrtout(std_out, message,'COLL') 649 650 read (tmp_unit,*,err=10,iomsg=errmsg) k11d,k22d,k33d 651 write(message, '(a,f12.7,a,f12.7,a,f12.7)' )' k11d=',k11d,'; k22d=',k22d,'; k33d=',k33d 652 call wrtout(ab_out,message,'COLL') 653 call wrtout(std_out, message,'COLL') 654 end if 655 656 if (lmax > 2) then 657 read (tmp_unit,*,err=10,iomsg=errmsg) rrf,h11f 658 write(message, '(a,f12.7,a,f12.7)' )' rrf =',rrf,'; h11f=',h11f 659 call wrtout(ab_out,message,'COLL') 660 call wrtout(std_out, message,'COLL') 661 662 read (tmp_unit,*,err=10,iomsg=errmsg) k11f 663 write(message, '(a,f12.7)' )' k11f=',k11f 664 call wrtout(ab_out,message,'COLL') 665 call wrtout(std_out, message,'COLL') 666 end if 667 668 if (abs(h11s)>1.d-08) nproj(1)=1 669 if (abs(h22s)>1.d-08) nproj(1)=2 670 if (abs(h33s)>1.d-08) nproj(1)=3 671 672 if (abs(h11p)>1.d-08) then 673 nproj(2)=1 674 if (lmax<1) then 675 write(message, '(a,a,a,a,i5,a,e12.4,a,a)' ) ch10,& 676 & ' psp3in : COMMENT -',ch10,& 677 & ' input lmax=',lmax,' does not agree with input h11p=',h11p,ch10,& 678 & ' setting lmax to 1' 679 call wrtout(ab_out,message,'COLL') 680 call wrtout(std_out, message,'COLL') 681 lmax=1 682 end if 683 end if 684 685 if (abs(h22p)>1.d-08) then 686 nproj(2)=2 687 if (lmax<1) then 688 write(message, '(a,a,a,a,i5,a,e12.4,a,a)' ) ch10,& 689 & ' psp3in : COMMENT -',ch10,& 690 & ' input lmax=',lmax,' does not agree with input h22p=',h22p,ch10,& 691 & ' setting lmax to 1' 692 call wrtout(ab_out,message,'COLL') 693 call wrtout(std_out, message,'COLL') 694 lmax=1 695 end if 696 end if 697 698 if (abs(h33p)>1.d-08) then 699 nproj(2)=3 700 if (lmax<1) then 701 write(message, '(a,a,a,a,i5,a,e12.4,a,a)' ) ch10,& 702 & ' psp3in : COMMENT -',ch10,& 703 & ' input lmax=',lmax,' does not agree with input h33p=',h33p,ch10,& 704 & ' setting lmax to 1' 705 call wrtout(ab_out,message,'COLL') 706 call wrtout(std_out, message,'COLL') 707 lmax=1 708 end if 709 end if 710 711 if (abs(h11d)>1.d-08) then 712 nproj(3)=1 713 if (lmax<2) then 714 write(message, '(a,a,a,a,i5,a,e12.4,a,a)' ) ch10,& 715 & ' psp3in : COMMENT -',ch10,& 716 & ' input lmax=',lmax,' does not agree with input h11d=',h11d,ch10,& 717 & ' setting lmax to 2' 718 call wrtout(ab_out,message,'COLL') 719 call wrtout(std_out, message,'COLL') 720 lmax=2 721 end if 722 end if 723 724 if (abs(h22d)>1.d-08) then 725 nproj(3)=2 726 if (lmax<2) then 727 write(message, '(a,a,a,a,i5,a,e12.4,a,a)' ) ch10,& 728 & ' psp3in : COMMENT -',ch10,& 729 & ' input lmax=',lmax,' does not agree with input h22d=',h22d,ch10,& 730 & ' setting lmax to 2' 731 call wrtout(ab_out,message,'COLL') 732 call wrtout(std_out, message,'COLL') 733 lmax=2 734 end if 735 end if 736 737 if (abs(h33d)>1.d-08) then 738 nproj(3)=3 739 if (lmax<2) then 740 write(message, '(a,a,a,a,i5,a,e12.4,a,a)' ) ch10,& 741 & ' psp3in : COMMENT -',ch10,& 742 & ' input lmax=',lmax,' does not agree with input h33d=',h33d,ch10,& 743 & ' setting lmax to 2' 744 call wrtout(ab_out,message,'COLL') 745 call wrtout(std_out, message,'COLL') 746 lmax=2 747 end if 748 end if 749 750 if (abs(h11f)>1.d-08) then 751 nproj(4)=1 752 if (lmax<3) then 753 write(message, '(a,a,a,a,i5,a,e12.4,a,a)' ) ch10,& 754 & ' psp3in : COMMENT -',ch10,& 755 & ' input lmax=',lmax,' does not agree with input h11f=',h11f,ch10,& 756 & ' setting lmax to 3' 757 call wrtout(ab_out,message,'COLL') 758 call wrtout(std_out, message,'COLL') 759 lmax=3 760 end if 761 end if 762 763 if(pspso/=0) then 764 765 if (abs(k11p)>1.d-08) then 766 nproj(psps%mpsang+1)=1 767 if (lmax<1) then 768 write(message, '(a,a,a,a,i5,a,e12.4,a,a)' ) ch10,& 769 & ' psp3in : COMMENT -',ch10,& 770 & ' input lmax=',lmax,' does not agree with input k11p=',k11p,ch10,& 771 & ' setting lmax to 1' 772 call wrtout(ab_out,message,'COLL') 773 call wrtout(std_out, message,'COLL') 774 lmax=1 775 end if 776 end if 777 778 if (abs(k22p)>1.d-08) then 779 nproj(psps%mpsang+1)=2 780 if (lmax<1) then 781 write(message, '(a,a,a,a,i5,a,e12.4,a,a)' ) ch10,& 782 & ' psp3in : COMMENT -',ch10,& 783 & ' input lmax=',lmax,' does not agree with input k22p=',k22p,ch10,& 784 & ' setting lmax to 1' 785 call wrtout(ab_out,message,'COLL') 786 call wrtout(std_out, message,'COLL') 787 lmax=1 788 end if 789 end if 790 791 792 if (abs(k33p)>1.d-08) then 793 nproj(psps%mpsang+1)=3 794 if (lmax<1) then 795 write(message, '(a,a,a,a,i5,a,e12.4,a,a)' ) ch10,& 796 & ' psp3in : COMMENT -',ch10,& 797 & ' input lmax=',lmax,' does not agree with input k33p=',k33p,ch10,& 798 & ' setting lmax to 1' 799 call wrtout(ab_out,message,'COLL') 800 call wrtout(std_out, message,'COLL') 801 lmax=1 802 end if 803 end if 804 805 if (abs(k11d)>1.d-08) then 806 nproj(psps%mpsang+2)=1 807 if (lmax<2) then 808 write(message, '(a,a,a,a,i5,a,e12.4,a,a)' ) ch10,& 809 & ' psp3in : COMMENT -',ch10,& 810 & ' input lmax=',lmax,' does not agree with input k11d=',k11d,ch10,& 811 & ' setting lmax to 2' 812 call wrtout(ab_out,message,'COLL') 813 call wrtout(std_out, message,'COLL') 814 lmax=2 815 end if 816 end if 817 818 if (abs(k22d)>1.d-08) then 819 nproj(psps%mpsang+2)=2 820 if (lmax<2) then 821 write(message, '(a,a,a,a,i5,a,e12.4,a,a)' ) ch10,& 822 & ' psp3in : COMMENT -',ch10,& 823 & ' input lmax=',lmax,' does not agree with input k22d=',k22d,ch10,& 824 & ' setting lmax to 2' 825 call wrtout(ab_out,message,'COLL') 826 call wrtout(std_out, message,'COLL') 827 lmax=2 828 end if 829 end if 830 831 if (abs(k33d)>1.d-08) then 832 nproj(psps%mpsang+2)=3 833 if (lmax<2) then 834 write(message, '(a,a,a,a,i5,a,e12.4,a,a)' ) ch10,& 835 & ' psp3in : COMMENT -',ch10,& 836 & ' input lmax=',lmax,' does not agree with input k33d=',k33d,ch10,& 837 & ' setting lmax to 2' 838 call wrtout(ab_out,message,'COLL') 839 call wrtout(std_out, message,'COLL') 840 lmax=2 841 end if 842 end if 843 844 if (abs(k11f)>1.d-08) then 845 nproj(psps%mpsang+3)=1 846 if (lmax<3) then 847 write(message, '(a,a,a,a,i5,a,e12.4,a,a)' ) ch10,& 848 & ' psp3in : COMMENT -',ch10,& 849 & ' input lmax=',lmax,' does not agree with input k11f=',k11f,ch10,& 850 & ' setting lmax to 3' 851 call wrtout(ab_out,message,'COLL') 852 call wrtout(std_out, message,'COLL') 853 lmax=3 854 end if 855 end if 856 857 end if 858 859 !Store the coefficients. 860 psps%gth_params%set(ipsp) = .true. 861 psps%gth_params%psppar(0, :, ipsp) = (/ rloc, cc1, cc2, cc3, cc4, zero, zero /) 862 psps%gth_params%psppar(1, :, ipsp) = (/ rrs, h11s, h22s, h33s, zero, zero, zero /) 863 psps%gth_params%psppar(2, :, ipsp) = (/ rrp, h11p, h22p, h33p, zero, zero, zero /) 864 psps%gth_params%psppar(3, :, ipsp) = (/ rrd, h11d, h22d, h33d, zero, zero, zero /) 865 psps%gth_params%psppar(4, :, ipsp) = (/ rrf, h11f, zero, zero, zero, zero, zero /) 866 867 !Store the k coefficients 868 psps%gth_params%psp_k_par(1, :, ipsp) = (/ zero, zero, zero /) 869 psps%gth_params%psp_k_par(2, :, ipsp) = (/ k11p, k22p, k33p /) 870 psps%gth_params%psp_k_par(3, :, ipsp) = (/ k11d, k22d, k33d /) 871 psps%gth_params%psp_k_par(4, :, ipsp) = (/ k11f, zero, zero /) 872 873 !Additionnal wavelet parameters 874 if (dtset%usewvl == 1) then 875 call wvl_descr_psp_fill(psps%gth_params, ipsp, 0, int(psps%zionpsp(ipsp)), int(psps%znuclpsp(ipsp)), tmp_unit) 876 end if 877 878 !Initialize array indlmn array giving l,m,n,ln,lm,s for i=lmn 879 nso=1;if(pspso/=0) nso=2 880 index=0;iln=0;indlmn(:,:)=0 881 do nn=1,nso 882 do ipsang=1+(nn-1)*(lmax+1),nn*lmax+1 883 if (nproj(ipsang)>0) then 884 ll=ipsang-(nn-1)*lmax-1 885 do kk=1,nproj(ipsang) 886 iln=iln+1 887 do mm=1,2*ll*psps%useylm+1 888 index=index+1 889 indlmn(1,index)=ll 890 indlmn(2,index)=mm-ll*psps%useylm-1 891 indlmn(3,index)=kk 892 indlmn(4,index)=ll*ll+(1-psps%useylm)*ll+mm 893 indlmn(5,index)=iln 894 indlmn(6,index)=nn 895 end do 896 end do 897 end if 898 end do 899 end do 900 901 ABI_MALLOC(dvlspl,(psps%mqgrid_ff)) 902 !First, the local potential -- compute on q grid and fit spline 903 call psp2lo(cc1,cc2,cc3,cc4,dvlspl,epsatm,psps%mqgrid_ff,psps%qgrid_ff,vlspl(:,1),rloc,.true.,yp1,ypn,zion) 904 ABI_FREE(dvlspl) 905 906 !DEBUG 907 !write(std_out,*)' psp3in : after psp2lo ' 908 !ENDDEBUG 909 910 !Fit spline to q^2 V(q) (Numerical Recipes subroutine) 911 ABI_MALLOC(work_space,(psps%mqgrid_ff)) 912 ABI_MALLOC(work_spl,(psps%mqgrid_ff)) 913 call spline (psps%qgrid_ff,vlspl(:,1),psps%mqgrid_ff,yp1,ypn,work_spl) 914 vlspl(:,2)=work_spl(:) 915 ABI_FREE(work_space) 916 ABI_FREE(work_spl) 917 918 !Second, compute KB energies and form factors and fit splines 919 ekb(:)=zero 920 921 !Check if any nonlocal projectors are being used 922 mproj=maxval(nproj) 923 if (mproj>0) then 924 925 ABI_MALLOC(ekb_sr,(psps%mpsang,mproj)) 926 ABI_MALLOC(ffspl_sr,(psps%mqgrid_ff,2,psps%mpsang,mproj)) 927 ABI_MALLOC(ekb_so,(psps%mpsang,mproj)) 928 ABI_MALLOC(ffspl_so,(psps%mqgrid_ff,2,psps%mpsang,mproj)) 929 930 call psp3nl(ekb_sr,ffspl_sr,h11s,h22s,h33s,h11p,h22p,h33p,h11d,h22d,& 931 & h33d,h11f,mproj,psps%mpsang,psps%mqgrid_ff,psps%qgrid_ff,rrd,rrf,rrp,rrs) 932 if(pspso/=0) then 933 call psp3nl(ekb_so,ffspl_so,zero,zero,zero,k11p,k22p,k33p,k11d,& 934 & k22d,k33d,k11f,mproj,psps%mpsang,psps%mqgrid_ff,psps%qgrid_ff,rrd,rrf,rrp,rrs) 935 end if 936 937 938 ! Convert ekb and ffspl 939 iln0=0 940 do jj=1,psps%lmnmax 941 iln=indlmn(5,jj) 942 if (iln>iln0) then 943 iln0=iln 944 if (indlmn(6,jj)<=1) then 945 ekb(iln)=ekb_sr(1+indlmn(1,jj),indlmn(3,jj)) 946 ffspl(:,:,iln)=ffspl_sr(:,:,1+indlmn(1,jj),indlmn(3,jj)) 947 else 948 ekb(iln)=ekb_so(1+indlmn(1,jj),indlmn(3,jj)) 949 ffspl(:,:,iln)=ffspl_so(:,:,1+indlmn(1,jj),indlmn(3,jj)) 950 end if 951 end if 952 end do 953 954 ABI_FREE(ekb_sr) 955 ABI_FREE(ffspl_sr) 956 ABI_FREE(ekb_so) 957 ABI_FREE(ffspl_so) 958 end if 959 960 return 961 962 ! Handle IO error 963 10 continue 964 ABI_ERROR(errmsg) 965 966 end subroutine psp3in
m_psp_hgh/psp10in [ Functions ]
[ Top ] [ m_psp_hgh ] [ Functions ]
NAME
psp10in
FUNCTION
Initialize pspcod=10 pseudopotentials (formalism is the same as in HGH psps PRB58,3641(1998) [[cite:Hartwigsen1998]], but the full h and k matrices are read, allowing for using also subsequent developments such as Theor. Chem. Acc. 114, 145 (2005) [[cite:Dolg2005]]: continue to read the file, then compute the corresponding local and non-local potentials.
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
SOURCE
1366 subroutine psp10in(dtset, ekb, epsatm, ffspl, indlmn, ipsp, lmax, nproj, psps, pspso, vlspl, zion) 1367 1368 !Arguments ------------------------------------ 1369 !scalars 1370 integer,intent(in) :: ipsp,pspso 1371 integer,intent(inout) :: lmax 1372 real(dp),intent(in) :: zion 1373 real(dp),intent(out) :: epsatm 1374 type(dataset_type),intent(in) :: dtset 1375 type(pseudopotential_type),intent(inout) :: psps 1376 !arrays 1377 integer,intent(out) :: indlmn(6,psps%lmnmax),nproj(psps%mpssoang) 1378 real(dp),intent(out) :: ekb(psps%lnmax) !vz_i 1379 real(dp),intent(inout) :: ffspl(psps%mqgrid_ff,2,psps%lnmax) !vz_i 1380 real(dp),intent(out) :: vlspl(psps%mqgrid_ff,2) 1381 1382 !Local variables------------------------------- 1383 !scalars 1384 integer :: ii,iln,iln0,index,ipsang,jj,kk,ll,mm,mproj,nn,nnonloc,nprl,nso 1385 real(dp) :: rloc,yp1,ypn 1386 character(len=500) :: message,errmsg 1387 !arrays 1388 integer,allocatable :: dummy_nproj(:) 1389 real(dp) :: cc(4) 1390 real(dp),allocatable :: dvlspl(:),ekb_so(:,:),ekb_sr(:,:),ffspl_so(:,:,:,:) 1391 real(dp),allocatable :: ffspl_sr(:,:,:,:),hij(:,:,:),kij(:,:,:),rr(:) 1392 real(dp),allocatable :: work_space(:),work_spl(:) 1393 1394 ! *************************************************************************** 1395 1396 !Set various terms to 0 in case not defined below 1397 !HGH values 1398 rloc=zero ; cc(:)=zero 1399 nproj(1:psps%mpssoang)=0 1400 1401 !Read and write different lines of the pseudopotential file 1402 1403 read (tmp_unit,*,err=10,iomsg=errmsg) rloc,nn,(cc(jj),jj=1,nn) 1404 write(message, '(a,f12.7)' ) ' rloc=',rloc 1405 call wrtout(ab_out,message,'COLL') 1406 call wrtout(std_out, message,'COLL') 1407 write(message, '(a,i1,a,4f12.7)' )' cc(1:',nn,')=',(cc(jj),jj=1,nn) 1408 call wrtout(ab_out,message,'COLL') 1409 call wrtout(std_out, message,'COLL') 1410 1411 !Read the number of the non-local projectors 1412 read (tmp_unit,*,err=10,iomsg=errmsg) nnonloc 1413 if (nnonloc/=lmax+1) then 1414 write(message, '(a,a,a,a,i5,a,i5,a,a,a,a,i5)' ) ch10,& 1415 & ' psp10in : COMMENT -',ch10,& 1416 & ' input lmax=',lmax,' does not agree with input nnonloc=',nnonloc,ch10,& 1417 & ' which has to be lmax+1.',ch10,& 1418 & ' Setting lmax to ',nnonloc-1 1419 call wrtout(ab_out,message,'COLL') 1420 call wrtout(std_out, message,'COLL') 1421 lmax=1 1422 end if 1423 ABI_MALLOC(rr,(0:lmax)) 1424 ABI_MALLOC(hij,(0:lmax,3,3)) 1425 ABI_MALLOC(kij,(0:lmax,3,3)) 1426 rr(:)=zero; hij(:,:,:)=zero; kij(:,:,:)=zero 1427 1428 !Read and echo the coefficients of non-local projectors 1429 prjloop: do ll=0,lmax 1430 read (tmp_unit,*,err=10,iomsg=errmsg) rr(ll),nprl,(hij(ll,1,jj),jj=1,nprl) 1431 do ii=2,nprl 1432 read (tmp_unit,*,err=10,iomsg=errmsg) (hij(ll,ii,jj),jj=ii,nprl) 1433 end do 1434 nproj(ll+1)=nprl 1435 write(message, '(a,i3,a,f12.7,2a,3f12.7,2a,12x,2f12.7,2a,24x,f12.7)' )& 1436 & ' for angular momentum l =',ll,' r(l) =',rr(ll),ch10,& 1437 & ' h11, h12, h13 =', (hij(ll,1,jj),jj=1,3),ch10,& 1438 & ' h22, h23 =', (hij(ll,2,jj),jj=2,3),ch10,& 1439 & ' h33 =', (hij(ll,3,jj),jj=3,3) 1440 call wrtout(ab_out,message,'COLL') 1441 call wrtout(std_out, message,'COLL') 1442 if (ll==0) cycle 1443 do ii=1,nprl 1444 read (tmp_unit,*,err=10,iomsg=errmsg) (kij(ll,ii,jj),jj=ii,nprl) 1445 end do 1446 write(message, '(a,3f12.7,2a,12x,2f12.7,2a,24x,f12.7)' )& 1447 & ' k11, k12, k13 =', (kij(ll,1,jj),jj=1,3),ch10,& 1448 & ' k22, k23 =', (kij(ll,2,jj),jj=2,3),ch10,& 1449 & ' k33 =', (kij(ll,3,jj),jj=3,3) 1450 call wrtout(ab_out,message,'COLL') 1451 call wrtout(std_out, message,'COLL') 1452 end do prjloop 1453 1454 if(pspso/=0) then 1455 1456 ! MJV 10/2008: is this correct? For the normal HGH psp there are cases 1457 ! where there are more SO projectors than SR ones! e.g. Pb with 12 electrons. 1458 do ll=1,lmax 1459 nproj(psps%mpsang+ll)=nproj(ll+1) 1460 end do 1461 1462 end if 1463 1464 !Store the coefficients. 1465 psps%gth_params%set(ipsp) = .true. 1466 psps%gth_params%psppar(0, :, ipsp) = (/ rloc, cc(1), cc(2), cc(3), cc(4), zero, zero /) 1467 do ii=1,4 1468 ll=ii-1 1469 if (ll>lmax) then 1470 psps%gth_params%psppar(ii,:,ipsp) = (/ zero, zero, zero, zero, zero, zero, zero /) 1471 else 1472 psps%gth_params%psppar(ii,:,ipsp) =& 1473 & (/ rr(ll), hij(ll,1,1), hij(ll,2,2), hij(ll,3,3), hij(ll,1,2), hij(ll,1,3), hij(ll,2,3) /) 1474 end if 1475 end do 1476 1477 !Additionnal wavelet parameters 1478 if (dtset%usewvl == 1) then 1479 call wvl_descr_psp_fill(psps%gth_params, ipsp, 0, int(psps%zionpsp(ipsp)), int(psps%znuclpsp(ipsp)), tmp_unit) 1480 end if 1481 1482 !Initialize array indlmn array giving l,m,n,ln,lm,s for i=lmn 1483 nso=1;if(pspso/=0) nso=2 1484 index=0;iln=0;indlmn(:,:)=0 1485 do nn=1,nso 1486 do ipsang=1+(nn-1)*(lmax+1),nn*lmax+1 1487 if (nproj(ipsang)>0) then 1488 ll=ipsang-(nn-1)*lmax-1 1489 do kk=1,nproj(ipsang) 1490 iln=iln+1 1491 do mm=1,2*ll*psps%useylm+1 1492 index=index+1 1493 indlmn(1,index)=ll 1494 indlmn(2,index)=mm-ll*psps%useylm-1 1495 indlmn(3,index)=kk 1496 indlmn(4,index)=ll*ll+(1-psps%useylm)*ll+mm 1497 indlmn(5,index)=iln 1498 indlmn(6,index)=nn 1499 end do 1500 end do 1501 end if 1502 end do 1503 end do 1504 1505 ABI_MALLOC(dvlspl,(psps%mqgrid_ff)) 1506 !First, the local potential -- compute on q grid and fit spline 1507 call psp2lo(cc(1),cc(2),cc(3),cc(4),dvlspl,epsatm,psps%mqgrid_ff,psps%qgrid_ff,& 1508 & vlspl(:,1),rloc,.true.,yp1,ypn,zion) 1509 ABI_FREE(dvlspl) 1510 1511 !Fit spline to q^2 V(q) (Numerical Recipes subroutine) 1512 ABI_MALLOC(work_space,(psps%mqgrid_ff)) 1513 ABI_MALLOC(work_spl,(psps%mqgrid_ff)) 1514 call spline (psps%qgrid_ff,vlspl(:,1),psps%mqgrid_ff,yp1,ypn,work_spl) 1515 vlspl(:,2)=work_spl(:) 1516 ABI_FREE(work_space) 1517 ABI_FREE(work_spl) 1518 1519 !Second, compute KB energies and form factors and fit splines 1520 ekb(:)=zero 1521 1522 !Check if any nonlocal projectors are being used 1523 mproj=maxval(nproj) 1524 1525 if (mproj>0) then 1526 1527 ABI_MALLOC(ekb_sr,(psps%mpsang,mproj)) 1528 ABI_MALLOC(ffspl_sr,(psps%mqgrid_ff,2,psps%mpsang,mproj)) 1529 ABI_MALLOC(ekb_so,(psps%mpsang,mproj)) 1530 ABI_MALLOC(ffspl_so,(psps%mqgrid_ff,2,psps%mpsang,mproj)) 1531 1532 call psp10nl(ekb_sr,ffspl_sr,hij,lmax,mproj,psps%mpsang,psps%mqgrid_ff,& 1533 & nproj,psps%qgrid_ff,rr) 1534 if(pspso/=0) then 1535 ABI_MALLOC(dummy_nproj,(psps%mpsang)) 1536 dummy_nproj(1)=0 1537 do ll=1,lmax 1538 dummy_nproj(ll+1)=nproj(psps%mpsang+ll) 1539 end do 1540 call psp10nl(ekb_so,ffspl_so,kij,lmax,mproj,psps%mpsang,psps%mqgrid_ff,& 1541 & dummy_nproj,psps%qgrid_ff,rr) 1542 ABI_FREE(dummy_nproj) 1543 end if 1544 1545 ! Convert ekb and ffspl 1546 iln0=0 1547 do jj=1,psps%lmnmax 1548 iln=indlmn(5,jj) 1549 if (iln>iln0) then 1550 iln0=iln 1551 if (indlmn(6,jj)<=1) then 1552 ekb(iln)=ekb_sr(1+indlmn(1,jj),indlmn(3,jj)) 1553 ffspl(:,:,iln)=ffspl_sr(:,:,1+indlmn(1,jj),indlmn(3,jj)) 1554 else 1555 ekb(iln)=ekb_so(1+indlmn(1,jj),indlmn(3,jj)) 1556 ffspl(:,:,iln)=ffspl_so(:,:,1+indlmn(1,jj),indlmn(3,jj)) 1557 end if 1558 end if 1559 end do 1560 1561 ABI_FREE(ekb_sr) 1562 ABI_FREE(ffspl_sr) 1563 ABI_FREE(ekb_so) 1564 ABI_FREE(ffspl_so) 1565 end if 1566 1567 ABI_FREE(rr) 1568 ABI_FREE(hij) 1569 ABI_FREE(kij) 1570 1571 return 1572 1573 ! Handle IO error 1574 10 continue 1575 ABI_ERROR(errmsg) 1576 1577 end subroutine psp10in
m_psp_hgh/psp10nl [ Functions ]
[ Top ] [ m_psp_hgh ] [ Functions ]
NAME
psp10nl
FUNCTION
Hartwigsen-Goedecker-Hutter nonlocal pseudopotential (from preprint of 1998). Uses Gaussians for fully nonlocal form, analytic expressions.
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
SOURCE
1607 subroutine psp10nl(ekb,ffspl,hij,lmax,mproj,mpsang,mqgrid,nproj,qgrid,rr) 1608 1609 !Arguments ------------------------------------ 1610 !scalars 1611 integer,intent(in) :: lmax,mproj,mpsang,mqgrid 1612 !arrays 1613 integer,intent(in) :: nproj(mpsang) 1614 real(dp),intent(in) :: hij(0:lmax,3,3),qgrid(mqgrid),rr(0:lmax) 1615 real(dp),intent(out) :: ekb(mpsang,mproj),ffspl(mqgrid,2,mpsang,mproj) 1616 1617 !Local variables------------------------------- 1618 !scalars 1619 integer :: info,ipack,iproj,iqgrid,jproj,ll,numproj 1620 real(dp) :: qmax,rrl 1621 character(len=500) :: message 1622 character :: jobz,uplo 1623 !arrays 1624 real(dp) :: ap(2,9),rwork1(9),work1(2,9),ww(3),yp1j(3),ypnj(3) 1625 real(dp),allocatable :: ppspl(:,:,:,:),uu(:,:),work(:),zz(:,:,:) 1626 1627 ! ************************************************************************* 1628 1629 ABI_MALLOC(ppspl,(mqgrid,2,mpsang,mproj)) 1630 ABI_MALLOC(work,(mqgrid)) 1631 1632 qmax=qgrid(mqgrid) 1633 jobz='v' 1634 uplo='u' 1635 ekb(:,:)=zero 1636 1637 lloop: do ll=0,lmax 1638 ap(:,:)=zero 1639 numproj=nproj(ll+1) 1640 1641 ! Fill up the matrix in packed storage 1642 prjloop: do jproj=1,numproj 1643 priloop: do iproj=1,jproj 1644 ipack=iproj+(jproj-1)*jproj/2 1645 if(mod((jproj-1)*jproj,2)/=0) then 1646 ABI_ERROR("odd") 1647 end if 1648 ap(1,ipack)=hij(ll,iproj,jproj) 1649 end do priloop 1650 end do prjloop 1651 1652 if(numproj/=0)then 1653 1654 ABI_MALLOC(uu,(numproj,numproj)) 1655 ABI_MALLOC(zz,(2,numproj,numproj)) 1656 1657 if (numproj > 1) then 1658 call ZHPEV(jobz,uplo,numproj,ap,ww,zz,numproj,work1,rwork1,info) 1659 uu(:,:)=zz(1,:,:) 1660 else 1661 ww(1)=hij(ll,1,1) 1662 uu(1,1)=one 1663 end if 1664 1665 ! Initialization of ekb, and spline fitting 1666 1667 if (ll==0) then ! s channel 1668 1669 rrl=rr(0) 1670 do iproj=1,numproj 1671 ekb(1,iproj)=ww(iproj)*32.d0*(rrl**3)*(pi**2.5d0)/(4.d0*pi)**2 1672 if(iproj==1)then 1673 do iqgrid=1,mqgrid 1674 ppspl(iqgrid,1,1,1)=exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrl)**2) 1675 end do 1676 yp1j(1)=zero 1677 ypnj(1)=-(two_pi*rrl)**2*qmax*exp(-0.5d0*(two_pi*qmax*rrl)**2) 1678 else if(iproj==2)then 1679 do iqgrid=1,mqgrid 1680 ppspl(iqgrid,1,1,2)=2.0d0/sqrt(15.0d0) & 1681 & *exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrl)**2) & 1682 & *( 3.d0-(two_pi*qgrid(iqgrid)*rrl)**2 ) 1683 end do 1684 yp1j(2)=zero 1685 ypnj(2)=2.0d0/sqrt(15.0d0)*(two_pi*rrl)**2*qmax & 1686 & *exp(-0.5d0*(two_pi*qmax*rrl)**2) * (-5.d0+(two_pi*qmax*rrl)**2) 1687 else if(iproj==3)then 1688 do iqgrid=1,mqgrid 1689 ppspl(iqgrid,1,1,3)=(4.0d0/3.0d0)/sqrt(105.0d0)*& 1690 & exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrl)**2) * & 1691 & (15.0d0-10.0d0*(two_pi*qgrid(iqgrid)*rrl)**2 + & 1692 & (two_pi*qgrid(iqgrid)*rrl)**4) 1693 end do 1694 yp1j(3)=zero 1695 ypnj(3)=(4.0d0/3.0d0)/sqrt(105.0d0)*exp(-0.5d0*(two_pi*qmax*rrl)**2) * & 1696 & (two_pi*rrl)**2*qmax*(-35.0d0+14d0*(two_pi*qmax*rrl)**2-(two_pi*qmax*rrl)**4) 1697 end if 1698 call spline(qgrid,ppspl(:,1,1,iproj),mqgrid,& 1699 & yp1j(iproj),ypnj(iproj),ppspl(:,2,1,iproj)) 1700 end do 1701 1702 else if (ll==1) then ! p channel 1703 1704 rrl=rr(1) 1705 do iproj=1,numproj 1706 ekb(2,iproj)=ww(iproj)*64.d0*(rrl**5)*(pi**2.5d0)/(4.d0*pi)**2 1707 if(iproj==1)then 1708 do iqgrid=1,mqgrid 1709 ppspl(iqgrid,1,2,1)=(1.0d0/sqrt(3.0d0))* & 1710 & exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrl)**2) * (two_pi*qgrid(iqgrid)) 1711 end do 1712 yp1j(1)=two_pi*(1.0d0/sqrt(3.0d0)) 1713 ypnj(1)=-two_pi*((two_pi*qmax*rrl)**2-1.d0)*exp(-0.5d0*(two_pi*qmax*rrl)**2)*& 1714 & (1.0d0/sqrt(3.0d0)) 1715 else if(iproj==2)then 1716 do iqgrid=1,mqgrid 1717 ppspl(iqgrid,1,2,2)=(2.0d0/sqrt(105.0d0))* & 1718 & exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrl)**2) * & 1719 & (two_pi*qgrid(iqgrid))*(5.0d0-(two_pi*qgrid(iqgrid)*rrl)**2) 1720 end do 1721 yp1j(2)=(5.0d0*two_pi)*(2.0d0/sqrt(105.0d0)) 1722 ypnj(2)=(2.0d0/sqrt(105.0d0))*two_pi*exp(-0.5d0*(two_pi*qmax*rrl)**2)* & 1723 & (-8*(two_pi*qmax*rrl)**2 + (two_pi*qmax*rrl)**4 + 5.0d0) 1724 else if(iproj==3)then 1725 do iqgrid=1,mqgrid 1726 ppspl(iqgrid,1,2,3)=(4.0d0/3.0d0)/sqrt(1155d0)*& 1727 & exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrl)**2) * & 1728 & (two_pi*qgrid(iqgrid))*& 1729 & (35.0d0-14.0d0*(two_pi*qgrid(iqgrid)*rrl)**2+(two_pi*qgrid(iqgrid)*rrl)**4) 1730 end do 1731 yp1j(3)=(35.0d0*two_pi)*(4.0d0/3.0d0)/sqrt(1155.0d0) 1732 ypnj(3)=(4.0d0/3.0d0)/sqrt(1155.0d0)*two_pi*exp(-0.5d0*(two_pi*qmax*rrl)**2)* & 1733 & (35.0d0-77.0d0*(two_pi*qmax*rrl)**2+19.0d0*(two_pi*qmax*rrl)**4 - & 1734 & (two_pi*qmax*rrl)**6) 1735 end if 1736 call spline(qgrid,ppspl(:,1,2,iproj),mqgrid,& 1737 & yp1j(iproj),ypnj(iproj),ppspl(:,2,2,iproj)) 1738 end do 1739 1740 else if (ll==2) then ! d channel 1741 1742 ! If there is a third projector. Warning : only two projectors are allowed. 1743 if ( numproj>2 ) then 1744 write(message, '(3a)' )& 1745 & ' only two d-projectors are allowed ',ch10,& 1746 & ' Action: check your pseudopotential file.' 1747 ABI_ERROR(message) 1748 end if 1749 1750 rrl=rr(2) 1751 do iproj=1,numproj 1752 ekb(3,iproj)=ww(iproj)*128.d0*(rrl**7)*(pi**2.5d0)/(4.d0*pi)**2 1753 if(iproj==1)then 1754 do iqgrid=1,mqgrid 1755 ppspl(iqgrid,1,3,1)=(1.0d0/sqrt(15.0d0))* & 1756 & exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrl)**2) * (two_pi*qgrid(iqgrid))**2 1757 end do 1758 yp1j(1)=zero 1759 ypnj(1)=(1.0d0/sqrt(15.0d0))*(two_pi**2)*& 1760 & exp(-0.5d0*(two_pi*qmax*rrl)**2)*qmax*(2d0-(two_pi*qmax*rrl)**2) 1761 else if(iproj==2)then 1762 do iqgrid=1,mqgrid 1763 ppspl(iqgrid,1,3,2)=(2.0d0/3.0d0)/sqrt(105.0d0)* & 1764 & exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrl)**2) * & 1765 & ((two_pi*qgrid(iqgrid))**2)*(7.0d0-(two_pi*qgrid(iqgrid)*rrl)**2) 1766 end do 1767 yp1j(2)=zero 1768 ypnj(2)=(2.0d0/3.0d0)/sqrt(105.0d0)*exp(-0.5d0*(two_pi*qmax*rrl)**2)* & 1769 & qmax*(two_pi**2)*( (two_pi*qmax*rrl)**4 - 11.0d0*(two_pi*qmax*rrl)**2 + 14.0d0) 1770 end if 1771 call spline(qgrid,ppspl(:,1,3,iproj),mqgrid,& 1772 & yp1j(iproj),ypnj(iproj),ppspl(:,2,3,iproj)) 1773 end do 1774 1775 else if (ll==3) then ! f channel 1776 1777 ! If there is a second projector. Warning : only one projector is allowed. 1778 if ( numproj>1 ) then 1779 write(message, '(a,a,a)' )& 1780 & 'only one f-projector is allowed ',ch10,& 1781 & 'Action: check your pseudopotential file.' 1782 ABI_ERROR(message) 1783 end if 1784 1785 rrl=rr(3) 1786 ekb(4,1)=ww(1)*(256.0d0/105.0d0)*(rrl**9)*(pi**2.5d0)/(4.d0*pi)**2 1787 do iqgrid=1,mqgrid 1788 ppspl(iqgrid,1,4,1)=(two_pi*qgrid(iqgrid))**3* & 1789 & exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrl)**2) 1790 end do 1791 ! Compute yp1,ypn=derivatives of f(q) at q=0, q=qgrid(mqgrid) 1792 yp1j(1)=zero 1793 ypnj(1)=(two_pi**3)*qmax**2*exp(-0.5d0*(two_pi*qmax*rrl)**2)*& 1794 & (3.0d0-(two_pi*qmax*rrl)**2) 1795 ! Fit spline to get second derivatives by spline fit 1796 call spline(qgrid,ppspl(:,1,4,1),mqgrid,& 1797 & yp1j(1),ypnj(1),ppspl(:,2,4,1)) 1798 1799 else 1800 ABI_ERROR("lmax>3?") 1801 end if 1802 1803 ! Linear combination using the eigenvectors 1804 ffspl(:,:,ll+1,:)=zero 1805 do jproj=1,numproj 1806 do iproj=1,numproj 1807 do iqgrid=1,mqgrid 1808 ffspl(iqgrid,1:2,ll+1,jproj)=ffspl(iqgrid,1:2,ll+1,jproj) & 1809 & +uu(iproj,jproj)*ppspl(iqgrid,1:2,ll+1,iproj) 1810 end do 1811 end do 1812 end do 1813 1814 ABI_FREE(uu) 1815 ABI_FREE(zz) 1816 1817 ! End condition on numproj(/=0) 1818 end if 1819 1820 end do lloop 1821 1822 ABI_FREE(ppspl) 1823 ABI_FREE(work) 1824 1825 end subroutine psp10nl
m_psp_hgh/psp2in [ Functions ]
[ Top ] [ m_psp_hgh ] [ Functions ]
NAME
psp2in
FUNCTION
Initialize pspcod=2 pseudopotentials (GTH format): continue to read the file, then compute the corresponding local and non-local potentials.
INPUTS
dtset <type(dataset_type)>=all input variables in this dataset ipsp=id in the array of the pseudo-potential. lmax=value of lmax mentioned at the second line of the psp file 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) epsatm=$ (4\pi)\int_0^\infty [r^2 (V(r)+Zv/r) dr]$ (hartree) ffspl(mqgrid,2,lnmax)=Kleinman-Bylander form factor f_l(q) and second derivative from spline fit for each angular momentum and each projector indlmn(6,i)= array giving l,m,n,lm,ln,s for i=ln (if useylm=0) or i=lmn (if useylm=1) nproj(mpsang)=number of projection functions for each angular momentum vlspl(mqgrid_vl,2)=q^2 Vloc(q) and second derivatives from spline fit dvlspl(mqgrid_vl,2)=dVloc(r)/dr and second derivatives from spline fit (only allocated if vlspl_recipSpace is false.
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 | mqgrid_ff(IN)=dimension of q (or G) grid for nl form factors (array ffspl) | mqgrid_vl(IN)=dimension of q (or G) grid or r grid (if vlspl_recipSpace = .false.) | qgrid_ff(mqgrid_ff)(IN)=values of q on grid from 0 to qmax (bohr^-1) for nl form factors | qgrid_vl(mqgrid_vl)(IN)=values of q on grid from 0 to qmax (bohr^-1) for Vloc | if vlspl_recipSpace is .true. else values of r on grid from | 0 to 2pi / qmax * mqgrid_ff (bohr). | useylm(IN)=governs the way the nonlocal operator is to be applied: | 1=using Ylm, 0=using Legendre polynomials | vlspl_recipSpace(IN)=.true. if pseudo are expressed in reciprocal space. | gth_params(OUT)=store GTH coefficients and parameters.
SOURCE
107 subroutine psp2in(dtset,ekb,epsatm,ffspl,indlmn,ipsp,lmax,nproj,psps,vlspl,dvlspl,zion) 108 109 !Arguments ------------------------------------ 110 !scalars 111 integer,intent(in) :: ipsp,lmax 112 real(dp),intent(in) :: zion 113 real(dp),intent(out) :: epsatm 114 type(dataset_type),intent(in) :: dtset 115 type(pseudopotential_type),intent(inout) :: psps 116 !arrays 117 integer,intent(out) :: indlmn(6,psps%lmnmax),nproj(psps%mpsang) 118 real(dp),intent(out) :: dvlspl(psps%mqgrid_vl,2),ekb(psps%lnmax) 119 real(dp),intent(inout) :: ffspl(psps%mqgrid_ff,2,psps%lnmax) !vz_i 120 real(dp),intent(out) :: vlspl(psps%mqgrid_vl,2) 121 122 !Local variables------------------------------- 123 !scalars 124 integer :: iln,index,ipsang,kk,ll,mm 125 real(dp) :: cc1,cc2,cc3,cc4,h1p,h1s,h2s,rloc,rrp,rrs 126 real(dp) :: yp1,ypn 127 character(len=500) :: message,errmsg 128 !arrays 129 real(dp),allocatable :: work_space(:),work_spl(:) 130 real(dp),allocatable :: dvloc(:) 131 132 ! *************************************************************************** 133 134 !Set various terms to 0 in case not defined below 135 !GTH values 136 rloc=0.d0 137 cc1=0.d0 138 cc2=0.d0 139 cc3=0.d0 140 cc4=0.d0 141 rrs=0.d0 142 h1s=0.d0 143 h2s=0.d0 144 rrp=0.d0 145 h1p=0.d0 146 nproj(1:psps%mpsang)=0 147 148 !Read and write different lines of the pseudopotential file 149 read (tmp_unit,*, err=10, iomsg=errmsg) rloc,cc1,cc2,cc3,cc4 150 write(message, '(a,f12.7)' ) ' rloc=',rloc 151 call wrtout(ab_out,message,'COLL') 152 call wrtout(std_out, message,'COLL') 153 write(message, '(a,f12.7,a,f12.7,a,f12.7,a,f12.7)' )' cc1=',cc1,'; cc2=',cc2,'; cc3=',cc3,'; cc4=',cc4 154 call wrtout(ab_out,message,'COLL') 155 call wrtout(std_out, message,'COLL') 156 157 read (tmp_unit,*, err=10, iomsg=errmsg) rrs,h1s,h2s 158 write(message, '(a,f12.7,a,f12.7,a,f12.7)' )' rrs=',rrs,'; h1s=',h1s,'; h2s=',h2s 159 call wrtout(ab_out,message,'COLL') 160 call wrtout(std_out, message,'COLL') 161 162 read (tmp_unit,*, err=10, iomsg=errmsg) rrp,h1p 163 write(message, '(a,f12.7,a,f12.7)' )' rrp=',rrp,'; h1p=',h1p 164 call wrtout(ab_out,message,'COLL') 165 call wrtout(std_out, message,'COLL') 166 167 !Store the coefficients. 168 psps%gth_params%set(ipsp) = .true. 169 psps%gth_params%psppar(0, :, ipsp) = (/ rloc, cc1, cc2, cc3, cc4, 0.d0, 0.d0 /) 170 psps%gth_params%psppar(1, :, ipsp) = (/ rrs, h1s, h2s, 0.d0, 0.d0, 0.d0, 0.d0 /) 171 psps%gth_params%psppar(2, :, ipsp) = (/ rrp, h1p, 0.d0, 0.d0, 0.d0, 0.d0, 0.d0 /) 172 if (dtset%usewvl == 1) then 173 call wvl_descr_psp_fill(psps%gth_params, ipsp, 0, int(psps%zionpsp(ipsp)), int(psps%znuclpsp(ipsp)), tmp_unit) 174 end if 175 176 if (abs(h1s)>1.d-08) nproj(1)=1 177 if (abs(h2s)>1.d-08) nproj(1)=2 178 179 if (abs(h1p)>1.d-08) then 180 if(psps%mpsang<2)then 181 write(message, '(a,es12.4,a,a,a,i2,a)' )& 182 & 'With non-zero h1p (=',h1p,'), mpsang should be at least 2,',ch10,& 183 & 'while mpsang=',psps%mpsang,'.' 184 ABI_ERROR(message) 185 end if 186 nproj(2)=1 187 if (lmax<1) then 188 write(message, '(a,i5,a,e12.4,a,a,a,a)' )& 189 & 'Input lmax=',lmax,' disagree with input h1p=',h1p,'.',& 190 & 'Your pseudopotential is incoherent.',ch10,& 191 & 'Action: correct your pseudopotential file.' 192 ABI_ERROR(message) 193 end if 194 end if 195 196 !Initialize array indlmn array giving l,m,n,lm,ln,s for i=lmn 197 index=0;iln=0;indlmn(:,:)=0 198 do ipsang=1,lmax+1 199 if(nproj(ipsang)>0)then 200 ll=ipsang-1 201 do kk=1,nproj(ipsang) 202 iln=iln+1 203 do mm=1,2*ll*psps%useylm+1 204 index=index+1 205 indlmn(1,index)=ll 206 indlmn(2,index)=mm-ll*psps%useylm-1 207 indlmn(3,index)=kk 208 indlmn(4,index)=ll*ll+(1-psps%useylm)*ll+mm 209 indlmn(5,index)=iln 210 indlmn(6,index)=1 211 end do 212 end do 213 end if 214 end do 215 216 !First, the local potential -- 217 !compute q^2V(q) or V(r) 218 !MJV NOTE: psp2lo should never be called with dvspl unallocated, which 219 !is possible unless .not.psps%vlspl_recipSpace 220 ABI_MALLOC(dvloc,(psps%mqgrid_vl)) 221 call psp2lo(cc1,cc2,cc3,cc4,dvloc,epsatm,psps%mqgrid_vl,psps%qgrid_vl,& 222 & vlspl(:,1),rloc,psps%vlspl_recipSpace,yp1,ypn,zion) 223 224 !Fit spline to (q^2)V(q) or V(r) 225 ABI_MALLOC(work_space,(psps%mqgrid_vl)) 226 ABI_MALLOC(work_spl,(psps%mqgrid_vl)) 227 call spline (psps%qgrid_vl,vlspl(:,1),psps%mqgrid_vl,yp1,ypn,work_spl) 228 vlspl(:,2)=work_spl(:) 229 if (.not.psps%vlspl_recipSpace) then 230 dvlspl(:,1) = dvloc 231 call spline (psps%qgrid_vl,dvlspl(:,1),psps%mqgrid_vl,yp1,ypn,work_spl) 232 dvlspl(:,2)=work_spl(:) 233 end if 234 235 ABI_FREE(work_space) 236 ABI_FREE(work_spl) 237 ABI_FREE(dvloc) 238 239 240 !Second, compute KB energies and form factors and fit splines 241 ekb(:)=0.0d0 242 !First check if any nonlocal projectors are being used 243 if (maxval(nproj(1:lmax+1))>0) then 244 call psp2nl(ekb,ffspl,h1p,h1s,h2s,psps%lnmax,psps%mqgrid_ff,psps%qgrid_ff,rrp,rrs) 245 end if 246 247 return 248 249 ! Handle IO error 250 10 continue 251 ABI_ERROR(errmsg) 252 253 end subroutine psp2in
m_psp_hgh/psp2nl [ Functions ]
[ Top ] [ m_psp_hgh ] [ Functions ]
NAME
psp2nl
FUNCTION
Goedecker-Teter-Hutter nonlocal pseudopotential (from preprint of 1996). Uses Gaussians for fully nonlocal form, analytic expressions.
INPUTS
h1p=factor defining strength of 1st projector for l=1 channel h1s=factor defining strength of 1st projector for l=0 channel h2s=factor defining strength of 2nd projector for l=0 channel lnmax=max. number of (l,n) components over all type of psps mqgrid=number of grid points for qgrid qgrid(mqgrid)=array of |G| values rrp=core radius for p channel (bohr) rrs=core radius for s channel (bohr)
OUTPUT
ekb(lnmax)=Kleinman-Bylander energy ffspl(mqgrid,2,lnmax)=Kleinman-Bylander form factor f_l(q) and second derivative from spline fit for each angular momentum and each projector
SOURCE
282 subroutine psp2nl(ekb,ffspl,h1p,h1s,h2s,lnmax,mqgrid,qgrid,rrp,rrs) 283 284 !Arguments ------------------------------------ 285 !scalars 286 integer,intent(in) :: lnmax,mqgrid 287 real(dp),intent(in) :: h1p,h1s,h2s,rrp,rrs 288 !arrays 289 real(dp),intent(in) :: qgrid(mqgrid) 290 real(dp),intent(inout) :: ekb(lnmax),ffspl(mqgrid,2,lnmax) !vz_i 291 292 !Local variables------------------------------- 293 !scalars 294 integer :: iln,iqgrid 295 real(dp) :: qmax,yp1,ypn 296 !arrays 297 real(dp),allocatable :: work(:) 298 299 ! ************************************************************************* 300 301 ABI_MALLOC(work,(mqgrid)) 302 303 !Kleinman-Bylander energies ekb were set to zero in calling program 304 305 !Compute KB energies 306 iln=0 307 if (abs(h1s)>1.d-12) then 308 iln=iln+1 309 ekb(iln)=h1s*32.d0*rrs**3*(pi**(2.5d0)/(4.d0*pi)**2) 310 end if 311 if (abs(h2s)>1.d-12) then 312 iln=iln+1 313 ekb(iln) =h2s*(128.d0/15.d0)*rrs**3*(pi**(2.5d0)/(4.d0*pi)**2) 314 end if 315 if (abs(h1p)>1.d-12) then 316 iln=iln+1 317 ekb(iln)=h1p*(64.d0/3.d0)*rrp**5*(pi**(2.5d0)/(4.d0*pi)**2) 318 end if 319 320 !Compute KB form factor 321 iln=0 322 323 !l=0 first projector 324 if (abs(h1s)>1.d-12) then 325 iln=iln+1 326 do iqgrid=1,mqgrid 327 ffspl(iqgrid,1,iln)=exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrs)**2) 328 end do 329 ! Compute yp1,ypn=derivatives of f(q) at q=0, q=qgrid(mqgrid) 330 yp1=0.d0 331 qmax=qgrid(mqgrid) 332 ypn=-4.d0*pi**2*qmax*rrs**2*exp(-0.5d0*(two_pi*qmax*rrs)**2) 333 ! Fit spline to get second derivatives by spline fit 334 call spline(qgrid,ffspl(:,1,iln),mqgrid,yp1,ypn,ffspl(:,2,iln)) 335 ! else 336 ! or else put first projector nonlocal correction at l=0 to 0 337 ! ffspl(:,:,iln)=0.0d0 338 end if 339 340 !l=0 second projector 341 if (abs(h2s)>1.d-12) then 342 iln=iln+1 343 do iqgrid=1,mqgrid 344 ffspl(iqgrid,1,iln)=exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrs)**2) * & 345 & (3.d0-(two_pi*qgrid(iqgrid)*rrs)**2) 346 end do 347 ! Compute yp1,ypn=derivatives of f(q) at q=0, q=qgrid(mqgrid) 348 yp1=0.d0 349 qmax=qgrid(mqgrid) 350 ypn=4.d0*pi**2*qmax*rrs**2*exp(-0.5d0*(two_pi*qmax*rrs)**2) * & 351 & (-5.d0+(two_pi*qmax*rrs)**2) 352 ! Fit spline to get second derivatives by spline fit 353 call spline(qgrid,ffspl(:,1,iln),mqgrid,yp1,ypn,ffspl(:,2,iln)) 354 ! else if(mproj>=2)then 355 ! or else put second projector nonlocal correction at l=0 to 0 356 ! ffspl(:,:,iln)=0.0d0 357 end if 358 359 !l=1 first projector 360 if (abs(h1p)>1.d-12) then 361 iln=iln+1 362 do iqgrid=1,mqgrid 363 ffspl(iqgrid,1,iln)=exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrp)**2) * & 364 & (two_pi*qgrid(iqgrid)) 365 end do 366 ! Compute yp1,ypn=derivatives of f(q) at q=0, q=qgrid(mqgrid) 367 yp1=two_pi 368 qmax=qgrid(mqgrid) 369 ypn=-two_pi*((two_pi*qmax*rrp)**2-1.d0) * exp(-0.5d0*(two_pi*qmax*rrp)**2) 370 ! Fit spline to get second derivatives by spline fit 371 call spline(qgrid,ffspl(:,1,iln),mqgrid,yp1,ypn,ffspl(:,2,iln)) 372 ! else if(mpsang>=2)then 373 ! or else put first projector l=1 nonlocal correction to 0 374 ! ffspl(:,:,iln)=0.0d0 375 end if 376 377 ABI_FREE(work) 378 379 end subroutine psp2nl
m_psp_hgh/psp3nl [ Functions ]
[ Top ] [ m_psp_hgh ] [ 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
SOURCE
1005 subroutine psp3nl(ekb,ffspl,h11s,h22s,h33s,h11p,h22p,h33p,h11d,h22d,& 1006 & h33d,h11f,mproj,mpsang,mqgrid,qgrid,rrd,rrf,rrp,rrs) 1007 1008 !Arguments ------------------------------------ 1009 !scalars 1010 integer,intent(in) :: mproj,mpsang,mqgrid 1011 real(dp),intent(in) :: h11d,h11f,h11p,h11s,h22d,h22p,h22s,h33d,h33p,h33s,rrd 1012 real(dp),intent(in) :: rrf,rrp,rrs 1013 !arrays 1014 real(dp),intent(in) :: qgrid(mqgrid) 1015 real(dp),intent(out) :: ekb(mpsang,mproj),ffspl(mqgrid,2,mpsang,mproj) 1016 1017 !Local variables------------------------------- 1018 !scalars 1019 integer :: info,iproj,iqgrid,ldz,mu,nproj,nu 1020 real(dp) :: qmax 1021 character(len=500) :: message 1022 character :: jobz,uplo 1023 !arrays 1024 real(dp) :: ap(2,9),rwork1(9),work1(2,9),ww(3),yp1j(3),ypnj(3) 1025 real(dp),allocatable :: ppspl(:,:,:,:),uu(:,:),work(:),zz(:,:,:) 1026 1027 ! ************************************************************************* 1028 1029 ABI_MALLOC(ppspl,(mqgrid,2,mpsang,mproj)) 1030 ABI_MALLOC(work,(mqgrid)) 1031 1032 qmax=qgrid(mqgrid) 1033 jobz='v' 1034 uplo='u' 1035 ekb(:,:)=0.0d0 1036 1037 !--------------------------------------------------------------- 1038 !Treat s channel 1039 1040 nproj=0 1041 ap(:,:)=0.0d0 1042 !If there is at least one s-projector 1043 if ( abs(h11s) >= 1.0d-8 ) then 1044 nproj=1 ; ldz=1 ; ap(1,1)=h11s 1045 end if 1046 nproj=1 1047 !If there is a second projector 1048 if ( abs(h22s) >= 1.0d-8 ) then 1049 nproj=2 ; ldz=2 ; ap(1,3)=h22s 1050 ap(1,2)=-0.5d0*sqrt(0.6d0)*h22s 1051 end if 1052 !If there is a third projector 1053 if ( abs(h33s) >= 1.0d-8 ) then 1054 nproj=3 ; ldz=3 ; ap(1,6)=h33s 1055 ap(1,4)=0.5d0*sqrt(5.d0/21.d0)*h33s 1056 ap(1,5)=-0.5d0*sqrt(100.d0/63.d0)*h33s 1057 end if 1058 1059 if(nproj/=0)then 1060 1061 ABI_MALLOC(uu,(nproj,nproj)) 1062 ABI_MALLOC(zz,(2,nproj,nproj)) 1063 1064 if (nproj > 1) then 1065 call ZHPEV(jobz,uplo,nproj,ap,ww,zz,ldz,work1,rwork1,info) 1066 uu(:,:)=zz(1,:,:) 1067 else 1068 ww(1)=h11s 1069 uu(1,1)=1.0d0 1070 end if 1071 1072 ! Initialization of ekb, and spline fitting 1073 do iproj=1,nproj 1074 ekb(1,iproj)=ww(iproj)*32.d0*(rrs**3)*(pi**2.5d0)/(4.d0*pi)**2 1075 if(iproj==1)then 1076 do iqgrid=1,mqgrid 1077 ppspl(iqgrid,1,1,1)=exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrs)**2) 1078 end do 1079 yp1j(1)=0.d0 1080 ypnj(1)=-(two_pi*rrs)**2*qmax*exp(-0.5d0*(two_pi*qmax*rrs)**2) 1081 else if(iproj==2)then 1082 do iqgrid=1,mqgrid 1083 ppspl(iqgrid,1,1,2)=2.0d0/sqrt(15.0d0) & 1084 & *exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrs)**2) & 1085 & *( 3.d0-(two_pi*qgrid(iqgrid)*rrs)**2 ) 1086 end do 1087 yp1j(2)=0.0d0 1088 ypnj(2)=2.0d0/sqrt(15.0d0)*(two_pi*rrs)**2*qmax & 1089 & *exp(-0.5d0*(two_pi*qmax*rrs)**2) * (-5.d0+(two_pi*qmax*rrs)**2) 1090 else if(iproj==3)then 1091 do iqgrid=1,mqgrid 1092 ppspl(iqgrid,1,1,3)=(4.0d0/3.0d0)/sqrt(105.0d0)*& 1093 & exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrs)**2) * & 1094 & (15.0d0-10.0d0*(two_pi*qgrid(iqgrid)*rrs)**2 + & 1095 & (two_pi*qgrid(iqgrid)*rrs)**4) 1096 end do 1097 yp1j(3)=0.0d0 1098 ypnj(3)=(4.0d0/3.0d0)/sqrt(105.0d0)*exp(-0.5d0*(two_pi*qmax*rrs)**2) * & 1099 & (two_pi*rrs)**2*qmax*(-35.0d0+14d0*(two_pi*qmax*rrs)**2-(two_pi*qmax*rrs)**4) 1100 end if 1101 call spline(qgrid,ppspl(:,1,1,iproj),mqgrid,& 1102 & yp1j(iproj),ypnj(iproj),ppspl(:,2,1,iproj)) 1103 end do 1104 1105 ! Linear combination using the eigenvectors 1106 ffspl(:,:,1,:)=0.0d0 1107 do mu=1,nproj 1108 do nu=1,nproj 1109 do iqgrid=1,mqgrid 1110 ffspl(iqgrid,1:2,1,mu)=ffspl(iqgrid,1:2,1,mu) & 1111 & +uu(nu,mu)*ppspl(iqgrid,1:2,1,nu) 1112 end do 1113 end do 1114 end do 1115 1116 ABI_FREE(uu) 1117 ABI_FREE(zz) 1118 end if ! End condition on nproj(/=0) 1119 1120 !-------------------------------------------------------------------- 1121 !Now treat p channel 1122 1123 nproj=0 1124 ap(:,:)=0.0d0 1125 !If there is at least one projector 1126 if ( abs(h11p) >= 1.0d-8 ) then 1127 nproj=1 ; ldz=1 ; ap(1,1)=h11p 1128 end if 1129 !If there is a second projector 1130 if ( abs(h22p) >= 1.0d-8 ) then 1131 nproj=2 ; ldz=2 ; ap(1,3)=h22p 1132 ap(1,2)=-0.5d0*sqrt(5.d0/7.d0)*h22p 1133 end if 1134 !If there is a third projector 1135 if ( abs(h33p) >= 1.0d-8 ) then 1136 nproj=3 ; ldz=3 ; ap(1,6)=h33p 1137 ap(1,4)= (1.d0/6.d0)*sqrt(35.d0/11.d0)*h33p 1138 ap(1,5)=-(1.d0/6.d0)*(14.d0/sqrt(11.d0))*h33p 1139 end if 1140 1141 if(nproj/=0)then 1142 1143 ABI_MALLOC(uu,(nproj,nproj)) 1144 ABI_MALLOC(zz,(2,nproj,nproj)) 1145 1146 if (nproj > 1) then 1147 call ZHPEV(jobz,uplo,nproj,ap,ww,zz,ldz,work1,rwork1,info) 1148 uu(:,:)=zz(1,:,:) 1149 else 1150 ww(1)=h11p 1151 uu(1,1)=1.0d0 1152 end if 1153 1154 ! Initialization of ekb, and spline fitting 1155 do iproj=1,nproj 1156 ekb(2,iproj)=ww(iproj)*64.d0*(rrp**5)*(pi**2.5d0)/(4.d0*pi)**2 1157 if(iproj==1)then 1158 do iqgrid=1,mqgrid 1159 ppspl(iqgrid,1,2,1)=(1.0d0/sqrt(3.0d0))* & 1160 & exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrp)**2) * (two_pi*qgrid(iqgrid)) 1161 end do 1162 yp1j(1)=two_pi*(1.0d0/sqrt(3.0d0)) 1163 ypnj(1)=-two_pi*((two_pi*qmax*rrp)**2-1.d0)*exp(-0.5d0*(two_pi*qmax*rrp)**2)*& 1164 & (1.0d0/sqrt(3.0d0)) 1165 else if(iproj==2)then 1166 do iqgrid=1,mqgrid 1167 ppspl(iqgrid,1,2,2)=(2.0d0/sqrt(105.0d0))* & 1168 & exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrp)**2) * & 1169 & (two_pi*qgrid(iqgrid))*(5.0d0-(two_pi*qgrid(iqgrid)*rrp)**2) 1170 end do 1171 yp1j(2)=(5.0d0*two_pi)*(2.0d0/sqrt(105.0d0)) 1172 ypnj(2)=(2.0d0/sqrt(105.0d0))*two_pi*exp(-0.5d0*(two_pi*qmax*rrp)**2)* & 1173 & (-8*(two_pi*qmax*rrp)**2 + (two_pi*qmax*rrp)**4 + 5.0d0) 1174 else if(iproj==3)then 1175 do iqgrid=1,mqgrid 1176 ppspl(iqgrid,1,2,3)=(4.0d0/3.0d0)/sqrt(1155d0)*& 1177 & exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrp)**2) * & 1178 & (two_pi*qgrid(iqgrid))*& 1179 & (35.0d0-14.0d0*(two_pi*qgrid(iqgrid)*rrp)**2+(two_pi*qgrid(iqgrid)*rrp)**4) 1180 end do 1181 yp1j(3)=(35.0d0*two_pi)*(4.0d0/3.0d0)/sqrt(1155.0d0) 1182 ypnj(3)=(4.0d0/3.0d0)/sqrt(1155.0d0)*two_pi*exp(-0.5d0*(two_pi*qmax*rrp)**2)* & 1183 & (35.0d0-77.0d0*(two_pi*qmax*rrp)**2+19.0d0*(two_pi*qmax*rrp)**4 - & 1184 & (two_pi*qmax*rrp)**6) 1185 end if 1186 call spline(qgrid,ppspl(:,1,2,iproj),mqgrid,& 1187 & yp1j(iproj),ypnj(iproj),ppspl(:,2,2,iproj)) 1188 end do 1189 1190 ! Linear combination using the eigenvectors 1191 ffspl(:,:,2,:)=0.0d0 1192 do mu=1,nproj 1193 do nu=1,nproj 1194 do iqgrid=1,mqgrid 1195 ffspl(iqgrid,1:2,2,mu)=ffspl(iqgrid,1:2,2,mu) & 1196 & +uu(nu,mu)*ppspl(iqgrid,1:2,2,nu) 1197 end do 1198 end do 1199 end do 1200 1201 ABI_FREE(uu) 1202 ABI_FREE(zz) 1203 end if ! End condition on nproj(/=0) 1204 1205 !----------------------------------------------------------------------- 1206 !Now treat d channel. 1207 1208 nproj=0 1209 ap(:,:)=0.0d0 1210 !If there is at least one projector 1211 if ( abs(h11d) >= 1.0d-8 ) then 1212 nproj=1 ; ldz=1 ; ap(1,1)=h11d 1213 end if 1214 !If there is a second projector 1215 if ( abs(h22d) >= 1.0d-8 ) then 1216 nproj=2 ; ldz=2 ; ap(1,3)=h22d 1217 ap(1,2)=-0.5d0*sqrt(7.d0/9.d0)*h22d 1218 end if 1219 !If there is a third projector. Warning : only two projectors are allowed. 1220 if ( abs(h33d) >= 1.0d-8 ) then 1221 write(message, '(a,a,a)' )& 1222 & ' only two d-projectors are allowed ',ch10,& 1223 & ' Action: check your pseudopotential file.' 1224 ABI_ERROR(message) 1225 ! nproj=3 ; ldz=3 ; ap(1,6)=h33d 1226 ! ap(1,4)= 0.5d0*sqrt(63.d0/143.d0)*h33d 1227 ! ap(1,5)= -0.5d0*(18.d0/sqrt(143.d0))*h33d 1228 end if 1229 1230 if(nproj/=0)then 1231 1232 ABI_MALLOC(uu,(nproj,nproj)) 1233 ABI_MALLOC(zz,(2,nproj,nproj)) 1234 1235 if (nproj > 1) then 1236 call ZHPEV(jobz,uplo,nproj,ap,ww,zz,ldz,work1,rwork1,info) 1237 uu(:,:)=zz(1,:,:) 1238 else 1239 ww(1)=h11d 1240 uu(1,1)=1.0d0 1241 end if 1242 1243 ! Initialization of ekb, and spline fitting 1244 do iproj=1,nproj 1245 ekb(3,iproj)=ww(iproj)*128.d0*(rrd**7)*(pi**2.5d0)/(4.d0*pi)**2 1246 if(iproj==1)then 1247 do iqgrid=1,mqgrid 1248 ppspl(iqgrid,1,3,1)=(1.0d0/sqrt(15.0d0))* & 1249 & exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrd)**2) * (two_pi*qgrid(iqgrid))**2 1250 end do 1251 yp1j(1)=0.0d0 1252 ypnj(1)=(1.0d0/sqrt(15.0d0))*(two_pi**2)*& 1253 & exp(-0.5d0*(two_pi*qmax*rrd)**2)*qmax*(2d0-(two_pi*qmax*rrd)**2) 1254 else if(iproj==2)then 1255 do iqgrid=1,mqgrid 1256 ppspl(iqgrid,1,3,2)=(2.0d0/3.0d0)/sqrt(105.0d0)* & 1257 & exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrd)**2) * & 1258 & ((two_pi*qgrid(iqgrid))**2)*(7.0d0-(two_pi*qgrid(iqgrid)*rrd)**2) 1259 end do 1260 yp1j(2)=0.0d0 1261 ypnj(2)=(2.0d0/3.0d0)/sqrt(105.0d0)*exp(-0.5d0*(two_pi*qmax*rrd)**2)* & 1262 & qmax*(two_pi**2)*( (two_pi*qmax*rrd)**4 - 11.0d0*(two_pi*qmax*rrd)**2 + 14.0d0) 1263 end if 1264 call spline(qgrid,ppspl(:,1,3,iproj),mqgrid,& 1265 & yp1j(iproj),ypnj(iproj),ppspl(:,2,3,iproj)) 1266 end do 1267 1268 ! Linear combination using the eigenvectors 1269 ffspl(:,:,3,:)=0.0d0 1270 do mu=1,nproj 1271 do nu=1,nproj 1272 do iqgrid=1,mqgrid 1273 ffspl(iqgrid,1:2,3,mu)=ffspl(iqgrid,1:2,3,mu) & 1274 & +uu(nu,mu)*ppspl(iqgrid,1:2,3,nu) 1275 end do 1276 end do 1277 end do 1278 1279 ABI_FREE(uu) 1280 ABI_FREE(zz) 1281 end if ! End condition on nproj(/=0) 1282 1283 !----------------------------------------------------------------------- 1284 !Treat now f channel (max one projector ! - so do not use ppspl) 1285 1286 !l=3 first projector 1287 if (abs(h11f)>1.d-12) then 1288 ekb(4,1)=h11f*(256.0d0/105.0d0)*(rrf**9)*(pi**2.5d0)/(4.d0*pi)**2 1289 do iqgrid=1,mqgrid 1290 ffspl(iqgrid,1,4,1)=(two_pi*qgrid(iqgrid))**3* & 1291 & exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrf)**2) 1292 end do 1293 ! Compute yp1,ypn=derivatives of f(q) at q=0, q=qgrid(mqgrid) 1294 yp1j(1)=0d0 1295 ypnj(1)=(two_pi**3)*qmax**2*exp(-0.5d0*(two_pi*qmax*rrf)**2)*& 1296 & (3.0d0-(two_pi*qmax*rrf)**2) 1297 ! Fit spline to get second derivatives by spline fit 1298 call spline(qgrid,ffspl(:,1,4,1),mqgrid,& 1299 & yp1j(1),ypnj(1),ffspl(:,2,4,1)) 1300 end if 1301 1302 !----------------------------------------------------------------------- 1303 1304 ABI_FREE(ppspl) 1305 ABI_FREE(work) 1306 1307 end subroutine psp3nl