orbitalset.F90

Go to the documentation of this file.

!! Copyright (C) 2016 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 orbitalset_oct_m
  use accel_oct_m
  use batch_oct_m
  use batch_ops_oct_m
  use blas_oct_m
  use comm_oct_m
  use debug_oct_m
  use distributed_oct_m
  use global_oct_m
  use hardware_oct_m
  use ions_oct_m
  use kpoints_oct_m
  use lalg_basic_oct_m
  use math_oct_m
  use mesh_oct_m
  use mesh_function_oct_m
  use messages_oct_m
  use poisson_oct_m
  use profiling_oct_m
  use species_oct_m
  use submesh_oct_m
  use types_oct_m
  use wfs_elec_oct_m
 
  implicit none
 
  private
 
  public ::                         &
    orbitalset_t,                   &
    orbitalset_init,                &
    orbitalset_end,                 &
    orbitalset_update_phase,        &
    orbitalset_update_phase_shift,  &
    dorbitalset_get_coefficients,   &
    zorbitalset_get_coefficients,   &
    dorbitalset_get_coeff_batch,    &
    zorbitalset_get_coeff_batch,    &
    dorbitalset_get_position_matrix_elem, &
    zorbitalset_get_position_matrix_elem, &
    dorbitalset_add_to_batch,       &
    zorbitalset_add_to_batch,       &
    orbitalset_set_jln
 
  type orbitalset_t
    ! Components are public by default
    integer             :: nn, ll, ii
    real(real64)        :: jj
    integer             :: norbs
    integer             :: ndim
    integer             :: iatom
    type(submesh_t)     :: sphere
    complex(real64), allocatable  :: phase(:,:)
    !                                         !< if the sphere cross the border of the box
    real(real64)        :: Ueff
    real(real64)        :: Ubar, Jbar
    real(real64)        :: alpha
    integer             :: nneighbors
    !                                         !< interaction is considered
    real(real64), allocatable  :: V_ij(:,:)
    real(real64), allocatable  :: coulomb_IIJJ(:,:,:,:,:)
    complex(real64), allocatable  :: zcoulomb_IIJJ(:,:,:,:,:,:,:)
    integer, allocatable:: map_os(:)
    complex(real64), allocatable  :: phase_shift(:,:)
 
    real(real64)        :: radius
    class(species_t), pointer :: spec => null()
    integer             :: spec_index
 
    real(real64), allocatable  :: dorb(:,:,:)
    complex(real64), allocatable  :: zorb(:,:,:)
    complex(real64), allocatable  :: eorb_submesh(:,:,:,:)
    complex(real64), allocatable  :: eorb_mesh(:,:,:,:)
    integer             :: ldorbs
    type(accel_mem_t)   :: dbuff_orb, zbuff_orb
    type(accel_mem_t), allocatable   :: buff_eorb (:)
 
    logical             :: use_submesh
 
    type(poisson_t)     :: poisson
  contains
    procedure :: set_species_index => orbitalset_set_species_index
  end type orbitalset_t
 
contains
 
  subroutine orbitalset_init(this)
    type(orbitalset_t),             intent(inout) :: this
 
    push_sub(orbitalset_init)
 
    this%iatom = -1
    this%nneighbors = 0
    this%nn = 0
    this%ll = 0
    this%jj = m_one
    this%ii = 0
    this%iatom = 0
    this%ndim = 1
    this%spec_index = 0
    this%norbs = 0
 
    this%Ueff = m_zero
    this%Ubar = m_zero
    this%Jbar = m_zero
    this%alpha = m_zero
    this%radius = m_zero
 
    pop_sub(orbitalset_init)
  end subroutine orbitalset_init
 
 
  subroutine orbitalset_end(this)
    type(orbitalset_t), intent(inout) :: this
 
    integer :: ik
 
    push_sub(orbitalset_end)
 
    safe_deallocate_a(this%phase)
    safe_deallocate_a(this%dorb)
    safe_deallocate_a(this%zorb)
    safe_deallocate_a(this%eorb_submesh)
    safe_deallocate_a(this%eorb_mesh)
    nullify(this%spec)
    call submesh_end(this%sphere)
 
    safe_deallocate_a(this%V_ij)
    safe_deallocate_a(this%coulomb_IIJJ)
    safe_deallocate_a(this%map_os)
    safe_deallocate_a(this%phase_shift)
 
    if(accel_is_enabled()) then
      call accel_release_buffer(this%dbuff_orb)
      call accel_release_buffer(this%zbuff_orb)
      if(allocated(this%buff_eorb)) then
        do ik = lbound(this%buff_eorb, dim=1), ubound(this%buff_eorb, dim=1)
          call accel_release_buffer(this%buff_eorb(ik))
        end do
        safe_deallocate_a(this%buff_eorb)
      end if
    end if
 
    pop_sub(orbitalset_end)
  end subroutine orbitalset_end
 
  subroutine orbitalset_set_jln(this, jj, ll, nn)
    type(orbitalset_t), intent(inout) :: this
    real(real64),       intent(in)    :: jj
    integer,            intent(in)    :: ll, nn
 
    push_sub(orbitalset_set_jln)
 
    this%jj = jj
    this%ll = ll
    this%nn = nn
 
    pop_sub(orbitalset_set_jln)
  end subroutine orbitalset_set_jln
 
 
  subroutine orbitalset_update_phase(os, dim, kpt, kpoints, spin_polarized, vec_pot, vec_pot_var, kpt_max)
    type(orbitalset_t),            intent(inout) :: os
    integer,                       intent(in)    :: dim
    type(distributed_t),           intent(in)    :: kpt
    type(kpoints_t),               intent(in)    :: kpoints
    logical,                       intent(in)    :: spin_polarized
    real(real64), optional,  allocatable, intent(in)    :: vec_pot(:)
    real(real64), optional,  allocatable, intent(in)    :: vec_pot_var(:, :)
    integer, optional,             intent(in)    :: kpt_max
 
    integer :: ns, iq, is, ikpoint, im, idim, kpt_end
    real(real64)   :: kr, kpoint(1:dim), dx(1:dim)
    integer :: iorb
 
    push_sub(orbitalset_update_phase)
 
    ns = os%sphere%np
 
    kpt_end = kpt%end
    if (present(kpt_max)) kpt_end = min(kpt_max, kpt_end)
 
    do iq = kpt%start, kpt_end
      !This is durty but avoids to refer to states_get_kpoint_index
      if (spin_polarized) then
        ikpoint = 1 + (iq - 1)/2
      else
        ikpoint = iq
      end if
 
      ! if this fails, it probably means that sb is not compatible with std
      assert(ikpoint <= kpoints_number(kpoints))
 
      kpoint(1:dim) = kpoints%get_point(ikpoint)
 
      !$omp parallel do private(dx, kr)
      do is = 1, ns
        ! this is only the correction to the global phase, that can
        ! appear if the sphere crossed the boundary of the cell.
        dx = os%sphere%rel_x(1:dim, is) - os%sphere%mesh%x(os%sphere%map(is), 1:dim) + os%sphere%center(1:dim)
        kr = dot_product(kpoint, dx)
        if (present(vec_pot)) then
          if (allocated(vec_pot)) kr = kr + dot_product(vec_pot, dx)
        end if
 
        if (present(vec_pot_var)) then
          if (allocated(vec_pot_var)) kr = kr + sum(vec_pot_var(1:dim, os%sphere%map(is)) &
            *(os%sphere%rel_x(1:dim, is)+os%sphere%center))
        end if
 
        os%phase(is, iq) = exp(m_zi*kr)
      end do
 
      if (.not. os%use_submesh) then
        ! We now compute the so-called Bloch sum of the localized orbitals
        do idim = 1, os%ndim
          do im = 1, os%norbs
            os%eorb_mesh(:, im, idim, iq) = m_z0
            do is = 1, ns
              os%eorb_mesh(os%sphere%map(is),im,idim,iq) = os%eorb_mesh(os%sphere%map(is),im,idim,iq) &
                + os%zorb(is, idim, im) * os%phase(is, iq)
            end do
          end do
        end do
      else ! In the case of the isolated system, we still use the submesh
        !$omp parallel private(idim, is)
        do im = 1, os%norbs
          do idim = 1, os%ndim
            !$omp do simd
            do is = 1, ns
              os%eorb_submesh(is,idim,im,iq) = os%zorb(is,idim,im)*os%phase(is, iq)
            end do
            !$omp end do simd
          end do
        end do
        !$omp end parallel
      end if
 
      if(accel_is_enabled() .and. os%ndim == 1) then
        if(os%use_submesh) then
          do iorb = 1, os%norbs
            call accel_write_buffer(os%buff_eorb(iq), ns, os%eorb_submesh(:, 1, iorb, iq), offset = (iorb - 1)*os%ldorbs)
          end do
        else
          do iorb = 1, os%norbs
            call accel_write_buffer(os%buff_eorb(iq), os%sphere%mesh%np, &
              os%eorb_mesh(:, iorb, 1, iq), offset = (iorb - 1)*os%ldorbs)
          end do
        end if
      end if
    end do
 
    pop_sub(orbitalset_update_phase)
  end subroutine orbitalset_update_phase
 
 
  subroutine orbitalset_update_phase_shift(os, dim, kpt, kpoints, spin_polarized, vec_pot, vec_pot_var, kpt_max)
    type(orbitalset_t),            intent(inout) :: os
    integer,                       intent(in)    :: dim
    type(distributed_t),           intent(in)    :: kpt
    type(kpoints_t),               intent(in)    :: kpoints
    logical,                       intent(in)    :: spin_polarized
    real(real64), optional,  allocatable, intent(in)    :: vec_pot(:)
    real(real64), optional,  allocatable, intent(in)    :: vec_pot_var(:, :)
    integer, optional,             intent(in)    :: kpt_max
 
    integer :: iq, ikpoint
    real(real64)   :: kr, kpoint(dim), dx(dim)
    integer :: inn, kpt_end
 
    push_sub(orbitalset_update_phase_shift)
 
    kpt_end = kpt%end
    if(present(kpt_max)) kpt_end = min(kpt_max, kpt_end)
 
    do iq = kpt%start, kpt_end
      !This is durty but avoids to refer to states_get_kpoint_index
      if(spin_polarized) then
        ikpoint = 1 + (iq - 1)/2
      else
        ikpoint = iq
      end if
 
      ! if this fails, it probably means that sb is not compatible with std
      assert(ikpoint <= kpoints_number(kpoints))
 
      kpoint(1:dim) = kpoints%get_point(ikpoint)
 
      if (os%nneighbors > 0) then
        do inn = 1, os%nneighbors
          dx(1:dim) = os%V_ij(inn,1:dim)
          kr = sum(kpoint(1:dim)*dx(1:dim))
          if (present(vec_pot)) then
            if (allocated(vec_pot)) kr = kr + sum(vec_pot(1:dim)*dx(1:dim))
          end if
 
          !At the moment the uniform vector potential is in vec_pot_var
          if (present(vec_pot_var)) then
            if (allocated(vec_pot_var))  kr = kr + sum(vec_pot_var(1:dim, 1)*dx(1:dim))
          end if
 
          !The sign is different as this is applied on the wavefunction and not the orbitals
          os%phase_shift(inn, iq) = exp(-m_zi*kr)
        end do
      end if
    end do
 
 
    pop_sub(orbitalset_update_phase_shift)
  end subroutine orbitalset_update_phase_shift
 
  ! ---------------------------------------------------------
  subroutine orbitalset_set_species_index(os, ions)
    class(orbitalset_t),     intent(inout) :: os
    type(ions_t),            intent(in)    :: ions
 
    integer :: ja
 
    push_sub(orbitalset_set_species_index)
 
    os%spec_index = os%spec%get_index()
    do ja = 1, ions%natoms
      if(ions%atom(ja)%species == os%spec) then
        os%spec_index = min(os%spec_index, ions%atom(ja)%species%get_index())
      end if
    end do
 
    pop_sub(orbitalset_set_species_index)
  end subroutine orbitalset_set_species_index
 
#include "undef.F90"
#include "real.F90"
#include "orbitalset_inc.F90"
 
#include "undef.F90"
#include "complex.F90"
#include "orbitalset_inc.F90"
 
end module orbitalset_oct_m