nonlocal_pseudopotential.F90

Go to the documentation of this file.

!! Copyright (C) 2009 X. 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 nonlocal_pseudopotential_oct_m
  use accel_oct_m
  use accel_blas_oct_m
  use batch_oct_m
  use batch_ops_oct_m
  use blas_oct_m
  use debug_oct_m
  use electron_space_oct_m
  use epot_oct_m
  use global_oct_m
  use hardware_oct_m
  use hgh_projector_oct_m
  use kb_projector_oct_m
  use lalg_basic_oct_m
  use math_oct_m
  use mesh_oct_m
  use messages_oct_m
  use mpi_oct_m
  use profiling_oct_m
  use projector_oct_m
  use projector_matrix_oct_m
  use ps_oct_m
  use rkb_projector_oct_m
  use space_oct_m
  use states_elec_oct_m
  use states_elec_dim_oct_m
  use submesh_oct_m
  use types_oct_m
  use wfs_elec_oct_m
 
  implicit none
 
  private
 
  public ::                                         &
    nonlocal_pseudopotential_t,                     &
    nonlocal_pseudopotential_accel_rebuild,         &
    nonlocal_pseudopotential_rebind_projectors,     &
    projection_t
 
  type nonlocal_pseudopotential_t
    private
    type(projector_matrix_t), allocatable, public :: projector_matrices(:)
    integer,                               public :: nprojector_matrices
    logical,                               public :: apply_projector_matrices
    logical,                               public :: has_non_local_potential
    integer                                       :: full_projection_size
    integer,                               public :: max_npoints
    integer,                               public :: total_points
    integer                                       :: max_nprojs
    logical                                       :: projector_mix
    complex(real64),          allocatable, public :: projector_phases(:, :, :, :)
    integer,                  allocatable, public :: projector_to_atom(:)
    integer                                       :: nregions
    integer,                  allocatable         :: regions(:)
    integer,                               public :: nphase
    type(accel_mem_t)                             :: buff_offsets
    type(accel_mem_t)                             :: buff_matrices
    type(accel_mem_t)                             :: buff_maps
    type(accel_mem_t)                             :: buff_sphere_weights
    type(accel_mem_t)                             :: buff_scals
    type(accel_mem_t)                             :: buff_position
    type(accel_mem_t)                             :: buff_pos
    type(accel_mem_t)                             :: buff_invmap
    type(accel_mem_t)                             :: buff_invmap_mat
    type(accel_mem_t)                             :: buff_point_to_mat
    type(accel_mem_t)                             :: buff_proj_to_mat
    type(accel_mem_t),                     public :: buff_projector_phases
    type(accel_mem_t),          pointer,   public :: buff_bra_phasepsi => null()
    type(accel_mem_t),          pointer,   public :: buff_bra_projection_temp => null()
    type(accel_mem_t)                             :: buff_mix
    logical                                       :: projector_self_overlap
    real(real64),             pointer,     public :: spin(:,:,:) => null()
  contains
 
    procedure :: init => nonlocal_pseudopotential_init
    
    procedure :: build => nonlocal_pseudopotential_build_proj
    
    procedure :: end => nonlocal_pseudopotential_destroy_proj
    
    procedure :: has_self_overlap => nonlocal_pseudopotential_self_overlap
    
    procedure :: dstart => dnonlocal_pseudopotential_start
    
    procedure :: zstart => znonlocal_pseudopotential_start
    
    procedure :: dfinish => dnonlocal_pseudopotential_finish
    
    procedure :: zfinish => znonlocal_pseudopotential_finish
    
    procedure :: dforce => dnonlocal_pseudopotential_force
    
    procedure :: zforce => znonlocal_pseudopotential_force
    
    procedure :: dposition_commutator => dnonlocal_pseudopotential_position_commutator
    
    procedure :: zposition_commutator => znonlocal_pseudopotential_position_commutator
    
    procedure :: dr_vn_local => dnonlocal_pseudopotential_r_vnlocal
    
    procedure :: zr_vn_local => znonlocal_pseudopotential_r_vnlocal
    
  end type nonlocal_pseudopotential_t
 
 
  !
  type projection_t
    private
    real(real64),      allocatable :: dprojection(:, :)
    complex(real64),   allocatable :: zprojection(:, :)
    type(accel_mem_t)              :: buff_projection
    type(accel_mem_t)              :: buff_spin_to_phase
    integer                        :: cuda_stream_projection_DtH
  end type projection_t
 
contains
 
  ! ---------------------------------------------------------
  subroutine nonlocal_pseudopotential_init(this)
    class(nonlocal_pseudopotential_t), intent(inout) :: this
 
    push_sub(nonlocal_pseudopotential_init)
 
    this%apply_projector_matrices = .false.
    this%has_non_local_potential = .false.
    this%nprojector_matrices = 0
 
    this%projector_self_overlap = .false.
 
    pop_sub(nonlocal_pseudopotential_init)
  end subroutine nonlocal_pseudopotential_init
 
  ! ---------------------------------------------------------
  subroutine nonlocal_pseudopotential_rebind_projectors(this, epot)
    class(nonlocal_pseudopotential_t), target, intent(inout) :: this
    type(epot_t),               target, intent(in)   :: epot
 
    integer :: imat, iatom
    type(projector_matrix_t), pointer :: pmat
 
    push_sub(nonlocal_pseudopotential_rebind_projectors)
 
    ! After the assignment these still point to the buffers of the source: drop them, do not free them
    nullify(this%buff_bra_phasepsi)
    nullify(this%buff_bra_projection_temp)
    safe_allocate(this%buff_bra_phasepsi)
    safe_allocate(this%buff_bra_projection_temp)
 
    if (.not. allocated(this%projector_matrices) .or. this%nprojector_matrices == 0) then
      pop_sub(nonlocal_pseudopotential_rebind_projectors)
      return
    end if
 
    assert(allocated(this%projector_to_atom))
    assert(size(this%projector_matrices) >= this%nprojector_matrices)
    assert(size(this%projector_to_atom) >= this%nprojector_matrices)
 
    do imat = 1, this%nprojector_matrices
      pmat => this%projector_matrices(imat)
      iatom = this%projector_to_atom(imat)
 
      assert(iatom >= 1 .and. iatom <= epot%natoms)
      assert(allocated(epot%proj))
      assert(allocated(epot%proj(iatom)%sphere%map))
      assert(allocated(epot%proj(iatom)%sphere%rel_x))
 
      pmat%map => epot%proj(iatom)%sphere%map
      pmat%position => epot%proj(iatom)%sphere%rel_x
    end do
 
    pop_sub(nonlocal_pseudopotential_rebind_projectors)
  end subroutine nonlocal_pseudopotential_rebind_projectors
 
 
 
  !--------------------------------------------------------
  subroutine nonlocal_pseudopotential_destroy_proj(this)
    class(nonlocal_pseudopotential_t), intent(inout) :: this
 
    integer :: iproj
 
    push_sub(nonlocal_pseudopotential_destroy_proj)
 
    if (allocated(this%projector_matrices)) then
 
      if (accel_is_enabled()) then
        call accel_free_buffer(this%buff_offsets)
        call accel_free_buffer(this%buff_matrices)
        call accel_free_buffer(this%buff_maps)
        if (accel_buffer_is_allocated(this%buff_sphere_weights)) call accel_free_buffer(this%buff_sphere_weights)
        call accel_free_buffer(this%buff_scals)
        call accel_free_buffer(this%buff_position)
        call accel_free_buffer(this%buff_pos)
        call accel_free_buffer(this%buff_invmap)
        call accel_free_buffer(this%buff_invmap_mat)
        call accel_free_buffer(this%buff_point_to_mat)
        call accel_free_buffer(this%buff_proj_to_mat)
        if (this%projector_mix) call accel_free_buffer(this%buff_mix)
        if (allocated(this%projector_phases)) call accel_free_buffer(this%buff_projector_phases)
      end if
 
      do iproj = 1, this%nprojector_matrices
        call projector_matrix_deallocate(this%projector_matrices(iproj))
      end do
      safe_deallocate_a(this%regions)
      safe_deallocate_a(this%projector_matrices)
      safe_deallocate_a(this%projector_phases)
      safe_deallocate_a(this%projector_to_atom)
    end if
 
    if (associated(this%buff_bra_phasepsi)) then
      if (accel_buffer_is_allocated(this%buff_bra_phasepsi)) call accel_free_buffer(this%buff_bra_phasepsi)
    end if
    if (associated(this%buff_bra_projection_temp)) then
      if (accel_buffer_is_allocated(this%buff_bra_projection_temp)) call accel_free_buffer(this%buff_bra_projection_temp)
    end if
    safe_deallocate_p(this%buff_bra_phasepsi)
    safe_deallocate_p(this%buff_bra_projection_temp)
 
    pop_sub(nonlocal_pseudopotential_destroy_proj)
  end subroutine nonlocal_pseudopotential_destroy_proj
 
  !-----------------------------------------------------------------
  !
  subroutine nonlocal_pseudopotential_build_proj(this, space, mesh, epot)
    class(nonlocal_pseudopotential_t), target, intent(inout) :: this
    class(space_t),                         intent(in)    :: space
    class(mesh_t),                    intent(in)    :: mesh
    type(epot_t),             target, intent(in)    :: epot
 
    integer :: iatom, iproj, ll, lmax, lloc, mm, ic, jc
    integer :: nmat, imat, ip, iorder
    integer :: nregion, jatom, katom, iregion
    integer, allocatable :: order(:), head(:), region_count(:)
    logical, allocatable :: atom_counted(:)
    logical :: overlap
    type(projector_matrix_t), pointer :: pmat
    type(kb_projector_t),     pointer :: kb_p
    type(rkb_projector_t),    pointer :: rkb_p
    type(hgh_projector_t),    pointer :: hgh_p
 
    push_sub(nonlocal_pseudopotential_build_proj)
 
    call profiling_in("ATOM_COLORING")
 
    ! this is most likely a very inefficient algorithm, O(natom**2) or
    ! O(natom**3), probably it should be replaced by something better.
 
    safe_allocate(order(1:epot%natoms))         ! order(iregion) = ?
    safe_allocate(head(1:epot%natoms + 1))      ! head(iregion) points to the first atom in region iregion
    safe_allocate(region_count(1:epot%natoms))  ! region_count(iregion): number of atoms in that region
    safe_allocate(atom_counted(1:epot%natoms))
 
    this%projector_self_overlap = .false.
    atom_counted = .false.
    order = -1
 
    head(1) = 1
    nregion = 0
    do
      nregion = nregion + 1
      assert(nregion <= epot%natoms)
 
      region_count(nregion) = 0
 
      do iatom = 1, epot%natoms
        if (atom_counted(iatom)) cycle
 
        overlap = .false.
 
        if (.not. projector_is(epot%proj(iatom), proj_none)) then
          assert(associated(epot%proj(iatom)%sphere%mesh))
          do jatom = 1, region_count(nregion)
            katom = order(head(nregion) + jatom - 1)
            if (projector_is(epot%proj(katom), proj_none)) cycle
            overlap = submesh_overlap(epot%proj(iatom)%sphere, epot%proj(katom)%sphere, space)
            if (overlap) exit
          end do
        end if
 
        if (.not. overlap) then
          ! iatom did not overlap with any previously counted atoms:
          ! iatom will be added to the current region
          region_count(nregion) = region_count(nregion) + 1
          order(head(nregion) - 1 + region_count(nregion)) = iatom
          atom_counted(iatom) = .true.
        end if
 
      end do
 
      head(nregion + 1) = head(nregion) + region_count(nregion)
 
      if (all(atom_counted)) exit
    end do
 
    safe_deallocate_a(atom_counted)
    safe_deallocate_a(region_count)
 
    call messages_write('The atoms can be separated in ')
    call messages_write(nregion)
    call messages_write(' non-overlapping groups.')
    call messages_info(debug_only=.true.)
 
    do iregion = 1, nregion
      do iatom = head(iregion), head(iregion + 1) - 1
        if (.not. projector_is(epot%proj(order(iatom)), proj_kb)) cycle
        do jatom = head(iregion), iatom - 1
          if (.not. projector_is(epot%proj(order(jatom)), proj_kb)) cycle
          assert(.not. submesh_overlap(epot%proj(order(iatom))%sphere, epot%proj(order(jatom))%sphere, space))
        end do
      end do
    end do
 
    call profiling_out("ATOM_COLORING")
 
    ! deallocate previous projectors
    call this%end()
 
    ! The two bra scratch buffers.  They are pointers so that
    ! X(nonlocal_pseudopotential_start), whose `this` is intent(in), can fill them.
    safe_allocate(this%buff_bra_phasepsi)
    safe_allocate(this%buff_bra_projection_temp)
 
    ! count projectors
    this%nprojector_matrices = 0
    this%apply_projector_matrices = .false.
    this%has_non_local_potential = .false.
    this%nregions = nregion
 
    !We determine if we have only local potential or not.
    do iorder = 1, epot%natoms
      iatom = order(iorder)
 
      if (.not. projector_is_null(epot%proj(iatom))) then
        this%has_non_local_potential = .true.
        exit
      end if
    end do
 
    do iorder = 1, epot%natoms
      iatom = order(iorder)
 
      if (.not. projector_is_null(epot%proj(iatom))) then
        this%nprojector_matrices = this%nprojector_matrices + 1
        this%apply_projector_matrices = .true.
      end if
    end do
 
    ! This is currently the only not supported case
    if (mesh%use_curvilinear)  this%apply_projector_matrices = .false.
 
    if (.not. this%apply_projector_matrices) then
      safe_deallocate_a(order)
      safe_deallocate_a(head)
 
      pop_sub(nonlocal_pseudopotential_build_proj)
      return
    end if
 
 
    safe_allocate(this%projector_matrices(1:this%nprojector_matrices))
    safe_allocate(this%regions(1:this%nregions + 1))
    safe_allocate(this%projector_to_atom(1:epot%natoms))
 
    this%full_projection_size = 0
    this%regions(this%nregions + 1) = this%nprojector_matrices + 1
 
    this%projector_mix = .false.
 
    iproj = 0
    do iregion = 1, this%nregions
      this%regions(iregion) = iproj + 1
      do iorder = head(iregion), head(iregion + 1) - 1
 
        iatom = order(iorder)
 
        if (projector_is(epot%proj(iatom), proj_none)) cycle
 
        iproj = iproj + 1
 
        pmat => this%projector_matrices(iproj)
 
        this%projector_to_atom(iproj) = iatom
 
        lmax = epot%proj(iatom)%lmax
        lloc = epot%proj(iatom)%lloc
 
        if (projector_is(epot%proj(iatom), proj_kb)) then
 
          ! count the number of projectors for this matrix
          nmat = 0
          do ll = 0, lmax
            if (ll == lloc) cycle
            do mm = -ll, ll
              nmat = nmat + epot%proj(iatom)%kb_p(ll, mm)%n_c
            end do
          end do
 
          call projector_matrix_allocate(pmat, nmat, epot%proj(iatom)%sphere, has_mix_matrix = .false.)
 
          ! generate the matrix
          pmat%dprojectors = m_zero
          imat = 1
          do ll = 0, lmax
            if (ll == lloc) cycle
            do mm = -ll, ll
              kb_p =>  epot%proj(iatom)%kb_p(ll, mm)
              do ic = 1, kb_p%n_c
                call lalg_copy(pmat%npoints, kb_p%p(:, ic), pmat%dprojectors(:, imat))
                ! the per-point part of the volume element is applied on the bra
                ! side (sphere_weights) for meshes with per-point volumes
                pmat%scal(imat) = kb_p%e(ic)*mesh%volume_element
                imat = imat + 1
              end do
            end do
          end do
 
          this%projector_self_overlap = this%projector_self_overlap .or. epot%proj(iatom)%sphere%overlap
 
        else if (projector_is(epot%proj(iatom), proj_hgh)) then
 
          this%projector_mix = .true.
 
          ! count the number of projectors for this matrix
          nmat = 0
          do ll = 0, lmax
            if (ll == lloc) cycle
            do mm = -ll, ll
              nmat = nmat + 3
            end do
          end do
 
          call projector_matrix_allocate(pmat, nmat, epot%proj(iatom)%sphere, &
            has_mix_matrix = .true., is_cmplx = (epot%proj_reltype == spin_orbit))
 
          ! generate the matrix
          if (epot%proj_reltype == spin_orbit) then
            pmat%zprojectors = m_zero
            pmat%zmix = m_zero
          else
            pmat%dprojectors = m_zero
            pmat%dmix = m_zero
          end if
 
          imat = 1
          do ll = 0, lmax
            if (ll == lloc) cycle
            do mm = -ll, ll
              hgh_p =>  epot%proj(iatom)%hgh_p(ll, mm)
 
              ! HGH pseudos mix different components, so we need to
              ! generate a matrix that mixes the projections
              if (epot%proj_reltype == spin_orbit) then
                do ic = 1, 3
                  do jc = 1, 3
                    pmat%zmix(imat - 1 + ic, imat - 1 + jc, 1) = hgh_p%h(ic, jc) + m_half*mm*hgh_p%k(ic, jc)
                    pmat%zmix(imat - 1 + ic, imat - 1 + jc, 2) = hgh_p%h(ic, jc) - m_half*mm*hgh_p%k(ic, jc)
 
                    if (mm < ll) then
                      pmat%zmix(imat - 1 + ic, imat + 3 - 1 + jc, 3) = m_half*hgh_p%k(ic, jc) * &
                        sqrt(real(ll*(ll+1)-mm*(mm+1), real64))
                    end if
 
                    if (-mm < ll) then
                      pmat%zmix(imat - 1 + ic, imat - 3 - 1 + jc, 4) = m_half*hgh_p%k(ic, jc) * &
                        sqrt(real(ll*(ll+1)-mm*(mm-1), real64))
                    end if
                  end do
                end do
              else
                do ic = 1, 3
                  do jc = 1, 3
                    pmat%dmix(imat - 1 + ic, imat - 1 + jc) = hgh_p%h(ic, jc)
                  end do
                end do
              end if
 
              do ic = 1, 3
                if (epot%proj_reltype == spin_orbit) then
                  call lalg_copy(pmat%npoints, hgh_p%zp(:, ic), pmat%zprojectors(:, imat))
                else
                  call lalg_copy(pmat%npoints, hgh_p%dp(:, ic), pmat%dprojectors(:, imat))
                end if
                pmat%scal(imat) = mesh%volume_element
                imat = imat + 1
              end do
 
            end do
          end do
 
          this%projector_self_overlap = this%projector_self_overlap .or. epot%proj(iatom)%sphere%overlap
 
        else if (projector_is(epot%proj(iatom), proj_rkb)) then
          assert(epot%proj_reltype == spin_orbit)
 
          this%projector_mix = .true.
 
          ! count the number of projectors for this matrix
          nmat = 0
          if (lloc /= 0) nmat = nmat + epot%proj(iatom)%kb_p(1, 1)%n_c
 
          do ll = 1, lmax
            if (ll == lloc) cycle
            do mm = -ll, ll
              nmat = nmat + epot%proj(iatom)%rkb_p(ll, mm)%n_c
            end do
          end do
 
          call projector_matrix_allocate(pmat, nmat, epot%proj(iatom)%sphere, &
            has_mix_matrix = .true., is_cmplx = .true.)
 
          pmat%zprojectors = m_zero
          pmat%zmix = m_zero
 
          imat = 1
          if (lloc /= 0) then
            kb_p => epot%proj(iatom)%kb_p(1, 1)
 
            do ic = 1, kb_p%n_c
              pmat%zmix(ic, ic, 1:2) = kb_p%e(ic)
              do ip = 1, pmat%npoints
                pmat%zprojectors(ip, ic) = kb_p%p(ip, ic)
              end do
              pmat%scal(ic) = mesh%volume_element
            end do
            imat = kb_p%n_c + 1
            nullify(kb_p)
          end if
 
          do ll = 1, lmax
            if (ll == lloc) cycle
            do mm = -ll, ll
              rkb_p =>  epot%proj(iatom)%rkb_p(ll, mm)
 
              ! See rkb_projector.F90 for understanding the indices
              do ic = 0, rkb_p%n_c/2-1
                pmat%zmix(imat + ic*2, imat + ic*2, 1) = rkb_p%f(ic*2+1, 1, 1)
                pmat%zmix(imat + ic*2, imat + ic*2, 2) = rkb_p%f(ic*2+1, 2, 2)
 
                pmat%zmix(imat + ic*2+1, imat + ic*2+1, 1) = rkb_p%f(ic*2+2, 1, 1)
                pmat%zmix(imat + ic*2+1, imat + ic*2+1, 2) = rkb_p%f(ic*2+2, 2, 2)
 
                if (mm < ll) then
                  pmat%zmix(imat + ic*2+rkb_p%n_c, imat + ic*2, 4) =  rkb_p%f(ic*2+1, 2, 1)
                  pmat%zmix(imat + ic*2+1+rkb_p%n_c, imat + ic*2+1, 4) =  rkb_p%f(ic*2+2, 2, 1)
                end if
 
                if (-mm < ll) then
                  pmat%zmix(imat + ic*2-rkb_p%n_c, imat + ic*2, 3) =  rkb_p%f(ic*2+1, 1, 2)
                  pmat%zmix(imat + ic*2+1-rkb_p%n_c, imat + ic*2+1, 3) =  rkb_p%f(ic*2+2, 1, 2)
                end if
              end do
 
              do ic = 1, rkb_p%n_c
                call lalg_copy(pmat%npoints, rkb_p%ket(:, ic, 1, 1), pmat%zprojectors(:, imat))
                pmat%scal(imat) = mesh%volume_element
                imat = imat + 1
              end do
            end do
 
            nullify(rkb_p)
          end do
 
          this%projector_self_overlap = this%projector_self_overlap .or. epot%proj(iatom)%sphere%overlap
 
        else
          cycle
        end if
 
        pmat%map => epot%proj(iatom)%sphere%map
        pmat%position => epot%proj(iatom)%sphere%rel_x
 
        pmat%regions = epot%proj(iatom)%sphere%regions
 
        this%full_projection_size = this%full_projection_size + pmat%nprojs
 
        ! the per-point volume weights of the sphere, applied to the bra side
        if (mesh%use_vol_pp) then
          safe_allocate(pmat%sphere_weights(1:max(pmat%npoints, 1)))
          do ip = 1, pmat%npoints
            pmat%sphere_weights(ip) = mesh%vol_pp(pmat%map(ip))
          end do
        end if
 
      end do
    end do
 
    if (mesh%parallel_in_domains) then
      call mesh%mpi_grp%allreduce_inplace(this%projector_self_overlap, 1, mpi_logical, mpi_lor)
    end if
 
    safe_deallocate_a(order)
    safe_deallocate_a(head)
 
    this%total_points = 0
    this%max_npoints = 0
    this%max_nprojs = 0
    do imat = 1, this%nprojector_matrices
      pmat => this%projector_matrices(imat)
 
      this%max_npoints = max(this%max_npoints, pmat%npoints)
      this%max_nprojs = max(this%max_nprojs, pmat%nprojs)
      this%total_points = this%total_points + pmat%npoints
    end do
 
    if (accel_is_enabled()) then
      call nonlocal_pseudopotential_build_accel_buffers(this, space, mesh)
      call nonlocal_pseudopotential_build_projector_phase_accel_buffer(this)
    end if
 
    pop_sub(nonlocal_pseudopotential_build_proj)
 
  end subroutine nonlocal_pseudopotential_build_proj
 
  ! ----------------------------------------------------------------------------------
  subroutine nonlocal_pseudopotential_accel_rebuild(this, space, mesh)
    class(nonlocal_pseudopotential_t), target, intent(inout) :: this
    class(space_t),                            intent(in)    :: space
    class(mesh_t),                             intent(in)    :: mesh
 
    push_sub(nonlocal_pseudopotential_accel_rebuild)
 
    if (.not. accel_is_enabled()) then
      pop_sub(nonlocal_pseudopotential_accel_rebuild)
      return
    end if
 
    call nonlocal_pseudopotential_detach_accel_buffers(this)
 
    if (.not. allocated(this%projector_matrices) .or. this%nprojector_matrices <= 0) then
      pop_sub(nonlocal_pseudopotential_accel_rebuild)
      return
    end if
 
    call nonlocal_pseudopotential_build_accel_buffers(this, space, mesh)
    call nonlocal_pseudopotential_build_projector_phase_accel_buffer(this)
 
    pop_sub(nonlocal_pseudopotential_accel_rebuild)
  end subroutine nonlocal_pseudopotential_accel_rebuild
 
  ! ----------------------------------------------------------------------------------
  subroutine nonlocal_pseudopotential_detach_accel_buffers(this)
    class(nonlocal_pseudopotential_t), intent(inout) :: this
 
    push_sub(nonlocal_pseudopotential_detach_accel_buffers)
 
    call accel_detach_buffer(this%buff_offsets)
    call accel_detach_buffer(this%buff_matrices)
    call accel_detach_buffer(this%buff_maps)
    call accel_detach_buffer(this%buff_scals)
    call accel_detach_buffer(this%buff_position)
    call accel_detach_buffer(this%buff_pos)
    call accel_detach_buffer(this%buff_invmap)
    call accel_detach_buffer(this%buff_invmap_mat)
    call accel_detach_buffer(this%buff_point_to_mat)
    call accel_detach_buffer(this%buff_proj_to_mat)
    call accel_detach_buffer(this%buff_projector_phases)
    call accel_detach_buffer(this%buff_mix)
 
    pop_sub(nonlocal_pseudopotential_detach_accel_buffers)
  end subroutine nonlocal_pseudopotential_detach_accel_buffers
 
  ! ----------------------------------------------------------------------------------
  subroutine nonlocal_pseudopotential_build_accel_buffers(this, space, mesh)
    class(nonlocal_pseudopotential_t), target, intent(inout) :: this
    class(space_t),                            intent(in)    :: space
    class(mesh_t),                             intent(in)    :: mesh
 
    integer, parameter   :: OFFSET_SIZE = 6
    integer, parameter   :: POINTS = 1, projs = 2, matrix = 3, map = 4, scal = 5, mix = 6
    integer              :: imat, matrix_size, scal_size
    integer              :: ip, is, ii, ipos, mix_offset
    integer, allocatable :: cnt(:), invmap(:, :), invmap2(:), pos(:)
    integer, allocatable :: invmap_mat(:, :), invmap_mat2(:)
    integer, allocatable :: offsets(:, :), point_to_mat(:), proj_to_mat(:)
    type(projector_matrix_t), pointer :: pmat
 
    push_sub(nonlocal_pseudopotential_build_accel_buffers)
 
    assert(allocated(this%projector_matrices))
    assert(this%nprojector_matrices > 0)
 
    safe_allocate(offsets(1:offset_size, 1:this%nprojector_matrices))
    safe_allocate(cnt(1:mesh%np))
 
    cnt = 0
 
    matrix_size = 0
    this%total_points = 0
    scal_size = 0
    this%max_npoints = 0
    this%max_nprojs = 0
    mix_offset = 0
    do imat = 1, this%nprojector_matrices
      pmat => this%projector_matrices(imat)
 
      this%max_npoints = max(this%max_npoints, pmat%npoints)
      this%max_nprojs = max(this%max_nprojs, pmat%nprojs)
 
      offsets(points, imat) = pmat%npoints
      offsets(projs, imat) = pmat%nprojs
 
      offsets(matrix, imat) = matrix_size
      matrix_size = matrix_size + pmat%npoints*pmat%nprojs
 
      offsets(map, imat) = this%total_points
      this%total_points = this%total_points + pmat%npoints
 
      offsets(scal, imat) = scal_size
      scal_size = scal_size + pmat%nprojs
 
      offsets(mix, imat) = mix_offset
      if (allocated(pmat%dmix)) then
        mix_offset = mix_offset + pmat%nprojs**2
      else if (allocated(pmat%zmix)) then
        mix_offset = mix_offset + 4*pmat%nprojs**2
      else
        offsets(mix, imat) = -1
      end if
 
      do is = 1, pmat%npoints
        ip = pmat%map(is)
        cnt(ip) = cnt(ip) + 1
      end do
    end do
 
    safe_allocate(invmap(1:max(maxval(cnt), 1), 1:mesh%np))
    safe_allocate(invmap2(1:max(maxval(cnt)*mesh%np, 1)))
    safe_allocate(invmap_mat(1:max(maxval(cnt), 1), 1:mesh%np))
    safe_allocate(invmap_mat2(1:max(maxval(cnt)*mesh%np, 1)))
    safe_allocate(pos(1:mesh%np + 1))
 
    cnt = 0
    ii = 0
    do imat = 1, this%nprojector_matrices
      pmat => this%projector_matrices(imat)
      do is = 1, pmat%npoints
        ip = pmat%map(is)
        cnt(ip) = cnt(ip) + 1
        invmap(cnt(ip), ip) = ii
        invmap_mat(cnt(ip), ip) = imat - 1
        ii = ii + 1
      end do
    end do
 
    ipos = 0
    pos(1) = 0
    do ip = 1, mesh%np
      do ii = 1, cnt(ip)
        ipos = ipos + 1
        invmap2(ipos) = invmap(ii, ip)
        invmap_mat2(ipos) = invmap_mat(ii, ip)
      end do
      pos(ip + 1) = ipos
    end do
 
    if (this%projector_matrices(1)%is_cmplx) then
      call accel_create_buffer(this%buff_matrices, accel_mem_read_only, type_cmplx, matrix_size)
    else
      call accel_create_buffer(this%buff_matrices, accel_mem_read_only, type_float, matrix_size)
    end if
    call accel_create_buffer(this%buff_maps, accel_mem_read_only, type_integer, this%total_points)
    if (mesh%use_vol_pp) then
      call accel_create_buffer(this%buff_sphere_weights, accel_mem_read_only, type_float, this%total_points)
    end if
    call accel_create_buffer(this%buff_position, accel_mem_read_only, type_float, 3*this%total_points)
    call accel_create_buffer(this%buff_scals, accel_mem_read_only, type_float, scal_size)
 
    if (mix_offset > 0) then
      if (allocated(this%projector_matrices(1)%zmix)) then
        call accel_create_buffer(this%buff_mix, accel_mem_read_only, type_cmplx, mix_offset)
      else
        call accel_create_buffer(this%buff_mix, accel_mem_read_only, type_float, mix_offset)
      end if
    end if
 
    do imat = 1, this%nprojector_matrices
      pmat => this%projector_matrices(imat)
      if (pmat%npoints > 0) then
        if (pmat%is_cmplx) then
          call accel_write_buffer(this%buff_matrices, pmat%npoints, pmat%nprojs, pmat%zprojectors, &
            offset = offsets(matrix, imat))
        else
          call accel_write_buffer(this%buff_matrices, pmat%npoints, pmat%nprojs, pmat%dprojectors, &
            offset = offsets(matrix, imat))
        end if
        call accel_write_buffer(this%buff_maps, pmat%npoints, pmat%map, offset = offsets(map, imat))
        if (mesh%use_vol_pp) then
          call accel_write_buffer(this%buff_sphere_weights, pmat%npoints, pmat%sphere_weights, &
            offset = offsets(map, imat))
        end if
        call accel_write_buffer(this%buff_position, space%dim, pmat%npoints, pmat%position, &
          offset = 3*offsets(map, imat))
      end if
      call accel_write_buffer(this%buff_scals, pmat%nprojs, pmat%scal, offset = offsets(scal, imat))
      if (offsets(mix, imat) /= -1) then
        if (allocated(pmat%zmix)) then
          call accel_write_buffer(this%buff_mix, pmat%nprojs, pmat%nprojs, 4, pmat%zmix, offset = offsets(mix, imat))
        else
          call accel_write_buffer(this%buff_mix, pmat%nprojs, pmat%nprojs, pmat%dmix, offset = offsets(mix, imat))
        end if
      end if
    end do
 
    call accel_create_buffer(this%buff_offsets, accel_mem_read_only, type_integer, offset_size*this%nprojector_matrices)
    call accel_write_buffer(this%buff_offsets, offset_size, this%nprojector_matrices, offsets)
 
    call accel_create_buffer(this%buff_pos, accel_mem_read_only, type_integer, mesh%np + 1)
    call accel_write_buffer(this%buff_pos, mesh%np + 1, pos)
 
    call accel_create_buffer(this%buff_invmap, accel_mem_read_only, type_integer, ipos)
    call accel_write_buffer(this%buff_invmap, ipos, invmap2)
 
    call accel_create_buffer(this%buff_invmap_mat, accel_mem_read_only, type_integer, ipos)
    call accel_write_buffer(this%buff_invmap_mat, ipos, invmap_mat2)
 
    ! 0-based matrix index of every packed projector point and every projector
    safe_allocate(point_to_mat(1:max(this%total_points, 1)))
    safe_allocate(proj_to_mat(1:max(scal_size, 1)))
    do imat = 1, this%nprojector_matrices
      pmat => this%projector_matrices(imat)
      point_to_mat(offsets(map, imat) + 1:offsets(map, imat) + pmat%npoints) = imat - 1
      proj_to_mat(offsets(scal, imat) + 1:offsets(scal, imat) + pmat%nprojs) = imat - 1
    end do
 
    call accel_create_buffer(this%buff_point_to_mat, accel_mem_read_only, type_integer, this%total_points)
    call accel_write_buffer(this%buff_point_to_mat, this%total_points, point_to_mat)
    call accel_create_buffer(this%buff_proj_to_mat, accel_mem_read_only, type_integer, scal_size)
    call accel_write_buffer(this%buff_proj_to_mat, scal_size, proj_to_mat)
 
    safe_deallocate_a(offsets)
    safe_deallocate_a(point_to_mat)
    safe_deallocate_a(proj_to_mat)
    safe_deallocate_a(cnt)
    safe_deallocate_a(invmap)
    safe_deallocate_a(invmap2)
    safe_deallocate_a(invmap_mat)
    safe_deallocate_a(invmap_mat2)
    safe_deallocate_a(pos)
 
    pop_sub(nonlocal_pseudopotential_build_accel_buffers)
  end subroutine nonlocal_pseudopotential_build_accel_buffers
 
  ! ----------------------------------------------------------------------------------
  subroutine nonlocal_pseudopotential_build_projector_phase_accel_buffer(this)
    class(nonlocal_pseudopotential_t), intent(inout) :: this
 
    integer :: ik, imat, iphase, nphase, offset, npoints
 
    push_sub(nonlocal_pseudopotential_build_projector_phase_accel_buffer)
 
    if (.not. allocated(this%projector_phases)) then
      pop_sub(nonlocal_pseudopotential_build_projector_phase_accel_buffer)
      return
    end if
 
    nphase = size(this%projector_phases, 2)
    this%nphase = nphase
 
    call accel_create_buffer(this%buff_projector_phases, accel_mem_read_only, type_cmplx, &
      this%total_points*nphase*size(this%projector_phases, 4))
 
    offset = 0
    do ik = lbound(this%projector_phases, 4), ubound(this%projector_phases, 4)
      do imat = 1, this%nprojector_matrices
        npoints = this%projector_matrices(imat)%npoints
        do iphase = 1, nphase
          if (npoints > 0) then
            call accel_write_buffer(this%buff_projector_phases, npoints, this%projector_phases(1:, iphase, imat, ik), &
              offset = offset, async=.true.)
          end if
          offset = offset + npoints
        end do
      end do
    end do
    call accel_finish()
 
    pop_sub(nonlocal_pseudopotential_build_projector_phase_accel_buffer)
  end subroutine nonlocal_pseudopotential_build_projector_phase_accel_buffer
 
  ! ----------------------------------------------------------------------------------
  !
  logical pure function nonlocal_pseudopotential_self_overlap(this) result(projector_self_overlap)
    class(nonlocal_pseudopotential_t), intent(in) :: this
 
    projector_self_overlap = this%projector_self_overlap
  end function nonlocal_pseudopotential_self_overlap
 
#include "undef.F90"
#include "real.F90"
#include "nonlocal_pseudopotential_inc.F90"
 
#include "undef.F90"
#include "complex.F90"
#include "nonlocal_pseudopotential_inc.F90"
 
end module nonlocal_pseudopotential_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: