TABLE OF CONTENTS


ABINIT/m_pspini [ Modules ]

[ Top ] [ Modules ]

NAME

  m_pspini

FUNCTION

  Initialize pseudopotential datastructures from files.

COPYRIGHT

  Copyright (C) 1998-2018 ABINIT group (DCA, XG, GMR, MT, FrD, AF, 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 .

PARENTS

CHILDREN

SOURCE

21 #if defined HAVE_CONFIG_H
22 #include "config.h"
23 #endif
24 
25 #include "abi_common.h"
26 
27 module m_pspini
28 
29  use defs_basis
30  use defs_datatypes
31  use defs_abitypes
32  use m_errors
33  use m_abicore
34  use m_xmpi
35  use m_psxml2ab
36  !use m_psps
37 
38  use m_time,      only : timab
39  use m_io_tools,  only : open_file
40  use m_pawrad,    only : pawrad_type
41  use m_pawtab,    only : pawtab_type, pawtab_set_flags
42  use m_psps,      only : psps_print, psps_ncwrite, nctab_init, nctab_free, nctab_mixalch, test_xml_xmlpaw_upf, &
43                          nctab_eval_tcorespl
44  use m_pawpsp,    only : pawpsp_bcast, pawpsp_read_pawheader, pawpsp_read_header_xml,&
45                          pawpsp_header_type, pawpsp_wvl, pawpsp_7in, pawpsp_17in
46  use m_pawxmlps,  only : paw_setup_free,paw_setuploc
47  use m_pspheads,  only : pawpsxml2ab
48 #if defined HAVE_BIGDFT
49  use BigDFT_API, only : dictionary, atomic_info, dict_init, dict_free, UNINITIALIZED
50 #endif
51 
52  use m_psp1,       only : psp1in
53  use m_psp5,       only : psp5in
54  use m_psp6,       only : psp6in
55  use m_psp8,       only : psp8in
56  use m_psp9,       only : psp9in
57  use m_upf2abinit, only : upf2abinit
58  use m_psp_hgh,    only : psp2in, psp3in, psp10in
59  use m_wvl_descr_psp,  only : wvl_descr_psp_fill
60 
61  implicit none
62 
63  private

ABINIT/psp_dump_outputs [ Functions ]

[ Top ] [ Functions ]

NAME

 psp_dump_outputs

FUNCTION

 (To be described ...)

COPYRIGHT

 Copyright (C) 2017-2018 ABINIT group (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 .
 For the initials of contributors, see ~abinit/doc/developers/contributors.txt

INPUTS

 (to be filled)

OUTPUT

 (to be filled)

SIDE EFFECTS

 (to be filled)

PARENTS

      pspatm

CHILDREN

SOURCE

1451 subroutine psp_dump_outputs(pfx,pspcod,lmnmax,lnmax,mpssoang, &
1452 &      mqgrid,n1xccc,mmax,maxrad,epsatm,qchrg,xcccrc,nctab, &
1453 &      indlmn,nproj,ekb,ffspl,vlspl,xccc1d)
1454 
1455  use defs_basis
1456  use m_errors
1457 
1458  use defs_datatypes, only : nctab_t
1459 
1460 !This section has been created automatically by the script Abilint (TD).
1461 !Do not modify the following lines by hand.
1462 #undef ABI_FUNC
1463 #define ABI_FUNC 'psp_dump_outputs'
1464 !End of the abilint section
1465 
1466  implicit none
1467 
1468 !Arguments ------------------------------------
1469 !scalars
1470  character(len=*), intent(in) :: pfx
1471  integer,intent(in) :: pspcod,lmnmax,lnmax,mpssoang,mqgrid,n1xccc
1472  integer,intent(in) :: mmax
1473  real(dp),intent(in) :: maxrad,epsatm,qchrg,xcccrc
1474  type(nctab_t),intent(in) :: nctab
1475 !arrays
1476  integer,intent(in) :: indlmn(6,lmnmax),nproj(mpssoang)
1477  real(dp),intent(in) :: ekb(lnmax),ffspl(mqgrid,2,lnmax),vlspl(mqgrid,2)
1478  real(dp),intent(in) :: xccc1d(n1xccc,6)
1479 
1480 !Local variables ------------------------------
1481 !scalars
1482  integer, parameter :: dump = 64
1483  integer :: ierr, i, j ,k
1484  character(len=500) :: msg
1485 
1486  ! *********************************************************************
1487 
1488  open(unit=dump, file=trim(pfx)//"_psp_info.yaml", status='REPLACE', err=10, iostat=ierr)
1489 
1490  write(dump,'(3a)') "%YAML 1.2", ch10, "---"
1491 
1492  write(dump, '(2a)') ch10, "# Pseudopotential info"
1493  write(dump, '(a,1x,i8)') "pspcod:", pspcod
1494 
1495  write(dump, '(2a)') ch10, "# Array dimensions"
1496  write(dump, '(a)') "dims:"
1497  write(dump, '(4x,a,1x,i8)') "lmnmax:", lmnmax
1498  write(dump, '(4x,a,1x,i8)') "lnmax:", lnmax
1499  write(dump, '(4x,a,1x,i8)') "mpssoang:", mpssoang
1500  write(dump, '(4x,a,1x,i8)') "mqgrid:", mqgrid
1501  write(dump, '(4x,a,1x,i8)') "n1xccc:", n1xccc
1502  write(dump, '(4x,a,1x,i8)') "mmax:", mmax
1503 
1504  write(dump, '(2a)') ch10, "# Quantities"
1505  write(dump, '(a,1x,e12.5)') "maxrad:", maxrad
1506  write(dump, '(a,1x,e12.5)') "epsatm:", epsatm
1507  write(dump, '(a,1x,e12.5)') "qchrg:", qchrg
1508  write(dump, '(a,1x,e12.5)') "xcccrc:", xcccrc
1509 
1510  write(dump, '(2a)') ch10, "# Structure: nctab"
1511  write(dump, '(a)') "nctab:"
1512  write(dump,'(4x,a,":",1x,i4)') "mqgrid_vl", nctab%mqgrid_vl
1513  write(dump,'(4x,a,":",1x,l4)') "has_tvale", nctab%has_tvale
1514  write(dump,'(4x,a,":",1x,l4)') "has_tcore", nctab%has_tcore
1515  write(dump,'(4x,a,":",1x,e12.5)') "dncdq0", nctab%dncdq0
1516  write(dump,'(4x,a,":",1x,e12.5)') "d2ncdq0", nctab%d2ncdq0
1517  write(dump,'(4x,a,":",1x,e12.5)') "dnvdq0", nctab%dnvdq0
1518 
1519  if ( nctab%has_tvale ) then
1520    write(dump, '(2a)') ch10, "# Array: nctab_tvalespl(mqgrid_vl,2)"
1521    write(dump, '(a)') "nctab_tvalespl:"
1522    do j=1,2
1523      do i=1,nctab%mqgrid_vl
1524        if ( i == 1 ) then
1525          write(dump,'(4x,a,1x,e12.5)') "- -", nctab%tvalespl(i,j)
1526        else
1527          write(dump,'(4x,a,1x,e12.5)') "  -", nctab%tvalespl(i,j)
1528        end if
1529      end do
1530    end do
1531  end if
1532 
1533  if ( nctab%has_tcore ) then
1534    write(dump, '(2a)') ch10, "# Array: nctab_tcorespl(mqgrid_vl,2)"
1535    write(dump, '(a)') "nctab_tcorespl:"
1536    do j=1,2
1537      do i=1,nctab%mqgrid_vl
1538        if ( i == 1 ) then
1539          write(dump,'(4x,a,1x,e12.5)') "- -", nctab%tcorespl(i,j)
1540        else
1541          write(dump,'(4x,a,1x,e12.5)') "  -", nctab%tcorespl(i,j)
1542        end if
1543      end do
1544    end do
1545  end if
1546 
1547  write(dump, '(2a)') ch10, "# Array: integer indlmn(6,lmnmax)"
1548  write(dump, '(a)') "indlmn:"
1549  do i=1,lmnmax
1550    write(dump,'(4x,a,i4,5(",",i4),a)') "- [", indlmn(:,i), "]"
1551  end do
1552 
1553  write(dump, '(2a)') ch10, "# Array: integer nproj(mpssoang)"
1554  write(dump, '(a)') "nproj:"
1555  do i=1,mpssoang
1556    write(dump,'(4x,"-",1x,i4)') nproj(i)
1557  end do
1558 
1559  write(dump, '(2a)') ch10, "# Array: double ekb(lnmax)"
1560  write(dump, '(a)') "ekb:"
1561  do i=1,lnmax
1562    write(dump,'(4x,"-",1x,e12.5)') ekb(i)
1563  end do
1564 
1565  write(dump, '(2a)') ch10, "# Array: ffspl(mqgrid,2,lnmax)"
1566  write(dump, '(a)') "ffspl:"
1567  do k=1,lnmax
1568    do j=1,2
1569      do i=1,mqgrid
1570        if ( (i == 1) .and. (j == 1) ) then
1571          write(dump,'(4x,a,1x,e12.5)') "- - -", ffspl(i,j,k)
1572        else if ( i == 1 ) then
1573          write(dump,'(4x,a,1x,e12.5)') "  - -", ffspl(i,j,k)
1574        else
1575          write(dump,'(4x,a,1x,e12.5)') "    -", ffspl(i,j,k)
1576        end if
1577      end do
1578    end do
1579  end do
1580 
1581  write(dump, '(2a)') ch10, "# Array: vlspl(mqgrid,2)"
1582  write(dump, '(a)') "vlspl:"
1583  do j=1,2
1584    do i=1,mqgrid
1585      if ( i == 1 ) then
1586        write(dump,'(4x,a,1x,e12.5)') "- -", vlspl(i,j)
1587      else
1588        write(dump,'(4x,a,1x,e12.5)') "  -", vlspl(i,j)
1589      end if
1590    end do
1591  end do
1592 
1593  write(dump, '(2a)') ch10, "# Array: xccc1d(n1xccc,6)"
1594  write(dump, '(a)') "xccc1d:"
1595  do j=1,6
1596    do i=1,mqgrid
1597      if ( i == 1 ) then
1598        write(dump,'(4x,a,1x,e12.5)') "- -", xccc1d(i,j)
1599      else
1600        write(dump,'(4x,a,1x,e12.5)') "  -", xccc1d(i,j)
1601      end if
1602    end do
1603  end do
1604 
1605  write (dump,'(2a)') ch10, "..."
1606 
1607  close(dump)
1608 
1609  return
1610  10 continue
1611 
1612  if ( ierr /= 0 ) then
1613    write(msg,'(a,a,a,i8)') "Error writing pseudopotential information", ch10, "IOSTAT=", ierr
1614    MSG_WARNING(msg)
1615  end if
1616 
1617 end subroutine psp_dump_outputs

ABINIT/pspatm [ Functions ]

[ Top ] [ 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 "UPF PWSCF format" (pspcod=11)

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

PARENTS

      pspini

CHILDREN

      nctab_eval_tcorespl,pawpsp_17in,pawpsp_7in,pawpsp_bcast
      pawpsp_read_header_xml,pawpsp_read_pawheader,pawpsp_wvl,psp10in,psp1in
      psp2in,psp3in,psp5in,psp6in,psp8in,psp9in,psp_dump_outputs
      psxml2abheader,test_xml_xmlpaw_upf,timab,upf2abinit,wrtout
      wvl_descr_psp_fill

SOURCE

 830 subroutine pspatm(dq,dtset,dtfil,ekb,epsatm,ffspl,indlmn,ipsp,pawrad,pawtab,&
 831 &  psps,vlspl,dvlspl,xcccrc,xccc1d,nctab,comm_mpi)
 832 
 833 
 834 !This section has been created automatically by the script Abilint (TD).
 835 !Do not modify the following lines by hand.
 836 #undef ABI_FUNC
 837 #define ABI_FUNC 'pspatm'
 838 !End of the abilint section
 839 
 840  implicit none
 841 
 842 !Arguments ---------------------------------------------
 843 !scalars
 844  integer,intent(in) :: ipsp
 845  integer, optional,intent(in) :: comm_mpi
 846  real(dp),intent(in) :: dq
 847  real(dp),intent(out) :: epsatm,xcccrc
 848  type(dataset_type),intent(in) :: dtset
 849  type(datafiles_type),intent(in) :: dtfil
 850  type(pawrad_type),intent(inout) :: pawrad
 851  type(pawtab_type),intent(inout) :: pawtab
 852  type(nctab_t),intent(inout) :: nctab
 853  type(pseudopotential_type),intent(inout) :: psps
 854 !arrays
 855  integer,intent(out) :: indlmn(6,psps%lmnmax)
 856  real(dp),intent(out) :: dvlspl(psps%mqgrid_vl,2)
 857  real(dp),intent(inout) :: ekb(psps%dimekb*(1-psps%usepaw))
 858  real(dp),intent(inout) :: ffspl(psps%mqgrid_ff,2,psps%lnmax)
 859  real(dp),intent(out) :: vlspl(psps%mqgrid_vl,2)
 860  real(dp),intent(inout) :: xccc1d(psps%n1xccc*(1-psps%usepaw),6)
 861 
 862 !Local variables ---------------------------------------
 863 !scalars
 864  integer :: ii,il,ilmn,iln,iln0,lloc,lmax,me,mmax
 865  integer :: paral_mode,pspcod,pspdat,pspxc,useupf,usexml,xmlpaw,unt
 866  real(dp) :: maxrad,qchrg,r2well,zion,znucl
 867  character(len=500) :: message,errmsg
 868  character(len=fnlen) :: title
 869  character(len=fnlen) :: filnam
 870  type(pawpsp_header_type):: pawpsp_header
 871  type(pspheader_type) :: pspheads_tmp
 872 !arrays
 873  integer,allocatable :: nproj(:)
 874  real(dp) :: tsec(2),ecut_tmp(3,2)
 875  real(dp),allocatable :: e990(:),e999(:),ekb1(:),ekb2(:),epspsp(:),rcpsp(:)
 876  real(dp),allocatable :: rms(:)
 877 #if defined HAVE_PSML
 878 !!  usexml= 0 for non xml ps format ; =1 for xml ps format
 879  character(len=3) :: atmsymb
 880  character(len=30) :: creator
 881  type(pspheader_type) :: psphead
 882 #endif
 883 
 884 ! ******************************************************************************
 885 
 886 !paral_mode defines how we access to the psp file
 887 !  paral_mode=0: all processes access to the file (sequentially)
 888 !  paral_mode=1: only proc. 0 access to the file and then broadcast
 889  paral_mode=0
 890  if (present(comm_mpi)) then
 891    if (psps%usepaw==1.and.xmpi_comm_size(comm_mpi)>1) paral_mode=1
 892  end if
 893  me=0;if (paral_mode==1) me=xmpi_comm_rank(comm_mpi)
 894 
 895  if (paral_mode == 1) then
 896    ABI_CHECK(psps%usepaw==1, "paral_mode==1 is only compatible with PAW, see call to pawpsp_bcast below")
 897  end if
 898 
 899  nctab%has_tvale = .False.; nctab%has_tcore = .False.
 900 
 901  if (me==0) then
 902 !  Dimensions of form factors and Vloc q grids must be the same in Norm-Conserving case
 903    if (psps%usepaw==0 .and. psps%mqgrid_ff/=psps%mqgrid_vl) then
 904      write(message, '(a,a,a,a,a)' )&
 905 &     'Dimension of q-grid for nl form factors (mqgrid_ff)',ch10,&
 906 &     'is different from dimension of q-grid for Vloc (mqgrid_vl) !',ch10,&
 907 &     'This is not allowed for norm-conserving psp.'
 908      MSG_ERROR(message)
 909    end if
 910 
 911    write(message, '(a,t38,a)' )'- pspatm: opening atomic psp file',trim(psps%filpsp(ipsp))
 912    call wrtout(ab_out,  message,'COLL')
 913    call wrtout(std_out,  message,'COLL')
 914    !write(message, "(2a)")"- md5: ",trim(psps%md5_pseudos(ipsps))
 915 
 916    !  Check if the file pseudopotential file is written in (XML| XML-PAW | UPF)
 917    call test_xml_xmlpaw_upf(psps%filpsp(ipsp), usexml, xmlpaw, useupf)
 918 
 919 !  ----------------------------------------------------------------------------
 920 !  allocate nproj here: can be read in now for UPF
 921    ABI_ALLOCATE(nproj,(psps%mpssoang))
 922    nproj(:)=0
 923 
 924    if (usexml /= 1 .and. useupf /= 1) then
 925 
 926 !    Open the atomic data file, and read the three first lines
 927 !    These three first lines have a similar format in all allowed psp files
 928 
 929 !    Open atomic data file (note: formatted input file)
 930      if (open_file(psps%filpsp(ipsp), message, unit=tmp_unit, form='formatted', status='old') /= 0) then
 931        MSG_ERROR(message)
 932      end if
 933      rewind (unit=tmp_unit,err=10,iomsg=errmsg)
 934 
 935 !    Read and write some description of file from first line (character data)
 936      read (tmp_unit,'(a)',err=10,iomsg=errmsg) title
 937      write(message, '(a,a)' ) '- ',trim(title)
 938      call wrtout(ab_out,message,'COLL')
 939      call wrtout(std_out,  message,'COLL')
 940 
 941 !    Read and write more data describing psp parameters
 942      read (tmp_unit,*,err=10,iomsg=errmsg) znucl,zion,pspdat
 943      write(message,'(a,f9.5,f10.5,2x,i8,t47,a)')'-',znucl,zion,pspdat,'znucl, zion, pspdat'
 944      call wrtout(ab_out,message,'COLL')
 945      call wrtout(std_out,message,'COLL')
 946 
 947      read (tmp_unit,*,err=10,iomsg=errmsg) pspcod,pspxc,lmax,lloc,mmax,r2well
 948      if(pspxc<0) then
 949        write(message, '(i5,i8,2i5,i10,f10.5,t47,a)' ) &
 950 &       pspcod,pspxc,lmax,lloc,mmax,r2well,'pspcod,pspxc,lmax,lloc,mmax,r2well'
 951      else
 952        write(message, '(4i5,i10,f10.5,t47,a)' ) &
 953 &       pspcod,pspxc,lmax,lloc,mmax,r2well,'pspcod,pspxc,lmax,lloc,mmax,r2well'
 954      end if
 955      call wrtout(ab_out,message,'COLL')
 956      call wrtout(std_out,  message,'COLL')
 957 
 958    else if (usexml == 1 .and. xmlpaw == 0) then
 959 
 960 ! the following is probably useless - already read in everything in inpspheads
 961 #if defined HAVE_PSML
 962 !     write(message,'(a,a)') &
 963 !&     '- pspatm: Reading pseudopotential header in XML form from ', trim(psps%filpsp(ipsp))
 964 !     call wrtout(ab_out,message,'COLL')
 965 !     call wrtout(std_out,  message,'COLL')
 966 
 967      call psxml2abheader( psps%filpsp(ipsp), psphead, atmsymb, creator, 0 )
 968      znucl = psphead%znuclpsp
 969      zion = psphead%zionpsp
 970      pspdat = psphead%pspdat
 971      pspcod = psphead%pspcod
 972      pspxc =  psphead%pspxc
 973      lmax = psphead%lmax
 974      !lloc   = 0 ! does this mean s? in psml case the local potential can be different from any l channel
 975      lloc = -1
 976      mmax = -1
 977      r2well = 0
 978 
 979      write(message,'(a,1x,a3,3x,a)') "-",atmsymb,trim(creator)
 980      call wrtout(ab_out,message,'COLL')
 981      call wrtout(std_out,message,'COLL')
 982      write(message,'(a,f9.5,f10.5,2x,i8,t47,a)')'-',znucl,zion,pspdat,'znucl, zion, pspdat'
 983      call wrtout(ab_out,message,'COLL')
 984      call wrtout(std_out,message,'COLL')
 985      if(pspxc<0) then
 986        write(message, '(i5,i8,2i5,i10,f10.5,t47,a)' ) &
 987 &       pspcod,pspxc,lmax,lloc,mmax,r2well,'pspcod,pspxc,lmax,lloc,mmax,r2well'
 988      else
 989        write(message, '(4i5,i10,f10.5,t47,a)' ) &
 990 &       pspcod,pspxc,lmax,lloc,mmax,r2well,'pspcod,pspxc,lmax,lloc,mmax,r2well'
 991      end if
 992      call wrtout(ab_out,message,'COLL')
 993      call wrtout(std_out,message,'COLL')
 994 #else
 995      write(message,'(a,a)')  &
 996 &     'ABINIT is not compiled with XML support for reading this type of pseudopotential ', &
 997 &     trim(psps%filpsp(ipsp))
 998      MSG_BUG(message)
 999 #endif
1000 ! END useless
1001    else if (usexml == 1 .and. xmlpaw == 1) then
1002      write(message,'(a,a)')  &
1003 &     '- pspatm : Reading pseudopotential header in XML form from ', trim(psps%filpsp(ipsp))
1004      call wrtout(ab_out,message,'COLL')
1005      call wrtout(std_out,  message,'COLL')
1006 
1007 !    Return header informations
1008      call pawpsxml2ab(psps%filpsp(ipsp),ecut_tmp, pspheads_tmp,0)
1009      lmax=pspheads_tmp%lmax
1010      pspxc=pspheads_tmp%pspxc
1011      znucl=pspheads_tmp%znuclpsp
1012      pawpsp_header%basis_size=pspheads_tmp%pawheader%basis_size
1013      pawpsp_header%l_size=pspheads_tmp%pawheader%l_size
1014      pawpsp_header%lmn_size=pspheads_tmp%pawheader%lmn_size
1015      pawpsp_header%mesh_size=pspheads_tmp%pawheader%mesh_size
1016      pawpsp_header%pawver=pspheads_tmp%pawheader%pawver
1017      pawpsp_header%shape_type=pspheads_tmp%pawheader%shape_type
1018      pawpsp_header%rpaw=pspheads_tmp%pawheader%rpaw
1019      pawpsp_header%rshp=pspheads_tmp%pawheader%rshp
1020      lloc=0; pspcod=17
1021 
1022    else if (useupf == 1) then
1023      if (psps%usepaw /= 0) then
1024        MSG_ERROR("UPF format not allowed with PAW (USPP part not read yet)")
1025      end if
1026 
1027      pspcod = 11
1028      r2well = 0
1029 
1030 !    should initialize znucl,zion,pspxc,lmax,lloc,mmax
1031      call upf2abinit (psps%filpsp(ipsp), znucl, zion, pspxc, lmax, lloc, mmax, &
1032 &     psps, epsatm, xcccrc, indlmn, ekb, ffspl, nproj, vlspl, xccc1d)
1033 
1034    else
1035      MSG_ERROR("You should not be here! erroneous type or pseudopotential input")
1036    end if
1037 
1038 !  ------------------------------------------------------------------------------
1039 !  Check data for consistency against main routine input
1040 
1041 !  Does required spin-orbit characteristics agree with format
1042 !  TODO: in case of pspcod 5 (phoney) and 8 (oncvpsp) this is not specific enough.
1043 !  they can be non-SOC as well.
1044 !  HGH is ok - can always turn SOC on or off.
1045 !  PAW is ok - can be used with or without SOC
1046 !  write(std_out,*) pspso
1047    if((pspcod/=3).and.(pspcod/=5).and.(pspcod/=8).and.(pspcod/=10).and. &
1048 &     (pspcod/=7).and.(pspcod/=17))then
1049 !    If pspso requires internal characteristics, set it to 1 for non-HGH psps
1050      if(psps%pspso(ipsp)==1) psps%pspso(ipsp)=0
1051      if(psps%pspso(ipsp)/=0)then
1052        write(message, '(3a,i0,3a)' )&
1053 &       'Pseudopotential file cannot give spin-orbit characteristics,',ch10,&
1054 &       'while pspso(itypat)= ',psps%pspso(ipsp),'.',ch10,&
1055 &       'Action: check your pseudopotential and input files for consistency.'
1056        MSG_ERROR(message)
1057      end if
1058    end if
1059 
1060 !  Does nuclear charge znuclpsp agree with psp input znucl
1061    !write(std_out,*)znucl,ipsp,psps%znuclpsp(ipsp)
1062    !MGNAG: v5[66] gives NAG in %znuclpsp if -nan
1063    if (abs(psps%znuclpsp(ipsp)-znucl)>tol8) then
1064      write(message, '(a,f10.5,2a,f10.5,5a)' )&
1065 &     'Pseudopotential file znucl=',znucl,ch10,&
1066 &     'does not equal input znuclpsp=',psps%znuclpsp(ipsp),' better than 1e-08 .',ch10,&
1067 &     'znucl is read from the psp file in pspatm, while',ch10,&
1068 &     'znuclpsp is read in iofn2.'
1069      MSG_BUG(message)
1070    end if
1071 
1072 !  Is the highest angular momentum within limits?
1073 !  Recall mpsang is 1+highest l for nonlocal correction.
1074 !  Nonlocal corrections for s, p, d, and f are supported.
1075    if (lmax+1>psps%mpsang) then
1076      write(message, '(a,i0,a,i0,a,a)' )&
1077 &     'input lmax+1= ',lmax+1,' exceeds mpsang= ',psps%mpsang,ch10,&
1078 &     'indicates input lmax too large for dimensions.'
1079      MSG_BUG(message)
1080    end if
1081 
1082 !  Check several choices for ixc against pspxc
1083 !  ixc is from ABINIT code; pspxc is from atomic psp file
1084    if (dtset%ixc==0) then
1085      MSG_WARNING('Note that input ixc=0 => no xc is being used.')
1086    else if(dtset%ixc/=pspxc) then
1087      write(message, '(a,i0,3a,i0,9a)' )&
1088 &     'Pseudopotential file pspxc= ',pspxc,',',ch10,&
1089 &     'not equal to input ixc= ',dtset%ixc,'.',ch10,&
1090 &     'These parameters must agree to get the same xc ',ch10,&
1091 &     'in ABINIT code as in psp construction.',ch10,&
1092 &     'Action: check psp design or input file.',ch10,&
1093 &     'Assume experienced user. Execution will continue.'
1094      MSG_WARNING(message)
1095    end if
1096 
1097    if (lloc>lmax .and. pspcod/=4 .and. pspcod/=8 .and. pspcod/=10) then
1098      write(message, '(a,2i12,a,a,a,a)' )&
1099 &     'lloc,lmax=',lloc,lmax,ch10,&
1100 &     'chosen l of local psp exceeds range from input data.',ch10,&
1101 &     'Action: check pseudopotential input file.'
1102      MSG_ERROR(message)
1103    end if
1104 
1105 !  Does the pspcod agree with type of calculation (paw or not)?
1106    if (((pspcod/=7.and.pspcod/=17).and.psps%usepaw==1).or.((pspcod==7.or.pspcod==17).and.psps%usepaw==0)) then
1107      write(message, '(a,i2,a,a,i0,a)' )&
1108 &     'In reading atomic psp file, finds pspcod= ',pspcod,ch10,&
1109 &     'This is not an allowed value with usepaw= ',psps%usepaw,'.'
1110      MSG_BUG(message)
1111    end if
1112 
1113    if (.not.psps%vlspl_recipSpace .and. &
1114 &   (pspcod /= 2 .and. pspcod /= 3 .and. pspcod /= 10 .and. pspcod /= 7)) then
1115 !    The following "if" statement can substitute the one just before once libBigDFT
1116 !    has been upgraded to include pspcod 10
1117 !    if (.not.psps%vlspl_recipSpace .and. (pspcod /= 2 .and. pspcod /= 3 .and. pspcod /= 10)) then
1118      write(message, '(a,i2,a,a)' )&
1119 &     'In reading atomic psp file, finds pspcod=',pspcod,ch10,&
1120 &     'This is not an allowed value with real space computation.'
1121      MSG_BUG(message)
1122    end if
1123 
1124 !  MJV 16/6/2009 added pspcod 11 for upf format
1125    if( pspcod<1 .or. (pspcod>11.and.pspcod/=17) ) then
1126      write(message, '(a,i0,4a)' )&
1127 &     'In reading atomic psp file, finds pspcod= ',pspcod,ch10,&
1128 &     'This is not an allowed value. Allowed values are 1 to 11 .',ch10,&
1129 &     'Action: check pseudopotential input file.'
1130      MSG_ERROR(message)
1131    end if
1132 
1133 !  -----------------------------------------------------------------------
1134 !  Set various terms to 0 in case not defined below
1135    ABI_ALLOCATE(e990,(psps%mpssoang))
1136    ABI_ALLOCATE(e999,(psps%mpssoang))
1137    ABI_ALLOCATE(rcpsp,(psps%mpssoang))
1138    ABI_ALLOCATE(rms,(psps%mpssoang))
1139    ABI_ALLOCATE(epspsp,(psps%mpssoang))
1140    ABI_ALLOCATE(ekb1,(psps%mpssoang))
1141    ABI_ALLOCATE(ekb2,(psps%mpssoang))
1142    e990(:)=zero ;e999(:)=zero
1143    rcpsp(:)=zero;rms(:)=zero
1144    ekb1(:)=zero ;ekb2(:)=zero
1145    epspsp(:)=zero
1146    qchrg=0
1147 
1148 !  ----------------------------------------------------------------------
1149    if(pspcod==1 .or. pspcod==4)then
1150 
1151 !    Teter pseudopotential (pspcod=1 or 4)
1152      call psp1in(dq,ekb,ekb1,ekb2,epsatm,epspsp,&
1153 &     e990,e999,ffspl,indlmn,lloc,lmax,psps%lmnmax,psps%lnmax,&
1154 &     mmax,psps%mpsang,psps%mqgrid_ff,nproj,psps%n1xccc,pspcod,qchrg,psps%qgrid_ff,&
1155 &     rcpsp,rms,psps%useylm,vlspl,xcccrc,xccc1d,zion,psps%znuclpsp(ipsp))
1156 
1157    else if (pspcod==2)then
1158 
1159 !    GTH pseudopotential
1160      call psp2in(dtset,ekb,epsatm,ffspl,indlmn,ipsp,lmax,nproj,psps,vlspl,dvlspl,zion)
1161      xccc1d(:,:)=0.0d0 ; qchrg=0.0d0 ; xcccrc=0.0d0
1162 
1163    else if (pspcod==3)then
1164 
1165 !    HGH pseudopotential
1166      call psp3in(dtset,ekb,epsatm,ffspl,indlmn,ipsp,lmax,nproj,psps, psps%pspso(ipsp), &
1167 &     vlspl,zion)
1168      xccc1d(:,:)=0.0d0 ; qchrg=0.0d0 ; xcccrc=0.0d0
1169 
1170    else if (pspcod==5)then
1171 
1172 !    Old phoney pseudopotentials
1173      call psp5in(ekb,ekb1,ekb2,epsatm,epspsp,&
1174 &     e990,e999,ffspl,indlmn,lloc,lmax,psps%lmnmax,psps%lnmax,&
1175 &     mmax,psps%mpsang,psps%mpssoang,psps%mqgrid_ff,nproj,psps%n1xccc,psps%pspso(ipsp),qchrg,psps%qgrid_ff,&
1176 &     rcpsp,rms,psps%useylm,vlspl,xcccrc,xccc1d,zion,psps%znuclpsp(ipsp))
1177 
1178    else if (pspcod==6)then
1179 
1180 !    FHI pseudopotentials
1181      call psp6in(ekb,epsatm,ffspl,indlmn,lloc,lmax,psps%lmnmax,psps%lnmax,mmax,&
1182 &     psps%mpsang,psps%mqgrid_ff,nproj,psps%n1xccc,psps%optnlxccc,psps%positron,qchrg,psps%qgrid_ff,psps%useylm,vlspl,&
1183 &     xcccrc,xccc1d,zion,psps%znuclpsp(ipsp))
1184 
1185    else if (pspcod==7)then
1186 !    PAW "pseudopotentials"
1187      call pawpsp_7in(epsatm,ffspl,dtset%icoulomb,dtset%ixc,&
1188 &     lmax,psps%lnmax,mmax,psps%mqgrid_ff,psps%mqgrid_vl,&
1189 &     pawrad,pawtab,dtset%pawxcdev,psps%qgrid_ff,psps%qgrid_vl,&
1190 &     dtset%usewvl,dtset%usexcnhat_orig,vlspl,xcccrc,dtset%xclevel,&
1191 &     dtset%xc_denpos,zion,psps%znuclpsp(ipsp))
1192 
1193    else if (pspcod==8)then
1194 
1195 !    DRH pseudopotentials
1196      call psp8in(ekb,epsatm,ffspl,indlmn,lloc,lmax,psps%lmnmax,psps%lnmax,mmax,&
1197 &     psps%mpsang,psps%mpssoang,psps%mqgrid_ff,psps%mqgrid_vl,nproj,psps%n1xccc,psps%pspso(ipsp),&
1198 &     qchrg,psps%qgrid_ff,psps%qgrid_vl,psps%useylm,vlspl,xcccrc,xccc1d,zion,psps%znuclpsp(ipsp),nctab,maxrad)
1199 
1200 #if defined DEV_YP_DEBUG_PSP
1201      call psp_dump_outputs("DBG",pspcod,psps%lmnmax,psps%lnmax,psps%mpssoang, &
1202 &     psps%mqgrid_ff,psps%n1xccc,mmax,maxrad,epsatm,qchrg,xcccrc,nctab, &
1203 &     indlmn,nproj,ekb,ffspl,vlspl,xccc1d)
1204 #endif
1205 
1206    else if (pspcod==9)then
1207 
1208 #if defined HAVE_PSML
1209      call psp9in(psps%filpsp(ipsp),ekb,epsatm,ffspl,indlmn,lloc,lmax,psps%lmnmax,psps%lnmax,mmax,&
1210 &     psps%mpsang,psps%mpssoang,psps%mqgrid_ff,psps%mqgrid_vl,nproj,psps%n1xccc, &
1211 &     psps%pspso(ipsp),qchrg,psps%qgrid_ff,psps%qgrid_vl,psps%useylm,vlspl,&
1212 &     xcccrc,xccc1d,zion,psps%znuclpsp(ipsp),nctab,maxrad)
1213 
1214 #if defined DEV_YP_DEBUG_PSP
1215      call psp_dump_outputs("DBG",pspcod,psps%lmnmax,psps%lnmax,psps%mpssoang, &
1216 &     psps%mqgrid_ff,psps%n1xccc,mmax,maxrad,epsatm,qchrg,xcccrc,nctab, &
1217 &     indlmn,nproj,ekb,ffspl,vlspl,xccc1d)
1218 #endif
1219 #else
1220      write(message,'(2a)')  &
1221 &     'ABINIT is not compiled with XML support for reading this type of pseudopotential ', &
1222 &     trim(psps%filpsp(ipsp))
1223      MSG_BUG(message)
1224 #endif
1225 
1226    else if (pspcod==10)then
1227 
1228 !    HGH pseudopotential, full h/k matrix read
1229      call psp10in(dtset,ekb,epsatm,ffspl,indlmn,ipsp,lmax,nproj,psps, psps%pspso(ipsp), &
1230 &     vlspl,zion)
1231      xccc1d(:,:)=0.0d0 ; qchrg=0.0d0 ; xcccrc=0.0d0
1232 
1233 !    NB for pspcod 11 the reading has already been done above.
1234    else if (pspcod==17)then
1235 !    PAW XML "pseudopotentials"
1236      call pawpsp_17in(epsatm,ffspl,dtset%icoulomb,ipsp,dtset%ixc,lmax,&
1237 &     psps%lnmax,mmax,psps%mqgrid_ff,psps%mqgrid_vl,pawpsp_header,pawrad,pawtab,&
1238 &     dtset%pawxcdev,psps%qgrid_ff,psps%qgrid_vl,dtset%usewvl,&
1239 &     dtset%usexcnhat_orig,vlspl,xcccrc,&
1240 &     dtset%xclevel,dtset%xc_denpos,pspheads_tmp%zionpsp,psps%znuclpsp(ipsp))
1241      call paw_setup_free(paw_setuploc)
1242    end if
1243 
1244    close (unit=tmp_unit)
1245 
1246 !  ----------------------------------------------------------------------
1247    if (pspcod==2 .or. pspcod==3 .or. pspcod==10)then
1248      write(message, '(a,a,a,a,a,a,a,a,a,a)' )ch10,&
1249 &     ' pspatm : COMMENT -',ch10,&
1250 &     '  the projectors are not normalized,',ch10,&
1251 &     '  so that the KB energies are not consistent with ',ch10,&
1252 &     '  definition in PRB44, 8503 (1991). ',ch10,& ! [[cite:Gonze1991]]
1253 &     '  However, this does not influence the results obtained hereafter.'
1254      call wrtout(ab_out,message,'COLL')
1255      call wrtout(std_out,message,'COLL')
1256 !    The following lines are added to keep backward compatibilty
1257      maxrad=zero
1258 #if defined HAVE_BIGDFT
1259      do ii=1,size(psps%gth_params%psppar,1)-1 ! psppar first dim begins at 0
1260        if (psps%gth_params%psppar(ii,0,ipsp)/=zero) maxrad=max(maxrad,psps%gth_params%psppar(ii,0,ipsp))
1261      end do
1262      if (abs(maxrad)<=tol12) then
1263        psps%gth_params%radii_cf(ipsp,3)=zero
1264      else
1265        psps%gth_params%radii_cf(ipsp,3)=max( &
1266 &       min(dtset%wvl_crmult*psps%gth_params%radii_cf(ipsp,1),15._dp*maxrad)/dtset%wvl_frmult, &
1267 &       psps%gth_params%radii_cf(ipsp,2))
1268      end if
1269 #endif
1270    end if
1271 
1272    if (pspcod/=7.and.pspcod/=17) then
1273      write(message, '(a,f14.8,a,a)' ) '  pspatm : epsatm=',epsatm,ch10,&
1274 &     '         --- l  ekb(1:nproj) -->'
1275      call wrtout(ab_out,message,'COLL')
1276      call wrtout(std_out,  message,'COLL')
1277      iln0=0
1278      do ilmn=1,psps%lmnmax
1279        iln=indlmn(5,ilmn)
1280        if (iln>iln0) then
1281          il=indlmn(1,ilmn)
1282          if (indlmn(6,ilmn)==1) then
1283            iln0=iln0+nproj(il+1)
1284            !if (dtset%optdriver == RUNL_SIGMA) then
1285            !  do ii=0,nproj(il+1)-1
1286            !    ekb(iln+ii) = zero
1287            !  end do
1288            !end if
1289            write(message, '(13x,i1,4f12.6)' ) il,(ekb(iln+ii),ii=0,nproj(il+1)-1)
1290          else
1291            iln0=iln0+nproj(il+psps%mpsang)
1292            write(message, '(2x,a,i1,4f12.6)' ) 'spin-orbit ',il,(ekb(iln+ii),ii=0,nproj(il+psps%mpsang)-1)
1293          end if
1294          call wrtout(ab_out,message,'COLL')
1295          call wrtout(std_out,message,'COLL')
1296        end if
1297      end do
1298    end if
1299 
1300    ! NC: Evalute spline-fit of the model core charge in reciprocal space.
1301    ! TODO: Be careful, because we will be using the PAW part in which tcore is always avaiable!
1302    ! Should add a test with 2 NC pseudos: one with NLCC and the other without!
1303    if (psps%usepaw == 0) then
1304      call nctab_eval_tcorespl(nctab, psps%n1xccc, xcccrc, xccc1d, psps%mqgrid_vl, psps%qgrid_vl)
1305    end if
1306 
1307    write(message,'(3a)') ' pspatm: atomic psp has been read ',' and splines computed',ch10
1308    call wrtout(ab_out,message,'COLL')
1309    call wrtout(std_out,message,'COLL')
1310 
1311    ABI_DEALLOCATE(e990)
1312    ABI_DEALLOCATE(e999)
1313    ABI_DEALLOCATE(rcpsp)
1314    ABI_DEALLOCATE(rms)
1315    ABI_DEALLOCATE(ekb1)
1316    ABI_DEALLOCATE(ekb2)
1317    ABI_DEALLOCATE(epspsp)
1318    ABI_DEALLOCATE(nproj)
1319 
1320    if (dtset%prtvol > 9 .and. psps%usepaw==0 .and. psps%lmnmax>3) then
1321 
1322      write (filnam, '(a,i0,a)') trim(dtfil%fnameabo_pspdata), ipsp, ".dat"
1323      if (open_file(filnam, message, newunit=unt) /= 0) then
1324        MSG_ERROR(message)
1325      end if
1326      write (unt,*) '# Pseudopotential data in reciprocal space as used by ABINIT'
1327      write (unt,'(a)', ADVANCE='NO') '# index       vlocal   '
1328      if (psps%lnmax > 0) &
1329 &     write (unt,'(a,I3)', ADVANCE='NO')   '           1st proj(l=', indlmn(1,1)
1330      if (psps%lnmax > 1) &
1331 &     write (unt,'(a,I3)', ADVANCE='NO')   ')            2nd(l=', indlmn(1,2)
1332      if (psps%lnmax > 2) &
1333 &     write (unt,'(a,I3,a)', ADVANCE='NO') ')            3rd(l=', indlmn(1,3), ')'
1334      write (unt,*)
1335 
1336      do ii = 1, psps%mqgrid_vl
1337        write(unt, '(I5,E24.16)', ADVANCE='NO') ii, vlspl(ii,1)
1338        if (psps%lnmax > 0) write(unt, '(E24.16)', ADVANCE='NO') ffspl(ii,1,1)
1339        if (psps%lnmax > 1) write(unt, '(E24.16)', ADVANCE='NO') ffspl(ii,1,2)
1340        if (psps%lnmax > 2) write(unt, '(E24.16)', ADVANCE='NO') ffspl(ii,1,3)
1341        write(unt, *)
1342      end do
1343      close(unt)
1344 
1345      write (filnam, '(a,i0,a)') trim(dtfil%fnameabo_nlcc_derivs), ipsp, ".dat"
1346      if (open_file(filnam, message, newunit=unt) /= 0) then
1347        MSG_ERROR(message)
1348      end if
1349      write (unt,*) '# Non-linear core corrections'
1350      write (unt,*) '#  r, pseudocharge, 1st, 2nd, 3rd, 4th, 5th derivatives'
1351      do ii = 1, psps%n1xccc
1352        write (unt,*) xcccrc*(ii-1)/(psps%n1xccc-1), xccc1d(ii,1), xccc1d(ii,2), &
1353 &       xccc1d(ii,3), xccc1d(ii,4), &
1354 &       xccc1d(ii,5), xccc1d(ii,6)
1355      end do
1356      close(unt)
1357    end if
1358 
1359  end if !me=0
1360 
1361  if (paral_mode==1) then
1362    call timab(48,1,tsec)
1363    call pawpsp_bcast(comm_mpi,epsatm,ffspl,pawrad,pawtab,vlspl,xcccrc)
1364    call timab(48,2,tsec)
1365  end if
1366 
1367  if (psps%usepaw==1) then
1368    indlmn(:,:)=0
1369    indlmn(1:6,1:pawtab%lmn_size)=pawtab%indlmn(1:6,1:pawtab%lmn_size)
1370  end if
1371 
1372 !--------------------------------------------------------------------
1373 !WVL+PAW:
1374  if (dtset%usepaw==1 .and. (dtset%icoulomb /= 0 .or. dtset%usewvl==1)) then
1375 #if defined HAVE_BIGDFT
1376    psps%gth_params%psppar(:,:,ipsp) = UNINITIALIZED(1._dp)
1377    psps%gth_params%radii_cf(ipsp,:) = UNINITIALIZED(1._dp)
1378    call wvl_descr_psp_fill(psps%gth_params, ipsp, psps%pspxc(1), int(psps%zionpsp(ipsp)), int(psps%znuclpsp(ipsp)), 0)
1379 #endif
1380 
1381 !  The following lines are added to keep backward compatibilty
1382    maxrad=zero
1383 #if defined HAVE_BIGDFT
1384    do ii=1,size(psps%gth_params%psppar,1)-1 ! psppar first dim begins at 0
1385      if (psps%gth_params%psppar(ii,0,ipsp)/=zero) maxrad=max(maxrad,psps%gth_params%psppar(ii,0,ipsp))
1386    end do
1387    if (abs(maxrad)<=tol12) then
1388 !== MT COMMENT
1389 !    Damien wants to activate this (in order to directly compare to bigDFT):
1390      psps%gth_params%radii_cf(ipsp,3)= psps%gth_params%radii_cf(ipsp,2)
1391 !    But, this changes strongly file references.
1392 !    So, I keep this, waiting for Tonatiuh s validation
1393      psps%gth_params%radii_cf(ipsp,3)= &
1394 &     (psps%gth_params%radii_cf(ipsp,1)+psps%gth_params%radii_cf(ipsp,2))*half
1395 !== MT COMMENT
1396    else
1397      psps%gth_params%radii_cf(ipsp,3)=max( &
1398 &     min(dtset%wvl_crmult*psps%gth_params%radii_cf(ipsp,1),15._dp*maxrad)/dtset%wvl_frmult, &
1399 &     psps%gth_params%radii_cf(ipsp,2))
1400    end if
1401    if(present(comm_mpi)) then
1402      call pawpsp_wvl(psps%filpsp(ipsp),pawrad,pawtab,dtset%usewvl,dtset%wvl_ngauss,comm_mpi)
1403    else
1404      call pawpsp_wvl(psps%filpsp(ipsp),pawrad,pawtab,dtset%usewvl,dtset%wvl_ngauss)
1405    end if
1406 #endif
1407  end if
1408 
1409 !end of WVL+PAW section
1410 !----------------------------------------------------
1411 
1412  return
1413 
1414  ! Handle IO error
1415  10 continue
1416  MSG_ERROR(errmsg)
1417 
1418 end subroutine pspatm

ABINIT/pspcor [ Functions ]

[ Top ] [ 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

PARENTS

      pspini

CHILDREN

SOURCE

682 subroutine pspcor(ecore,epsatm,natom,ntypat,typat,zion)
683 
684 
685 !This section has been created automatically by the script Abilint (TD).
686 !Do not modify the following lines by hand.
687 #undef ABI_FUNC
688 #define ABI_FUNC 'pspcor'
689 !End of the abilint section
690 
691  implicit none
692 
693 !Arguments ------------------------------------
694 !scalars
695  integer,intent(in) :: natom,ntypat
696  real(dp),intent(out) :: ecore
697 !arrays
698  integer,intent(in) :: typat(natom)
699  real(dp),intent(in) :: epsatm(ntypat),zion(ntypat)
700 
701 !Local variables-------------------------------
702 !scalars
703  integer :: ia
704  real(dp) :: charge,esum
705 
706 ! *************************************************************************
707 
708  charge = 0.d0
709  esum = 0.d0
710  do ia=1,natom
711 !  compute pseudocharge:
712    charge=charge+zion(typat(ia))
713 !  add pseudocore energies together:
714    esum = esum + epsatm(typat(ia))
715  end do
716 
717  ecore=charge*esum
718 
719 end subroutine pspcor

ABINIT/pspini [ Functions ]

[ Top ] [ 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-2018 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.

PARENTS

      bethe_salpeter,eph,gstate,nonlinear,respfn,screening,sigma,wfk_analyze

CHILDREN

      nctab_free,nctab_init,nctab_mixalch,pawtab_set_flags,pspatm,pspcor
      psps_ncwrite,psps_print,timab,wrtout,xmpi_sum

SOURCE

135 subroutine pspini(dtset,dtfil,ecore,gencond,gsqcut,gsqcutdg,pawrad,pawtab,psps,rprimd,comm_mpi)
136 
137 
138 !This section has been created automatically by the script Abilint (TD).
139 !Do not modify the following lines by hand.
140 #undef ABI_FUNC
141 #define ABI_FUNC 'pspini'
142 !End of the abilint section
143 
144  implicit none
145 
146 !Arguments ------------------------------------
147 !scalars
148  integer, optional,intent(in) :: comm_mpi
149  integer,intent(out) :: gencond
150  real(dp),intent(in) :: gsqcut,gsqcutdg
151  real(dp),intent(out) :: ecore
152  type(dataset_type),intent(in) :: dtset
153  type(datafiles_type),intent(in) :: dtfil
154 !arrays
155  real(dp),intent(in) :: rprimd(3,3)
156 !no_abirules
157  type(pseudopotential_type), target,intent(inout) :: psps
158  type(pawrad_type), intent(inout) :: pawrad(psps%ntypat*psps%usepaw)
159  type(pawtab_type), intent(inout) :: pawtab(psps%ntypat*psps%usepaw)
160 
161 !Local variables-------------------------------
162 !scalars
163  integer,parameter :: npspmax=50
164  integer,save :: dimekb_old=0,ifirst=1,ixc_old=-1,lmnmax_old=0,lnmax_old=0
165  integer,save :: mpssoang_old=0,mqgridff_old=0,mqgridvl_old=0,optnlxccc_old=-1
166  integer,save :: paw_size_old=-1,pawxcdev_old=-1,positron_old=-2,usepaw_old=-1
167  integer,save :: usexcnhat_old=-1,usewvl_old=-1,useylm_old=-1
168  integer :: comm_mpi_,ierr,ii,ilang,ilmn,ilmn0,iproj,ipsp,ipspalch
169  integer :: ispin,itypalch,itypat,mtypalch,npsp,npspalch,ntypalch
170  integer :: ntypat,ntyppure,paw_size
171  logical :: has_kij,has_tproj,has_tvale,has_nabla,has_shapefncg,has_vminushalf,has_wvl
172  real(dp),save :: ecore_old=zero,gsqcut_old=zero,gsqcutdg_old=zero
173  real(dp) :: dq,epsatm_psp,qmax,rmax,xcccrc
174  character(len=500) :: message
175  type(pawrad_type) :: pawrad_dum
176  type(pawtab_type) :: pawtab_dum
177  type(nctab_t) :: nctab_dum
178  type(nctab_t),pointer :: nctab_ptr
179 !arrays
180  integer :: paw_options(9)
181  integer,save :: paw_options_old(9)=(/-1,-1,-1,-1,-1,-1,-1,-1,-1/)
182  integer,save :: pspso_old(npspmax),pspso_zero(npspmax)
183  integer,allocatable :: indlmn_alch(:,:,:),new_pspso(:)
184  integer,pointer :: indlmn(:,:)
185  real(dp),save :: epsatm(npspmax)
186  real(dp) :: tsec(2)
187  real(dp),allocatable :: dvlspl(:,:),dvlspl_alch(:,:,:),ekb(:),ekb_alch(:,:)
188  real(dp),allocatable :: epsatm_alch(:),ffspl(:,:,:),ffspl_alch(:,:,:,:)
189  real(dp),allocatable :: vlspl(:,:),vlspl_alch(:,:,:),xccc1d(:,:)
190  real(dp),allocatable :: xccc1d_alch(:,:,:),xcccrc_alch(:)
191  type(nctab_t),target,allocatable :: nctab_alch(:)
192 
193 ! *************************************************************************
194 
195  DBG_ENTER("COLL")
196 
197 !Keep track of time spent in this subroutine
198  call timab(15,1,tsec)
199 
200 !-------------------------------------------------------------
201 ! Some initializations
202 !-------------------------------------------------------------
203 
204 !Useful sizes
205  ntypat=psps%ntypat
206  mtypalch=psps%mtypalch
207  npsp=psps%npsp
208  if (npsp>npspmax) then
209    MSG_BUG("npsp>npspmax in pspini !")
210  end if
211 
212 !Size of grids for atomic data represented in reciprocal space
213 
214 !Set up q grids, make qmax 20% larger than largest expected:
215  qmax=1.2d0 * sqrt(gsqcut)
216 !ffnl is always computed in reciprocal space
217  dq=qmax/(one*(psps%mqgrid_ff-1))
218  do ii=1,psps%mqgrid_ff
219    psps%qgrid_ff(ii)=(ii-1)*dq
220  end do
221  if (psps%usepaw==1) qmax=1.2d0 * sqrt(gsqcutdg)
222 !If vlspl is computed in real space, qgrid contains a real space mesh
223 !the max is taken as the biggest distance in the box.
224  if (psps%vlspl_recipSpace) then
225    dq=qmax/(one*(psps%mqgrid_vl-1))
226  else
227    rmax = (rprimd(1, 1) + rprimd(1, 2) + rprimd(1, 3)) ** 2
228    rmax = rmax + (rprimd(2, 1) + rprimd(2, 2) + rprimd(2, 3)) ** 2
229    rmax = rmax + (rprimd(3, 1) + rprimd(3, 2) + rprimd(3, 3)) ** 2
230    rmax = sqrt(rmax)
231    dq = rmax /(one*(psps%mqgrid_vl-1))
232  end if
233  do ii=1,psps%mqgrid_vl
234    psps%qgrid_vl(ii)=(ii-1)*dq
235  end do
236 
237 !Determine whether new optional data requests have changed
238  paw_options=0;paw_size=0
239  if (psps%usepaw==1) then
240    paw_size=size(pawtab)
241    has_kij=(dtset%positron/=0)
242    has_tvale=.true. ! Will be modified later (depending on PAW dataset format)
243    has_nabla=.false.
244    has_shapefncg=(dtset%optdriver==RUNL_GSTATE.and.((dtset%iprcel>=20.and.dtset%iprcel<70).or.dtset%iprcel>=80))
245    has_wvl=(dtset%usewvl==1.or.dtset%icoulomb/=0)
246    has_tproj=(dtset%usewvl==1) ! projectors will be free at the end of the psp reading
247    has_vminushalf=(maxval(dtset%ldaminushalf)==1)
248    if (has_kij)       paw_options(1)=1
249    if (has_tvale)     paw_options(2)=1
250    if (has_nabla)     paw_options(5)=1
251    if (has_shapefncg) paw_options(6)=1
252    if (has_wvl)       paw_options(7)=1
253    if (has_tproj)     paw_options(8)=1
254    if (has_vminushalf)paw_options(9)=1
255    !if (dtset%prtvclmb /= 0) then
256    paw_options(3) = 1
257    paw_options(4) = 1
258    !end if
259  end if
260 
261 !Determine whether the spin-orbit characteristic has changed
262 !Do not forget that the SO is not consistent with alchemy presently
263  ABI_ALLOCATE(new_pspso,(npsp))
264  if (ifirst==1) pspso_old(:)=-1
265  if (ifirst==1) pspso_zero(:)=-1
266  do ipsp=1,npsp
267    new_pspso(ipsp)=1
268 !  No new characteristics if it is equal to the old one,
269 !  or, if it is one, the old one is equal to the intrinsic characteristic one.
270    if (psps%pspso(ipsp)==pspso_old(ipsp).or. &
271 &   (psps%pspso(ipsp)==1.and.pspso_old(ipsp)==pspso_zero(ipsp))) then
272      new_pspso(ipsp)=0
273    end if
274 !  No new characteristics if PAW
275    if (psps%usepaw==1) new_pspso(ipsp)=0
276 !  Prepare the saving of the intrinsic pseudopotential characteristics
277    if(psps%pspso(ipsp)==1) pspso_zero(ipsp)=0
278  end do
279 
280 !Compute the general condition for new computation of pseudopotentials
281  gencond=0
282  if(   ixc_old /= dtset%ixc                &
283 & .or. mqgridff_old /= psps%mqgrid_ff      &
284 & .or. mqgridvl_old /= psps%mqgrid_vl      &
285 & .or. mpssoang_old /= psps%mpssoang       &
286 & .or. abs(gsqcut_old-gsqcut)>1.0d-10      &
287 & .or. (psps%usepaw==1.and.abs(gsqcutdg_old-gsqcutdg)>1.0d-10) &
288 & .or. dimekb_old /= psps%dimekb           &
289 & .or. lmnmax_old /= psps%lmnmax           &
290 & .or. lnmax_old  /= psps%lnmax            &
291 & .or. optnlxccc_old /= psps%optnlxccc     &
292 & .or. usepaw_old /= psps%usepaw           &
293 & .or. useylm_old /= psps%useylm           &
294 & .or. pawxcdev_old /= dtset%pawxcdev      &
295 & .or. positron_old /= dtset%positron      &
296 & .or. usewvl_old /= dtset%usewvl          &
297 & .or. paw_size_old /= paw_size            &
298 & .or. usexcnhat_old/=dtset%usexcnhat_orig &
299 & .or. any(paw_options_old(:)/=paw_options(:)) &
300 & .or. sum(new_pspso(:))/=0                &
301 & .or. mtypalch>0                          &
302 & .or. (dtset%usewvl==1.and.psps%usepaw==1)&
303 & ) gencond=1
304 
305  if (present(comm_mpi).and.psps%usepaw==1) then
306    if(xmpi_comm_size(comm_mpi)>1)then
307      call xmpi_sum(gencond,comm_mpi,ierr)
308    end if
309    if (gencond/=0) gencond=1
310  end if
311  ABI_DEALLOCATE(new_pspso)
312 
313 !-------------------------------------------------------------
314 ! Following section is only reached when new computation
315 ! of pseudopotentials is needed
316 !-------------------------------------------------------------
317 
318  if (gencond==1) then
319 
320    write(message, '(a,a)' ) ch10,&
321 &   '--- Pseudopotential description ------------------------------------------------'
322    call wrtout(ab_out,message,'COLL')
323 
324    ABI_ALLOCATE(ekb,(psps%dimekb*(1-psps%usepaw)))
325    ABI_ALLOCATE(xccc1d,(psps%n1xccc*(1-psps%usepaw),6))
326    ABI_ALLOCATE(ffspl,(psps%mqgrid_ff,2,psps%lnmax))
327    ABI_ALLOCATE(vlspl,(psps%mqgrid_vl,2))
328    if (.not.psps%vlspl_recipSpace) then
329      ABI_ALLOCATE(dvlspl,(psps%mqgrid_vl,2))
330    else
331      ABI_ALLOCATE(dvlspl,(0,0))
332    end if
333 
334 !  PAW: reset flags for optional data
335    if (psps%usepaw==1) then
336      call pawtab_set_flags(pawtab,has_kij=paw_options(1),has_tvale=paw_options(2),&
337 &     has_vhnzc=paw_options(3),has_vhtnzc=paw_options(4),&
338 &     has_nabla=paw_options(5),has_shapefncg=paw_options(6),&
339 &     has_wvl=paw_options(7),has_tproj=paw_options(8))
340 ! the following have to be included in pawtab_set_flags
341      do ipsp=1,psps%ntypat
342        pawtab(ipsp)%has_vminushalf=dtset%ldaminushalf(ipsp)
343      end do
344    end if
345 
346 !  Read atomic pseudopotential data and get transforms
347 !  for each atom type : two cases, alchemy or not.
348 
349 !  No alchemical pseudoatom, in all datasets, npsp=ntypat
350    if(mtypalch==0)then
351 
352      do ipsp=1,npsp
353 
354        xcccrc=zero
355        ekb(:)=zero;ffspl(:,:,:)=zero;vlspl(:,:)=zero
356        if (.not.psps%vlspl_recipSpace) dvlspl(:, :)=zero
357        if (psps%usepaw==0) xccc1d(:,:)=zero
358        indlmn=>psps%indlmn(:,:,ipsp)
359        indlmn(:,:)=0
360 
361        write(message, '(a,i4,a,t38,a)' ) &
362 &       '- pspini: atom type',ipsp,'  psp file is',trim(psps%filpsp(ipsp))
363        call wrtout(ab_out,message,'COLL')
364        call wrtout(std_out,message,'COLL')
365 
366 !      Read atomic psp V(r) and wf(r) to get local and nonlocal psp:
367 !      Cannot use the same call in case of bound checking, because of pawrad/pawtab
368        if(psps%usepaw==0)then
369          call pspatm(dq,dtset,dtfil,ekb,epsatm(ipsp),ffspl,indlmn,ipsp,&
370 &         pawrad_dum,pawtab_dum,psps,vlspl,dvlspl,xcccrc,xccc1d,psps%nctab(ipsp))
371          psps%ekb(:,ipsp)=ekb(:)
372          psps%xccc1d(:,:,ipsp)=xccc1d(:,:)
373        else
374          comm_mpi_=xmpi_comm_self;if (present(comm_mpi)) comm_mpi_=comm_mpi
375          call pspatm(dq,dtset,dtfil,ekb,epsatm(ipsp),ffspl,indlmn,ipsp,&
376 &         pawrad(ipsp),pawtab(ipsp),psps,vlspl,dvlspl,xcccrc,xccc1d,nctab_dum,&
377 &         comm_mpi=comm_mpi_)
378        end if
379 
380        ! Copy data to psps datastructure.
381        psps%xcccrc(ipsp)=xcccrc
382        psps%znucltypat(ipsp)=psps%znuclpsp(ipsp)
383        psps%ffspl(:,:,:,ipsp)=ffspl(:,:,:)
384        psps%vlspl(:,:,ipsp)=vlspl(:,:)
385        if (.not.psps%vlspl_recipSpace) psps%dvlspl(:, :, ipsp) = dvlspl(:, :)
386      end do ! ipsp
387 
388    else ! if mtypalch/=0
389 
390      npspalch=psps%npspalch
391      ntyppure=npsp-npspalch
392      ntypalch=psps%ntypalch
393      ABI_ALLOCATE(epsatm_alch,(npspalch))
394      ABI_ALLOCATE(ekb_alch,(psps%dimekb,npspalch*(1-psps%usepaw)))
395      ABI_ALLOCATE(ffspl_alch,(psps%mqgrid_ff,2,psps%lnmax,npspalch))
396      ABI_ALLOCATE(xccc1d_alch,(psps%n1xccc*(1-psps%usepaw),6,npspalch))
397      ABI_ALLOCATE(xcccrc_alch,(npspalch))
398      ABI_ALLOCATE(vlspl_alch,(psps%mqgrid_vl,2,npspalch))
399      if (.not.psps%vlspl_recipSpace) then
400        ABI_ALLOCATE(dvlspl_alch,(psps%mqgrid_vl,2,npspalch))
401      end if
402      ABI_ALLOCATE(indlmn,(6,psps%lmnmax))
403      ABI_ALLOCATE(indlmn_alch,(6,psps%lmnmax,npspalch))
404 
405      ! Allocate NC tables used for mixing.
406      if (psps%usepaw == 0) then
407        ABI_DT_MALLOC(nctab_alch, (npspalch))
408        do ipspalch=1,npspalch
409          call nctab_init(nctab_alch(ipspalch), psps%mqgrid_vl, .False., .False.)
410        end do
411      end if
412 
413      do ipsp=1,npsp
414        write(message, '(a,i4,a,t38,a)' ) &
415 &       '- pspini: atom type',ipsp,'  psp file is',trim(psps%filpsp(ipsp))
416        call wrtout(ab_out,message,'COLL')
417 
418        xcccrc=zero
419        ekb(:)=zero;ffspl(:,:,:)=zero;vlspl(:,:)=zero
420        if (.not.psps%vlspl_recipSpace) dvlspl(:, :)=zero
421        if (psps%usepaw==0) xccc1d(:,:)=zero
422        indlmn(:,:)=0
423 
424 !      Read atomic psp V(r) and wf(r) to get local and nonlocal psp:
425        if (psps%usepaw==0) then
426          if (ipsp <= ntyppure) then
427            ! Store data in nctab if pure atom.
428            nctab_ptr => psps%nctab(ipsp)
429          else
430            ! Store data in nctab_alch (to be mixed afterwards).
431            nctab_ptr => nctab_alch(ipsp-ntyppure)
432          end if
433 
434          call pspatm(dq,dtset,dtfil,ekb,epsatm_psp,ffspl,indlmn,ipsp,&
435 &         pawrad_dum,pawtab_dum,psps,vlspl,dvlspl,xcccrc,xccc1d,nctab_ptr)
436 
437        else if (psps%usepaw==1) then
438          comm_mpi_=xmpi_comm_self;if (present(comm_mpi)) comm_mpi_=comm_mpi
439          call pspatm(dq,dtset,dtfil,ekb,epsatm_psp,ffspl,indlmn,ipsp,&
440 &         pawrad(ipsp),pawtab(ipsp),psps,vlspl,dvlspl,xcccrc,xccc1d,nctab_dum,&
441 &         comm_mpi=comm_mpi_)
442        end if
443 
444        if (ipsp<=ntyppure) then
445 !        Pure pseudopotentials, leading to pure pseudoatoms
446          epsatm(ipsp)=epsatm_psp
447          psps%znucltypat(ipsp)=psps%znuclpsp(ipsp)
448          if (psps%usepaw==0) psps%ekb(:,ipsp)=ekb(:)
449          psps%ffspl(:,:,:,ipsp)=ffspl(:,:,:)
450          psps%vlspl(:,:,ipsp)=vlspl(:,:)
451          if (.not.psps%vlspl_recipSpace) psps%dvlspl(:, :, ipsp)=dvlspl(:, :)
452          if (psps%usepaw==0) psps%xccc1d(:,:,ipsp)=xccc1d(:,:)
453          psps%xcccrc(ipsp)=xcccrc
454          psps%indlmn(:,:,ipsp)=indlmn(:,:)
455 
456        else
457 !        Pseudopotentials for alchemical generation
458          ipspalch=ipsp-ntyppure
459          epsatm_alch(ipspalch)=epsatm_psp
460          ffspl_alch(:,:,:,ipspalch)=ffspl(:,:,:)
461          vlspl_alch(:,:,ipspalch)=vlspl(:,:)
462          if (.not.psps%vlspl_recipSpace) dvlspl_alch(:,:,ipspalch)=dvlspl(:,:)
463          if (psps%usepaw==0) then
464            ekb_alch(:,ipspalch)=ekb(:)
465            xccc1d_alch(:,:,ipspalch)=xccc1d(:,:)
466          end if
467          xcccrc_alch(ipspalch)=xcccrc
468          indlmn_alch(:,:,ipspalch)=indlmn(:,:)
469 !        write(std_out,'(a,6i4)' )' pspini : indlmn_alch(:,1,ipspalch)=',indlmn_alch(:,1,ipspalch)
470 !        write(std_out,'(a,6i4)' )' pspini : indlmn_alch(:,2,ipspalch)=',indlmn_alch(:,2,ipspalch)
471        end if
472 
473      end do ! ipsp
474 
475      ! Generate data for alchemical pseudos.
476      do itypalch=1,ntypalch
477        itypat=itypalch+ntyppure
478        psps%znucltypat(itypat)=200.0+itypalch    ! Convention for alchemical pseudoatoms
479        vlspl(:,:)=zero
480        if (.not.psps%vlspl_recipSpace) dvlspl(:, :) = zero
481        epsatm(itypat)=zero
482        xcccrc=zero
483        if (psps%usepaw==0) xccc1d(:,:)=zero
484 
485 !      Here, linear combination of the quantities
486 !      MG: FIXME I think that the mixing of xcccrc is wrong when the xxccrc are different!
487 !      but this is minor bug since alchemical pseudos should not have XCCC (?)
488        do ipspalch=1,npspalch
489          epsatm(itypat) = epsatm(itypat) + epsatm_alch(ipspalch) * psps%mixalch(ipspalch,itypalch)
490          vlspl(:,:) = vlspl(:,:) + vlspl_alch(:,:,ipspalch) * psps%mixalch(ipspalch,itypalch)
491          if (.not.psps%vlspl_recipSpace) then
492            dvlspl(:,:) = dvlspl(:,:) + dvlspl_alch(:,:,ipspalch) * psps%mixalch(ipspalch,itypalch)
493          end if
494          xcccrc = xcccrc + xcccrc_alch(ipspalch) * psps%mixalch(ipspalch,itypalch)
495          if (psps%usepaw==0) then
496            xccc1d(:,:) = xccc1d(:,:) + xccc1d_alch(:,:,ipspalch) * psps%mixalch(ipspalch,itypalch)
497          end if
498        end do ! ipspalch
499 
500        psps%vlspl(:,:,itypat)=vlspl(:,:)
501        if (.not.psps%vlspl_recipSpace) psps%dvlspl(:, :, itypat) = dvlspl(:, :)
502        if (psps%usepaw==0) then
503          psps%xccc1d(:,:,itypat)=xccc1d(:,:)
504        end if
505        psps%xcccrc(itypat)=xcccrc
506 
507        if (abs(xcccrc) > tol6) then
508          write(std_out, *)"xcccrc", xcccrc
509          MSG_WARNING("Alchemical pseudopotential with nlcc!")
510        end if
511 
512 !      Combine the different non-local projectors : for the scalar part then
513 !      the spin-orbit part, treat the different angular momenta
514 !      WARNING : this coding does not work for PAW
515        ilmn=0
516        psps%indlmn(:,:,itypat)=0
517        do ispin=1,2
518          do ilang=0,3
519            if(ispin==2 .and. ilang==0)cycle
520            iproj=0
521            do ipspalch=1,npspalch
522              if(abs(psps%mixalch(ipspalch,itypalch))>tol10)then
523                do ilmn0=1,psps%lmnmax
524                  if(indlmn_alch(5,ilmn0,ipspalch)/=0)then
525                    if(indlmn_alch(6,ilmn0,ipspalch)==ispin)then
526                      if(indlmn_alch(1,ilmn0,ipspalch)==ilang)then
527                        ilmn=ilmn+1         ! increment the counter
528                        iproj=iproj+1       ! increment the counter, this does not work for PAW
529                        if(ilmn>psps%lmnmax)then
530                          MSG_BUG('Problem with the alchemical pseudopotentials : ilmn>lmnmax.')
531                        end if
532                        psps%indlmn(1,ilmn,itypat)=ilang
533                        psps%indlmn(2,ilmn,itypat)=indlmn_alch(2,ilmn0,ipspalch)
534                        psps%indlmn(3,ilmn,itypat)=iproj                       ! This does not work for PAW
535                        psps%indlmn(4,ilmn,itypat)=ilmn                        ! This does not work for PAW
536                        psps%indlmn(5,ilmn,itypat)=ilmn
537                        psps%indlmn(6,ilmn,itypat)=ispin
538                        ! The two lines below do not work for PAW
539                        if (psps%usepaw==0) then
540                          psps%ekb(ilmn,itypat)=psps%mixalch(ipspalch,itypalch) *ekb_alch(ilmn0,ipspalch)
541                        end if
542                        psps%ffspl(:,:,ilmn,itypat)=ffspl_alch(:,:,ilmn0,ipspalch)
543                      end if ! ilang is OK
544                    end if ! ispin is OK
545                  end if ! ilmn0 exist
546                end do ! ilmn0
547              end if ! mixalch>tol10
548            end do ! ipspalch
549          end do ! ilang
550        end do ! ispin
551 
552      end do ! itypalch
553 
554      ABI_DEALLOCATE(epsatm_alch)
555      ABI_DEALLOCATE(ekb_alch)
556      ABI_DEALLOCATE(ffspl_alch)
557      ABI_DEALLOCATE(xccc1d_alch)
558      ABI_DEALLOCATE(xcccrc_alch)
559      ABI_DEALLOCATE(vlspl_alch)
560      if (.not.psps%vlspl_recipSpace) then
561        ABI_DEALLOCATE(dvlspl_alch)
562      end if
563      ABI_DEALLOCATE(indlmn_alch)
564      ABI_DEALLOCATE(indlmn)
565 
566      ! Mix NC tables.
567      if (psps%usepaw == 0) then
568        call nctab_mixalch(nctab_alch, npspalch, ntypalch, psps%algalch, psps%mixalch, psps%nctab(ntyppure+1:))
569        do ipspalch=1,npspalch
570          call nctab_free(nctab_alch(ipspalch))
571        end do
572        ABI_DT_FREE(nctab_alch)
573      end if
574    end if ! mtypalch
575 
576    ABI_DEALLOCATE(ekb)
577    ABI_DEALLOCATE(ffspl)
578    ABI_DEALLOCATE(vlspl)
579    ABI_DEALLOCATE(xccc1d)
580 
581    if (.not.psps%vlspl_recipSpace) then
582      ABI_DEALLOCATE(dvlspl)
583    end if
584 
585  end if !  End condition of new computation needed
586 
587 !-------------------------------------------------------------
588 ! Following section is always executed
589 !-------------------------------------------------------------
590 !One should move this section of code outside of pspini,
591 !but epsatm is needed, so should be in the psp datastructure.
592 !Compute pseudo correction energy. Will differ from an already
593 !computed one if the number of atom differ ...
594  call pspcor(ecore,epsatm,dtset%natom,ntypat,dtset%typat,psps%ziontypat)
595  if(abs(ecore_old-ecore)>tol8*abs(ecore_old+ecore))then
596    write(message, '(2x,es15.8,t50,a)' ) ecore,'ecore*ucvol(ha*bohr**3)'
597 !  ecore is useless if iscf<=0, but at least it has been initialized
598    if(dtset%iscf>=0)then
599      call wrtout(ab_out,message,'COLL')
600    end if
601    call wrtout(std_out,message,'COLL')
602  end if
603 
604 !End of pseudopotential output section
605  write(message, '(2a)' )&
606 & '--------------------------------------------------------------------------------',ch10
607  call wrtout(ab_out,message,'COLL')
608 
609 !-------------------------------------------------------------
610 ! Keep track of this call to the routine
611 !-------------------------------------------------------------
612 
613  if (ifirst==1) ifirst=0
614 
615  mqgridff_old=psps%mqgrid_ff
616  mqgridvl_old=psps%mqgrid_vl
617  mpssoang_old=psps%mpssoang
618  ixc_old=dtset%ixc
619  gsqcut_old=gsqcut;if (psps%usepaw==1) gsqcutdg_old=gsqcutdg
620  lmnmax_old=psps%lmnmax
621  lnmax_old=psps%lnmax
622  optnlxccc_old=psps%optnlxccc
623  usepaw_old=psps%usepaw
624  dimekb_old=psps%dimekb
625  useylm_old=psps%useylm
626  pawxcdev_old=dtset%pawxcdev
627  positron_old=dtset%positron
628  usewvl_old = dtset%usewvl
629  usexcnhat_old=dtset%usexcnhat_orig
630  paw_size_old=paw_size
631  ecore_old=ecore
632  paw_options_old(:)=paw_options(:)
633 
634  do ipsp=1,npsp
635    pspso_old(ipsp)=psps%pspso(ipsp)
636    if(pspso_zero(ipsp)==0)pspso_zero(ipsp)=psps%pspso(ipsp)
637  end do
638  psps%mproj = maxval(psps%indlmn(3,:,:))
639 
640  if (gencond == 1) call psps_print(psps,std_out,dtset%prtvol)
641 
642  ! Write the PSPS.nc file and exit here if requested by the user.
643  if (abs(dtset%prtpsps) == 1) then
644    if (xmpi_comm_rank(xmpi_world) == 0) call psps_ncwrite(psps, trim(dtfil%filnam_ds(4))//"_PSPS.nc")
645    if (dtset%prtpsps == -1) then
646      MSG_ERROR_NODUMP("prtpsps == -1 ==> aborting now")
647    end if
648  end if
649 
650  call timab(15,2,tsec)
651 
652  DBG_EXIT("COLL")
653 
654 end subroutine pspini