ps_hgh.F90

Go to the documentation of this file.

!! Copyright (C) 2002-2006 M. Marques, A. Castro, A. Rubio, G. Bertsch
!! Copyright (C) 2020 N. Tancogne-Dejean
!!
!! This program is free software; you can redistribute it and/or modify
!! it under the terms of the GNU General Public License as published by
!! the Free Software Foundation; either version 2, or (at your option)
!! any later version.
!!
!! This program is distributed in the hope that it will be useful,
!! but WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!! GNU General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with this program; if not, write to the Free Software
!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
!! 02110-1301, USA.
!!
 
#include "global.h"
 
module ps_hgh_oct_m
  
  use atomic_oct_m
  use debug_oct_m
  use global_oct_m
  use io_oct_m
  use, intrinsic :: iso_fortran_env
  use loct_math_oct_m
  use logrid_oct_m
  use messages_oct_m
  use namespace_oct_m
  use profiling_oct_m
 
  implicit none
 
  private
  public  ::      &
    ps_hgh_t,     &
    hgh_init,     &
    hgh_process,  &
    hgh_debug,    &
    hgh_end,      &
    hgh_get_eigen
 
  type ps_hgh_t
    ! Components are public by default
 
    character(len=256), private :: title
    integer,            private :: pspdat ! date of creation of PP (DDMMYY)
    integer,            private :: pspcod ! code for the pseudopotential (3 for .hgh)
    integer,            private :: pspxc  ! exchange-correlation used to generate the psp
    integer                     :: lmax   ! Maximum l to use
    integer,            private :: mmax   ! Maximum number of points in real space grid
    real(real64),       private :: r2well ! ??
 
    
    character(len=5), private :: atom_name
    integer          :: z_val
    real(real64)     :: rlocal
    real(real64), private   :: rc(0:3)
    real(real64), private   :: c(1:4)
    real(real64)     :: h(0:3, 1:3, 1:3)
    real(real64)     :: k(0:3, 1:3, 1:3)
 
    type(valconf_t)  :: conf
    integer          :: l_max
 
    real(real64), allocatable   :: vlocal(:)
    real(real64), allocatable   :: kb(:,:,:)
    real(real64), allocatable   :: kbr(:)
    real(real64), allocatable   :: rphi(:,:)
    real(real64), allocatable, private :: eigen(:)
 
    type(logrid_t) :: g
  end type ps_hgh_t
 
  real(real64), parameter :: eps = 1.0e-8_real64
 
