states_elec_calc.F90

Go to the documentation of this file.

!! Copyright (C) 2002-2006 M. Marques, A. Castro, A. Rubio, G. Bertsch
!!
!! 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 states_elec_calc_oct_m
#ifndef NDEBUG
  use, intrinsic :: ieee_exceptions
#endif
  use accel_oct_m
  use accel_blas_oct_m
  use batch_oct_m
  use batch_ops_oct_m
  use blas_oct_m
  use blacs_oct_m
  use iso_c_binding
  use comm_oct_m
  use debug_oct_m
  use derivatives_oct_m
  use electron_space_oct_m
  use fourier_space_oct_m
  use global_oct_m
  use grid_oct_m
  use hardware_oct_m
  use kpoints_oct_m
  use lalg_adv_oct_m
  use lalg_basic_oct_m
  use math_oct_m
  use messages_oct_m
  use mesh_oct_m
  use mesh_batch_oct_m
  use mesh_function_oct_m
  use mpi_oct_m
  use mpi_lib_oct_m
  use namespace_oct_m
  use par_vec_oct_m
  use pblas_oct_m
  use poisson_oct_m
  use profiling_oct_m
  use scalapack_oct_m
  use singularity_oct_m
  use space_oct_m
  use states_abst_oct_m
  use states_elec_oct_m
  use states_elec_dim_oct_m
  use states_elec_parallel_oct_m
  use types_oct_m
  use wfs_elec_oct_m
 
  implicit none
 
  private
 
  public ::                              &
    states_elec_calc_norms,              &
    states_elec_orthogonalize,           &
    states_elec_rotate,                  &
    dstates_elec_calc_orth_test,         &
    zstates_elec_calc_orth_test,         &
    dstates_elec_orthogonalization,      &
    zstates_elec_orthogonalization,      &
    dstates_elec_orthogonalize_single,   &
    zstates_elec_orthogonalize_single,   &
    dstates_elec_orthogonalize_single_batch,   &
    zstates_elec_orthogonalize_single_batch,   &
    dstates_elec_orthogonalization_full, &
    zstates_elec_orthogonalization_full, &
    dstates_elec_residue,                &
    zstates_elec_residue,                &
    dstates_elec_matrix,                 &
    zstates_elec_matrix,                 &
    dstates_elec_calc_overlap,           &
    zstates_elec_calc_overlap,           &
    dstates_elec_calc_projections,       &
    zstates_elec_calc_projections,       &
    dstates_elec_rrqr_decomposition,     &
    zstates_elec_rrqr_decomposition,     &
    states_elec_calc_orbital_difference
 
  interface states_elec_rotate
    module procedure dstates_elec_rotate, zstates_elec_rotate
  end interface states_elec_rotate
 
contains
 
  ! ---------------------------------------------------------
  subroutine states_elec_orthogonalize(st, namespace, mesh)
    type(states_elec_t),  intent(inout) :: st
    type(namespace_t),    intent(in)    :: namespace
    class(mesh_t),        intent(in)    :: mesh
 
    integer :: ik
 
    push_sub(states_elec_orthogonalize)
 
    do ik = st%d%kpt%start, st%d%kpt%end
      if (states_are_real(st)) then
        call dstates_elec_orthogonalization_full(st, namespace, mesh, ik)
      else
        call zstates_elec_orthogonalization_full(st, namespace, mesh, ik)
      end if
    end do
 
    pop_sub(states_elec_orthogonalize)
  end subroutine states_elec_orthogonalize
 
  ! ---------------------------------------------------------
  !
  subroutine states_elec_calc_norms(grid, kpoints, st, norm_ks)
    type(grid_t),             intent(in)    :: grid
    type(kpoints_t),          intent(in)    :: kpoints
    type(states_elec_t),      intent(in)    :: st
    real(real64), contiguous,        intent(out)   :: norm_ks(:, :)
 
    integer :: ik_ispin, iblock, minst, maxst
 
    push_sub(states_elec_calc_norms)
 
    ! If state or k is parallel, norm_ks is only partially touched, then gets summed over processes,
    ! hence initialise to zero
    norm_ks = m_zero
 
    do ik_ispin = st%d%kpt%start, st%d%kpt%end
      do iblock = st%group%block_start, st%group%block_end
        minst = states_elec_block_min(st, iblock)
        maxst = states_elec_block_max(st, iblock)
        ! Compute norm, defering grid communication until outside the loops
        call mesh_batch_nrm2(grid, st%group%psib(iblock, ik_ispin), norm_ks(minst:maxst, ik_ispin), reduce=.false.)
      enddo
    enddo
 
    ! Reduction over domains
    if (grid%parallel_in_domains) then
      norm_ks = norm_ks**2
      call grid%allreduce(norm_ks)
      norm_ks = sqrt(norm_ks)
    end if
 
    ! Reduction over states and k-points
    call comm_allreduce(st%st_kpt_mpi_grp, norm_ks)
 
    pop_sub(states_elec_calc_norms)
 
  end subroutine states_elec_calc_norms
 
  function states_elec_calc_orbital_difference(mesh, st, psib1, psib2) result(diff)
    class(mesh_t),         intent(in) :: mesh
    type(states_elec_t),   intent(in) :: st
    type(wfs_elec_t),      intent(in) :: psib1(:, st%d%kpt%start:)
    type(wfs_elec_t),      intent(in) :: psib2(:, st%d%kpt%start:)
 
    real(real64) :: diff
 
    integer :: ik, ib, minst, maxst
    type(wfs_elec_t) :: delta_psi
    real(real64), allocatable :: norm_k(:)
 
    push_sub(states_elec_calc_orbital_difference)
 
    safe_allocate(norm_k(1:maxval(st%group%block_size)))
 
    diff = m_zero
    do ik = st%d%kpt%start, st%d%kpt%end
      do ib = st%group%block_start, st%group%block_end
        minst = states_elec_block_min(st, ib)
        maxst = states_elec_block_max(st, ib)
        call psib2(ib, ik)%copy_to(delta_psi, copy_data=.true.)
        call batch_axpy(mesh%np, -m_one, psib1(ib, ik), delta_psi)
        call mesh_batch_nrm2(mesh, delta_psi, norm_k, reduce=.false.)
        diff = diff + st%kweights(ik) * sum(st%occ(minst:maxst, ik) * norm_k(1:psib1(ib, ik)%nst)**2)
        call delta_psi%end()
      enddo
    enddo
    safe_deallocate_a(norm_k)
 
    call comm_allreduce(st%dom_st_kpt_mpi_grp, diff)
    diff = sqrt(diff / st%qtot)
 
    pop_sub(states_elec_calc_orbital_difference)
 
  end function states_elec_calc_orbital_difference
 
 
#include "undef.F90"
#include "real.F90"
#include "states_elec_calc_inc.F90"
 
#include "undef.F90"
#include "complex.F90"
#include "states_elec_calc_inc.F90"
#include "undef.F90"
 
end module states_elec_calc_oct_m
 
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: