ps_xml.F90

Go to the documentation of this file.

!! Copyright (C) 2015-2018 Xavier Andrade
!!
!! 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_xml_oct_m
  use debug_oct_m
  use global_oct_m
  use, intrinsic :: iso_fortran_env
  use messages_oct_m
  use namespace_oct_m
  use profiling_oct_m
  use pseudo_oct_m
  use string_oct_m
 
  implicit none
 
  private
  public ::     &
    ps_xml_t,       &
    ps_xml_init,    &
    ps_xml_end
 
  type ps_xml_t
    ! Components are public by default
    logical            :: initialized
    logical            :: kleinman_bylander
    logical            :: nlcc
    logical            :: has_density
    integer            :: atomic_number
    real(real64)       :: mass
    real(real64)       :: valence_charge
    integer            :: lmax
    integer            :: llocal
    integer            :: nchannels
    integer            :: grid_size
    integer            :: nwavefunctions
    real(real64), allocatable :: grid(:)
    real(real64), allocatable :: weights(:)
    real(real64), allocatable :: potential(:, :)
    real(real64), allocatable :: wavefunction(:, :)
    real(real64), allocatable :: projector(:, :, :)
    real(real64), allocatable :: dij(:, :, :)
    real(real64), allocatable :: qij(:, :, :)
    logical :: has_qij_integrals
    logical :: has_radial_augmentation
    logical :: radial_augmentation_has_angular_channels
    integer :: radial_augmentation_lmax
    real(real64), allocatable :: radial_augmentation(:, :, :, :)
    real(real64), allocatable :: nlcc_density(:)
    real(real64), allocatable :: density(:)
    integer, allocatable :: wf_n(:)
    integer, allocatable :: wf_l(:)
    real(real64), allocatable :: wf_occ(:)
    type(pseudo_t)     :: pseudo
    integer, allocatable :: kb_radius(:,:)
  end type ps_xml_t
 
contains
 
  ! ---------------------------------------------------------
  subroutine ps_xml_init(this, namespace, filename, fmt, ierr)
    type(ps_xml_t),    intent(inout) :: this
    type(namespace_t), intent(in)    :: namespace
    character(len=*),  intent(in)    :: filename
    integer,           intent(in)    :: fmt
    integer,           intent(out)   :: ierr
 
    integer :: ll, ii, jj, ic, jc, ip, nproj
    type(pseudo_t) :: pseudo
 
    push_sub(ps_xml_init)
 
    this%initialized = .false.
 
    call pseudo_init(pseudo, string_f_to_c(filename), fmt, ierr)
 
    if (ierr == pseudo_status_file_not_found) then
      call messages_write("Pseudopotential file '" // trim(filename) // "' not found")
      call messages_fatal(namespace=namespace)
    end if
 
    if (ierr == pseudo_status_unknown_format) then
      call messages_write("Cannot determine the format for pseudopotential file '" // trim(filename) // "'")
      call messages_fatal(namespace=namespace)
    end if
 
    if (ierr == pseudo_status_unsupported_type_ultrasoft) then
      call messages_write("Ultrasoft pseudopotential file '" // trim(filename) // "' not supported")
      call messages_fatal(namespace=namespace)
    end if
 
    if (ierr == pseudo_status_unsupported_type_paw) then
      call messages_write("PAW pseudopotential file '" // trim(filename) // "' not supported")
      call messages_fatal(namespace=namespace)
    end if
 
    if (ierr == pseudo_status_unsupported_type) then
      call messages_write("Pseudopotential file '" // trim(filename) // "' not supported")
      call messages_fatal(namespace=namespace)
    end if
 
    if (ierr == pseudo_status_format_not_supported) then
      call messages_write("Pseudopotential file '" // trim(filename) // "' not supported")
      call messages_fatal(namespace=namespace)
    end if
 
    this%initialized = .true.
    this%valence_charge = pseudo_valence_charge(pseudo)
    this%mass = pseudo_mass(pseudo)
    this%lmax = pseudo_lmax(pseudo)
    this%llocal = pseudo_llocal(pseudo)
    this%nchannels = pseudo_nchannels(pseudo)
 
    this%kleinman_bylander = (pseudo_type(pseudo) == pseudo_type_kleinman_bylander .or. &
      pseudo_type(pseudo) == pseudo_type_ultrasoft)
    this%has_qij_integrals = pseudo_has_qij_integrals(pseudo)
    if (pseudo_type(pseudo) == pseudo_type_ultrasoft .and. .not. this%has_qij_integrals) then
      call messages_write("Ultrasoft pseudopotential file '" // trim(filename) // "' does not contain PP_Q or PP_QIJ")
      call messages_fatal(namespace=namespace)
    end if
    this%has_radial_augmentation = pseudo_has_radial_augmentation(pseudo)
    this%radial_augmentation_has_angular_channels = pseudo_radial_augmentation_has_angular_channels(pseudo)
    this%radial_augmentation_lmax = pseudo_radial_augmentation_lmax(pseudo)
 
    this%grid_size = pseudo_mesh_size(pseudo)
    nproj = pseudo_nprojectors(pseudo)
 
    safe_allocate(this%grid(1:this%grid_size))
    safe_allocate(this%weights(1:this%grid_size))
 
    call pseudo_grid(pseudo, this%grid(1))
    call pseudo_grid_weights(pseudo, this%weights(1))
 
    if (.not. this%kleinman_bylander) then
 
      safe_allocate(this%potential(1:this%grid_size, 0:this%lmax))
      safe_allocate(this%wavefunction(1:this%grid_size, 0:this%lmax))
 
      do ll = 0, this%lmax
        if (.not. pseudo_has_radial_function(pseudo, ll)) then
          call messages_write("The pseudopotential file '"//trim(filename)//"' does not contain")
          call messages_new_line()
          call messages_write("the wave functions. Octopus cannot use it.")
          call messages_fatal(namespace=namespace)
        end if
 
        call pseudo_radial_function(pseudo, ll, this%wavefunction(1, ll))
        call pseudo_radial_potential(pseudo, ll, this%potential(1, ll))
      end do
 
    else
 
      safe_allocate(this%potential(1:this%grid_size, this%llocal:this%llocal))
 
      call pseudo_local_potential(pseudo, this%potential(1, this%llocal))
 
      safe_allocate(this%projector(1:this%grid_size, 0:this%lmax, 1:this%nchannels))
 
      this%projector = m_zero
 
      safe_allocate(this%kb_radius(0:this%lmax, 1:this%nchannels))
      this%kb_radius = this%grid_size
 
      do ll = 0, this%lmax
        if (this%llocal == ll) cycle
        do ic = 1,  pseudo_nprojectors_per_l(pseudo, ll)
          call pseudo_projector(pseudo, ll, ic, this%projector(1, ll, ic))
 
          do ip = this%grid_size, 1, -1
            if (abs(this%projector(ip, ll, ic)) > m_epsilon) then
              this%kb_radius(ll, ic) = ip
              exit
            end if
          end do
        end do
      end do
 
      if (this%has_radial_augmentation) then
        safe_allocate(this%radial_augmentation(1:this%grid_size, 1:nproj, 1:nproj, 0:this%radial_augmentation_lmax))
 
        do ii = 1, nproj
          do jj = 1, nproj
            do ll = 0, this%radial_augmentation_lmax
              call pseudo_radial_augmentation(pseudo, ii, jj, ll, this%radial_augmentation(1, ii, jj, ll))
            end do
          end do
        end do
      end if
 
 
      safe_allocate(this%dij(0:this%lmax, 1:this%nchannels, 1:this%nchannels))
      safe_allocate(this%qij(0:this%lmax, 1:this%nchannels, 1:this%nchannels))
 
      this%dij = 0.0_real64
      this%qij = 0.0_real64
      do ll = 0, this%lmax
        do ic = 1, this%nchannels
          do jc = 1, this%nchannels
            this%dij(ll, ic, jc) = pseudo_dij(pseudo, ll, ic, jc)
            if (this%has_qij_integrals) this%qij(ll, ic, jc) = pseudo_qij(pseudo, ll, ic, jc)
          end do
        end do
      end do
 
      this%nwavefunctions = pseudo_nwavefunctions(pseudo)
 
      safe_allocate(this%wavefunction(1:this%grid_size, 1:this%nwavefunctions))
      safe_allocate(this%wf_n(1:this%nwavefunctions))
      safe_allocate(this%wf_l(1:this%nwavefunctions))
      safe_allocate(this%wf_occ(1:this%nwavefunctions))
 
      do ii = 1, this%nwavefunctions
        call pseudo_wavefunction(pseudo, ii, this%wf_n(ii), this%wf_l(ii), this%wf_occ(ii), this%wavefunction(1, ii))
      end do
 
    end if
 
    this%has_density = pseudo_has_density(pseudo)
 
    if (this%has_density) then
      safe_allocate(this%density(1:this%grid_size))
      call pseudo_density(pseudo, this%density(1))
    end if
 
    this%nlcc = pseudo_has_nlcc(pseudo)
    if (this%nlcc) then
      safe_allocate(this%nlcc_density(1:this%grid_size))
      call pseudo_nlcc_density(pseudo, this%nlcc_density(1))
    end if
 
    if (.not. this%kleinman_bylander) call ps_xml_check_normalization(this, namespace)
 
    this%pseudo = pseudo
 
    pop_sub(ps_xml_init)
  end subroutine ps_xml_init
 
  ! ---------------------------------------------------------
  subroutine ps_xml_check_normalization(this, namespace)
    type(ps_xml_t),    intent(in) :: this
    type(namespace_t), intent(in) :: namespace
 
    integer :: ll, ip
    real(real64)   :: nrm, rr
 
    push_sub(ps_xml_check_normalization)
 
    !  checking normalization of the wavefunctions
    do ll = 0, this%lmax
      nrm = m_zero
      do ip = 1, this%grid_size
        rr = this%grid(ip)
        nrm = nrm + this%wavefunction(ip, ll)**2*this%weights(ip)*rr**2
      end do
      nrm = sqrt(nrm)
 
      nrm = abs(nrm - m_one)
      if (nrm > 1.0e-5_real64) then
        write(message(1), '(a,i2,a)') "Eigenstate for l = ", ll, ' is not normalized'
        write(message(2), '(a, f12.6,a)') '(abs(1 - norm) = ', nrm, ')'
        call messages_warning(2, namespace=namespace)
      end if
 
    end do
 
    pop_sub(ps_xml_check_normalization)
  end subroutine ps_xml_check_normalization
 
  ! ---------------------------------------------------------
  subroutine ps_xml_end(this)
    type(ps_xml_t), intent(inout) :: this
 
    push_sub(ps_xml_end)
 
    call pseudo_end(this%pseudo)
 
    safe_deallocate_a(this%grid)
    safe_deallocate_a(this%weights)
    safe_deallocate_a(this%potential)
    safe_deallocate_a(this%wavefunction)
    safe_deallocate_a(this%projector)
    safe_deallocate_a(this%dij)
    safe_deallocate_a(this%qij)
    safe_deallocate_a(this%radial_augmentation)
    safe_deallocate_a(this%nlcc_density)
    safe_deallocate_a(this%density)
    safe_deallocate_a(this%wf_n)
    safe_deallocate_a(this%wf_l)
    safe_deallocate_a(this%wf_occ)
    safe_deallocate_a(this%kb_radius)
 
    this%initialized = .false.
 
    pop_sub(ps_xml_end)
  end subroutine ps_xml_end
 
end module ps_xml_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: