TABLE OF CONTENTS
ABINIT/psden [ Functions ]
NAME
psden
FUNCTION
Calculate a pseudo-density from an original density on a radial grid (regular or logarithmic)
INPUTS
ilog=1 if grid is logarithmic, else 0 mesh= dimension of nc nc(mesh)= density to be pseudized rc= cut-off radius rad(mesh) = radial mesh
OUTPUT
ff(mesh)= pseudized density
Optional:
ff1(mesh)= 1st derivative of pseudo density (only r<rc modified)
ff2(mesh)= 2nd derivative of pseudo density (only r<rc modified)
NOTES
ff=exp(-(a+b.r^2+c.r^4))
PARENTS
psp6cc
CHILDREN
ctrap
SOURCE
588 subroutine psden(ilog,ff,mesh,nc,rc,rad,ff1,ff2) 589 590 use defs_basis 591 use m_profiling_abi 592 use m_errors 593 594 use m_numeric_tools, only : ctrap 595 596 !This section has been created automatically by the script Abilint (TD). 597 !Do not modify the following lines by hand. 598 #undef ABI_FUNC 599 #define ABI_FUNC 'psden' 600 !End of the abilint section 601 602 implicit none 603 604 !Arguments ------------------------------------ 605 !scalars 606 integer,intent(in) :: ilog,mesh 607 real(dp),intent(in) :: rc 608 !arrays 609 real(dp),intent(in) :: nc(mesh),rad(mesh) 610 real(dp),intent(out) :: ff(mesh) 611 real(dp),intent(inout),optional :: ff1(mesh),ff2(mesh) 612 613 !Local variables------------------------------- 614 !scalars 615 integer :: ii,nc1 616 real(dp) :: aa,aa1,aa2,bb,cc,c1,c3,f0,f0p,norm1,norm2,rc1,step 617 !arrays 618 real(dp),allocatable :: fpir(:),gg(:) 619 620 ! ************************************************************************* 621 622 rc1=rc/four 623 624 ABI_ALLOCATE(fpir,(mesh)) 625 fpir(1:mesh)=four_pi*rad(1:mesh)**2 626 if (ilog==1) fpir(1:mesh)=fpir(1:mesh)*rad(1:mesh) 627 628 if (ilog==0) then 629 step=rad(2)-rad(1) 630 nc1=int(rc1/step)+1 631 rc1=(nc1-1)*step 632 else if (ilog==1) then 633 step=log(rad(2)/rad(1)) 634 nc1=int(log(rc1/rad(1))/step)+1 635 rc1=rad(nc1) 636 end if 637 ff(1:nc1)=nc(1:nc1)*fpir(1:nc1) 638 call ctrap(nc1,ff(1:nc1),step,c3) 639 if (ilog==1) c3=c3+half*ff(1) 640 f0=nc(nc1);c1=-log(f0) 641 f0p=half*(nc(nc1+1)-nc(nc1-1))/step 642 643 ii=0;aa1=zero;norm1=c3+one 644 do while (norm1>c3.and.ii<100) 645 ii=ii+1;aa1=aa1+one 646 aa=c1-aa1*rc1**4+rc1*(f0p/f0+four*aa1*rc1**3)*half 647 bb=-half*(f0p/f0+four*aa1*rc1**3)/rc1 648 ff(1:nc1)=fpir(1:nc1)*exp(-aa-bb*rad(1:nc1)**2-aa1*rad(1:nc1)**4) 649 call ctrap(nc1,ff(1:nc1),step,norm1) 650 if (ilog==1) norm1=norm1+half*ff(1) 651 end do 652 if (ii==100) then 653 MSG_ERROR('Big pb 1 in psden !') 654 end if 655 656 ii=0;aa2=zero;norm2=c3-one 657 do while (norm2<c3.and.ii<100) 658 ii=ii+1;aa2=aa2-one 659 aa=c1-aa2*rc1**4+rc1*(f0p/f0+four*aa2*rc1**3)*half 660 bb=-half*(f0p/f0+four*aa2*rc1**3)/rc1 661 ff(1:nc1)=fpir(1:nc1)*exp(-aa-bb*rad(1:nc1)**2-aa2*rad(1:nc1)**4) 662 call ctrap(nc1,ff(1:nc1),step,norm2) 663 if (ilog==1) norm2=norm2+half*ff(1) 664 end do 665 if (ii==100) then 666 MSG_ERROR('Big pb 2 in psden !') 667 end if 668 669 do while (abs(norm2-c3)>tol10) 670 671 cc=(aa1+aa2)*half 672 aa=c1-cc*rc1**4+rc1*(f0p/f0+four*cc*rc1**3)*half 673 bb=-half*(f0p/f0+four*cc*rc1**3)/rc1 674 ff(1:nc1)=fpir(1:nc1)*exp(-aa-bb*rad(1:nc1)**2-cc*rad(1:nc1)**4) 675 call ctrap (nc1,ff(1:nc1),step,norm2) 676 if (ilog==1) norm2=norm2+half*ff(1) 677 if ((norm1-c3)*(norm2-c3)>zero) then 678 aa1=cc 679 norm1=norm2 680 else 681 aa2=cc 682 end if 683 684 end do ! while 685 686 ff(1)=exp(-aa);if (ilog==1) ff(1)=ff(1)*exp(-bb*rad(1)**2-cc*rad(1)**4) 687 ff(2:nc1)=ff(2:nc1)/fpir(2:nc1) 688 if (nc1<mesh) ff(nc1+1:mesh)=nc(nc1+1:mesh) 689 if (present(ff1)) ff1(1:nc1)=-(two*bb*rad(1:nc1)+four*cc*rad(1:nc1)**3)*ff(1:nc1) 690 if (present(ff2)) ff2(1:nc1)=-(two*bb+12.0_dp*cc*rad(1:nc1)**2)*ff(1:nc1) & 691 & +(two*bb*rad(1:nc1)+four*cc*rad(1:nc1)**3)**2*ff(1:nc1) 692 693 ABI_ALLOCATE(gg,(mesh)) 694 gg(1:mesh)=fpir(1:mesh)*ff(1:mesh) 695 call ctrap(mesh,gg(1:mesh),step,norm1) 696 if (ilog==1) norm1=norm1+half*gg(1) 697 write(std_out,*) 'psden: tild_nc integral= ',norm1 698 ABI_DEALLOCATE(gg) 699 700 ABI_DEALLOCATE(fpir) 701 702 end subroutine psden
ABINIT/psp6cc [ Functions ]
NAME
psp6cc
FUNCTION
Compute the core charge density, for use in the XC core correction, following the function definition valid for the format 6 of pseudopotentials.
INPUTS
mmax=maximum number of points in real space grid in the psp file n1xccc=dimension of xccc1d ; 0 if no XC core correction is used rchrg=cut-off radius for the core density znucl=nuclear number of atom as specified in psp file
OUTPUT
xccc1d(n1xccc,6)= 1D core charge function and its five first derivatives
Optional output:
vh_tnzc(mmax) = Hartree potential induced by density tild_[n_Z+n_core]
(pseudized [n_Z+n_core], where n_Z=ions, n_core=core electrons)
using a simple pseudization scheme
PARENTS
psp6in
CHILDREN
psden,smooth,spline,splint,vhtnzc
SOURCE
373 subroutine psp6cc(mmax,n1xccc,rchrg,xccc1d,znucl,& 374 & vh_tnzc) ! optional argument 375 376 use defs_basis 377 use m_splines 378 use m_profiling_abi 379 380 use m_numeric_tools, only : smooth 381 382 !This section has been created automatically by the script Abilint (TD). 383 !Do not modify the following lines by hand. 384 #undef ABI_FUNC 385 #define ABI_FUNC 'psp6cc' 386 use interfaces_64_psp, except_this_one => psp6cc 387 !End of the abilint section 388 389 implicit none 390 391 !Arguments ------------------------------------ 392 !scalars 393 integer,intent(in) :: mmax,n1xccc 394 real(dp),intent(in) :: rchrg,znucl 395 !arrays 396 real(dp),intent(inout) :: xccc1d(n1xccc,6) !vz_i 397 real(dp),intent(out),optional :: vh_tnzc(mmax) 398 399 !Local variables------------------------------- 400 !scalars 401 integer :: i1xccc,irad 402 real(dp) :: der1,dern 403 character(len=500) :: errmsg 404 !arrays 405 real(dp),allocatable :: ff(:),ff1(:),ff2(:),ff3(:),gg(:),gg1(:),gg2(:),gg3(:) 406 real(dp),allocatable :: gg4(:),nc(:),rad(:),work(:),xx(:) 407 408 !********************************************************************** 409 410 ABI_ALLOCATE(ff,(mmax)) 411 ABI_ALLOCATE(ff1,(mmax)) 412 ABI_ALLOCATE(ff2,(mmax)) 413 ABI_ALLOCATE(ff3,(mmax)) 414 ABI_ALLOCATE(rad,(mmax)) 415 ABI_ALLOCATE(gg,(n1xccc)) 416 ABI_ALLOCATE(gg1,(n1xccc)) 417 ABI_ALLOCATE(gg2,(n1xccc)) 418 ABI_ALLOCATE(gg3,(n1xccc)) 419 ABI_ALLOCATE(gg4,(n1xccc)) 420 ABI_ALLOCATE(work,(n1xccc)) 421 ABI_ALLOCATE(xx,(n1xccc)) 422 423 ! 424 !read from pp file the model core charge (ff) and first (ff1) and 425 !second (ff2) derivative on logarithmic mesh mmax; rad is the radial grid 426 !the input functions contain the 4pi factor, it must be rescaled. 427 428 do irad=1,mmax 429 read(tmp_unit,*, err=10, iomsg=errmsg) rad(irad),ff(irad),ff1(irad),ff2(irad) 430 ff(irad)=ff(irad)/four_pi 431 ff1(irad)=ff1(irad)/four_pi 432 ff2(irad)=ff2(irad)/four_pi 433 end do 434 435 !Optional output: VHartree(tild_[n_Z+n_core]) 436 if (present(vh_tnzc)) then 437 ABI_ALLOCATE(nc,(mmax)) 438 nc=ff ! n_core 439 call psden(1,ff,mmax,nc,rchrg,rad,ff1=ff1,ff2=ff2) 440 call vhtnzc(ff,rchrg,vh_tnzc,mmax,rad,znucl) 441 ABI_DEALLOCATE(nc) 442 end if 443 444 rad(1)=zero 445 446 !calculate third derivative ff3 on logarithmic grid 447 der1=ff2(1) 448 dern=ff2(mmax) 449 call spline(rad,ff1,mmax,der1,dern,ff3) 450 451 !generate uniform mesh xx in the box cut by rchrg: 452 453 do i1xccc=1,n1xccc 454 xx(i1xccc)=(i1xccc-1)* rchrg/dble(n1xccc-1) 455 end do 456 457 ! 458 !now interpolate core charge and derivatives on the uniform grid 459 ! 460 !core charge, input=ff, output=gg 461 call splint(mmax,rad,ff,ff2,n1xccc,xx,gg) 462 463 !first derivative input=ff1, output=gg1 464 call splint(mmax,rad,ff1,ff3,n1xccc,xx,gg1) 465 466 !normalize gg1 467 gg1(:)=gg1(:)*rchrg 468 469 !now calculate second to fourth derivative by forward differences 470 !to avoid numerical noise uses a smoothing function 471 472 call smooth(gg1,n1xccc,10) 473 474 gg2(n1xccc)=zero 475 do i1xccc=1,n1xccc-1 476 gg2(i1xccc)=(gg1(i1xccc+1)-gg1(i1xccc))*dble(n1xccc-1) 477 end do 478 479 call smooth(gg2,n1xccc,10) 480 481 gg3(n1xccc)=zero 482 do i1xccc=1,n1xccc-1 483 gg3(i1xccc)=(gg2(i1xccc+1)-gg2(i1xccc))*dble(n1xccc-1) 484 end do 485 486 call smooth(gg3,n1xccc,10) 487 488 gg4(n1xccc)=zero 489 do i1xccc=1,n1xccc-1 490 gg4(i1xccc)=(gg3(i1xccc+1)-gg3(i1xccc))*dble(n1xccc-1) 491 end do 492 493 call smooth(gg4,n1xccc,10) 494 495 !write on xcc1d 496 xccc1d(:,1)=gg(:) 497 xccc1d(:,2)=gg1(:) 498 xccc1d(:,3)=gg2(:) 499 xccc1d(:,4)=gg3(:) 500 xccc1d(:,5)=gg4(:) 501 502 !WARNING : fifth derivative not yet computed 503 xccc1d(:,6)=zero 504 505 !DEBUG 506 !note: the normalization condition is the following: 507 !4pi rchrg /dble(n1xccc-1) sum xx^2 xccc1d(:,1) = qchrg 508 ! 509 !norm=zero 510 !do i1xccc=1,n1xccc 511 !norm = norm + four_pi*rchrg/dble(n1xccc-1)*& 512 !& xx(i1xccc)**2*xccc1d(i1xccc,1) 513 !end do 514 !write(std_out,*) ' norm=',norm 515 ! 516 !write(std_out,*)' psp1cc : output of core charge density and derivatives ' 517 !write(std_out,*)' xx gg gg1 ' 518 !do i1xccc=1,n1xccc 519 !write(10, '(3es14.6)' ) xx(i1xccc),xccc1d(i1xccc,1),xccc1d(i1xccc,2) 520 !end do 521 !write(std_out,*)' xx gg2 gg3 ' 522 !do i1xccc=1,n1xccc 523 !write(11, '(3es14.6)' ) xx(i1xccc),xccc1d(i1xccc,3),xccc1d(i1xccc,4) 524 !end do 525 !write(std_out,*)' xx gg4 gg5 ' 526 !do i1xccc=1,n1xccc 527 !write(12, '(3es14.6)' ) xx(i1xccc),xccc1d(i1xccc,5),xccc1d(i1xccc,6) 528 !end do 529 !write(std_out,*)' psp1cc : debug done, stop ' 530 !stop 531 !ENDDEBUG 532 533 ABI_DEALLOCATE(ff) 534 ABI_DEALLOCATE(ff1) 535 ABI_DEALLOCATE(ff2) 536 ABI_DEALLOCATE(ff3) 537 ABI_DEALLOCATE(gg) 538 ABI_DEALLOCATE(gg1) 539 ABI_DEALLOCATE(gg2) 540 ABI_DEALLOCATE(gg3) 541 ABI_DEALLOCATE(gg4) 542 ABI_DEALLOCATE(rad) 543 ABI_DEALLOCATE(work) 544 ABI_DEALLOCATE(xx) 545 546 return 547 548 ! Handle IO error 549 10 continue 550 MSG_ERROR(errmsg) 551 552 end subroutine psp6cc
ABINIT/psp6cc_drh [ Functions ]
NAME
psp6cc_drh
FUNCTION
Compute the core charge density, for use in the XC core correction, following the function definition valid for the format 6 of pseudopotentials. Version modified by DHamann, with consistent treatment of the derivatives in this routine and the remaining of the code.
INPUTS
mmax=maximum number of points in real space grid in the psp file n1xccc=dimension of xccc1d ; 0 if no XC core correction is used rchrg=cut-off radius for the core density
OUTPUT
xccc1d(n1xccc,6)= 1D core charge function and its five first derivatives
PARENTS
psp6in
CHILDREN
cc_derivatives
NOTES
Test version by DRH - requires very smooth model core charge
SOURCE
856 subroutine psp6cc_drh(mmax,n1xccc,rchrg,xccc1d) 857 858 use defs_basis 859 use m_profiling_abi 860 use m_errors 861 862 !This section has been created automatically by the script Abilint (TD). 863 !Do not modify the following lines by hand. 864 #undef ABI_FUNC 865 #define ABI_FUNC 'psp6cc_drh' 866 use interfaces_64_psp, except_this_one => psp6cc_drh 867 !End of the abilint section 868 869 implicit none 870 871 !Arguments ------------------------------------ 872 !scalars 873 integer,intent(in) :: mmax,n1xccc 874 real(dp),intent(in) :: rchrg 875 !arrays 876 real(dp),intent(inout) :: xccc1d(n1xccc,6) !vz_i 877 878 !Local variables------------------------------- 879 !scalars 880 integer :: irad 881 character(len=500) :: errmsg 882 !arrays 883 real(dp),allocatable :: ff(:),ff1(:),ff2(:),rad(:) 884 885 !********************************************************************** 886 887 ABI_ALLOCATE(ff,(mmax)) 888 ABI_ALLOCATE(ff1,(mmax)) 889 ABI_ALLOCATE(ff2,(mmax)) 890 ABI_ALLOCATE(rad,(mmax)) 891 892 ! 893 !read from pp file the model core charge (ff) and first (ff1) and 894 !second (ff2) derivative on logarithmic mesh mmax; rad is the radial grid 895 !the input functions contain the 4pi factor, it must be rescaled. 896 897 !***drh test 898 write(std_out,'(a,2i6)') 'drh:psp6cc_drh - mmax,n1xccc',mmax,n1xccc 899 !***end drh test 900 do irad=1,mmax 901 read(tmp_unit,*,err=10,iomsg=errmsg) rad(irad),ff(irad),ff1(irad),ff2(irad) 902 ff(irad)=ff(irad)/4.d0/pi 903 ff1(irad)=ff1(irad)/4.d0/pi 904 ff2(irad)=ff2(irad)/4.d0/pi 905 end do 906 rad(1)=0.d0 907 908 call cc_derivatives(rad,ff,ff1,ff2,mmax,n1xccc,rchrg,xccc1d) 909 910 ABI_DEALLOCATE(ff) 911 ABI_DEALLOCATE(ff1) 912 ABI_DEALLOCATE(ff2) 913 ABI_DEALLOCATE(rad) 914 915 return 916 917 ! Handle IO error 918 10 continue 919 MSG_ERROR(errmsg) 920 921 end subroutine psp6cc_drh
ABINIT/psp6in [ Functions ]
NAME
psp6in
FUNCTION
Initialize pspcod=6 (Pseudopotentials from the fhi98pp code): continue to read the corresponding file, then compute the local and non-local potentials.
COPYRIGHT
Copyright (C) 1999-2018 ABINIT group (XG, AF, GJ,FJ,MT, DRH) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt . For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
INPUTS
lloc=angular momentum choice of local pseudopotential
lmax=value of lmax mentioned at the second line of the psp file
lmnmax=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=max. number of (l,n) components over all type of psps
mmax=maximum number of points in real space grid in the psp file
angular momentum of nonlocal pseudopotential
mpsang= 1+maximum angular momentum for nonlocal pseudopotentials
mqgrid=dimension of q (or G) grid for arrays.
n1xccc=dimension of xccc1d ; 0 if no XC core correction is used
optnlxccc=option for nl XC core correction (input variable)
positron=0 if electron GS calculation
1 if positron GS calculation
2 if electron GS calculation in presence of the positron
qgrid(mqgrid)=values of q (or |G|) on grid from 0 to qmax
useylm=governs the way the nonlocal operator is to be applied:
1=using Ylm, 0=using Legendre polynomials
zion=nominal valence of atom as specified in psp file
znucl=nuclear number 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)+\frac{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 qchrg is not used, and could be suppressed later vlspl(mqgrid,2)=q^2 Vloc(q) and second derivatives from spline fit xcccrc=XC core correction cutoff radius (bohr) xccc1d(n1xccc,6)=1D core charge function and five derivatives, from psp file
PARENTS
pspatm
CHILDREN
psp5lo,psp5nl,psp6cc,psp6cc_drh,spline,wrtout
SOURCE
66 #if defined HAVE_CONFIG_H 67 #include "config.h" 68 #endif 69 70 #include "abi_common.h" 71 72 subroutine psp6in(ekb,epsatm,ffspl,indlmn,lloc,lmax,lmnmax,lnmax,& 73 & mmax,mpsang,mqgrid,nproj,n1xccc,optnlxccc,positron,qchrg,qgrid,& 74 & useylm,vlspl,xcccrc,xccc1d,zion,znucl) 75 76 use defs_basis 77 use m_splines 78 use m_errors 79 use m_profiling_abi 80 81 !This section has been created automatically by the script Abilint (TD). 82 !Do not modify the following lines by hand. 83 #undef ABI_FUNC 84 #define ABI_FUNC 'psp6in' 85 use interfaces_14_hidewrite 86 use interfaces_64_psp, except_this_one => psp6in 87 !End of the abilint section 88 89 implicit none 90 91 !Arguments ------------------------------------ 92 !scalars 93 integer,intent(in) :: lloc,lmax,lmnmax,lnmax,mmax,mpsang,mqgrid,n1xccc 94 integer,intent(in) :: optnlxccc,positron,useylm 95 real(dp),intent(in) :: zion,znucl 96 real(dp),intent(out) :: epsatm,qchrg,xcccrc 97 !arrays 98 integer,intent(out) :: indlmn(6,lmnmax),nproj(mpsang) 99 real(dp),intent(in) :: qgrid(mqgrid) 100 real(dp),intent(out) :: ekb(lnmax),vlspl(mqgrid,2) !vz_i 101 real(dp),intent(inout) :: ffspl(mqgrid,2,lnmax) !vz_i 102 real(dp),intent(inout) :: xccc1d(n1xccc,6) !vz_i 103 104 !Local variables------------------------------- 105 !scalars 106 integer :: ii,index,ipsang,irad,jj,jpsang,mm,mmax2 107 real(dp) :: al,al_announced,amesh,fchrg,ratio,rchrg,yp1,ypn 108 character(len=3) :: testxc 109 character(len=500) :: message,errmsg 110 !arrays 111 real(dp),allocatable :: ekb_tmp(:),ffspl_tmp(:,:,:),rad(:),vloc(:) 112 !real(dp),allocatable :: radbis 113 real(dp),allocatable :: vpspll(:,:),wfll(:,:),work_space(:),work_spl(:) 114 115 ! *************************************************************************** 116 117 !File format of formatted fhi psp input, as adapted for use 118 !by the ABINIT code (the 3 first lines 119 !have already been read in calling -pspatm- routine) : 120 121 !(1) title (character) line 122 !(2) znucl,zion,pspdat 123 !(3) pspcod,pspxc,lmax,lloc,mmax,r2well 124 !(4) rchrg,fchrg,qchrg 125 !Note : prior to version 2.2, this 4th line started with 4-- , 126 !and no core-correction was available. 127 !(5)-(18) -empty- 128 !(19) mmax, amesh ( mesh increment r(m+1)/r(m) ) 129 !Then, for ll=0,lmax : 130 !for irad=1,mmax : irad, r(irad), upsp(irad,ll), vpsp(irad,ll) 131 132 read (tmp_unit, '(a3)', err=10, iomsg=errmsg) testxc 133 if(testxc/='4--')then 134 backspace(tmp_unit) 135 read (tmp_unit,*, err=10, iomsg=errmsg) rchrg,fchrg,qchrg 136 write(message, '(3f20.14,t64,a)' ) rchrg,fchrg,qchrg,'rchrg,fchrg,qchrg' 137 call wrtout(ab_out,message,'COLL') 138 call wrtout(std_out, message,'COLL') 139 else 140 write(message, '(a)' ) ' No XC core correction.' 141 call wrtout(ab_out,message,'COLL') 142 call wrtout(std_out, message,'COLL') 143 rchrg=zero ; fchrg=zero ; qchrg=zero 144 end if 145 do ii=5,18 146 read(tmp_unit,*, err=10, iomsg=errmsg) 147 end do 148 149 if (positron==1.and.abs(fchrg)<=tol14) then 150 write(message,'(5a)')& 151 & 'You can only perform positronic ground-state calculations (positron=1)',ch10,& 152 & 'using fhi pseudopotentials with a core density (fchrg>0)',ch10,& 153 & 'Action: change your psp file (add fchrg>0).' 154 MSG_ERROR(message) 155 end if 156 !-------------------------------------------------------------------- 157 !Will now proceed at the reading of pots and wfs 158 159 !rad(:)=radial grid r(i) 160 !vpspll(:,1),...,vpspll(:,4)=nonlocal pseudopotentials 161 !wfll(:,1),...,wfll(:,4)=reference config. wavefunctions 162 163 ABI_ALLOCATE(rad,(mmax)) 164 ABI_ALLOCATE(vpspll,(mmax,mpsang)) 165 ABI_ALLOCATE(wfll,(mmax,mpsang)) 166 167 !Read atomic pseudopotential for each l, filling up arrays vpspll 168 !and wfll. Also set up rad array (actually read more than once) 169 !Note: put each l into vpspll(:,l+1) 170 do ipsang=1,lmax+1 171 nproj(ipsang)=1 172 read(tmp_unit,*, err=10, iomsg=errmsg)mmax2,amesh 173 if(ipsang==1)then 174 write(message, '(f10.6,t20,a)' ) amesh,' amesh (Hamman grid)' 175 al_announced=log(amesh) 176 call wrtout(ab_out,message,'COLL') 177 call wrtout(std_out, message,'COLL') 178 end if 179 do irad=1,mmax 180 read(tmp_unit,*, err=10, iomsg=errmsg)jj,rad(irad),wfll(irad,ipsang),vpspll(irad,ipsang) 181 ! DEBUG 182 ! Maybe the normalization is different 183 ! wfll(irad,ipsang)=wfll(irad,ipsang)/rad(irad) 184 ! ENDDEBUG 185 end do 186 end do 187 188 189 !Generate core charge function and derivatives, if needed 190 if(fchrg>tol14)then 191 192 if (positron==1) then 193 call psp6cc(mmax,n1xccc,rchrg,xccc1d,znucl,vh_tnzc=vpspll(:,lloc+1)) 194 else if(optnlxccc==1)then 195 call psp6cc(mmax,n1xccc,rchrg,xccc1d,znucl) 196 else if(optnlxccc==2)then 197 call psp6cc_drh(mmax,n1xccc,rchrg,xccc1d) 198 end if 199 200 ! The core charge function for pspcod=6 201 ! becomes zero beyond rchrg. Thus xcccrc must be set 202 ! equal to rchrg . 203 xcccrc=rchrg 204 else 205 xccc1d(:,:)=zero 206 xcccrc=zero 207 end if 208 209 !Compute in real(dp) al : the announced amesh is inaccurate. 210 ratio=rad(mmax)/rad(1) 211 al=log(ratio)/dble(mmax-1) 212 213 !DEBUG 214 !write(std_out,*)' psp6in : al ; al_announced =',al,al_announced 215 !allocate(radbis(mmax)) 216 !write(std_out,*)' psp6in : lloc ',lloc 217 !do ipsang=1,lmax+1 218 !write(std_out,*)' psp6in : ipsang ',ipsang 219 !do irad=1,mmax 220 !write(std_out,*)irad,rad(irad),wfll(irad,ipsang),vpspll(irad,ipsang) 221 !end do 222 !end do 223 !deallocate(radbis) 224 !ENDDEBUG 225 226 !vloc(:)=Vlocal(r), lloc=0, 1, or 2 or -1 for avg. 227 ABI_ALLOCATE(vloc,(mmax)) 228 !Copy appropriate nonlocal psp for use as local one 229 vloc( 1:mmax ) = vpspll( 1:mmax , lloc+1 ) 230 231 !-------------------------------------------------------------------- 232 !Carry out calculations for local (lloc) pseudopotential. 233 !Obtain Fourier transform (1-d sine transform) 234 !to get q^2 V(q). 235 236 call psp5lo(al,epsatm,mmax,mqgrid,qgrid,& 237 & vlspl(:,1),rad,vloc,yp1,ypn,zion) 238 239 !Fit spline to q^2 V(q) (Numerical Recipes subroutine) 240 ABI_ALLOCATE(work_space,(mqgrid)) 241 ABI_ALLOCATE(work_spl,(mqgrid)) 242 call spline (qgrid,vlspl(:,1),mqgrid,yp1,ypn,work_spl) 243 vlspl(:,2)=work_spl(:) 244 ABI_DEALLOCATE(work_space) 245 ABI_DEALLOCATE(work_spl) 246 247 !-------------------------------------------------------------------- 248 !Take care of non-local part 249 250 ABI_ALLOCATE(ekb_tmp,(mpsang)) 251 ABI_ALLOCATE(ffspl_tmp,(mqgrid,2,mpsang)) 252 253 !Zero out all Kleinman-Bylander energies to initialize 254 ekb_tmp(:)=zero 255 ekb(:)=zero 256 257 !Allow for option of no nonlocal corrections (lloc=lmax=0) 258 if (lloc==0.and.lmax==0) then 259 write(message, '(a,f5.1)' ) ' Note: local psp for atom with Z=',znucl 260 call wrtout(ab_out,message,'COLL') 261 call wrtout(std_out, message,'COLL') 262 else 263 264 ! ---------------------------------------------------------------------- 265 ! Compute KB form factors and fit splines 266 267 call psp5nl(al,ekb_tmp,ffspl_tmp,lmax,mmax,mpsang,mqgrid,qgrid,rad,vloc,& 268 & vpspll,wfll) 269 270 end if 271 272 jj=0;index=0;indlmn(:,:)=0 273 do ipsang=1,lmax+1 274 ! nproj had been set at 1, by default 275 if(abs(ekb_tmp(ipsang))<tol10)then 276 nproj(ipsang)=0 277 end if 278 ! Possible values for nproj in this routine : 0 or 1. 279 if(nproj(ipsang)==1)then 280 if (useylm==1) then 281 jj=jj+1 282 do mm=1,2*ipsang-1 283 index=index+1 284 indlmn(1,index)=ipsang-1 285 indlmn(2,index)=mm-ipsang 286 indlmn(3,index)=1 287 indlmn(4,index)=mm+(ipsang-1)*(ipsang-1) 288 indlmn(5,index)=jj 289 indlmn(6,index)=1 290 end do 291 else 292 jj=jj+1 293 index=index+1 294 indlmn(1,index)=ipsang-1 295 indlmn(2,index)=0 296 indlmn(3,index)=1 297 indlmn(4,index)=ipsang+(ipsang-1)*(ipsang-1) 298 indlmn(5,index)=jj 299 indlmn(6,index)=1 300 end if 301 end if 302 end do 303 !Transfer ekb and ffspl to their definitive location 304 jpsang=1 305 do ipsang=1,lmax+1 306 if(nproj(ipsang)/=0)then 307 ekb(jpsang)=ekb_tmp(ipsang) 308 ffspl(:,:,jpsang)=ffspl_tmp(:,:,ipsang) 309 jpsang=jpsang+1 310 if(jpsang>lnmax)then 311 write(message,'(3a,2i6)')& 312 & 'Problem with the dimension of the ekb and ffspl arrays.',ch10,& 313 & 'ipsang,lnmax=',ipsang,lnmax 314 end if 315 end if 316 end do 317 318 ABI_DEALLOCATE(ekb_tmp) 319 ABI_DEALLOCATE(ffspl_tmp) 320 321 !DEBUG 322 !write(std_out,*)' psp6in : exit ' 323 !write(std_out,*)' psp6in : indlmn(1:6,jj)' 324 !do jj=1,lmnmax 325 !write(std_out,*)indlmn(1:6,jj) 326 !end do 327 !ENDDEBUG 328 329 ABI_DEALLOCATE(vpspll) 330 ABI_DEALLOCATE(rad) 331 ABI_DEALLOCATE(vloc) 332 ABI_DEALLOCATE(wfll) 333 334 return 335 336 ! Handle IO error 337 10 continue 338 MSG_ERROR(errmsg) 339 340 end subroutine psp6in
ABINIT/vhtnzc [ Functions ]
NAME
vhtnzc
FUNCTION
Compute VHartree(tild[n_Z+n_core]) from input ncore
INPUTS
mesh=dimension of radial mesh nc= core density (to be pseudized) rad(mesh)=radial mesh rc=cut-off radius znucl=nuclear number of atom as specified in psp file
OUTPUT
vhtnzc(mesh) = hartree potential induced by density tild[n_Z+n_core] (pseudo core density + nucleus)
PARENTS
psp6cc
CHILDREN
ctrap
SOURCE
730 subroutine vhtnzc(nc,rc,vh_tnzc,mesh,rad,znucl) 731 732 use defs_basis 733 use m_profiling_abi 734 735 use m_numeric_tools, only : ctrap 736 737 !This section has been created automatically by the script Abilint (TD). 738 !Do not modify the following lines by hand. 739 #undef ABI_FUNC 740 #define ABI_FUNC 'vhtnzc' 741 !End of the abilint section 742 743 implicit none 744 745 !Arguments ------------------------------------ 746 !scalars 747 integer,intent(in) :: mesh 748 real(dp),intent(in) :: znucl 749 real(dp),intent(in) :: rc 750 !arrays 751 real(dp),intent(in) :: nc(mesh),rad(mesh) 752 real(dp),intent(out) :: vh_tnzc(mesh) 753 754 !Local variables------------------------------- 755 !scalars 756 integer :: ir,nc1 757 real(dp) :: gnorm,rc1,step,yp1,yp2,yp3 758 !arrays 759 real(dp),allocatable :: den1(:),den2(:),den3(:),den4(:),nzc(:),rvhn(:),shapefunc(:) 760 761 ! ************************************************************************* 762 763 rc1=rc/four 764 765 step=log(rad(2)/rad(1)) 766 nc1=int(log(rc1/rad(1))/step)+1 767 rc1=rad(nc1) 768 769 ABI_ALLOCATE(shapefunc,(mesh)) 770 shapefunc(1)=one 771 shapefunc(2:nc1)=(sin(pi*rad(2:nc1)/rc1)/(pi*rad(2:nc1)/rc1))**2 772 if (nc1<mesh) shapefunc(nc1+1:mesh)=zero 773 774 ABI_ALLOCATE(den1,(mesh)) 775 den1(1:mesh)=four_pi*shapefunc(1:mesh)*rad(1:mesh)**3 776 call ctrap(mesh,den1,step,gnorm) 777 gnorm =one/gnorm 778 ABI_DEALLOCATE(den1) 779 780 ABI_ALLOCATE(nzc,(mesh)) 781 nzc(1:mesh)=four*pi*nc(1:mesh)*rad(1:mesh)**2-four_pi*shapefunc(1:mesh)*rad(1:mesh)**2*znucl*gnorm 782 ABI_DEALLOCATE(shapefunc) 783 784 ABI_ALLOCATE(rvhn,(mesh)) 785 rvhn(1)=zero 786 787 ABI_ALLOCATE(den1,(mesh)) 788 ABI_ALLOCATE(den2,(mesh)) 789 ABI_ALLOCATE(den3,(mesh)) 790 ABI_ALLOCATE(den4,(mesh)) 791 792 den1(1)=zero;den2(1)=zero 793 do ir=2,mesh 794 den1(ir)= rad(ir)*nzc(ir) 795 den2(ir)= den1(ir)/rad(ir) 796 end do 797 798 !For first few points do stupid integral 799 den3(1)=zero;den4(1)=zero 800 do ir=2,mesh 801 call ctrap(ir,den1(1:ir),step,den3(ir)) 802 call ctrap(ir,den2(1:ir),step,den4(ir)) 803 end do 804 805 do ir=1,mesh 806 rvhn(ir)=den3(ir)+rad(ir)*(den4(mesh)-den4(ir)) 807 end do 808 809 ABI_DEALLOCATE(den1) 810 ABI_DEALLOCATE(den2) 811 ABI_DEALLOCATE(den3) 812 ABI_DEALLOCATE(den4) 813 814 vh_tnzc(2:mesh)=rvhn(2:mesh)/rad(2:mesh) 815 yp2=(vh_tnzc(3)-vh_tnzc(2))/(rad(3)-rad(2)) 816 yp3=(vh_tnzc(4)-vh_tnzc(3))/(rad(4)-rad(3)) 817 yp1=yp2+(yp2-yp3)*rad(2)/(rad(3)-rad(2)) 818 vh_tnzc(1)=vh_tnzc(2)-(yp1+yp2)*rad(2) 819 820 ABI_DEALLOCATE(nzc) 821 ABI_DEALLOCATE(rvhn) 822 823 end subroutine vhtnzc