TABLE OF CONTENTS
ABINIT/m_psp6 [ Modules ]
NAME
m_psp6
FUNCTION
Initialize pspcod=6 (Pseudopotentials from the fhi98pp code):
COPYRIGHT
Copyright (C) 1999-2024 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 .
SOURCE
16 #if defined HAVE_CONFIG_H 17 #include "config.h" 18 #endif 19 20 #include "abi_common.h" 21 22 module m_psp6 23 24 use defs_basis 25 use m_splines 26 use m_errors 27 use m_abicore 28 29 use m_numeric_tools, only : smooth, ctrap 30 use m_psptk, only : psp5lo, psp5nl, cc_derivatives 31 32 implicit none 33 34 private
m_psp6/psden [ Functions ]
[ Top ] [ m_psp6 ] [ 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))
SOURCE
535 subroutine psden(ilog,ff,mesh,nc,rc,rad,ff1,ff2) 536 537 !Arguments ------------------------------------ 538 !scalars 539 integer,intent(in) :: ilog,mesh 540 real(dp),intent(in) :: rc 541 !arrays 542 real(dp),intent(in) :: nc(mesh),rad(mesh) 543 real(dp),intent(out) :: ff(mesh) 544 real(dp),intent(inout),optional :: ff1(mesh),ff2(mesh) 545 546 !Local variables------------------------------- 547 !scalars 548 integer :: ii,nc1 549 real(dp) :: aa,aa1,aa2,bb,cc,c1,c3,f0,f0p,norm1,norm2,rc1,step 550 !arrays 551 real(dp),allocatable :: fpir(:),gg(:) 552 553 ! ************************************************************************* 554 555 rc1=rc/four 556 557 ABI_MALLOC(fpir,(mesh)) 558 fpir(1:mesh)=four_pi*rad(1:mesh)**2 559 if (ilog==1) fpir(1:mesh)=fpir(1:mesh)*rad(1:mesh) 560 561 if (ilog==0) then 562 step=rad(2)-rad(1) 563 nc1=int(rc1/step)+1 564 rc1=(nc1-1)*step 565 else if (ilog==1) then 566 step=log(rad(2)/rad(1)) 567 nc1=int(log(rc1/rad(1))/step)+1 568 rc1=rad(nc1) 569 end if 570 ff(1:nc1)=nc(1:nc1)*fpir(1:nc1) 571 call ctrap(nc1,ff(1:nc1),step,c3) 572 if (ilog==1) c3=c3+half*ff(1) 573 f0=nc(nc1);c1=-log(f0) 574 f0p=half*(nc(nc1+1)-nc(nc1-1))/step 575 576 ii=0;aa1=zero;norm1=c3+one 577 do while (norm1>c3.and.ii<100) 578 ii=ii+1;aa1=aa1+one 579 aa=c1-aa1*rc1**4+rc1*(f0p/f0+four*aa1*rc1**3)*half 580 bb=-half*(f0p/f0+four*aa1*rc1**3)/rc1 581 ff(1:nc1)=fpir(1:nc1)*exp(-aa-bb*rad(1:nc1)**2-aa1*rad(1:nc1)**4) 582 call ctrap(nc1,ff(1:nc1),step,norm1) 583 if (ilog==1) norm1=norm1+half*ff(1) 584 end do 585 if (ii==100) then 586 ABI_ERROR('Big pb 1 in psden !') 587 end if 588 589 ii=0;aa2=zero;norm2=c3-one 590 do while (norm2<c3.and.ii<100) 591 ii=ii+1;aa2=aa2-one 592 aa=c1-aa2*rc1**4+rc1*(f0p/f0+four*aa2*rc1**3)*half 593 bb=-half*(f0p/f0+four*aa2*rc1**3)/rc1 594 ff(1:nc1)=fpir(1:nc1)*exp(-aa-bb*rad(1:nc1)**2-aa2*rad(1:nc1)**4) 595 call ctrap(nc1,ff(1:nc1),step,norm2) 596 if (ilog==1) norm2=norm2+half*ff(1) 597 end do 598 if (ii==100) then 599 ABI_ERROR('Big pb 2 in psden !') 600 end if 601 602 do while (abs(norm2-c3)>tol10) 603 604 cc=(aa1+aa2)*half 605 aa=c1-cc*rc1**4+rc1*(f0p/f0+four*cc*rc1**3)*half 606 bb=-half*(f0p/f0+four*cc*rc1**3)/rc1 607 ff(1:nc1)=fpir(1:nc1)*exp(-aa-bb*rad(1:nc1)**2-cc*rad(1:nc1)**4) 608 call ctrap (nc1,ff(1:nc1),step,norm2) 609 if (ilog==1) norm2=norm2+half*ff(1) 610 if ((norm1-c3)*(norm2-c3)>zero) then 611 aa1=cc 612 norm1=norm2 613 else 614 aa2=cc 615 end if 616 617 end do ! while 618 619 ff(1)=exp(-aa);if (ilog==1) ff(1)=ff(1)*exp(-bb*rad(1)**2-cc*rad(1)**4) 620 ff(2:nc1)=ff(2:nc1)/fpir(2:nc1) 621 if (nc1<mesh) ff(nc1+1:mesh)=nc(nc1+1:mesh) 622 if (present(ff1)) ff1(1:nc1)=-(two*bb*rad(1:nc1)+four*cc*rad(1:nc1)**3)*ff(1:nc1) 623 if (present(ff2)) ff2(1:nc1)=-(two*bb+12.0_dp*cc*rad(1:nc1)**2)*ff(1:nc1) & 624 & +(two*bb*rad(1:nc1)+four*cc*rad(1:nc1)**3)**2*ff(1:nc1) 625 626 ABI_MALLOC(gg,(mesh)) 627 gg(1:mesh)=fpir(1:mesh)*ff(1:mesh) 628 call ctrap(mesh,gg(1:mesh),step,norm1) 629 if (ilog==1) norm1=norm1+half*gg(1) 630 !write(std_out,*) 'psden: tild_nc integral= ',norm1 631 ABI_FREE(gg) 632 633 ABI_FREE(fpir) 634 635 end subroutine psden
m_psp6/psp6cc [ Functions ]
[ Top ] [ m_psp6 ] [ 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
SOURCE
345 subroutine psp6cc(mmax,n1xccc,rchrg,xccc1d,znucl,& 346 & vh_tnzc) ! optional argument 347 348 !Arguments ------------------------------------ 349 !scalars 350 integer,intent(in) :: mmax,n1xccc 351 real(dp),intent(in) :: rchrg,znucl 352 !arrays 353 real(dp),intent(inout) :: xccc1d(n1xccc,6) !vz_i 354 real(dp),intent(out),optional :: vh_tnzc(mmax) 355 356 !Local variables------------------------------- 357 !scalars 358 integer :: i1xccc,irad 359 real(dp) :: der1,dern 360 character(len=500) :: errmsg 361 !arrays 362 real(dp),allocatable :: ff(:),ff1(:),ff2(:),ff3(:),gg(:),gg1(:),gg2(:),gg3(:) 363 real(dp),allocatable :: gg4(:),nc(:),rad(:),work(:),xx(:) 364 365 !********************************************************************** 366 367 ABI_MALLOC(ff,(mmax)) 368 ABI_MALLOC(ff1,(mmax)) 369 ABI_MALLOC(ff2,(mmax)) 370 ABI_MALLOC(ff3,(mmax)) 371 ABI_MALLOC(rad,(mmax)) 372 ABI_MALLOC(gg,(n1xccc)) 373 ABI_MALLOC(gg1,(n1xccc)) 374 ABI_MALLOC(gg2,(n1xccc)) 375 ABI_MALLOC(gg3,(n1xccc)) 376 ABI_MALLOC(gg4,(n1xccc)) 377 ABI_MALLOC(work,(n1xccc)) 378 ABI_MALLOC(xx,(n1xccc)) 379 380 !read from pp file the model core charge (ff) and first (ff1) and 381 !second (ff2) derivative on logarithmic mesh mmax; rad is the radial grid 382 !the input functions contain the 4pi factor, it must be rescaled. 383 384 do irad=1,mmax 385 read(tmp_unit,*, err=10, iomsg=errmsg) rad(irad),ff(irad),ff1(irad),ff2(irad) 386 ff(irad)=ff(irad)/four_pi 387 ff1(irad)=ff1(irad)/four_pi 388 ff2(irad)=ff2(irad)/four_pi 389 end do 390 391 !Optional output: VHartree(tild_[n_Z+n_core]) 392 if (present(vh_tnzc)) then 393 ABI_MALLOC(nc,(mmax)) 394 nc=ff ! n_core 395 call psden(1,ff,mmax,nc,rchrg,rad,ff1=ff1,ff2=ff2) 396 call vhtnzc(ff,rchrg,vh_tnzc,mmax,rad,znucl) 397 ABI_FREE(nc) 398 end if 399 400 rad(1)=zero 401 402 !calculate third derivative ff3 on logarithmic grid 403 der1=ff2(1) 404 dern=ff2(mmax) 405 call spline(rad,ff1,mmax,der1,dern,ff3) 406 407 !generate uniform mesh xx in the box cut by rchrg: 408 409 do i1xccc=1,n1xccc 410 xx(i1xccc)=(i1xccc-1)* rchrg/dble(n1xccc-1) 411 end do 412 413 !now interpolate core charge and derivatives on the uniform grid 414 !core charge, input=ff, output=gg 415 call splint(mmax,rad,ff,ff2,n1xccc,xx,gg) 416 417 !first derivative input=ff1, output=gg1 418 call splint(mmax,rad,ff1,ff3,n1xccc,xx,gg1) 419 420 !normalize gg1 421 gg1(:)=gg1(:)*rchrg 422 423 !now calculate second to fourth derivative by forward differences 424 !to avoid numerical noise uses a smoothing function 425 426 call smooth(gg1,n1xccc,10) 427 428 gg2(n1xccc)=zero 429 do i1xccc=1,n1xccc-1 430 gg2(i1xccc)=(gg1(i1xccc+1)-gg1(i1xccc))*dble(n1xccc-1) 431 end do 432 433 call smooth(gg2,n1xccc,10) 434 435 gg3(n1xccc)=zero 436 do i1xccc=1,n1xccc-1 437 gg3(i1xccc)=(gg2(i1xccc+1)-gg2(i1xccc))*dble(n1xccc-1) 438 end do 439 440 call smooth(gg3,n1xccc,10) 441 442 gg4(n1xccc)=zero 443 do i1xccc=1,n1xccc-1 444 gg4(i1xccc)=(gg3(i1xccc+1)-gg3(i1xccc))*dble(n1xccc-1) 445 end do 446 447 call smooth(gg4,n1xccc,10) 448 449 !write on xcc1d 450 xccc1d(:,1)=gg(:) 451 xccc1d(:,2)=gg1(:) 452 xccc1d(:,3)=gg2(:) 453 xccc1d(:,4)=gg3(:) 454 xccc1d(:,5)=gg4(:) 455 456 !WARNING : fifth derivative not yet computed 457 xccc1d(:,6)=zero 458 459 !note: the normalization condition is the following: 460 !4pi rchrg /dble(n1xccc-1) sum xx^2 xccc1d(:,1) = qchrg 461 ! 462 !norm=zero 463 !do i1xccc=1,n1xccc 464 !norm = norm + four_pi*rchrg/dble(n1xccc-1)*& 465 !& xx(i1xccc)**2*xccc1d(i1xccc,1) 466 !end do 467 !write(std_out,*) ' norm=',norm 468 ! 469 !write(std_out,*)' psp1cc : output of core charge density and derivatives ' 470 !write(std_out,*)' xx gg gg1 ' 471 !do i1xccc=1,n1xccc 472 !write(10, '(3es14.6)' ) xx(i1xccc),xccc1d(i1xccc,1),xccc1d(i1xccc,2) 473 !end do 474 !write(std_out,*)' xx gg2 gg3 ' 475 !do i1xccc=1,n1xccc 476 !write(11, '(3es14.6)' ) xx(i1xccc),xccc1d(i1xccc,3),xccc1d(i1xccc,4) 477 !end do 478 !write(std_out,*)' xx gg4 gg5 ' 479 !do i1xccc=1,n1xccc 480 !write(12, '(3es14.6)' ) xx(i1xccc),xccc1d(i1xccc,5),xccc1d(i1xccc,6) 481 !end do 482 !write(std_out,*)' psp1cc : debug done, stop ' 483 !stop 484 !ENDDEBUG 485 486 ABI_FREE(ff) 487 ABI_FREE(ff1) 488 ABI_FREE(ff2) 489 ABI_FREE(ff3) 490 ABI_FREE(gg) 491 ABI_FREE(gg1) 492 ABI_FREE(gg2) 493 ABI_FREE(gg3) 494 ABI_FREE(gg4) 495 ABI_FREE(rad) 496 ABI_FREE(work) 497 ABI_FREE(xx) 498 499 return 500 501 ! Handle IO error 502 10 continue 503 ABI_ERROR(errmsg) 504 505 end subroutine psp6cc
m_psp6/psp6cc_drh [ Functions ]
[ Top ] [ m_psp6 ] [ 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
NOTES
Test version by DRH - requires very smooth model core charge
SOURCE
764 subroutine psp6cc_drh(mmax,n1xccc,rchrg,xccc1d) 765 766 !Arguments ------------------------------------ 767 !scalars 768 integer,intent(in) :: mmax,n1xccc 769 real(dp),intent(in) :: rchrg 770 !arrays 771 real(dp),intent(inout) :: xccc1d(n1xccc,6) !vz_i 772 773 !Local variables------------------------------- 774 !scalars 775 integer :: irad 776 character(len=500) :: errmsg 777 !arrays 778 real(dp),allocatable :: ff(:),ff1(:),ff2(:),rad(:) 779 780 !********************************************************************** 781 782 ABI_MALLOC(ff,(mmax)) 783 ABI_MALLOC(ff1,(mmax)) 784 ABI_MALLOC(ff2,(mmax)) 785 ABI_MALLOC(rad,(mmax)) 786 787 ! 788 !read from pp file the model core charge (ff) and first (ff1) and 789 !second (ff2) derivative on logarithmic mesh mmax; rad is the radial grid 790 !the input functions contain the 4pi factor, it must be rescaled. 791 792 !write(std_out,'(a,2i6)') 'drh:psp6cc_drh - mmax,n1xccc',mmax,n1xccc 793 do irad=1,mmax 794 read(tmp_unit,*,err=10,iomsg=errmsg) rad(irad),ff(irad),ff1(irad),ff2(irad) 795 ff(irad)=ff(irad)/4.d0/pi 796 ff1(irad)=ff1(irad)/4.d0/pi 797 ff2(irad)=ff2(irad)/4.d0/pi 798 end do 799 rad(1)=0.d0 800 801 call cc_derivatives(rad,ff,ff1,ff2,mmax,n1xccc,rchrg,xccc1d) 802 803 ABI_FREE(ff) 804 ABI_FREE(ff1) 805 ABI_FREE(ff2) 806 ABI_FREE(rad) 807 808 return 809 810 ! Handle IO error 811 10 continue 812 ABI_ERROR(errmsg) 813 814 end subroutine psp6cc_drh
m_psp6/psp6in [ Functions ]
[ Top ] [ m_psp6 ] [ 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.
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
SOURCE
94 subroutine psp6in(ekb,epsatm,ffspl,indlmn,lloc,lmax,lmnmax,lnmax,& 95 & mmax,mpsang,mqgrid,nproj,n1xccc,optnlxccc,positron,qchrg,qgrid,& 96 & useylm,vlspl,xcccrc,xccc1d,zion,znucl) 97 98 !Arguments ------------------------------------ 99 !scalars 100 integer,intent(in) :: lloc,lmax,lmnmax,lnmax,mmax,mpsang,mqgrid,n1xccc 101 integer,intent(in) :: optnlxccc,positron,useylm 102 real(dp),intent(in) :: zion,znucl 103 real(dp),intent(out) :: epsatm,qchrg,xcccrc 104 !arrays 105 integer,intent(out) :: indlmn(6,lmnmax),nproj(mpsang) 106 real(dp),intent(in) :: qgrid(mqgrid) 107 real(dp),intent(out) :: ekb(lnmax),vlspl(mqgrid,2) !vz_i 108 real(dp),intent(inout) :: ffspl(mqgrid,2,lnmax) !vz_i 109 real(dp),intent(inout) :: xccc1d(n1xccc,6) !vz_i 110 111 !Local variables------------------------------- 112 !scalars 113 integer :: ii,index,ipsang,irad,jj,jpsang,mm,mmax2 114 real(dp) :: al,al_announced,amesh,fchrg,ratio,rchrg,yp1,ypn 115 character(len=3) :: testxc 116 character(len=500) :: message,errmsg 117 !arrays 118 real(dp),allocatable :: ekb_tmp(:),ffspl_tmp(:,:,:),rad(:),vloc(:) 119 !real(dp),allocatable :: radbis 120 real(dp),allocatable :: vpspll(:,:),wfll(:,:),work_space(:),work_spl(:) 121 122 ! *************************************************************************** 123 124 !File format of formatted fhi psp input, as adapted for use 125 !by the ABINIT code (the 3 first lines 126 !have already been read in calling -pspatm- routine) : 127 128 !(1) title (character) line 129 !(2) znucl,zion,pspdat 130 !(3) pspcod,pspxc,lmax,lloc,mmax,r2well 131 !(4) rchrg,fchrg,qchrg 132 !Note : prior to version 2.2, this 4th line started with 4-- , 133 !and no core-correction was available. 134 !(5)-(18) -empty- 135 !(19) mmax, amesh ( mesh increment r(m+1)/r(m) ) 136 !Then, for ll=0,lmax : 137 !for irad=1,mmax : irad, r(irad), upsp(irad,ll), vpsp(irad,ll) 138 139 read (tmp_unit, '(a3)', err=10, iomsg=errmsg) testxc 140 if(testxc/='4--')then 141 backspace(tmp_unit) 142 read (tmp_unit,*, err=10, iomsg=errmsg) rchrg,fchrg,qchrg 143 write(message, '(3f20.14,t64,a)' ) rchrg,fchrg,qchrg,'rchrg,fchrg,qchrg' 144 call wrtout(ab_out,message,'COLL') 145 call wrtout(std_out, message,'COLL') 146 else 147 write(message, '(a)' ) ' No XC core correction.' 148 call wrtout(ab_out,message,'COLL') 149 call wrtout(std_out, message,'COLL') 150 rchrg=zero ; fchrg=zero ; qchrg=zero 151 end if 152 do ii=5,18 153 read(tmp_unit,*, err=10, iomsg=errmsg) 154 end do 155 156 if (positron==1.and.abs(fchrg)<=tol14) then 157 write(message,'(5a)')& 158 & 'You can only perform positronic ground-state calculations (positron=1)',ch10,& 159 & 'using fhi pseudopotentials with a core density (fchrg>0)',ch10,& 160 & 'Action: change your psp file (add fchrg>0).' 161 ABI_ERROR(message) 162 end if 163 !-------------------------------------------------------------------- 164 !Will now proceed at the reading of pots and wfs 165 166 !rad(:)=radial grid r(i) 167 !vpspll(:,1),...,vpspll(:,4)=nonlocal pseudopotentials 168 !wfll(:,1),...,wfll(:,4)=reference config. wavefunctions 169 170 ABI_MALLOC(rad,(mmax)) 171 ABI_MALLOC(vpspll,(mmax,mpsang)) 172 ABI_MALLOC(wfll,(mmax,mpsang)) 173 174 !Read atomic pseudopotential for each l, filling up arrays vpspll 175 !and wfll. Also set up rad array (actually read more than once) 176 !Note: put each l into vpspll(:,l+1) 177 do ipsang=1,lmax+1 178 nproj(ipsang)=1 179 read(tmp_unit,*, err=10, iomsg=errmsg)mmax2,amesh 180 if(ipsang==1)then 181 write(message, '(f10.6,t20,a)' ) amesh,' amesh (Hamman grid)' 182 al_announced=log(amesh) 183 call wrtout(ab_out,message,'COLL') 184 call wrtout(std_out, message,'COLL') 185 end if 186 do irad=1,mmax 187 read(tmp_unit,*, err=10, iomsg=errmsg)jj,rad(irad),wfll(irad,ipsang),vpspll(irad,ipsang) 188 end do 189 end do 190 191 192 !Generate core charge function and derivatives, if needed 193 if(fchrg>tol14)then 194 195 if (positron==1) then 196 call psp6cc(mmax,n1xccc,rchrg,xccc1d,znucl,vh_tnzc=vpspll(:,lloc+1)) 197 else if(optnlxccc==1)then 198 call psp6cc(mmax,n1xccc,rchrg,xccc1d,znucl) 199 else if(optnlxccc==2)then 200 call psp6cc_drh(mmax,n1xccc,rchrg,xccc1d) 201 end if 202 203 ! The core charge function for pspcod=6 becomes zero beyond rchrg. 204 ! Thus xcccrc must be set equal to rchrg. 205 xcccrc=rchrg 206 else 207 xccc1d(:,:)=zero 208 xcccrc=zero 209 end if 210 211 !Compute in real(dp) al : the announced amesh is inaccurate. 212 ratio=rad(mmax)/rad(1) 213 al=log(ratio)/dble(mmax-1) 214 215 !vloc(:)=Vlocal(r), lloc=0, 1, or 2 or -1 for avg. 216 ABI_MALLOC(vloc,(mmax)) 217 !Copy appropriate nonlocal psp for use as local one 218 vloc( 1:mmax ) = vpspll( 1:mmax , lloc+1 ) 219 220 !-------------------------------------------------------------------- 221 !Carry out calculations for local (lloc) pseudopotential. 222 !Obtain Fourier transform (1-d sine transform) to get q^2 V(q). 223 224 call psp5lo(al,epsatm,mmax,mqgrid,qgrid,& 225 & vlspl(:,1),rad,vloc,yp1,ypn,zion) 226 227 !Fit spline to q^2 V(q) (Numerical Recipes subroutine) 228 ABI_MALLOC(work_space,(mqgrid)) 229 ABI_MALLOC(work_spl,(mqgrid)) 230 call spline (qgrid,vlspl(:,1),mqgrid,yp1,ypn,work_spl) 231 vlspl(:,2)=work_spl(:) 232 ABI_FREE(work_space) 233 ABI_FREE(work_spl) 234 235 !-------------------------------------------------------------------- 236 !Take care of non-local part 237 238 ABI_MALLOC(ekb_tmp,(mpsang)) 239 ABI_MALLOC(ffspl_tmp,(mqgrid,2,mpsang)) 240 241 !Zero out all Kleinman-Bylander energies to initialize 242 ekb_tmp(:)=zero 243 ekb(:)=zero 244 245 !Allow for option of no nonlocal corrections (lloc=lmax=0) 246 if (lloc==0.and.lmax==0) then 247 write(message, '(a,f5.1)' ) ' Note: local psp for atom with Z=',znucl 248 call wrtout(ab_out,message,'COLL') 249 call wrtout(std_out, message,'COLL') 250 else 251 252 ! ---------------------------------------------------------------------- 253 ! Compute KB form factors and fit splines 254 255 call psp5nl(al,ekb_tmp,ffspl_tmp,lmax,mmax,mpsang,mqgrid,qgrid,rad,vloc, vpspll,wfll) 256 257 end if 258 259 jj=0;index=0;indlmn(:,:)=0 260 do ipsang=1,lmax+1 261 ! nproj had been set at 1, by default 262 if(abs(ekb_tmp(ipsang))<tol10)then 263 nproj(ipsang)=0 264 end if 265 ! Possible values for nproj in this routine : 0 or 1. 266 if(nproj(ipsang)==1)then 267 if (useylm==1) then 268 jj=jj+1 269 do mm=1,2*ipsang-1 270 index=index+1 271 indlmn(1,index)=ipsang-1 272 indlmn(2,index)=mm-ipsang 273 indlmn(3,index)=1 274 indlmn(4,index)=mm+(ipsang-1)*(ipsang-1) 275 indlmn(5,index)=jj 276 indlmn(6,index)=1 277 end do 278 else 279 jj=jj+1 280 index=index+1 281 indlmn(1,index)=ipsang-1 282 indlmn(2,index)=0 283 indlmn(3,index)=1 284 indlmn(4,index)=ipsang+(ipsang-1)*(ipsang-1) 285 indlmn(5,index)=jj 286 indlmn(6,index)=1 287 end if 288 end if 289 end do 290 !Transfer ekb and ffspl to their definitive location 291 jpsang=1 292 do ipsang=1,lmax+1 293 if(nproj(ipsang)/=0)then 294 ekb(jpsang)=ekb_tmp(ipsang) 295 ffspl(:,:,jpsang)=ffspl_tmp(:,:,ipsang) 296 jpsang=jpsang+1 297 if(jpsang>lnmax)then 298 write(message,'(3a,2i6)')& 299 & 'Problem with the dimension of the ekb and ffspl arrays.',ch10,& 300 & 'ipsang,lnmax=',ipsang,lnmax 301 end if 302 end if 303 end do 304 305 ABI_FREE(ekb_tmp) 306 ABI_FREE(ffspl_tmp) 307 ABI_FREE(vpspll) 308 ABI_FREE(rad) 309 ABI_FREE(vloc) 310 ABI_FREE(wfll) 311 312 return 313 314 ! Handle IO error 315 10 continue 316 ABI_ERROR(errmsg) 317 318 end subroutine psp6in
m_psp6/vhtnzc [ Functions ]
[ Top ] [ m_psp6 ] [ 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)
SOURCE
657 subroutine vhtnzc(nc,rc,vh_tnzc,mesh,rad,znucl) 658 659 !Arguments ------------------------------------ 660 !scalars 661 integer,intent(in) :: mesh 662 real(dp),intent(in) :: znucl 663 real(dp),intent(in) :: rc 664 !arrays 665 real(dp),intent(in) :: nc(mesh),rad(mesh) 666 real(dp),intent(out) :: vh_tnzc(mesh) 667 668 !Local variables------------------------------- 669 !scalars 670 integer :: ir,nc1 671 real(dp) :: gnorm,rc1,step,yp1,yp2,yp3 672 !arrays 673 real(dp),allocatable :: den1(:),den2(:),den3(:),den4(:),nzc(:),rvhn(:),shapefunc(:) 674 675 ! ************************************************************************* 676 677 rc1=rc/four 678 679 step=log(rad(2)/rad(1)) 680 nc1=int(log(rc1/rad(1))/step)+1 681 rc1=rad(nc1) 682 683 ABI_MALLOC(shapefunc,(mesh)) 684 shapefunc(1)=one 685 shapefunc(2:nc1)=(sin(pi*rad(2:nc1)/rc1)/(pi*rad(2:nc1)/rc1))**2 686 if (nc1<mesh) shapefunc(nc1+1:mesh)=zero 687 688 ABI_MALLOC(den1,(mesh)) 689 den1(1:mesh)=four_pi*shapefunc(1:mesh)*rad(1:mesh)**3 690 call ctrap(mesh,den1,step,gnorm) 691 gnorm =one/gnorm 692 ABI_FREE(den1) 693 694 ABI_MALLOC(nzc,(mesh)) 695 nzc(1:mesh)=four*pi*nc(1:mesh)*rad(1:mesh)**2-four_pi*shapefunc(1:mesh)*rad(1:mesh)**2*znucl*gnorm 696 ABI_FREE(shapefunc) 697 698 ABI_MALLOC(rvhn,(mesh)) 699 rvhn(1)=zero 700 701 ABI_MALLOC(den1,(mesh)) 702 ABI_MALLOC(den2,(mesh)) 703 ABI_MALLOC(den3,(mesh)) 704 ABI_MALLOC(den4,(mesh)) 705 706 den1(1)=zero;den2(1)=zero 707 do ir=2,mesh 708 den1(ir)= rad(ir)*nzc(ir) 709 den2(ir)= den1(ir)/rad(ir) 710 end do 711 712 !For first few points do stupid integral 713 den3(1)=zero;den4(1)=zero 714 do ir=2,mesh 715 call ctrap(ir,den1(1:ir),step,den3(ir)) 716 call ctrap(ir,den2(1:ir),step,den4(ir)) 717 end do 718 719 do ir=1,mesh 720 rvhn(ir)=den3(ir)+rad(ir)*(den4(mesh)-den4(ir)) 721 end do 722 723 ABI_FREE(den1) 724 ABI_FREE(den2) 725 ABI_FREE(den3) 726 ABI_FREE(den4) 727 728 vh_tnzc(2:mesh)=rvhn(2:mesh)/rad(2:mesh) 729 yp2=(vh_tnzc(3)-vh_tnzc(2))/(rad(3)-rad(2)) 730 yp3=(vh_tnzc(4)-vh_tnzc(3))/(rad(4)-rad(3)) 731 yp1=yp2+(yp2-yp3)*rad(2)/(rad(3)-rad(2)) 732 vh_tnzc(1)=vh_tnzc(2)-(yp1+yp2)*rad(2) 733 734 ABI_FREE(nzc) 735 ABI_FREE(rvhn) 736 737 end subroutine vhtnzc