TABLE OF CONTENTS
ABINIT/m_pspini [ Modules ]
NAME
m_pspini
FUNCTION
Initialize pseudopotential datastructures from files.
COPYRIGHT
Copyright (C) 1998-2024 ABINIT group (DCA, XG, GMR, MT, FrD, AF, DRH, YP) 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_pspini 23 24 use defs_basis 25 use m_errors 26 use m_abicore 27 use m_xmpi 28 use m_psxml2ab 29 use m_dtset 30 use m_dtfil 31 !use m_psps 32 33 use defs_datatypes, only : pseudopotential_type, nctab_t, pspheader_type 34 use m_fstrings, only : ftoa, sjoin 35 use m_time, only : timab 36 use m_io_tools, only : open_file 37 use m_pawrad, only : pawrad_type 38 use m_pawtab, only : pawtab_type, pawtab_set_flags 39 use m_psps, only : psps_print, psps_ncwrite_path, nctab_init, nctab_free, nctab_mixalch, test_xml_xmlpaw_upf, & 40 nctab_eval_tcorespl 41 use m_pawpsp, only : pawpsp_bcast, pawpsp_read_pawheader, pawpsp_read_header_xml,& 42 pawpsp_header_type, pawpsp_wvl, pawpsp_7in, pawpsp_17in 43 use m_pawxmlps, only : paw_setup_free,paw_setuploc 44 use m_pspheads, only : pawpsxml2ab 45 #if defined HAVE_BIGDFT 46 use BigDFT_API, only : dictionary, atomic_info, dict_init, dict_free, UNINITIALIZED 47 #endif 48 49 use m_psp1, only : psp1in 50 use m_psp5, only : psp5in 51 use m_psp6, only : psp6in 52 use m_psp8, only : psp8in 53 use m_psp9, only : psp9in 54 use m_upf2abinit, only : upf1_to_abinit, upf2_to_abinit 55 use m_psp_hgh, only : psp2in, psp3in, psp10in 56 use m_wvl_descr_psp, only : wvl_descr_psp_fill 57 58 implicit none 59 60 private
ABINIT/psp_dump_outputs [ Functions ]
NAME
psp_dump_outputs
FUNCTION
Debugging routines used to dumo PSP data in Yaml format.
SOURCE
1359 subroutine psp_dump_outputs(pfx,pspcod,lmnmax,lnmax,mpssoang, & 1360 mqgrid,n1xccc,mmax,maxrad,epsatm,qchrg,xcccrc,nctab, & 1361 indlmn,nproj,ekb,ffspl,vlspl,xccc1d) 1362 1363 !Arguments ------------------------------------ 1364 !scalars 1365 character(len=*), intent(in) :: pfx 1366 integer,intent(in) :: pspcod,lmnmax,lnmax,mpssoang,mqgrid,n1xccc 1367 integer,intent(in) :: mmax 1368 real(dp),intent(in) :: maxrad,epsatm,qchrg,xcccrc 1369 type(nctab_t),intent(in) :: nctab 1370 !arrays 1371 integer,intent(in) :: indlmn(6,lmnmax),nproj(mpssoang) 1372 real(dp),intent(in) :: ekb(lnmax),ffspl(mqgrid,2,lnmax),vlspl(mqgrid,2) 1373 real(dp),intent(in) :: xccc1d(n1xccc,6) 1374 1375 !Local variables ------------------------------ 1376 !scalars 1377 integer, parameter :: dump = 64 1378 integer :: ierr, i, j ,k 1379 character(len=500) :: msg 1380 1381 ! ********************************************************************* 1382 1383 open(unit=dump, file=trim(pfx)//"_psp_info.yaml", status='REPLACE', err=10, iostat=ierr) 1384 1385 write(dump,'(3a)') "%YAML 1.2", ch10, "---" 1386 1387 write(dump, '(2a)') ch10, "# Pseudopotential info" 1388 write(dump, '(a,1x,i8)') "pspcod:", pspcod 1389 1390 write(dump, '(2a)') ch10, "# Array dimensions" 1391 write(dump, '(a)') "dims:" 1392 write(dump, '(4x,a,1x,i8)') "lmnmax:", lmnmax 1393 write(dump, '(4x,a,1x,i8)') "lnmax:", lnmax 1394 write(dump, '(4x,a,1x,i8)') "mpssoang:", mpssoang 1395 write(dump, '(4x,a,1x,i8)') "mqgrid:", mqgrid 1396 write(dump, '(4x,a,1x,i8)') "n1xccc:", n1xccc 1397 write(dump, '(4x,a,1x,i8)') "mmax:", mmax 1398 1399 write(dump, '(2a)') ch10, "# Quantities" 1400 write(dump, '(a,1x,e12.5)') "maxrad:", maxrad 1401 write(dump, '(a,1x,e12.5)') "epsatm:", epsatm 1402 write(dump, '(a,1x,e12.5)') "qchrg:", qchrg 1403 write(dump, '(a,1x,e12.5)') "xcccrc:", xcccrc 1404 1405 write(dump, '(2a)') ch10, "# Structure: nctab" 1406 write(dump, '(a)') "nctab:" 1407 write(dump,'(4x,a,":",1x,i4)') "mqgrid_vl", nctab%mqgrid_vl 1408 write(dump,'(4x,a,":",1x,l4)') "has_tvale", nctab%has_tvale 1409 write(dump,'(4x,a,":",1x,l4)') "has_tcore", nctab%has_tcore 1410 write(dump,'(4x,a,":",1x,e12.5)') "dncdq0", nctab%dncdq0 1411 write(dump,'(4x,a,":",1x,e12.5)') "d2ncdq0", nctab%d2ncdq0 1412 write(dump,'(4x,a,":",1x,e12.5)') "dnvdq0", nctab%dnvdq0 1413 1414 if ( nctab%has_tvale ) then 1415 write(dump, '(2a)') ch10, "# Array: nctab_tvalespl(mqgrid_vl,2)" 1416 write(dump, '(a)') "nctab_tvalespl:" 1417 do j=1,2 1418 do i=1,nctab%mqgrid_vl 1419 if ( i == 1 ) then 1420 write(dump,'(4x,a,1x,e12.5)') "- -", nctab%tvalespl(i,j) 1421 else 1422 write(dump,'(4x,a,1x,e12.5)') " -", nctab%tvalespl(i,j) 1423 end if 1424 end do 1425 end do 1426 end if 1427 1428 if ( nctab%has_tcore ) then 1429 write(dump, '(2a)') ch10, "# Array: nctab_tcorespl(mqgrid_vl,2)" 1430 write(dump, '(a)') "nctab_tcorespl:" 1431 do j=1,2 1432 do i=1,nctab%mqgrid_vl 1433 if ( i == 1 ) then 1434 write(dump,'(4x,a,1x,e12.5)') "- -", nctab%tcorespl(i,j) 1435 else 1436 write(dump,'(4x,a,1x,e12.5)') " -", nctab%tcorespl(i,j) 1437 end if 1438 end do 1439 end do 1440 end if 1441 1442 write(dump, '(2a)') ch10, "# Array: integer indlmn(6,lmnmax)" 1443 write(dump, '(a)') "indlmn:" 1444 do i=1,lmnmax 1445 write(dump,'(4x,a,i4,5(",",i4),a)') "- [", indlmn(:,i), "]" 1446 end do 1447 1448 write(dump, '(2a)') ch10, "# Array: integer nproj(mpssoang)" 1449 write(dump, '(a)') "nproj:" 1450 do i=1,mpssoang 1451 write(dump,'(4x,"-",1x,i4)') nproj(i) 1452 end do 1453 1454 write(dump, '(2a)') ch10, "# Array: double ekb(lnmax)" 1455 write(dump, '(a)') "ekb:" 1456 do i=1,lnmax 1457 write(dump,'(4x,"-",1x,e12.5)') ekb(i) 1458 end do 1459 1460 write(dump, '(2a)') ch10, "# Array: ffspl(mqgrid,2,lnmax)" 1461 write(dump, '(a)') "ffspl:" 1462 do k=1,lnmax 1463 do j=1,2 1464 do i=1,mqgrid 1465 if ( (i == 1) .and. (j == 1) ) then 1466 write(dump,'(4x,a,1x,e12.5)') "- - -", ffspl(i,j,k) 1467 else if ( i == 1 ) then 1468 write(dump,'(4x,a,1x,e12.5)') " - -", ffspl(i,j,k) 1469 else 1470 write(dump,'(4x,a,1x,e12.5)') " -", ffspl(i,j,k) 1471 end if 1472 end do 1473 end do 1474 end do 1475 1476 write(dump, '(2a)') ch10, "# Array: vlspl(mqgrid,2)" 1477 write(dump, '(a)') "vlspl:" 1478 do j=1,2 1479 do i=1,mqgrid 1480 if ( i == 1 ) then 1481 write(dump,'(4x,a,1x,e12.5)') "- -", vlspl(i,j) 1482 else 1483 write(dump,'(4x,a,1x,e12.5)') " -", vlspl(i,j) 1484 end if 1485 end do 1486 end do 1487 1488 write(dump, '(2a)') ch10, "# Array: xccc1d(n1xccc,6)" 1489 write(dump, '(a)') "xccc1d:" 1490 do j=1,6 1491 do i=1,mqgrid 1492 if ( i == 1 ) then 1493 write(dump,'(4x,a,1x,e12.5)') "- -", xccc1d(i,j) 1494 else 1495 write(dump,'(4x,a,1x,e12.5)') " -", xccc1d(i,j) 1496 end if 1497 end do 1498 end do 1499 1500 write (dump,'(2a)') ch10, "..." 1501 1502 close(dump) 1503 1504 return 1505 10 continue 1506 1507 if ( ierr /= 0 ) then 1508 write(msg,'(a,a,a,i8)') "Error writing pseudopotential information", ch10, "IOSTAT=", ierr 1509 ABI_WARNING(msg) 1510 end if 1511 1512 end subroutine psp_dump_outputs
ABINIT/pspatm [ Functions ]
NAME
pspatm
FUNCTION
Open atomic pseudopotential data file for a given atom, read the three first lines, make some checks, then call appropriate subroutine for the reading of V(r) and wf R(r) data for each angular momentum, and subsequent Fourier and Bessel function transforms for local and nonlocal potentials. Close psp file at end. Handles pseudopotential files produced by (pspcod=1 or 4) Teter code, or from the Goedecker-Teter-Hutter paper (pspcod=2), or from the Hartwigsen-Goedecker-Hutter paper (pspcod=3 or 10) or "Phoney pseudopotentials" (Hamman grid in real space) (pspcod=5) or "Troullier-Martins pseudopotentials" from the FHI (pspcod=6) or "XML format" (pspcod=9) or "UPF1 PWSCF format" (pspcod=11) or "UPF2 PWSCF format" (pspcod=12)
INPUTS
dq= spacing of the q-grid dtset <type(dataset_type)>=all input variables in this dataset | ixc=exchange-correlation choice from main routine data file | pawxcdev=choice of XC development in PAW formalism | usexcnhat_orig=choice for use of nhat in Vxc in PAW formalism | xclevel= XC functional level ipsp=id in the array of the currently read pseudo.
OUTPUT
ekb(dimekb)=
->NORM-CONSERVING PSPS ONLY (pspcod/=7):
(Real) Kleinman-Bylander energies (hartree)
{{\ \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 number of basis functions (l,n) (dimekb=lnmax)
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)
indlmn(6,i)= array giving l,m,n,lm,ln,s for i=ln (if useylm=0)
or i=lmn (if useylm=1)
pawrad <type(pawrad_type)>=paw radial mesh and related data
pawtab <type(pawtab_type)>=paw tabulated starting data
vlspl(mqgrid_vl,2)=q^2 Vloc(q) and second derivatives from spline fit
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
xcccrc=XC core correction cutoff radius (bohr) from psp file
xccc1d(n1xccc*(1-usepaw),6)=1D core charge function and five derivatives, from psp file (used in NC only)
nctab=<nctab_t>
has_tvale=True if the pseudo provides the valence density (used in NC only)
tvalespl(mqgrid_vl(1-usepaw),2)=the pseudo valence density and 2nd derivative in reciprocal space on a regular grid
(used in NC only)
SIDE EFFECTS
Input/Output : psps <type(pseudopotential_type)>=at output, values depending on the read pseudo are set. | dimekb(IN)=dimension of ekb (see module defs_psp.f) | filpsp(IN)=name of formatted external file containing atomic psp data. | 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 nl form factors (array ffspl) | mqgrid_vl(IN)=dimension of q (or G) grid for Vloc (array vlspl) | n1xccc(IN)=dimension of xccc1d ; 0 if no XC core correction is used | optnlxccc(IN)=option for nl XC core correction | positron(IN)=0 if electron GS calculation | 1 if positron GS calculation | 2 if electron GS calculation in presence of the positron | pspso(INOUT)=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 | if =1 : this input will be changed at output to 1, 2, 3, according | to the intrinsic characteristics of the psp file | 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 | usepaw(IN)= 0 for non paw calculation; =1 for paw calculation | 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. | znuclpsp(IN)=atomic number of atom as specified in input file to main routine
NOTES
Format expected for the three first lines of pseudopotentials (1) title (character) line (2) znucl,zion,pspdat (3) pspcod,pspxc,lmax,lloc,mmax,r2well Dimensions of form factors and Vloc q grids must be the same in Norm-Conserving case
SOURCE
786 subroutine pspatm(dq,dtset,dtfil,ekb,epsatm,ffspl,indlmn,ipsp,pawrad,pawtab,& 787 psps,vlspl,dvlspl,xcccrc,xccc1d,nctab,comm_mpi) 788 789 !Arguments --------------------------------------------- 790 !scalars 791 integer,intent(in) :: ipsp 792 integer, optional,intent(in) :: comm_mpi 793 real(dp),intent(in) :: dq 794 real(dp),intent(out) :: epsatm,xcccrc 795 type(dataset_type),intent(in) :: dtset 796 type(datafiles_type),intent(in) :: dtfil 797 type(pawrad_type),intent(inout) :: pawrad 798 type(pawtab_type),intent(inout) :: pawtab 799 type(nctab_t),intent(inout) :: nctab 800 type(pseudopotential_type),intent(inout) :: psps 801 !arrays 802 integer,intent(out) :: indlmn(6,psps%lmnmax) 803 real(dp),intent(out) :: dvlspl(psps%mqgrid_vl,2) 804 real(dp),intent(inout) :: ekb(psps%dimekb*(1-psps%usepaw)) 805 real(dp),intent(inout) :: ffspl(psps%mqgrid_ff,2,psps%lnmax) 806 real(dp),intent(out) :: vlspl(psps%mqgrid_vl,2) 807 real(dp),intent(inout) :: xccc1d(psps%n1xccc*(1-psps%usepaw),6) 808 809 !Local variables --------------------------------------- 810 !scalars 811 integer :: ii,il,ilmn,iln,iln0,lloc,lmax,me,mmax 812 integer :: paral_mode,pspcod,pspdat,pspxc,useupf,usexml,xmlpaw,unt 813 real(dp) :: maxrad,qchrg,r2well,zion,znucl 814 logical,parameter :: nc_debug = .False. 815 !logical,parameter :: nc_debug = .True. 816 character(len=500) :: msg,errmsg 817 character(len=fnlen) :: title, filnam 818 type(pawpsp_header_type):: pawpsp_header 819 type(pspheader_type) :: pspheads_tmp 820 !arrays 821 integer,allocatable :: nproj(:) 822 real(dp) :: tsec(2),ecut_tmp(3,2) 823 real(dp),allocatable :: e990(:),e999(:),ekb1(:),ekb2(:),epspsp(:),rcpsp(:) 824 real(dp),allocatable :: rms(:) 825 #if defined HAVE_LIBPSML 826 !! usexml= 0 for non xml ps format ; =1 for xml ps format 827 character(len=3) :: atmsymb 828 character(len=30) :: creator 829 type(pspheader_type) :: psphead 830 #endif 831 832 ! ****************************************************************************** 833 834 !paral_mode defines how we access to the psp file 835 ! paral_mode=0: all processes access to the file (sequentially) 836 ! paral_mode=1: only proc. 0 access to the file and then broadcast 837 paral_mode=0 838 if (present(comm_mpi)) then 839 if (psps%usepaw==1.and.xmpi_comm_size(comm_mpi)>1) paral_mode=1 840 end if 841 me=0;if (paral_mode==1) me=xmpi_comm_rank(comm_mpi) 842 843 if (paral_mode == 1) then 844 ABI_CHECK(psps%usepaw==1, "paral_mode==1 is only compatible with PAW, see call to pawpsp_bcast below") 845 end if 846 847 nctab%has_tvale = .False.; nctab%has_tcore = .False. 848 pspcod = -1 849 850 if (me==0) then 851 ! Dimensions of form factors and Vloc q grids must be the same in Norm-Conserving case 852 if (psps%usepaw==0 .and. psps%mqgrid_ff/=psps%mqgrid_vl) then 853 write(msg, '(a,a,a,a,a)' )& 854 'Dimension of q-grid for nl form factors (mqgrid_ff)',ch10,& 855 'is different from dimension of q-grid for Vloc (mqgrid_vl) !',ch10,& 856 'This is not allowed for norm-conserving psp.' 857 ABI_ERROR(msg) 858 end if 859 860 write(msg, '(a,t38,a)' )'- pspatm: opening atomic psp file',trim(psps%filpsp(ipsp)) 861 call wrtout([std_out, ab_out], msg) 862 863 ! Check if the file pseudopotential file is written in (XML| XML-PAW | UPF) 864 call test_xml_xmlpaw_upf(psps%filpsp(ipsp), usexml, xmlpaw, useupf) 865 866 ! ---------------------------------------------------------------------------- 867 ! allocate nproj here: can be read in now for UPF 868 ABI_MALLOC(nproj,(psps%mpssoang)) 869 nproj(:)=0 870 871 if (usexml /= 1 .and. useupf == 0) then 872 873 ! Open the atomic data file, and read the three first lines 874 ! These three first lines have a similar format in all allowed psp files 875 ! Open atomic data file (note: formatted input file) 876 if (open_file(psps%filpsp(ipsp), msg, unit=tmp_unit, form='formatted', status='old') /= 0) then 877 ABI_ERROR(msg) 878 end if 879 rewind (unit=tmp_unit,err=10,iomsg=errmsg) 880 881 ! Read and write some description of file from first line (character data) 882 read (tmp_unit,'(a)',err=10,iomsg=errmsg) title 883 write(msg, '(a,a)' ) '- ',trim(title) 884 call wrtout([std_out, ab_out], msg) 885 886 ! Read and write more data describing psp parameters 887 read (tmp_unit,*,err=10,iomsg=errmsg) znucl,zion,pspdat 888 write(msg,'(a,f9.5,f10.5,2x,i8,t47,a)')'-',znucl,zion,pspdat,'znucl, zion, pspdat' 889 call wrtout([std_out, ab_out], msg) 890 891 read (tmp_unit,*,err=10,iomsg=errmsg) pspcod,pspxc,lmax,lloc,mmax,r2well 892 if(pspxc<0) then 893 write(msg, '(i5,i8,2i5,i10,f10.5,t47,a)' ) & 894 pspcod,pspxc,lmax,lloc,mmax,r2well,'pspcod,pspxc,lmax,lloc,mmax,r2well' 895 else 896 write(msg, '(4i5,i10,f10.5,t47,a)' ) & 897 pspcod,pspxc,lmax,lloc,mmax,r2well,'pspcod,pspxc,lmax,lloc,mmax,r2well' 898 end if 899 call wrtout([std_out, ab_out], msg) 900 901 else if (usexml == 1 .and. xmlpaw == 0) then 902 903 ! the following is probably useless - already read in everything in inpspheads 904 #if defined HAVE_LIBPSML 905 call psxml2abheader( psps%filpsp(ipsp), psphead, atmsymb, creator, 0 ) 906 znucl = psphead%znuclpsp 907 zion = psphead%zionpsp 908 pspdat = psphead%pspdat 909 pspcod = psphead%pspcod 910 pspxc = psphead%pspxc 911 lmax = psphead%lmax 912 !lloc = 0 ! does this mean s? in psml case the local potential can be different from any l channel 913 lloc = -1 914 mmax = -1 915 r2well = 0 916 917 write(msg,'(a,1x,a3,3x,a)') "-",atmsymb,trim(creator) 918 call wrtout([std_out, ab_out], msg) 919 write(msg,'(a,f9.5,f10.5,2x,i8,t47,a)')'-',znucl,zion,pspdat,'znucl, zion, pspdat' 920 call wrtout([std_out, ab_out], msg) 921 if(pspxc<0) then 922 write(msg, '(i5,i8,2i5,i10,f10.5,t47,a)' ) & 923 pspcod,pspxc,lmax,lloc,mmax,r2well,'pspcod,pspxc,lmax,lloc,mmax,r2well' 924 else 925 write(msg, '(4i5,i10,f10.5,t47,a)' ) & 926 pspcod,pspxc,lmax,lloc,mmax,r2well,'pspcod,pspxc,lmax,lloc,mmax,r2well' 927 end if 928 call wrtout([std_out, ab_out], msg) 929 #else 930 write(msg,'(a,a)') & 931 'ABINIT is not compiled with XML support for reading this type of pseudopotential ', trim(psps%filpsp(ipsp)) 932 ABI_BUG(msg) 933 #endif 934 ! END useless 935 else if (usexml == 1 .and. xmlpaw == 1) then 936 write(msg,'(a,a)') & 937 '- pspatm : Reading pseudopotential header in XML form from ', trim(psps%filpsp(ipsp)) 938 call wrtout([std_out, ab_out], msg) 939 940 ! Return header information 941 call pawpsxml2ab(psps%filpsp(ipsp),ecut_tmp, pspheads_tmp,0) 942 lmax=pspheads_tmp%lmax 943 pspxc=pspheads_tmp%pspxc 944 znucl=pspheads_tmp%znuclpsp 945 pawpsp_header%basis_size=pspheads_tmp%pawheader%basis_size 946 pawpsp_header%l_size=pspheads_tmp%pawheader%l_size 947 pawpsp_header%lmn_size=pspheads_tmp%pawheader%lmn_size 948 pawpsp_header%mesh_size=pspheads_tmp%pawheader%mesh_size 949 pawpsp_header%pawver=pspheads_tmp%pawheader%pawver 950 pawpsp_header%shape_type=pspheads_tmp%pawheader%shape_type 951 pawpsp_header%rpaw=pspheads_tmp%pawheader%rpaw 952 pawpsp_header%rshp=pspheads_tmp%pawheader%rshp 953 lloc=0; pspcod=17 954 955 else if (useupf /= 0) then 956 if (psps%usepaw /= 0) then 957 ABI_ERROR("UPF format not allowed with PAW (USPP part not read yet)") 958 end if 959 960 r2well = 0; qchrg=zero 961 962 ! should initialize znucl,zion,pspxc,lmax,lloc,mmax 963 if (useupf == 1) then 964 pspcod = 11 965 call upf1_to_abinit(psps%filpsp(ipsp), znucl, zion, pspxc, lmax, lloc, mmax, & 966 psps, epsatm, xcccrc, indlmn, ekb, ffspl, nproj, vlspl, xccc1d) 967 else 968 pspcod = 12 969 call upf2_to_abinit(ipsp, psps%filpsp(ipsp), dtset%vloc_rcut, znucl, zion, pspxc, lmax, lloc, mmax, & 970 psps, epsatm, xcccrc, indlmn, ekb, ffspl, nproj, vlspl, xccc1d, nctab, maxrad) 971 972 if (nc_debug) then 973 call psp_dump_outputs("UPF2", pspcod, psps%lmnmax, psps%lnmax, psps%mpssoang, & 974 psps%mqgrid_ff, psps%n1xccc, mmax, maxrad, epsatm, qchrg, xcccrc, nctab, & 975 indlmn, nproj, ekb, ffspl, vlspl, xccc1d) 976 end if 977 end if 978 979 else 980 ABI_ERROR("You should not be here! erroneous type or pseudopotential input") 981 end if 982 983 ! ------------------------------------------------------------------------------ 984 ! Check data for consistency against main routine input 985 986 ! Does required spin-orbit characteristics agree with format 987 ! TODO: in case of pspcod 5 (phoney) and 8 (oncvpsp) this is not specific enough. 988 ! they can be non-SOC as well. 989 ! HGH is ok - can always turn SOC on or off. 990 ! PAW is ok - can be used with or without SOC 991 ! write(std_out,*) pspso 992 if((pspcod/=3).and.(pspcod/=5).and.(pspcod/=8).and.(pspcod/=10).and. pspcod /= 12 .and. & 993 (pspcod/=7).and.(pspcod/=17))then 994 ! If pspso requires internal characteristics, set it to 1 for non-HGH psps 995 if(psps%pspso(ipsp)==1) psps%pspso(ipsp)=0 996 if(psps%pspso(ipsp) /= 0)then 997 write(msg, '(3a,i0,3a)' )& 998 'Pseudopotential file cannot give spin-orbit characteristics,',ch10,& 999 'while pspso(itypat)= ',psps%pspso(ipsp),'.',ch10,& 1000 'Action: check your pseudopotential and input files for consistency.' 1001 ABI_ERROR(msg) 1002 end if 1003 end if 1004 1005 ! Does nuclear charge znuclpsp agree with psp input znucl 1006 if (abs(psps%znuclpsp(ipsp)-znucl) > tol8) then 1007 write(msg, '(a,f10.5,2a,f10.5,5a)' )& 1008 'Pseudopotential file znucl: ',znucl,ch10,& 1009 'does not equal input znuclpsp: ',psps%znuclpsp(ipsp),' better than 1e-08 .',ch10,& 1010 'znucl is read from the psp file in pspatm, while',ch10,& 1011 'znuclpsp is read in iofn2.' 1012 ABI_BUG(msg) 1013 end if 1014 1015 ! Is the highest angular momentum within limits? 1016 ! Recall mpsang is 1+highest l for nonlocal correction. 1017 ! Nonlocal corrections for s, p, d, and f are supported. 1018 if (lmax+1>psps%mpsang) then 1019 write(msg, '(a,i0,a,i0,a,a)' )& 1020 'input lmax+1: ',lmax+1,' exceeds mpsang: ',psps%mpsang,ch10,& 1021 'indicates input lmax too large for dimensions.' 1022 ABI_BUG(msg) 1023 end if 1024 1025 ! Check several choices for ixc against pspxc 1026 ! ixc is from ABINIT code; pspxc is from atomic psp file 1027 if (dtset%ixc==0) then 1028 ABI_WARNING('Note that input ixc=0 => no xc is being used.') 1029 else if(dtset%ixc/=pspxc) then 1030 write(msg, '(a,i0,a,i0,8a)' )& 1031 'Pseudopotential file pspxc: ',pspxc,', not equal to input ixc: ',dtset%ixc,'.',ch10,& 1032 'These parameters must agree to get the same xc in ABINIT code as in psp construction.',ch10,& 1033 'Action: check psp design or input file.',ch10,& 1034 'Assume experienced user. Execution will continue.' 1035 ABI_WARNING(msg) 1036 end if 1037 1038 if (lloc>lmax .and. pspcod/=4 .and. pspcod/=8 .and. pspcod/=10) then 1039 write(msg, '(a,2i12,a,a,a,a)' )& 1040 'lloc,lmax=',lloc,lmax,ch10,& 1041 'chosen l of local psp exceeds range from input data.',ch10,& 1042 'Action: check pseudopotential input file.' 1043 ABI_ERROR(msg) 1044 end if 1045 1046 ! Does the pspcod agree with type of calculation (paw or not)? 1047 if (((pspcod/=7.and.pspcod/=17).and.psps%usepaw==1).or.((pspcod==7.or.pspcod==17).and.psps%usepaw==0)) then 1048 write(msg, '(a,i0,a,a,i0,a)' )& 1049 'In reading atomic psp file, finds pspcod= ',pspcod,ch10,& 1050 'This is not an allowed value with usepaw= ',psps%usepaw,'.' 1051 ABI_BUG(msg) 1052 end if 1053 1054 if (.not.psps%vlspl_recipSpace .and. (pspcod /= 2 .and. pspcod /= 3 .and. pspcod /= 10 .and. pspcod /= 7)) then 1055 ! The following "if" statement can substitute the one just before once libBigDFT 1056 ! has been upgraded to include pspcod 10 1057 ! if (.not.psps%vlspl_recipSpace .and. (pspcod /= 2 .and. pspcod /= 3 .and. pspcod /= 10)) then 1058 write(msg, '(a,i0,2a)' )& 1059 'In reading atomic psp file, finds pspcod=',pspcod,ch10,& 1060 'This is not an allowed value with real space computation.' 1061 ABI_BUG(msg) 1062 end if 1063 1064 ! MJV 16/6/2009 added pspcod 11 for upf format 1065 !if( pspcod<1 .or. (pspcod>11.and.pspcod/=17) ) then 1066 if( pspcod<1) then 1067 write(msg, '(a,i0,4a)' )& 1068 'In reading atomic psp file, finds pspcod= ',pspcod,ch10,& 1069 'This is not an allowed value. Allowed values are 1-12 or 17 .',ch10,& 1070 'Action: check pseudopotential input file.' 1071 ABI_ERROR(msg) 1072 end if 1073 1074 ! ----------------------------------------------------------------------- 1075 ! Set various terms to 0 in case not defined below 1076 ABI_MALLOC(e990,(psps%mpssoang)) 1077 ABI_MALLOC(e999,(psps%mpssoang)) 1078 ABI_MALLOC(rcpsp,(psps%mpssoang)) 1079 ABI_MALLOC(rms,(psps%mpssoang)) 1080 ABI_MALLOC(epspsp,(psps%mpssoang)) 1081 ABI_MALLOC(ekb1,(psps%mpssoang)) 1082 ABI_MALLOC(ekb2,(psps%mpssoang)) 1083 e990(:)=zero ;e999(:)=zero 1084 rcpsp(:)=zero;rms(:)=zero 1085 ekb1(:)=zero ;ekb2(:)=zero 1086 epspsp(:)=zero 1087 qchrg=zero 1088 1089 ! ---------------------------------------------------------------------- 1090 if (pspcod==1 .or. pspcod==4)then 1091 1092 ! Teter pseudopotential (pspcod=1 or 4) 1093 call psp1in(dq,ekb,ekb1,ekb2,epsatm,epspsp,& 1094 e990,e999,ffspl,indlmn,lloc,lmax,psps%lmnmax,psps%lnmax,& 1095 mmax,psps%mpsang,psps%mqgrid_ff,nproj,psps%n1xccc,pspcod,qchrg,psps%qgrid_ff,& 1096 rcpsp,rms,psps%useylm,vlspl,xcccrc,xccc1d,zion,psps%znuclpsp(ipsp)) 1097 1098 else if (pspcod==2)then 1099 1100 ! GTH pseudopotential 1101 call psp2in(dtset,ekb,epsatm,ffspl,indlmn,ipsp,lmax,nproj,psps,vlspl,dvlspl,zion) 1102 xccc1d(:,:)=0.0d0 ; qchrg=0.0d0 ; xcccrc=0.0d0 1103 1104 else if (pspcod==3)then 1105 1106 ! HGH pseudopotential 1107 call psp3in(dtset,ekb,epsatm,ffspl,indlmn,ipsp,lmax,nproj,psps, psps%pspso(ipsp), vlspl,zion) 1108 xccc1d(:,:)=0.0d0 ; qchrg=0.0d0 ; xcccrc=0.0d0 1109 1110 else if (pspcod==5)then 1111 1112 ! Old phoney pseudopotentials 1113 call psp5in(ekb,ekb1,ekb2,epsatm,epspsp,& 1114 e990,e999,ffspl,indlmn,lloc,lmax,psps%lmnmax,psps%lnmax,& 1115 mmax,psps%mpsang,psps%mpssoang,psps%mqgrid_ff,nproj,psps%n1xccc,psps%pspso(ipsp),qchrg,psps%qgrid_ff,& 1116 rcpsp,rms,psps%useylm,vlspl,xcccrc,xccc1d,zion,psps%znuclpsp(ipsp)) 1117 1118 else if (pspcod==6)then 1119 ! FHI pseudopotentials 1120 call psp6in(ekb,epsatm,ffspl,indlmn,lloc,lmax,psps%lmnmax,psps%lnmax,mmax,& 1121 psps%mpsang,psps%mqgrid_ff,nproj,psps%n1xccc,psps%optnlxccc,psps%positron,qchrg,psps%qgrid_ff,psps%useylm,vlspl,& 1122 xcccrc,xccc1d,zion,psps%znuclpsp(ipsp)) 1123 1124 else if (pspcod==7)then 1125 ! PAW "pseudopotentials" 1126 call pawpsp_7in(epsatm,ffspl,dtset%icoulomb,ABS(dtset%hyb_mixing),dtset%ixc,& 1127 & lmax,psps%lnmax,mmax,psps%mqgrid_ff,psps%mqgrid_vl,& 1128 & pawrad,pawtab,dtset%pawxcdev,psps%qgrid_ff,psps%qgrid_vl,& 1129 & dtset%usewvl,dtset%usexcnhat_orig,vlspl,xcccrc,dtset%xclevel,& 1130 & dtset%xc_denpos,zion,psps%znuclpsp(ipsp),xc_taupos=dtset%xc_taupos) 1131 1132 else if (pspcod==8)then 1133 1134 ! DRH pseudopotentials 1135 call psp8in(ekb,epsatm,ffspl,indlmn,lloc,lmax,psps%lmnmax,psps%lnmax,mmax,& 1136 psps%mpsang,psps%mpssoang,psps%mqgrid_ff,psps%mqgrid_vl,nproj,psps%n1xccc,psps%pspso(ipsp),& 1137 qchrg,psps%qgrid_ff,psps%qgrid_vl,psps%useylm,vlspl,xcccrc,xccc1d,zion,psps%znuclpsp(ipsp),nctab,maxrad) 1138 1139 if (nc_debug) then 1140 call psp_dump_outputs("PSP8",pspcod,psps%lmnmax,psps%lnmax,psps%mpssoang, & 1141 psps%mqgrid_ff,psps%n1xccc,mmax,maxrad,epsatm,qchrg,xcccrc,nctab, & 1142 indlmn,nproj,ekb,ffspl,vlspl,xccc1d) 1143 end if 1144 1145 else if (pspcod==9)then 1146 1147 #if defined HAVE_LIBPSML 1148 call psp9in(psps%filpsp(ipsp),ekb,epsatm,ffspl,indlmn,lloc,lmax,psps%lmnmax,psps%lnmax,mmax,& 1149 psps%mpsang,psps%mpssoang,psps%mqgrid_ff,psps%mqgrid_vl,nproj,psps%n1xccc, & 1150 psps%pspso(ipsp),qchrg,psps%qgrid_ff,psps%qgrid_vl,psps%useylm,vlspl,& 1151 xcccrc,xccc1d,zion,psps%znuclpsp(ipsp),nctab,maxrad) 1152 1153 if (nc_debug) then 1154 call psp_dump_outputs("PSML",pspcod,psps%lmnmax,psps%lnmax,psps%mpssoang, & 1155 psps%mqgrid_ff,psps%n1xccc,mmax,maxrad,epsatm,qchrg,xcccrc,nctab, & 1156 indlmn,nproj,ekb,ffspl,vlspl,xccc1d) 1157 end if 1158 #else 1159 write(msg,'(2a)') & 1160 'ABINIT is not compiled with XML support for reading this type of pseudopotential ', trim(psps%filpsp(ipsp)) 1161 ABI_ERROR(msg) 1162 #endif 1163 1164 else if (pspcod==10)then 1165 1166 ! HGH pseudopotential, full h/k matrix read 1167 call psp10in(dtset,ekb,epsatm,ffspl,indlmn,ipsp,lmax,nproj,psps, psps%pspso(ipsp), vlspl,zion) 1168 xccc1d(:,:)=0.0d0 ; qchrg=0.0d0 ; xcccrc=0.0d0 1169 1170 ! NB for pspcod 11 the reading has already been done above. 1171 else if (pspcod==17)then 1172 ! PAW XML pseudopotentials 1173 call pawpsp_17in(epsatm,ffspl,dtset%icoulomb,ipsp,ABS(dtset%hyb_mixing),dtset%ixc,lmax,& 1174 & psps%lnmax,mmax,psps%mqgrid_ff,psps%mqgrid_vl,pawpsp_header,pawrad,pawtab,& 1175 & dtset%pawxcdev,psps%qgrid_ff,psps%qgrid_vl,dtset%usewvl,& 1176 & dtset%usexcnhat_orig,vlspl,xcccrc,& 1177 & dtset%xclevel,dtset%xc_denpos,pspheads_tmp%zionpsp,psps%znuclpsp(ipsp),& 1178 & xc_taupos=dtset%xc_taupos) 1179 call paw_setup_free(paw_setuploc) 1180 end if 1181 1182 close (unit=tmp_unit) 1183 1184 ! ---------------------------------------------------------------------- 1185 if (pspcod==2 .or. pspcod==3 .or. pspcod==10)then 1186 write(msg, '(a,a,a,a,a,a,a,a,a,a)' )ch10,& 1187 ' pspatm : COMMENT -',ch10,& 1188 ' the projectors are not normalized,',ch10,& 1189 ' so that the KB energies are not consistent with ',ch10,& 1190 ' definition in PRB44, 8503 (1991). ',ch10,& ! [[cite:Gonze1991]] 1191 ' However, this does not influence the results obtained hereafter.' 1192 call wrtout([std_out, ab_out], msg) 1193 ! The following lines are added to keep backward compatibilty 1194 maxrad=zero 1195 #if defined HAVE_BIGDFT 1196 do ii=1,size(psps%gth_params%psppar,1)-1 ! psppar first dim begins at 0 1197 if (psps%gth_params%psppar(ii,0,ipsp)/=zero) maxrad=max(maxrad,psps%gth_params%psppar(ii,0,ipsp)) 1198 end do 1199 if (abs(maxrad)<=tol12) then 1200 psps%gth_params%radii_cf(ipsp,3)=zero 1201 else 1202 psps%gth_params%radii_cf(ipsp,3)=max( & 1203 & min(dtset%wvl_crmult*psps%gth_params%radii_cf(ipsp,1),15._dp*maxrad)/dtset%wvl_frmult, & 1204 & psps%gth_params%radii_cf(ipsp,2)) 1205 end if 1206 #endif 1207 end if 1208 1209 if (pspcod/=7.and.pspcod/=17) then 1210 write(msg, '(a,f14.8,a,a)' ) ' pspatm : epsatm=',epsatm,ch10,' --- l ekb(1:nproj) -->' 1211 call wrtout([std_out, ab_out], msg) 1212 !print *, "nproj", nproj 1213 iln0=0 1214 do ilmn=1,psps%lmnmax 1215 iln=indlmn(5,ilmn) 1216 if (iln>iln0) then 1217 il=indlmn(1,ilmn) 1218 if (indlmn(6,ilmn)==1) then 1219 iln0=iln0+nproj(il+1) 1220 write(msg, '(13x,i1,4f12.6)' ) il,(ekb(iln+ii),ii=0,nproj(il+1)-1) 1221 else 1222 ! Note il = ll here i.e. the s channel in the SOC part is not included in nproj. 1223 iln0=iln0+nproj(il+psps%mpsang) 1224 if (dtset%spnorbscl /= one) then 1225 call wrtout([std_out, ab_out], & 1226 sjoin(" Rescaling spin-orbit KB energies using spnorbscl.", ftoa(dtset%spnorbscl))) 1227 ekb(iln:iln + nproj(il + psps%mpsang) - 1) = ekb(iln:iln + nproj(il + psps%mpsang) - 1) * dtset%spnorbscl 1228 end if 1229 write(msg, '(2x,a,i1,4f12.6)' ) 'spin-orbit ',il,(ekb(iln+ii),ii=0,nproj(il+psps%mpsang)-1) 1230 end if 1231 call wrtout([std_out, ab_out], msg) 1232 end if 1233 end do 1234 end if 1235 1236 ! NC: Evalute spline-fit of the model core charge in reciprocal space. 1237 ! TODO: Be careful, because we will be using the PAW routines in which tcore is always available! 1238 ! Should add a test with 2 NC pseudos: one with NLCC and the other without! 1239 if (psps%usepaw == 0) then 1240 call nctab_eval_tcorespl(nctab, psps%n1xccc, xcccrc, xccc1d, psps%mqgrid_vl, psps%qgrid_vl) 1241 end if 1242 1243 write(msg,'(3a)') ' pspatm: atomic psp has been read ',' and splines computed',ch10 1244 call wrtout([std_out, ab_out], msg) 1245 1246 ABI_FREE(e990) 1247 ABI_FREE(e999) 1248 ABI_FREE(rcpsp) 1249 ABI_FREE(rms) 1250 ABI_FREE(ekb1) 1251 ABI_FREE(ekb2) 1252 ABI_FREE(epspsp) 1253 ABI_FREE(nproj) 1254 1255 if (dtset%prtvol > 9 .and. psps%usepaw==0 .and. psps%lmnmax>3) then 1256 write (filnam, '(a,i0,a)') trim(dtfil%fnameabo_pspdata), ipsp, ".dat" 1257 if (open_file(filnam, msg, newunit=unt) /= 0) then 1258 ABI_ERROR(msg) 1259 end if 1260 write (unt,*) '# Pseudopotential data in reciprocal space as used by ABINIT' 1261 write (unt,'(a)', ADVANCE='NO') '# index vlocal ' 1262 if (psps%lnmax > 0) write (unt,'(a,I3)', ADVANCE='NO') ' 1st proj(l=', indlmn(1,1) 1263 if (psps%lnmax > 1) write (unt,'(a,I3)', ADVANCE='NO') ') 2nd(l=', indlmn(1,2) 1264 if (psps%lnmax > 2) write (unt,'(a,I3,a)', ADVANCE='NO') ') 3rd(l=', indlmn(1,3), ')' 1265 write (unt,*) 1266 1267 do ii = 1, psps%mqgrid_vl 1268 write(unt, '(I5,E24.16)', ADVANCE='NO') ii, vlspl(ii,1) 1269 if (psps%lnmax > 0) write(unt, '(E24.16)', ADVANCE='NO') ffspl(ii,1,1) 1270 if (psps%lnmax > 1) write(unt, '(E24.16)', ADVANCE='NO') ffspl(ii,1,2) 1271 if (psps%lnmax > 2) write(unt, '(E24.16)', ADVANCE='NO') ffspl(ii,1,3) 1272 write(unt, *) 1273 end do 1274 close(unt) 1275 1276 write (filnam, '(a,i0,a)') trim(dtfil%fnameabo_nlcc_derivs), ipsp, ".dat" 1277 if (open_file(filnam, msg, newunit=unt) /= 0) then 1278 ABI_ERROR(msg) 1279 end if 1280 write (unt,*) '# Non-linear core corrections' 1281 write (unt,*) '# r, pseudocharge, 1st, 2nd, 3rd, 4th, 5th derivatives' 1282 do ii = 1, psps%n1xccc 1283 write (unt,*) xcccrc*(ii-1)/(psps%n1xccc-1), xccc1d(ii,1), xccc1d(ii,2), & 1284 xccc1d(ii,3), xccc1d(ii,4), xccc1d(ii,5), xccc1d(ii,6) 1285 end do 1286 close(unt) 1287 end if 1288 1289 end if ! me=0 1290 1291 if (paral_mode==1) then 1292 call timab(48,1,tsec) 1293 call pawpsp_bcast(comm_mpi,epsatm,ffspl,pawrad,pawtab,vlspl,xcccrc) 1294 call timab(48,2,tsec) 1295 end if 1296 1297 if (psps%usepaw==1) then 1298 indlmn(:,:)=0 1299 indlmn(1:6,1:pawtab%lmn_size)=pawtab%indlmn(1:6,1:pawtab%lmn_size) 1300 end if 1301 1302 !-------------------------------------------------------------------- 1303 !WVL+PAW: 1304 if (dtset%usepaw==1 .and. (dtset%icoulomb /= 0 .or. dtset%usewvl==1)) then 1305 #if defined HAVE_BIGDFT 1306 psps%gth_params%psppar(:,:,ipsp) = UNINITIALIZED(1._dp) 1307 psps%gth_params%radii_cf(ipsp,:) = UNINITIALIZED(1._dp) 1308 call wvl_descr_psp_fill(psps%gth_params, ipsp, psps%pspxc(1), int(psps%zionpsp(ipsp)), int(psps%znuclpsp(ipsp)), 0) 1309 #endif 1310 1311 ! The following lines are added to keep backward compatibilty 1312 maxrad=zero 1313 #if defined HAVE_BIGDFT 1314 do ii=1,size(psps%gth_params%psppar,1)-1 ! psppar first dim begins at 0 1315 if (psps%gth_params%psppar(ii,0,ipsp)/=zero) maxrad=max(maxrad,psps%gth_params%psppar(ii,0,ipsp)) 1316 end do 1317 if (abs(maxrad)<=tol12) then 1318 !== MT COMMENT 1319 ! Damien wants to activate this (in order to directly compare to bigDFT): 1320 psps%gth_params%radii_cf(ipsp,3)= psps%gth_params%radii_cf(ipsp,2) 1321 ! But, this changes strongly file references. 1322 ! So, I keep this, waiting for Tonatiuh s validation 1323 psps%gth_params%radii_cf(ipsp,3) = (psps%gth_params%radii_cf(ipsp,1)+psps%gth_params%radii_cf(ipsp,2))*half 1324 !== MT COMMENT 1325 else 1326 psps%gth_params%radii_cf(ipsp,3)=max( & 1327 min(dtset%wvl_crmult*psps%gth_params%radii_cf(ipsp,1),15._dp*maxrad)/dtset%wvl_frmult, & 1328 psps%gth_params%radii_cf(ipsp,2)) 1329 end if 1330 if(present(comm_mpi)) then 1331 call pawpsp_wvl(psps%filpsp(ipsp),pawrad,pawtab,dtset%usewvl,dtset%wvl_ngauss,comm_mpi) 1332 else 1333 call pawpsp_wvl(psps%filpsp(ipsp),pawrad,pawtab,dtset%usewvl,dtset%wvl_ngauss) 1334 end if 1335 #endif 1336 end if 1337 1338 !end of WVL+PAW section 1339 !---------------------------------------------------- 1340 1341 return 1342 1343 ! Handle IO error 1344 10 continue 1345 ABI_ERROR(errmsg) 1346 1347 end subroutine pspatm
ABINIT/pspcor [ Functions ]
NAME
pspcor
FUNCTION
Compute ecore pseudoion-pseudoion correction energy from epsatm for different types of atoms in unit cell.
INPUTS
natom=number of atoms in cell ntypat=number of types of atoms typat(natom)=integer label of 'typat' for each atom in cell epsatm(ntypat)=pseudoatom energy for each type of atom zion(ntypat)=valence charge on each type of atom in cell
OUTPUT
ecore=resulting psion-psion energy in Hartrees
SOURCE
656 subroutine pspcor(ecore,epsatm,natom,ntypat,typat,zion) 657 658 !Arguments ------------------------------------ 659 !scalars 660 integer,intent(in) :: natom,ntypat 661 real(dp),intent(out) :: ecore 662 !arrays 663 integer,intent(in) :: typat(natom) 664 real(dp),intent(in) :: epsatm(ntypat),zion(ntypat) 665 666 !Local variables------------------------------- 667 !scalars 668 integer :: ia 669 real(dp) :: charge,esum 670 671 ! ************************************************************************* 672 673 charge = 0.d0 674 esum = 0.d0 675 do ia=1,natom 676 ! compute pseudocharge: 677 charge=charge+zion(typat(ia)) 678 ! add pseudocore energies together: 679 esum = esum + epsatm(typat(ia)) 680 end do 681 682 ecore=charge*esum 683 684 end subroutine pspcor
ABINIT/pspini [ Functions ]
NAME
pspini
FUNCTION
Looping over atom types 1 ... ntypat, read pseudopotential data filename, then call pspatm for each psp. Might combine the psps to generate pseudoatoms, thanks to alchemy. Also compute ecore=[Sum(i) zion(i)] * [Sum(i) epsatm(i)] by calling pspcor.
COPYRIGHT
Copyright (C) 1998-2024 ABINIT group (DCA, XG, GMR, MT) 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
dtset <type(dataset_type)>=all input variables in this dataset | iscf=parameter controlling scf or non-scf calculations | ixc=exchange-correlation choice as input to main routine | natom=number of atoms in unit cell | pawxcdev=choice of XC development in PAW formalism | prtvol= control output volume | typat(natom)=type (integer) for each atom | main routine, for each type of atom gsqcut=cutoff for G^2 based on ecut for basis sphere (bohr^-2) gsqcutdg=PAW only - cutoff for G^2 based on ecutdg (fine grid) for basis sphere (bohr^-2) rprimd(3,3)=dimensional primitive translations in real space (bohr) used to estimate real space mesh (if necessary)
OUTPUT
ecore=total psp core correction energy*ucvol (hartree*bohr^3) pawrad(ntypat*usepaw) <type(pawrad_type)>=paw radial mesh and related data pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data gencond=general condition for new computation of pseudopotentials (if gencond=1, new psps have been re-computed)
SIDE EFFECTS
psps <type(pseudopotential_type)>=at output, psps is completely initialized At the input, it is already partially or completely initialized.
NOTES
The interplay with the multi-dataset mode is interesting: the pseudopotentials are independent of the dataset, but the largest q vector, the spin-orbit characteristics, the use of Ylm as well as ixc play a role in the set up of pseudopotentials (ixc plays a very minor role, however). So, the pseudopotential data ought not be recomputed when gsqcut, gsqcutdg, mqgrid_ff, mqgrid_vl, npspso, ixc, dimekb and useylm do not change. In many cases, this routine is also called just to write the psp line of the header, without reading again the psp. This psp line is constant throughout run.
SOURCE
124 subroutine pspini(dtset,dtfil,ecore,gencond,gsqcut,gsqcutdg,pawrad,pawtab,psps,rprimd,comm_mpi) 125 126 !Arguments ------------------------------------ 127 !scalars 128 integer, optional,intent(in) :: comm_mpi 129 integer,intent(out) :: gencond 130 real(dp),intent(in) :: gsqcut,gsqcutdg 131 real(dp),intent(out) :: ecore 132 type(dataset_type),intent(in) :: dtset 133 type(datafiles_type),intent(in) :: dtfil 134 !arrays 135 real(dp),intent(in) :: rprimd(3,3) 136 type(pseudopotential_type), target,intent(inout) :: psps 137 type(pawrad_type), intent(inout) :: pawrad(psps%ntypat*psps%usepaw) 138 type(pawtab_type), intent(inout) :: pawtab(psps%ntypat*psps%usepaw) 139 140 !Local variables------------------------------- 141 !scalars 142 integer,parameter :: npspmax=50 143 integer,save :: dimekb_old=0,ifirst=1,ixc_old=-1,lmnmax_old=0,lnmax_old=0 144 integer,save :: mpssoang_old=0,mqgridff_old=0,mqgridvl_old=0,optnlxccc_old=-1 145 integer,save :: paw_size_old=-1,pawxcdev_old=-1,positron_old=-2,usekden_old=-1,usepaw_old=-1 146 integer,save :: usexcnhat_old=-1,usewvl_old=-1,useylm_old=-1 147 integer :: comm_mpi_,ierr,ii,ilang,ilmn,ilmn0,iproj,ipsp,ipspalch 148 integer :: ispin,itypalch,itypat,mtypalch,npsp,npspalch,ntypalch 149 integer :: ntypat,ntyppure,paw_size 150 logical :: has_coretau,has_kij,has_tproj,has_tvale,has_nabla,has_shapefncg,has_vminushalf,has_wvl 151 real(dp),save :: ecore_old=zero,gsqcut_old=zero,gsqcutdg_old=zero, spnorbscl_old=-one,hyb_mixing_old=-999.0_dp 152 real(dp) :: dq,epsatm_psp,qmax,rmax,xcccrc 153 character(len=500) :: msg 154 type(pawrad_type) :: pawrad_dum 155 type(pawtab_type) :: pawtab_dum 156 type(nctab_t) :: nctab_dum 157 type(nctab_t),pointer :: nctab_ptr 158 !arrays 159 integer :: paw_options(10) 160 integer,save :: paw_options_old(10)=(/-1,-1,-1,-1,-1,-1,-1,-1,-1,-1/) 161 integer,save :: pspso_old(npspmax),pspso_zero(npspmax) 162 integer,allocatable :: indlmn_alch(:,:,:),new_pspso(:) 163 integer,pointer :: indlmn(:,:) 164 real(dp),save :: epsatm(npspmax) 165 real(dp) :: tsec(2) 166 real(dp),allocatable :: dvlspl(:,:),dvlspl_alch(:,:,:),ekb(:),ekb_alch(:,:) 167 real(dp),allocatable :: epsatm_alch(:),ffspl(:,:,:),ffspl_alch(:,:,:,:) 168 real(dp),allocatable :: vlspl(:,:),vlspl_alch(:,:,:),xccc1d(:,:) 169 real(dp),allocatable :: xccc1d_alch(:,:,:),xcccrc_alch(:) 170 type(nctab_t),target,allocatable :: nctab_alch(:) 171 172 ! ************************************************************************* 173 174 DBG_ENTER("COLL") 175 176 ! Keep track of time spent in this subroutine 177 call timab(15,1,tsec) 178 179 !------------------------------------------------------------- 180 ! Some initializations 181 !------------------------------------------------------------- 182 183 !Useful sizes 184 ntypat=psps%ntypat 185 mtypalch=psps%mtypalch 186 npsp=psps%npsp 187 if (npsp>npspmax) then 188 ABI_BUG("npsp>npspmax in pspini !") 189 end if 190 191 ! Set up q grids for atomic data represented in reciprocal space 192 ! make qmax 20% larger than largest expected: 193 qmax=1.2d0 * sqrt(gsqcut) 194 !ffnl is always computed in reciprocal space 195 dq=qmax/(one*(psps%mqgrid_ff-1)) 196 do ii=1,psps%mqgrid_ff 197 psps%qgrid_ff(ii)=(ii-1)*dq 198 end do 199 if (psps%usepaw==1) qmax=1.2d0 * sqrt(gsqcutdg) 200 !If vlspl is computed in real space, qgrid contains a real space mesh 201 !the max is taken as the biggest distance in the box. 202 if (psps%vlspl_recipSpace) then 203 dq=qmax/(one*(psps%mqgrid_vl-1)) 204 else 205 rmax = (rprimd(1, 1) + rprimd(1, 2) + rprimd(1, 3)) ** 2 206 rmax = rmax + (rprimd(2, 1) + rprimd(2, 2) + rprimd(2, 3)) ** 2 207 rmax = rmax + (rprimd(3, 1) + rprimd(3, 2) + rprimd(3, 3)) ** 2 208 rmax = sqrt(rmax) 209 dq = rmax /(one*(psps%mqgrid_vl-1)) 210 end if 211 do ii=1,psps%mqgrid_vl 212 psps%qgrid_vl(ii)=(ii-1)*dq 213 end do 214 215 !Determine whether new optional data requests have changed 216 paw_options=0;paw_size=0 217 if (psps%usepaw==1) then 218 paw_size=size(pawtab) 219 has_kij=(dtset%positron/=0.or.abs(dtset%effmass_free-one)>tol8.or.dtset%orbmag>0) 220 has_tvale=.true. ! Will be modified later (depending on PAW dataset format) 221 has_nabla=.false. 222 has_shapefncg=(dtset%optdriver==RUNL_GSTATE.and.((dtset%iprcel>=20.and.dtset%iprcel<70).or.dtset%iprcel>=80)) 223 has_wvl=(dtset%usewvl==1.or.dtset%icoulomb/=0) 224 has_tproj=(dtset%usewvl==1) ! projectors will be free at the end of the psp reading 225 has_vminushalf=(maxval(dtset%ldaminushalf)==1) 226 has_coretau=(dtset%usekden>=1) 227 if (has_kij) paw_options(1)=1 228 if (has_tvale) paw_options(2)=1 229 if (has_nabla) paw_options(5)=1 230 if (has_shapefncg) paw_options(6)=1 231 if (has_wvl) paw_options(7)=1 232 if (has_tproj) paw_options(8)=1 233 if (has_vminushalf)paw_options(9)=1 234 if (has_coretau) paw_options(10)=1 235 !if (dtset%prtvclmb /= 0) then 236 paw_options(3) = 1 237 paw_options(4) = 1 238 !end if 239 end if 240 241 !Determine whether the spin-orbit characteristic has changed 242 !Do not forget that the SO is not consistent with alchemy presently 243 ABI_MALLOC(new_pspso,(npsp)) 244 if (ifirst==1) pspso_old(:)=-1 245 if (ifirst==1) pspso_zero(:)=-1 246 do ipsp=1,npsp 247 new_pspso(ipsp)=1 248 ! No new characteristics if it is equal to the old one, 249 ! or, if it is one, the old one is equal to the intrinsic characteristic one. 250 if (psps%pspso(ipsp)==pspso_old(ipsp).or. & 251 & (psps%pspso(ipsp)==1.and.pspso_old(ipsp)==pspso_zero(ipsp))) then 252 new_pspso(ipsp)=0 253 end if 254 ! No new characteristics if PAW 255 if (psps%usepaw==1) new_pspso(ipsp)=0 256 ! Prepare the saving of the intrinsic pseudopotential characteristics 257 if(psps%pspso(ipsp)==1) pspso_zero(ipsp)=0 258 end do 259 260 !Compute the general condition for new computation of pseudopotentials 261 gencond=0 262 if( ixc_old /= dtset%ixc & 263 & .or. ABS(hyb_mixing_old - dtset%hyb_mixing)>tol8 & 264 & .or. mqgridff_old /= psps%mqgrid_ff & 265 & .or. mqgridvl_old /= psps%mqgrid_vl & 266 & .or. mpssoang_old /= psps%mpssoang & 267 & .or. abs(gsqcut_old-gsqcut)>1.0d-10 & 268 & .or. (psps%usepaw==1.and.abs(gsqcutdg_old-gsqcutdg)>1.0d-10) & 269 & .or. dimekb_old /= psps%dimekb & 270 & .or. lmnmax_old /= psps%lmnmax & 271 & .or. lnmax_old /= psps%lnmax & 272 & .or. optnlxccc_old /= psps%optnlxccc & 273 & .or. usepaw_old /= psps%usepaw & 274 & .or. useylm_old /= psps%useylm & 275 & .or. pawxcdev_old /= dtset%pawxcdev & 276 & .or. (psps%usepaw == 0 .and. spnorbscl_old /= dtset%spnorbscl) & 277 & .or. positron_old /= dtset%positron & 278 & .or. usewvl_old /= dtset%usewvl & 279 & .or. paw_size_old /= paw_size & 280 & .or. usekden_old/=dtset%usekden & 281 & .or. usexcnhat_old/=dtset%usexcnhat_orig & 282 & .or. any(paw_options_old(:)/=paw_options(:)) & 283 & .or. sum(new_pspso(:))/=0 & 284 & .or. mtypalch>0 & 285 & .or. (dtset%usewvl==1.and.psps%usepaw==1)& 286 & ) gencond=1 287 288 if (present(comm_mpi).and.psps%usepaw==1) then 289 if(xmpi_comm_size(comm_mpi)>1)then 290 call xmpi_sum(gencond,comm_mpi,ierr) 291 end if 292 if (gencond/=0) gencond=1 293 end if 294 ABI_FREE(new_pspso) 295 296 !------------------------------------------------------------- 297 ! Following section is only reached when new computation 298 ! of pseudopotentials is needed 299 !------------------------------------------------------------- 300 301 if (gencond==1) then 302 303 write(msg, '(a,a)' ) ch10,'--- Pseudopotential description ------------------------------------------------' 304 call wrtout(ab_out,msg) 305 306 ABI_MALLOC(ekb,(psps%dimekb*(1-psps%usepaw))) 307 ABI_MALLOC(xccc1d,(psps%n1xccc*(1-psps%usepaw),6)) 308 ABI_MALLOC(ffspl,(psps%mqgrid_ff,2,psps%lnmax)) 309 ABI_MALLOC(vlspl,(psps%mqgrid_vl,2)) 310 if (.not.psps%vlspl_recipSpace) then 311 ABI_MALLOC(dvlspl,(psps%mqgrid_vl,2)) 312 else 313 ABI_MALLOC(dvlspl,(0,0)) 314 end if 315 316 ! PAW: reset flags for optional data 317 if (psps%usepaw==1) then 318 call pawtab_set_flags(pawtab,has_kij=paw_options(1),has_tvale=paw_options(2),& 319 & has_vhnzc=paw_options(3),has_vhtnzc=paw_options(4),& 320 & has_nabla=paw_options(5),has_shapefncg=paw_options(6),& 321 & has_wvl=paw_options(7),has_tproj=paw_options(8),& 322 & has_vminushalf=paw_options(9),has_coretau=paw_options(10)) 323 end if 324 325 ! Read atomic pseudopotential data and get transforms 326 ! for each atom type: two cases, alchemy or not. 327 328 if(mtypalch==0)then 329 ! No alchemical pseudoatom, in all datasets, npsp=ntypat 330 331 do ipsp=1,npsp 332 333 xcccrc=zero 334 ekb(:)=zero;ffspl(:,:,:)=zero;vlspl(:,:)=zero 335 if (.not.psps%vlspl_recipSpace) dvlspl(:, :)=zero 336 if (psps%usepaw==0) xccc1d(:,:)=zero 337 indlmn=>psps%indlmn(:,:,ipsp) 338 indlmn(:,:)=0 339 340 write(msg, '(a,i4,a,t38,a)' )'- pspini: atom type',ipsp,' psp file is',trim(psps%filpsp(ipsp)) 341 call wrtout([std_out, ab_out], msg) 342 343 ! Read atomic psp V(r) and wf(r) to get local and nonlocal psp: 344 ! Cannot use the same call in case of bound checking, because of pawrad/pawtab 345 if(psps%usepaw==0)then 346 call pspatm(dq,dtset,dtfil,ekb,epsatm(ipsp),ffspl,indlmn,ipsp,& 347 pawrad_dum,pawtab_dum,psps,vlspl,dvlspl,xcccrc,xccc1d,psps%nctab(ipsp)) 348 psps%ekb(:,ipsp)=ekb(:) 349 psps%xccc1d(:,:,ipsp)=xccc1d(:,:) 350 else 351 comm_mpi_=xmpi_comm_self;if (present(comm_mpi)) comm_mpi_=comm_mpi 352 call pspatm(dq,dtset,dtfil,ekb,epsatm(ipsp),ffspl,indlmn,ipsp,& 353 pawrad(ipsp),pawtab(ipsp),psps,vlspl,dvlspl,xcccrc,xccc1d,nctab_dum,comm_mpi=comm_mpi_) 354 if (dtset%usefock==1.and.pawtab(ipsp)%has_fock==0) then 355 ABI_BUG('The PAW data file does not contain Fock information. Change the PAW data file!') 356 end if 357 end if 358 359 ! Copy data to psps datastructure. 360 psps%xcccrc(ipsp)=xcccrc 361 psps%znucltypat(ipsp)=psps%znuclpsp(ipsp) 362 psps%ffspl(:,:,:,ipsp)=ffspl(:,:,:) 363 psps%vlspl(:,:,ipsp)=vlspl(:,:) 364 if (.not.psps%vlspl_recipSpace) psps%dvlspl(:, :, ipsp) = dvlspl(:, :) 365 end do ! ipsp 366 367 else 368 ! mtypalch/=0 369 370 npspalch=psps%npspalch 371 ntyppure=npsp-npspalch 372 ntypalch=psps%ntypalch 373 ABI_MALLOC(epsatm_alch,(npspalch)) 374 ABI_MALLOC(ekb_alch,(psps%dimekb,npspalch*(1-psps%usepaw))) 375 ABI_MALLOC(ffspl_alch,(psps%mqgrid_ff,2,psps%lnmax,npspalch)) 376 ABI_MALLOC(xccc1d_alch,(psps%n1xccc*(1-psps%usepaw),6,npspalch)) 377 ABI_MALLOC(xcccrc_alch,(npspalch)) 378 ABI_MALLOC(vlspl_alch,(psps%mqgrid_vl,2,npspalch)) 379 if (.not.psps%vlspl_recipSpace) then 380 ABI_MALLOC(dvlspl_alch,(psps%mqgrid_vl,2,npspalch)) 381 end if 382 ABI_MALLOC(indlmn,(6,psps%lmnmax)) 383 ABI_MALLOC(indlmn_alch,(6,psps%lmnmax,npspalch)) 384 385 ! Allocate NC tables used for mixing. 386 if (psps%usepaw == 0) then 387 ABI_MALLOC(nctab_alch, (npspalch)) 388 do ipspalch=1,npspalch 389 call nctab_init(nctab_alch(ipspalch), psps%mqgrid_vl, .False., .False.) 390 end do 391 end if 392 393 do ipsp=1,npsp 394 write(msg, '(a,i4,a,t38,a)' ) '- pspini: atom type',ipsp,' psp file is',trim(psps%filpsp(ipsp)) 395 call wrtout(ab_out,msg) 396 397 xcccrc=zero 398 ekb(:)=zero;ffspl(:,:,:)=zero;vlspl(:,:)=zero 399 if (.not.psps%vlspl_recipSpace) dvlspl(:, :)=zero 400 if (psps%usepaw==0) xccc1d(:,:)=zero 401 indlmn(:,:)=0 402 403 ! Read atomic psp V(r) and wf(r) to get local and nonlocal psp: 404 if (psps%usepaw==0) then 405 if (ipsp <= ntyppure) then 406 ! Store data in nctab if pure atom. 407 nctab_ptr => psps%nctab(ipsp) 408 else 409 ! Store data in nctab_alch (to be mixed afterwards). 410 nctab_ptr => nctab_alch(ipsp-ntyppure) 411 end if 412 413 call pspatm(dq,dtset,dtfil,ekb,epsatm_psp,ffspl,indlmn,ipsp,& 414 & pawrad_dum,pawtab_dum,psps,vlspl,dvlspl,xcccrc,xccc1d,nctab_ptr) 415 416 else if (psps%usepaw==1) then 417 comm_mpi_=xmpi_comm_self;if (present(comm_mpi)) comm_mpi_=comm_mpi 418 call pspatm(dq,dtset,dtfil,ekb,epsatm_psp,ffspl,indlmn,ipsp,& 419 & pawrad(ipsp),pawtab(ipsp),psps,vlspl,dvlspl,xcccrc,xccc1d,nctab_dum,& 420 & comm_mpi=comm_mpi_) 421 end if 422 423 if (ipsp<=ntyppure) then 424 ! Pure pseudopotentials, leading to pure pseudoatoms 425 epsatm(ipsp)=epsatm_psp 426 psps%znucltypat(ipsp)=psps%znuclpsp(ipsp) 427 if (psps%usepaw==0) psps%ekb(:,ipsp)=ekb(:) 428 psps%ffspl(:,:,:,ipsp)=ffspl(:,:,:) 429 psps%vlspl(:,:,ipsp)=vlspl(:,:) 430 if (.not.psps%vlspl_recipSpace) psps%dvlspl(:, :, ipsp)=dvlspl(:, :) 431 if (psps%usepaw==0) psps%xccc1d(:,:,ipsp)=xccc1d(:,:) 432 psps%xcccrc(ipsp)=xcccrc 433 psps%indlmn(:,:,ipsp)=indlmn(:,:) 434 435 else 436 ! Pseudopotentials for alchemical generation 437 ipspalch=ipsp-ntyppure 438 epsatm_alch(ipspalch)=epsatm_psp 439 ffspl_alch(:,:,:,ipspalch)=ffspl(:,:,:) 440 vlspl_alch(:,:,ipspalch)=vlspl(:,:) 441 if (.not.psps%vlspl_recipSpace) dvlspl_alch(:,:,ipspalch)=dvlspl(:,:) 442 if (psps%usepaw==0) then 443 ekb_alch(:,ipspalch)=ekb(:) 444 xccc1d_alch(:,:,ipspalch)=xccc1d(:,:) 445 end if 446 xcccrc_alch(ipspalch)=xcccrc 447 indlmn_alch(:,:,ipspalch)=indlmn(:,:) 448 ! write(std_out,'(a,6i4)' )' pspini : indlmn_alch(:,1,ipspalch)=',indlmn_alch(:,1,ipspalch) 449 ! write(std_out,'(a,6i4)' )' pspini : indlmn_alch(:,2,ipspalch)=',indlmn_alch(:,2,ipspalch) 450 end if 451 452 end do ! ipsp 453 454 ! Generate data for alchemical pseudos. 455 do itypalch=1,ntypalch 456 itypat=itypalch+ntyppure 457 psps%znucltypat(itypat)=200.0+itypalch ! Convention for alchemical pseudoatoms 458 vlspl(:,:)=zero 459 if (.not.psps%vlspl_recipSpace) dvlspl(:, :) = zero 460 epsatm(itypat)=zero 461 xcccrc=zero 462 if (psps%usepaw==0) xccc1d(:,:)=zero 463 464 ! Here, linear combination of the quantities 465 ! MG: FIXME I think that the mixing of xcccrc is wrong when the xxccrc are different! 466 ! but this is minor bug since alchemical pseudos should not have XCCC (?) 467 do ipspalch=1,npspalch 468 epsatm(itypat) = epsatm(itypat) + epsatm_alch(ipspalch) * psps%mixalch(ipspalch,itypalch) 469 vlspl(:,:) = vlspl(:,:) + vlspl_alch(:,:,ipspalch) * psps%mixalch(ipspalch,itypalch) 470 if (.not.psps%vlspl_recipSpace) then 471 dvlspl(:,:) = dvlspl(:,:) + dvlspl_alch(:,:,ipspalch) * psps%mixalch(ipspalch,itypalch) 472 end if 473 xcccrc = xcccrc + xcccrc_alch(ipspalch) * psps%mixalch(ipspalch,itypalch) 474 if (psps%usepaw==0) then 475 xccc1d(:,:) = xccc1d(:,:) + xccc1d_alch(:,:,ipspalch) * psps%mixalch(ipspalch,itypalch) 476 end if 477 end do ! ipspalch 478 479 psps%vlspl(:,:,itypat)=vlspl(:,:) 480 if (.not.psps%vlspl_recipSpace) psps%dvlspl(:, :, itypat) = dvlspl(:, :) 481 if (psps%usepaw==0) then 482 psps%xccc1d(:,:,itypat)=xccc1d(:,:) 483 end if 484 psps%xcccrc(itypat)=xcccrc 485 486 if (abs(xcccrc) > tol6) then 487 write(std_out, *)"xcccrc", xcccrc 488 ABI_WARNING("Alchemical pseudopotential with nlcc!") 489 end if 490 491 ! Combine the different non-local projectors : for the scalar part then 492 ! the spin-orbit part, treat the different angular momenta 493 ! WARNING : this coding does not work for PAW 494 ilmn=0 495 psps%indlmn(:,:,itypat)=0 496 do ispin=1,2 497 do ilang=0,3 498 if(ispin==2 .and. ilang==0)cycle 499 iproj=0 500 do ipspalch=1,npspalch 501 if(abs(psps%mixalch(ipspalch,itypalch))>tol10)then 502 do ilmn0=1,psps%lmnmax 503 if(indlmn_alch(5,ilmn0,ipspalch)/=0)then 504 if(indlmn_alch(6,ilmn0,ipspalch)==ispin)then 505 if(indlmn_alch(1,ilmn0,ipspalch)==ilang)then 506 ilmn=ilmn+1 ! increment the counter 507 iproj=iproj+1 ! increment the counter, this does not work for PAW 508 if(ilmn>psps%lmnmax)then 509 ABI_BUG('Problem with the alchemical pseudopotentials : ilmn>lmnmax.') 510 end if 511 psps%indlmn(1,ilmn,itypat)=ilang 512 psps%indlmn(2,ilmn,itypat)=indlmn_alch(2,ilmn0,ipspalch) 513 psps%indlmn(3,ilmn,itypat)=iproj ! This does not work for PAW 514 psps%indlmn(4,ilmn,itypat)=ilmn ! This does not work for PAW 515 psps%indlmn(5,ilmn,itypat)=ilmn 516 psps%indlmn(6,ilmn,itypat)=ispin 517 ! The two lines below do not work for PAW 518 if (psps%usepaw==0) then 519 psps%ekb(ilmn,itypat)=psps%mixalch(ipspalch,itypalch) *ekb_alch(ilmn0,ipspalch) 520 end if 521 psps%ffspl(:,:,ilmn,itypat)=ffspl_alch(:,:,ilmn0,ipspalch) 522 end if ! ilang is OK 523 end if ! ispin is OK 524 end if ! ilmn0 exist 525 end do ! ilmn0 526 end if ! mixalch>tol10 527 end do ! ipspalch 528 end do ! ilang 529 end do ! ispin 530 531 end do ! itypalch 532 533 ABI_FREE(epsatm_alch) 534 ABI_FREE(ekb_alch) 535 ABI_FREE(ffspl_alch) 536 ABI_FREE(xccc1d_alch) 537 ABI_FREE(xcccrc_alch) 538 ABI_FREE(vlspl_alch) 539 if (.not.psps%vlspl_recipSpace) then 540 ABI_FREE(dvlspl_alch) 541 end if 542 ABI_FREE(indlmn_alch) 543 ABI_FREE(indlmn) 544 545 ! Mix NC tables. 546 if (psps%usepaw == 0) then 547 call nctab_mixalch(nctab_alch, npspalch, ntypalch, psps%algalch, psps%mixalch, psps%nctab(ntyppure+1:)) 548 do ipspalch=1,npspalch 549 call nctab_free(nctab_alch(ipspalch)) 550 end do 551 ABI_FREE(nctab_alch) 552 end if 553 end if ! mtypalch 554 555 ABI_FREE(ekb) 556 ABI_FREE(ffspl) 557 ABI_FREE(vlspl) 558 ABI_FREE(xccc1d) 559 560 if (.not.psps%vlspl_recipSpace) then 561 ABI_FREE(dvlspl) 562 end if 563 564 end if ! End condition of new computation needed 565 566 !------------------------------------------------------------- 567 ! Following section is always executed 568 !------------------------------------------------------------- 569 !One should move this section of code outside of pspini, 570 !but epsatm is needed, so should be in the psp datastructure. 571 !Compute pseudo correction energy. Will differ from an already 572 !computed one if the number of atom differ ... 573 call pspcor(ecore,epsatm,dtset%natom,ntypat,dtset%typat,psps%ziontypat) 574 if(abs(ecore_old-ecore)>tol8*abs(ecore_old+ecore))then 575 write(msg, '(2x,es15.8,t50,a)' ) ecore,'ecore*ucvol(ha*bohr**3)' 576 ! ecore is useless if iscf<=0, but at least it has been initialized 577 if(dtset%iscf>=0) call wrtout(ab_out,msg) 578 call wrtout(std_out,msg) 579 end if 580 581 !End of pseudopotential output section 582 write(msg, '(2a)' )'--------------------------------------------------------------------------------',ch10 583 call wrtout(ab_out,msg) 584 585 !------------------------------------------------------------- 586 ! Keep track of this call to the routine 587 !------------------------------------------------------------- 588 589 if (ifirst==1) ifirst=0 590 591 mqgridff_old=psps%mqgrid_ff 592 mqgridvl_old=psps%mqgrid_vl 593 mpssoang_old=psps%mpssoang 594 spnorbscl_old = dtset%spnorbscl 595 ixc_old=dtset%ixc 596 hyb_mixing_old=dtset%hyb_mixing 597 gsqcut_old=gsqcut;if (psps%usepaw==1) gsqcutdg_old=gsqcutdg 598 lmnmax_old=psps%lmnmax 599 lnmax_old=psps%lnmax 600 optnlxccc_old=psps%optnlxccc 601 usepaw_old=psps%usepaw 602 dimekb_old=psps%dimekb 603 useylm_old=psps%useylm 604 pawxcdev_old=dtset%pawxcdev 605 positron_old=dtset%positron 606 usewvl_old = dtset%usewvl 607 usekden_old = dtset%usekden 608 usexcnhat_old=dtset%usexcnhat_orig 609 paw_size_old=paw_size 610 ecore_old=ecore 611 paw_options_old(:)=paw_options(:) 612 613 do ipsp=1,npsp 614 pspso_old(ipsp)=psps%pspso(ipsp) 615 if(pspso_zero(ipsp)==0)pspso_zero(ipsp)=psps%pspso(ipsp) 616 end do 617 psps%mproj = maxval(psps%indlmn(3,:,:)) 618 619 if (gencond == 1) call psps_print(psps,std_out,dtset%prtvol) 620 621 ! Write the PSPS.nc file and exit here if requested by the user. 622 if (abs(dtset%prtpsps) == 1) then 623 if (xmpi_comm_rank(xmpi_world) == 0) call psps_ncwrite_path(psps, trim(dtfil%filnam_ds(4))//"_PSPS.nc") 624 if (dtset%prtpsps == -1) then 625 ABI_ERROR_NODUMP("prtpsps == -1 ==> aborting now") 626 end if 627 end if 628 629 call timab(15,2,tsec) 630 631 DBG_EXIT("COLL") 632 633 end subroutine pspini