contains
 
  ! ---------------------------------------------------------
  subroutine hgh_init(psp, filename, namespace)
    type(ps_hgh_t),    intent(inout) :: psp
    character(len=*),  intent(in)    :: filename
    type(namespace_t), intent(in)    :: namespace
 
    integer :: iunit, i, np
    real(real64) :: aa, bb
 
    push_sub(hgh_init)
 
    iunit = io_open(trim(filename), action='read', form='formatted', status='old')
    i = load_params(iunit, psp, namespace)
    if (i /= 0) then
      call messages_write('Error reading hgh file')
      call messages_fatal(namespace=namespace)
    end if
    call io_close(iunit)
 
    ! Finds out psp%l_max. The most special cases are H, He, Li_sc and Be_sc, where psp%l_max = -1.
    psp%l_max = 0
    do while(psp%rc(psp%l_max) > 0.01_real64)
      psp%l_max = psp%l_max + 1
      if (psp%l_max > 3) exit
    end do
    psp%l_max = psp%l_max - 1
 
    call logrid_find_parameters(namespace, psp%conf%z, aa, bb, np)
 
    call logrid_init(psp%g, logrid_psf, aa, bb, np)
 
    ! Allocation of stuff.
    safe_allocate(psp%vlocal(1:psp%g%nrval))
    psp%vlocal = m_zero
    if (psp%l_max >= 0) then
      safe_allocate(psp%kbr(0:psp%l_max))
      safe_allocate(psp%kb(1:psp%g%nrval, 0:psp%l_max, 1:3))
      psp%kbr = m_zero
      psp%kb = m_zero
    end if
 
    pop_sub(hgh_init)
  end subroutine hgh_init
 
 
  ! ---------------------------------------------------------
  subroutine hgh_end(psp)
    type(ps_hgh_t), intent(inout) :: psp
 
    push_sub(hgh_end)
 
    if (psp%l_max >= 0) then
      safe_deallocate_a(psp%kbr)
      safe_deallocate_a(psp%kb)
    end if
    safe_deallocate_a(psp%vlocal)
    safe_deallocate_a(psp%rphi)
    safe_deallocate_a(psp%eigen)
    call logrid_end(psp%g)
 
    pop_sub(hgh_end)
  end subroutine hgh_end
 
 
  ! ---------------------------------------------------------
  subroutine hgh_process(psp, namespace)
    type(ps_hgh_t),    intent(inout) :: psp
    type(namespace_t), intent(in)    :: namespace
 
    integer :: l, i, ierr
 
    push_sub(hgh_process)
 
    safe_allocate(psp%rphi(1:psp%g%nrval, 1:psp%conf%p))
    safe_allocate(psp%eigen(1:psp%conf%p))
    psp%rphi = m_zero
    psp%eigen = m_zero
 
 
    ! Fixes the local potential
    call vlocalr_scalar(psp%g%rofi, psp%g%nrval, psp, psp%vlocal)
 
    ! And the projectors
    do l = 0, psp%l_max
      do i = 1, 3
        call projectorr_scalar(psp%g%rofi, psp%g%nrval, psp, i, l, psp%kb(:, l, i))
      end do
    end do
 
    ! get the pseudoatomic eigenfunctions (WARNING: This is not correctly done yet: "some" wavefunctions
    ! are obtained, but not the real ones!!!
    call solve_schroedinger(psp, ierr, namespace)
    if (ierr /= 0) then ! If the wavefunctions could not be found, we set its number to zero.
      write(message(1),'(a)') 'The algorithm that calculates atomic wavefunctions could not'
      write(message(2),'(a)') 'do its job. The program will continue, but expect poor'
      write(message(3),'(a)') 'convergence properties.'
      call messages_warning(3, namespace=namespace)
      psp%conf%p = 0
    end if
 
    ! Define the KB-projector cut-off radii
    call get_cutoff_radii(psp)
 
    pop_sub(hgh_process)
  end subroutine hgh_process
 
  ! ---------------------------------------------------------
  subroutine hgh_get_eigen(psp, eigen)
    type(ps_hgh_t), intent(in)  :: psp
    real(real64),   intent(out) :: eigen(:,:)
 
    integer :: i
 
    push_sub(hgh_get_eigen)
 
    do i = 1, psp%conf%p
      eigen(i, :) = psp%eigen(i)
    end do
 
    pop_sub(hgh_get_eigen)
  end subroutine hgh_get_eigen
 
  ! ---------------------------------------------------------
  function load_params(unit, params, namespace)
    integer,             intent(in)  :: unit        ! where to read from
    type(ps_hgh_t),      intent(out) :: params      ! obvious
    type(namespace_t),   intent(in)  :: namespace
 
    integer                     :: load_params ! 0 if success,
    ! 1 otherwise.
 
    integer :: i, iostat, j, k, nn, nnonloc, lloc
    character(len=VALCONF_STRING_LENGTH) :: line
 
    push_sub(load_params)
 
    ! Set initially everything to zero.
    params%c(1:4) = m_zero
    params%rlocal = m_zero
    params%rc = m_zero
    params%h = m_zero
    params%k = m_zero
 
    read(unit, *) params%title
    read(unit, *) params%conf%z, params%z_val, params%pspdat
    read(unit, *, iostat=iostat) params%pspcod, params%pspxc, params%l_max, lloc, params%mmax, params%r2well
 
    select case (params%pspcod)
    case (3)
 
      ! Reads the file in a hopefully smart way
      iostat = 1
      j = 5
      read(unit,'(a)') line
      do while((iostat /= 0) .and. (j > 0))
        j = j - 1
        read(line, *, iostat=iostat) params%rlocal, params%c(1:j)
      end do
      if (j<1) read(line, *, iostat=iostat) params%atom_name, params%z_val, params%rlocal
      if (iostat /= 0) then
        load_params = 1
        pop_sub(load_params)
        return
      end if
 
      read(unit,'(a)', iostat = iostat) line
      if (iostat /= 0) then
        load_params = 0
        pop_sub(load_params)
        return
      end if
      iostat = 1
      j = 4
      do while((iostat /= 0) .and. (j > 0))
        j = j - 1
        read(line, *, iostat=iostat) params%rc(0), (params%h(0, i, i), i = 1, j)
      end do
      if (j < 0) then
        load_params = 2
        pop_sub(load_params)
        return
      end if
 
      kloop: do k = 1, 3
        read(unit, '(a)', iostat = iostat) line
        if (iostat /= 0) exit kloop
        iostat = 1
        j = 4
        do while((iostat /= 0) .and. (j > 0))
          j = j - 1
          read(line, *, iostat = iostat) params%rc(k), (params%h(k, i, i), i = 1, j)
        end do
        if (abs(params%rc(k)) <= m_epsilon) exit kloop
        read(unit, '(a)') line
        iostat = 1
        j = 4
        do while((iostat /= 0) .and. (j>0))
          j = j - 1
          read(line, *, iostat = iostat) (params%k(k, i, i), i = 1, 3)
        end do
      end do kloop
 
      ! Fill in the rest of the parameter matrices...
      params%h(0, 1, 2) = -m_half      * sqrt(m_three/m_five)            * params%h(0, 2, 2)
      params%h(0, 1, 3) =  m_half      * sqrt(m_five/21.0_real64)         * params%h(0, 3, 3)
      params%h(0, 2, 3) = -m_half      * sqrt(100.0_real64/63.0_real64)    * params%h(0, 3, 3)
      params%h(1, 1, 2) = -m_half      * sqrt(m_five/7.0_real64)            * params%h(1, 2, 2)
      params%h(1, 1, 3) =  m_one/6.0_real64 * sqrt(35.0_real64/11.0_real64)     * params%h(1, 3, 3)
      params%h(1, 2, 3) = -m_one/6.0_real64 * (14.0_real64 / sqrt(11.0_real64)) * params%h(1, 3, 3)
      params%h(2, 1, 2) = -m_half      * sqrt(7.0_real64/9.0_real64)            * params%h(2, 2, 2)
      params%h(2, 1, 3) =  m_half      * sqrt(63.0_real64/143.0_real64)    * params%h(2, 3, 3)
      params%h(2, 2, 3) = -m_half      * (18.0_real64/sqrt(143.0_real64))  * params%h(2, 3, 3)
 
      params%k(0, 1, 2) = -m_half      * sqrt(m_three/m_five)            * params%k(0, 2, 2)
      params%k(0, 1, 3) =  m_half      * sqrt(m_five/21.0_real64)         * params%k(0, 3, 3)
      params%k(0, 2, 3) = -m_half      * sqrt(100.0_real64/63.0_real64)    * params%k(0, 3, 3)
      params%k(1, 1, 2) = -m_half      * sqrt(m_five/7.0_real64)            * params%k(1, 2, 2)
      params%k(1, 1, 3) =  m_one/6.0_real64 * sqrt(35.0_real64/11.0_real64)     * params%k(1, 3, 3)
      params%k(1, 2, 3) = -m_one/6.0_real64 * (14.0_real64 / sqrt(11.0_real64)) * params%k(1, 3, 3)
      params%k(2, 1, 2) = -m_half      * sqrt(7.0_real64/9.0_real64)            * params%k(2, 2, 2)
      params%k(2, 1, 3) =  m_half      * sqrt(63.0_real64/143.0_real64)    * params%k(2, 3, 3)
      params%k(2, 2, 3) = -m_half      * (18.0_real64/sqrt(143.0_real64))  * params%k(2, 3, 3)
 
 
    case (10)
 
      read(unit, *) params%rlocal, nn, params%c(1)
      read(unit, *) nnonloc
      read(unit, '(a)', iostat = iostat) line
      kloop2 : do k = 0, 3
        if (iostat /= 0) exit kloop2
        iostat = 1
        j = 4
        do while((iostat /= 0) .and. (j > 0))
          j = j - 1
          read(line, *, iostat = iostat) params%rc(k), nn, (params%h(k, 1, i), i = 1, j)
        end do
        if (abs(params%rc(k)) <= m_epsilon) exit kloop2
        read(unit, '(a)') line
        select case (nn)
        case (3)
          read(line, *) (params%h(1, 2, i), i = 2, 3)
          read(unit, '(a)') line
          read(line, *) params%h(1, 3, 3)
          !Reading the k matrix, if possible
          read(unit, '(a)', iostat=iostat) line
          if (iostat /= 0) exit kloop2
          read(line, *, iostat=iostat) (params%k(1, 1, i), i = 1, 3)
          if (iostat /= 0) continue !No k matrix
          read(unit, '(a)') line
          read(line, *) (params%k(1, 2, i), i = 2, 3)
          read(unit, '(a)') line
          read(line, *) params%k(1, 3, 3)
          read(unit, '(a)', iostat = iostat) line
        case (2)
          read(line, *) params%h(0, 2, 2)
          !Reading the k matrix, if possible
          read(unit, '(a)', iostat=iostat) line
          if (iostat /= 0) exit kloop2
          read(line, *, iostat=iostat) (params%k(1, 1, i), i = 1, 2)
          if (iostat /= 0) continue !No k matrix
          read(unit, '(a)') line
          read(line, *) params%k(1, 2, 2)
          read(unit, '(a)', iostat = iostat) line
        case default
          message(1) = "Error parsing the pseudopotential"
          call messages_fatal(1, namespace=namespace)
        end select
      end do kloop2
 
    case default
      message(1) = "Inconsistency in pseudopotential file:"
      write(message(2),'(a,i2)') "  expecting pspcod = 3, but found ", params%pspcod
      call messages_fatal(2, namespace=namespace)
    end select
 
 
    ! Parameters are symmetric.
    do k = 0, 3
      do i = 1, 3
        do j = i + 1, 3
          params%h(k, j, i) = params%h(k, i, j)
          params%k(k, j, i) = params%k(k, i, j)
        end do
      end do
    end do
 
    load_params = 0
    pop_sub(load_params)
  end function load_params
 
 
  ! ---------------------------------------------------------
  subroutine get_cutoff_radii(psp)
    type(ps_hgh_t), intent(inout)     :: psp
 
    integer  :: ir, l, i
    real(real64) :: dincv, tmp
    real(real64), parameter :: threshold = 1.0e-4_real64
 
    push_sub(get_cutoff_radii)
 
    do l = 0, psp%l_max
      tmp = m_zero
      do i = 1, 3
        do ir = psp%g%nrval, 2, -1
          dincv = abs(psp%kb(ir, l, i))
          if (dincv > threshold) exit
        end do
        tmp = psp%g%rofi(ir + 1)
        psp%kbr(l) = max(tmp, psp%kbr(l))
      end do
    end do
 
    pop_sub(get_cutoff_radii)
  end subroutine get_cutoff_radii
 
 
  ! ---------------------------------------------------------
  ! Local pseudopotential, both in real and reciprocal space.
  ! See Eq. (1)
  subroutine vlocalr_scalar(r, np, p, vloc)
    type(ps_hgh_t), intent(in)    :: p
    real(real64),   intent(in)    :: r(:)
    integer,        intent(in)    :: np
    real(real64),   intent(inout) :: vloc(:)
 
    integer :: ip
    real(real64) :: r1, r2
 
    ! The tolerance is set such that |erf(x/sqrt(2))-2*x/sqrt(2*pi)| is below 2e-16
    ! This leads to a tolerance of 1e-5.
    ! However, this values produces an arithmetic error for GSL so it is increase to 5e-5
    real(real64), parameter :: tol = 5e-5_real64
 
    push_sub(vlocalr_scalar)
 
    do ip = 1, np
      r1 = r(ip)/p%rlocal
 
      if (r1 < tol) then
        vloc(ip) = -(m_two * p%z_val) / (sqrt(m_two * m_pi) * p%rlocal) + p%c(1)
      else
        r2 = r1**2
 
        vloc(ip) = -(p%z_val / r(ip)) * loct_erf(r1 / sqrt(m_two))
        if(r2 < m_two * m_max_exp_arg) then !Else we are getting an underflow
          vloc(ip) = vloc(ip) + exp(-m_half*r2) * (p%c(1) + r2*(p%c(2) + r2*(p%c(3) + p%c(4)*r2)))
        end if
      end if
    end do
 
    pop_sub(vlocalr_scalar)
  end subroutine vlocalr_scalar
 
 
  ! ---------------------------------------------------------
  ! Projector in real space, see Eq. 3
  subroutine projectorr_scalar(r, np, p, i, l, proj)
    type(ps_hgh_t), intent(in)    :: p
    real(real64),   intent(in)    :: r(:)
    integer,        intent(in)    :: np
    integer,        intent(in)    :: i
    integer,        intent(in)    :: l
    real(real64),   intent(inout) :: proj(:)
 
    integer :: ip
    real(real64) :: x, y, rr, arg
 
    push_sub(projectorr_scalar)
 
    x = l + real(4*i-1, real64)  / m_two
    y = loct_gamma(x)
    x = sqrt(y)
    rr = m_one
 
    do ip = 1, np
 
      if (l /= 0 .or. i /= 1) then
        rr = r(ip) ** (l + 2*(i-1))
      end if
 
      arg = -r(ip)**2/(m_two*p%rc(l)**2)
      if(arg > m_min_exp_arg) then
        proj(ip) = sqrt(m_two) * rr * exp(-r(ip)**2/(m_two*p%rc(l)**2)) / &
          (p%rc(l)**(l + real(4*i-1, real64)  / m_two) * x)
      else
        proj(ip) = m_zero
      end if
    end do
 
    pop_sub(projectorr_scalar)
  end subroutine projectorr_scalar
 
 
  ! ---------------------------------------------------------
  subroutine solve_schroedinger(psp, ierr, namespace)
    type(ps_hgh_t),    intent(inout) :: psp
    integer,           intent(out)   :: ierr
    type(namespace_t), intent(in)    :: namespace
 
    integer :: iter, ir, l, nnode, nprin, i, j, irr, n, k
    real(real64) :: vtot, a2b4, diff, nonl
    real(real64), allocatable :: prev(:, :), rho(:, :), ve(:, :)
    real(real64), parameter :: tol = 1.0e-4_real64
    real(real64) :: e, z, dr, rmax
    real(real64), allocatable :: s(:), hato(:), g(:), y(:)
 
    push_sub(solve_schroedinger)
 
    ierr = 0
 
    ! Allocations...
    safe_allocate(   s(1:psp%g%nrval))
    safe_allocate(hato(1:psp%g%nrval))
    safe_allocate(   g(1:psp%g%nrval))
    safe_allocate(   y(1:psp%g%nrval))
    safe_allocate(prev(1:psp%g%nrval, 1))
    safe_allocate( rho(1:psp%g%nrval, 1))
    safe_allocate(  ve(1:psp%g%nrval, 1))
    hato = m_zero
    g = m_zero
    y = m_zero
    rho = m_zero
    ve = m_zero
 
    ! These numerical parameters have to be fixed for egofv to work.
    s(2:psp%g%nrval) = psp%g%drdi(2:psp%g%nrval)*psp%g%drdi(2:psp%g%nrval)
    s(1) = s(2)
    a2b4 = m_fourth*psp%g%a**2
 
    ! "Double" self-consistent loop: nonlocal and XC parts have to be calculated
    ! self-consistently.
    diff = 1e5_real64
    iter = 0
    self_consistent: do while(diff > tol)
      prev = rho
      iter = iter + 1
      do n = 1, psp%conf%p
        l = psp%conf%l(n)
        do ir = 2, psp%g%nrval
          vtot = 2*psp%vlocal(ir) + ve(ir, 1) + real(l*(l + 1), real64)/(psp%g%rofi(ir)**2)
          nonl = m_zero
          if (iter>2 .and. psp%l_max >= 0 .and. psp%rphi(ir, n) > 1.0e-7_real64) then
            do i = 1, 3
              do j = 1, 3
                do irr = 2, psp%g%nrval
                  if(psp%kb(irr,l,j) < 1e-100_real64 .or. psp%kb(ir, l, i) < 1e-100_real64) cycle !To avoid FPE
                  nonl = nonl + psp%h(l, i, j)*psp%kb(ir, l, i)* &
                    psp%g%drdi(irr)*psp%g%rofi(irr)*psp%rphi(irr, n)*psp%kb(irr,l,j)
                end do
              end do
            end do
            nonl = m_two * nonl / psp%rphi(ir, n) * psp%g%rofi(ir)
          end if
          vtot = vtot + nonl
          hato(ir) = vtot*s(ir) + a2b4
        end do
        hato(1) = hato(2)
        ! We will assume there is only the possibility of two equal l numbers...
        nnode = 1
        do k = 1, n - 1
          if (psp%conf%l(k) == psp%conf%l(n)) nnode = 2
        end do
        nprin = l + 1
        if (iter == 1) then
          e = -((psp%z_val/real(nprin, real64) )**2)
          z = psp%z_val
        else
          e = psp%eigen(n)
          z = psp%z_val
        end if
        dr = -1.0e5_real64
        rmax = psp%g%rofi(psp%g%nrval)
        call egofv(namespace, hato, s, psp%g%nrval, e, g, y, l, z, &
          real(psp%g%a, real64), real(psp%g%b, real64), rmax, nprin, nnode, dr, ierr)
        if (ierr /= 0) exit self_consistent
        psp%eigen(n) = e
 
        psp%rphi(2:psp%g%nrval, n) = g(2:psp%g%nrval) * sqrt(psp%g%drdi(2:psp%g%nrval))
        psp%rphi(1, n) = psp%rphi(2, n)
      end do
      rho = m_zero
      do n = 1, psp%conf%p
        rho(1:psp%g%nrval, 1) = rho(1:psp%g%nrval, 1) + psp%conf%occ(n,1)*psp%rphi(1:psp%g%nrval, n)**2
      end do
      if (iter>1) rho = m_half*(rho + prev)
      diff = sqrt(sum(psp%g%drdi(2:psp%g%nrval)*(rho(2:psp%g%nrval, 1)-prev(2:psp%g%nrval, 1))**2))
      call atomhxc(namespace, 'LDA', psp%g, 1, rho, ve)
 
    end do self_consistent
 
    if (ierr == 0) then
      !  checking normalization of the calculated wavefunctions
      !do l = 0, psp%l_max_occ
      do n = 1, psp%conf%p
        e = sqrt(sum(psp%g%drdi(2:psp%g%nrval)*psp%rphi(2:psp%g%nrval, n)**2))
        e = abs(e - m_one)
        if (e > 1.0e-5_real64) then
          write(message(1), '(a,i2,a)') "Eigenstate for n = ", n , ' is not normalized'
          write(message(2), '(a, f12.6,a)') '(abs(1-norm) = ', e, ')'
          call messages_warning(2, namespace=namespace)
        end if
      end do
 
      ! Output in Ha and not in stupid Rydbergs.
      psp%eigen = psp%eigen / m_two
    end if
 
    ! Deallocations.
    safe_deallocate_a(s)
    safe_deallocate_a(hato)
    safe_deallocate_a(g)
    safe_deallocate_a(y)
    safe_deallocate_a(prev)
    safe_deallocate_a(rho)
    safe_deallocate_a(ve)
 
    pop_sub(solve_schroedinger)
  end subroutine solve_schroedinger
 
 
  ! ---------------------------------------------------------
  subroutine hgh_debug(psp, dir, namespace)
    type(ps_hgh_t),    intent(in) :: psp
    character(len=*),  intent(in) :: dir
    type(namespace_t), intent(in) :: namespace
 
    integer :: hgh_unit, loc_unit, dat_unit, kbp_unit, wav_unit, i, l, k
    character(len=256) :: dirname
 
    push_sub(hgh_debug)
 
    ! Open files.
    dirname = trim(dir)//'/hgh.'//trim(psp%atom_name)
    call io_mkdir(trim(dirname), namespace)
    hgh_unit = io_open(trim(dirname)//'/hgh', namespace, action='write')
    loc_unit = io_open(trim(dirname)//'/local', namespace, action='write')
    dat_unit = io_open(trim(dirname)//'/info', namespace, action='write')
    kbp_unit = io_open(trim(dirname)//'/nonlocal', namespace, action='write')
    wav_unit = io_open(trim(dirname)//'/wave', namespace, action='write')
 
    ! Writes down the input file, to be checked against SHARE_DIR/pseudopotentials/HGH/ATOM_NAME.hgh
    write(hgh_unit,'(a5,i6,5f12.6)') psp%atom_name, psp%z_val, psp%rlocal, psp%c(1:4)
    write(hgh_unit,'(  11x,4f12.6)') psp%rc(0), (psp%h(0,i,i), i = 1, 3)
    do k = 1, 3
      write(hgh_unit,'(  11x,4f12.6)') psp%rc(k), (psp%h(k, i, i), i = 1, 3)
      write(hgh_unit,'(  23x,4f12.6)')            (psp%k(k, i, i), i = 1, 3)
    end do
 
    ! Writes down some info.
    write(dat_unit,'(a,i3)')        'lmax  = ', psp%l_max
    if (psp%l_max >= 0) then
      write(dat_unit,'(a,4f14.6)') 'kbr   = ', psp%kbr
    end if
    write(dat_unit,'(a,5f14.6)')    'eigen = ', psp%eigen
    write(dat_unit,'(a,5f14.6)')    'occ   = ', psp%conf%occ(1:psp%conf%p, 1)
    ! Writes down local part.
    do i = 1, psp%g%nrval
      write(loc_unit, *) psp%g%rofi(i), psp%vlocal(i)
    end do
 
    ! Writes down nonlocal part.
    if (psp%l_max >= 0) then
      do i = 1, psp%g%nrval
        write(kbp_unit, '(10es14.4)') psp%g%rofi(i), ( (psp%kb(i, l, k) ,k = 1, 3),l = 0, psp%l_max)
      end do
    end if
 
    ! And the pseudo-wavefunctions.
    do i = 1, psp%g%nrval
      write(wav_unit, *) psp%g%rofi(i), (psp%rphi(i, l), l = 1, psp%conf%p)
    end do
 
    ! Closes files and exits
    call io_close(hgh_unit)
    call io_close(loc_unit)
    call io_close(wav_unit)
    call io_close(dat_unit)
    call io_close(kbp_unit)
 
    pop_sub(hgh_debug)
  end subroutine hgh_debug
 
end module ps_hgh_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: