stress.F90

Go to the documentation of this file.

!! Copyright (C) 2002-2016 M. Marques, A. Castro, A. Rubio, G. Bertsch
!! Copyright (C) 2023 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 stress_oct_m
  use batch_ops_oct_m
  use boundaries_oct_m
  use comm_oct_m
  use debug_oct_m
  use density_oct_m
  use derivatives_oct_m
  use electron_space_oct_m
  use energy_oct_m
  use energy_calc_oct_m
  use epot_oct_m
  use global_oct_m
  use grid_oct_m
  use hamiltonian_elec_oct_m
  use io_oct_m
  use interaction_partner_oct_m
  use ions_oct_m
  use ion_interaction_oct_m
  use, intrinsic :: iso_fortran_env
  use kpoints_oct_m
  use ks_potential_oct_m
  use lalg_basic_oct_m
  use lda_u_oct_m
  use loct_math_oct_m
  use magnetic_constrain_oct_m
  use math_oct_m
  use mesh_oct_m
  use mesh_function_oct_m
  use mesh_batch_oct_m
  use messages_oct_m
  use mpi_oct_m
  use namespace_oct_m
  use nonlocal_pseudopotential_oct_m
  use poisson_oct_m
  use profiling_oct_m
  use ps_oct_m
  use pseudopotential_oct_m
  use space_oct_m
  use species_oct_m
  use species_pot_oct_m
  use species_factory_oct_m
  use splines_oct_m
  use states_elec_oct_m
  use states_elec_dim_oct_m
  use submesh_oct_m
  use symmetries_oct_m
  use symmetrizer_oct_m
  use types_oct_m
  use unit_oct_m
  use unit_system_oct_m
  use v_ks_oct_m
  use wfs_elec_oct_m
  use xc_oct_m
  use xc_f03_lib_m
  use xc_vdw_oct_m, only: d3_lib_options
  implicit none
 
  private
  public ::                    &
    stress_calculate,          &
    output_stress,             &
    output_pressure
 
contains
 
  ! ---------------------------------------------------------
  subroutine stress_calculate(namespace, gr, hm, st, ions, ks, ext_partners)
    type(namespace_t),           intent(in)    :: namespace
    type(grid_t),                intent(inout) :: gr
    type(hamiltonian_elec_t),    intent(inout) :: hm
    type(states_elec_t), target, intent(inout) :: st
    type(ions_t),                intent(inout) :: ions
    type(v_ks_t),                intent(in)    :: ks
    type(partner_list_t),        intent(in)    :: ext_partners
 
    real(real64), allocatable    :: rho_total(:)
    real(real64) :: stress(3,3) ! stress tensor in Cartesian coordinate
    real(real64) :: stress_kin(3,3), stress_Hartree(3,3), stress_xc(3,3), stress_xc_nlcc(3,3)
    real(real64) :: stress_ps(3,3), stress_ps_nl(3,3), stress_ps_local(3,3), stress_ii(3,3)
    real(real64) :: stress_hubbard(3,3)
    integer :: ip
    real(real64), allocatable :: vh(:)
    real(real64), allocatable :: grad_vh(:,:)
    real(real64) :: ehartree
    real(real64), contiguous, pointer :: rho(:)
 
    call profiling_in("STRESS_CALCULATE")
    push_sub(stress_calculate)
 
    if (st%wfs_type /= type_cmplx) then
      write(message(1),'(a)') 'The stress tensors for real wavefunctions has not been implemented!'
 
      if (hm%kpoints%full%npoints == 1) then
        write(message(2),'(a)') 'For testing this feature, you can add ForceComplex=yes to the input file'
        call messages_fatal(2, namespace=namespace)
      end if
 
      call messages_fatal(1, namespace=namespace)
    end if
 
    if (ions%space%periodic_dim == 1) then
      call messages_not_implemented("Stress tensor for 1D periodic systems", namespace=namespace)
    end if
 
    if (.not. ions%space%is_periodic()) then
      write(message(1),'(a)') 'The stress tensor cannot be computed for isolated systems'
      call messages_fatal(1, namespace=namespace)
    end if
 
    if (ks%vdw%vdw_correction /= option__vdwcorrection__none .and. .not. any(ks%vdw%vdw_correction == d3_lib_options)) then
      write(message(1),'(a)') 'The stress tensor is currently only implemented with DFT-D3 vdW correction'
      call messages_fatal(1, namespace=namespace)
    end if
 
    if (hm%pcm%run_pcm) then
      call messages_not_implemented('Stress tensor with PCM')
    end if
 
    if (allocated(hm%v_static)) then
      call messages_not_implemented('Stress tensor with static electric fields')
    end if
 
    if (ks%has_photons) then
      call messages_not_implemented('Stress tensor with photon modes')
    end if
 
    if (.not. hm%vnl%apply_projector_matrices) then
      call messages_not_implemented('Stress tensor with relativistic Kleinman-Bylander pseudopotential')
    end if
 
    if (hm%ep%reltype == scalar_relativistic_zora .or. hm%ep%reltype == fully_relativistic_zora) then
      call messages_not_implemented('Stress tensor with ZORA')
    end if
 
    ! Checks for the xc part of KS-DFT and GKS-DFT
    if (ks%theory_level == kohn_sham_dft .or. ks%theory_level == generalized_kohn_sham_dft) then
      if (.not. xc_is_energy_functional(hm%xc)) then
        call messages_not_implemented("Stress tensor with xc functionals that are not energy functionals")
      end if
 
      if ( .not. in_family(hm%xc%family, [xc_family_lda, xc_family_gga])) then
        write(message(1),'(a)') 'The stress tensor computation is currently only possible at the Kohn-Sham DFT level'
        write(message(2),'(a)') 'with LDA and GGA functionals or for independent particles.'
        call messages_fatal(2, namespace=namespace)
      end if
 
      if (in_family(hm%xc%family, [xc_family_gga]) .and. st%d%ispin == spinors) then
        call messages_not_implemented("Stress tensor for GGAs with spinors", namespace=namespace)
      end if
    end if
 
    if (hm%magnetic_constrain%level /= constrain_none) then
      call messages_not_implemented("Stress tensor with MagneticConstrain /= constrain_none")
    end if
 
    stress(:,:) = m_zero
 
    safe_allocate(rho_total(1:gr%np_part))
    do ip = 1, gr%np
      rho_total(ip) = sum(st%rho(ip, 1:st%d%nspin))
    end do
 
    ! As we rely on some of the full energy components, we need to recompute it first
    ! TODO: We should restrict the components of the energy needed to be computed
    call energy_calc_total(namespace, ions%space, hm, gr, st, ext_partners, iunit = -1, full = .true.)
 
    ! In order to get the electrostatic part (Hartree and local pseudopotential part),
    ! we need to get the Hartree potential and its gradient
    safe_allocate(vh(1:gr%np_part))
    safe_allocate(grad_vh(1:gr%np, 1:gr%der%dim))
    if (ks%theory_level /= independent_particles) then
      call lalg_copy(gr%np, hm%ks_pot%vhartree, vh)
    else
      if (hm%d%spin_channels > 1) then
        safe_allocate(rho(1:gr%np_part))
        call lalg_copy(gr%np, st%rho(:,1), rho)
        call lalg_axpy(gr%np, m_one, st%rho(:,2), rho)
      else
        rho => st%rho(:,1)
      end if
      ! In the case of independent particles, we use the electron density without NLCC
      call dpoisson_solve(hm%psolver, ions%namespace, vh, rho, all_nodes = .true.)
      if (hm%d%spin_channels > 1) then
        safe_deallocate_p(rho)
      else
        nullify(rho)
      end if
    end if
    ehartree = hm%energy%hartree
    ! We also compute the gradient here
    call dderivatives_grad(gr%der, vh, grad_vh)
 
    ! We now compute the various contributions to the stress tensor
 
    ! Stress from kinetic energy of electrons
    call stress_from_kinetic(gr, ions%space, hm, st, gr%symm, ions%latt%rcell_volume, stress_kin)
    stress = stress + stress_kin
 
    if (ks%theory_level == independent_particles) then
      stress_hartree = m_zero
      stress_xc = m_zero
    else
      call stress_from_hartree(gr, ions%space, ions%latt%rcell_volume, vh, grad_vh, ehartree, stress_hartree)
      stress = stress + stress_hartree
 
      call stress_from_xc(hm%energy, ions%latt%rcell_volume, ions%space%periodic_dim, stress_xc)
 
      ! Nonlinear core correction contribution
      if (allocated(st%rho_core)) then
        call stress_from_xc_nlcc(ions%latt%rcell_volume, gr, st, ions, hm%ks_pot%vxc, stress_xc_nlcc)
        stress_xc = stress_xc + stress_xc_nlcc
      end if
      ! Adds the beyond LDA contribution to the stress tensor
      stress_xc = stress_xc + ks%stress_xc_gga / ions%latt%rcell_volume
      stress = stress + stress_xc
    end if
 
    call stress_from_pseudo_local(gr, st, hm, ions, rho_total, vh, grad_vh, stress_ps_local)
    stress_ps = stress_ps_local
    stress = stress + stress_ps_local
 
    safe_deallocate_a(vh)
    safe_deallocate_a(grad_vh)
 
    call stress_from_pseudo_nonloc(gr, st, hm, ions, stress_ps_nl)
    stress_ps = stress_ps + stress_ps_nl
    stress = stress + stress_ps_nl
 
    call stress_from_hubbard(namespace, gr, st, hm, ions%space, ions%latt%rcell_volume, stress_hubbard)
    stress = stress + stress_hubbard
 
    call ion_interaction_stress(ions%ion_interaction, ions%space, ions%latt, ions%atom, ions%natoms, ions%pos, stress_ii)
    stress = stress + stress_ii
    ! Stress from kinetic energy of ion
    ! Stress from ion-field interaction
 
    ! Sign changed to fit conventional definition
    stress = -stress
 
    st%stress_tensors%kinetic = stress_kin
    st%stress_tensors%Hartree = stress_hartree
    st%stress_tensors%xc = stress_xc
    st%stress_tensors%ps_local = stress_ps_local
    st%stress_tensors%ps_nl = stress_ps_nl
    st%stress_tensors%hubbard = stress_hubbard
    st%stress_tensors%ion_ion = stress_ii
 
    ! Stress contribution from vdW D3
    if (ks%vdw%vdw_correction /= option__vdwcorrection__none) then
      st%stress_tensors%vdw = hm%ep%vdw_stress
    else
      st%stress_tensors%vdw = m_zero
    end if
    stress = stress + st%stress_tensors%vdw
 
    ! Symmetrize the stress tensor if we use k-point symmetries
    if (hm%kpoints%use_symmetries) then
      call dsymmetrize_tensor_cart(gr%symm, stress, use_non_symmorphic=.true.)
    end if
    ! We guarantee that the matrix is truely symmetric. There could be small numerical assymetries after symmetrization
    call dsymmetrize_matrix(ions%space%periodic_dim, stress)
 
    st%stress_tensors%total = stress
 
    ! Some sumrule for validation
    ! Sumrule is -3P_{kin}\Omega = 2 E_{kin}
    st%stress_tensors%kinetic_sumrule = m_zero
    ! Sumrule is -3P_{Hartree}\Omega = E_{Hartree}
    st%stress_tensors%Hartree_sumrule = m_zero
    if(ions%space%periodic_dim == 3) then
      st%stress_tensors%kinetic_sumrule = (stress_kin(1,1) + stress_kin(2,2) + stress_kin(3,3))*ions%latt%rcell_volume
      st%stress_tensors%kinetic_sumrule = st%stress_tensors%kinetic_sumrule - m_two * hm%energy%kinetic
 
      st%stress_tensors%hartree_sumrule = (stress_hartree(1,1) + stress_hartree(2,2) + stress_hartree(3,3))*ions%latt%rcell_volume
      st%stress_tensors%hartree_sumrule = st%stress_tensors%hartree_sumrule - hm%energy%hartree
    end if
 
    safe_deallocate_a(rho_total)
 
    pop_sub(stress_calculate)
    call profiling_out("STRESS_CALCULATE")
  end subroutine stress_calculate
 
  ! -------------------------------------------------------
  subroutine stress_from_kinetic(gr, space, hm, st, symm, rcell_volume, stress_kin)
    type(grid_t),                   intent(in)    :: gr
    class(space_t),                 intent(in)    :: space
    type(hamiltonian_elec_t),       intent(in)    :: hm
    type(states_elec_t),            intent(inout) :: st
    type(symmetries_t),             intent(in)    :: symm
    real(real64),                   intent(in)    :: rcell_volume
    real(real64),                   intent(out)   :: stress_kin(3, 3)
 
    integer :: ik, ist, idir, jdir, ib, minst, maxst
    complex(real64), allocatable :: stress_l_block(:)
    type(wfs_elec_t) :: psib, gpsib(space%dim)
 
    call profiling_in("STRESS_FROM_KINETIC")
    push_sub(stress_from_kinetic)
 
    stress_kin(:,:) = m_zero
 
    safe_allocate(stress_l_block(1:st%block_size))
 
    do ik = st%d%kpt%start, st%d%kpt%end
      if (st%kweights(ik) <= m_epsilon) cycle
 
      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 hamiltonian_elec_copy_and_set_phase(hm, gr, st%d%kpt, st%group%psib(ib, ik), psib)
 
        ! calculate the gradient
        call zderivatives_batch_grad(gr%der, psib, gpsib, set_bc=.false.)
 
        ! Accumulate the result
        do idir = 1, space%periodic_dim
          do jdir = idir, space%periodic_dim
            call zmesh_batch_dotp_vector(gr, gpsib(idir), gpsib(jdir), stress_l_block)
 
            do ist = minst, maxst
              stress_kin(idir,jdir) = stress_kin(idir,jdir) &
                + st%kweights(ik) * st%occ(ist, ik) &
                * real(stress_l_block(ist - minst + 1), real64)
            end do
          end do
        end do
 
        do idir = 1, space%dim
          call gpsib(idir)%end()
        end do
        call psib%end()
 
      end do
    end do
 
    if (st%parallel_in_states .or. st%d%kpt%parallel) then
      call comm_allreduce(st%st_kpt_mpi_grp, stress_kin)
    end if
 
 
    ! Symmetrize the kinetic stress tensor
    call upper_triangular_to_hermitian(space%periodic_dim, stress_kin)
 
    ! Symmetrize the stress tensor if we use k-point symmetries
    if (hm%kpoints%use_symmetries) then
      call dsymmetrize_tensor_cart(symm, stress_kin, use_non_symmorphic=.true.)
    end if
 
    stress_kin = stress_kin / rcell_volume
 
    call profiling_out("STRESS_FROM_KINETIC")
    pop_sub(stress_from_kinetic)
  end subroutine stress_from_kinetic
 
  ! -------------------------------------------------------
  subroutine stress_from_hartree(gr, space, volume, vh, grad_vh, ehartree, stress_Hartree)
    type(grid_t),       intent(in)    :: gr
    class(space_t),     intent(in)    :: space
    real(real64),       intent(in)    :: volume
    real(real64),       intent(in)    :: vh(:)
    real(real64),       intent(in)    :: grad_vh(:,:)
    real(real64),       intent(in)    :: ehartree
    real(real64),       intent(out)   :: stress_Hartree(3, 3)
 
    integer :: idir, jdir
 
    call profiling_in("STRESS_FROM_HARTREE")
    push_sub(stress_from_hartree)
 
    stress_hartree(:,:) = m_zero
 
    do idir = 1, space%periodic_dim
      do jdir = idir, space%periodic_dim
        stress_hartree(idir, jdir) = -dmf_dotp(gr, grad_vh(:,idir), grad_vh(:, jdir))/m_four/m_pi
      end do
      stress_hartree(idir, idir) = stress_hartree(idir, idir) + ehartree
    end do
 
    call upper_triangular_to_hermitian(space%periodic_dim, stress_hartree)
 
    stress_hartree =  stress_hartree/volume
 
    call profiling_out("STRESS_FROM_HARTREE")
    pop_sub(stress_from_hartree)
  end subroutine stress_from_hartree
 
 
  ! -------------------------------------------------------
  !
  ! Note: We assume hm%energy%echange, correlation, and intnvxc
  ! have already been calculated somewhere else.
  subroutine stress_from_xc(energy, rcell_volume, periodic_dim, stress_xc)
    type(energy_t),           intent(in)    :: energy
    real(real64),             intent(in)    :: rcell_volume
    integer,                  intent(in)    :: periodic_dim
    real(real64),             intent(out)   :: stress_xc(3, 3)
 
    integer :: idir
 
    call profiling_in("STRESS_FROM_XC")
    push_sub(stress_from_xc)
 
    stress_xc = m_zero
    do idir = 1, periodic_dim
      stress_xc(idir, idir) = - energy%exchange - energy%correlation + energy%intnvxc
    end do
    stress_xc(:,:) = stress_xc(:,:) / rcell_volume
 
    call profiling_out("STRESS_FROM_XC")
    pop_sub(stress_from_xc)
  end subroutine stress_from_xc
 
 
  ! -------------------------------------------------------
  subroutine stress_from_xc_nlcc(rcell_volume, gr, st, ions, vxc, stress_xc_nlcc)
    real(real64),             intent(in)    :: rcell_volume
    type(grid_t),             intent(in)    :: gr
    type(states_elec_t),      intent(in)    :: st
    type(ions_t),             intent(in)    :: ions
    real(real64),             intent(in)    :: vxc(:,:)
    real(real64),             intent(out)   :: stress_xc_nlcc(3, 3)
 
    integer :: idir, jdir, iat
    real(real64), allocatable :: gnlcc(:,:), gnlcc_x(:,:,:), vxc_tot(:)
 
    call profiling_in("STRESS_FROM_XC_NLCC")
    push_sub(stress_from_xc_nlcc)
 
    assert(allocated(st%rho_core))
 
    stress_xc_nlcc = m_zero
 
    ! We first accumulate the contribution from all the pseudo-ions
    safe_allocate(gnlcc(gr%np, gr%der%dim))
    safe_allocate(gnlcc_x(gr%np, gr%der%dim, gr%der%dim))
    gnlcc_x = m_zero
    do iat = ions%atoms_dist%start, ions%atoms_dist%end
      assert(ions%atom(iat)%species%is_ps())
      call species_get_nlcc_grad(ions%atom(iat)%species, ions%space, ions%latt, &
        ions%pos(:,iat), gr, gnlcc, gnlcc_x)
    end do
    safe_deallocate_a(gnlcc)
 
    if (ions%atoms_dist%parallel) then
      call comm_allreduce(ions%atoms_dist%mpi_grp, gnlcc_x)
    end if
 
    ! Sum over spin of the xc potential
    safe_allocate(vxc_tot(1:gr%np))
    call lalg_copy(gr%np, vxc(:, 1), vxc_tot)
    if(st%d%nspin > 1) call lalg_axpy(gr%np, m_one, vxc(:, 2), vxc_tot)
 
    do idir = 1, ions%space%periodic_dim
      do jdir = idir, ions%space%periodic_dim
        stress_xc_nlcc(idir, jdir) = dmf_dotp(gr, vxc_tot, gnlcc_x(:,idir, jdir))
      end do
    end do
    safe_deallocate_a(vxc_tot)
    safe_deallocate_a(gnlcc_x)
 
    call upper_triangular_to_hermitian(ions%space%periodic_dim, stress_xc_nlcc)
 
    stress_xc_nlcc(:,:) = stress_xc_nlcc(:,:) / rcell_volume
 
    call profiling_out("STRESS_FROM_XC_NLCC")
    pop_sub(stress_from_xc_nlcc)
  end subroutine stress_from_xc_nlcc
 
  ! -------------------------------------------------------
  subroutine stress_from_pseudo_nonloc(gr, st, hm, ions, stress_ps_nl)
    type(grid_t),      target,           intent(in) :: gr
    type(states_elec_t),              intent(inout) :: st
    type(hamiltonian_elec_t),            intent(in) :: hm
    type(ions_t),                        intent(in) :: ions
    real(real64),                       intent(out) :: stress_ps_nl(3, 3)
 
    integer :: ik, ist, idir, jdir
    integer :: ib, minst, maxst
    type(wfs_elec_t) :: psib, rvnl_psib(3), gpsib(3)
    complex(real64), allocatable :: stress_tmp(:)
 
    call profiling_in("STRESS_FROM_PSEUDO_NL")
    push_sub(stress_from_pseudo_nonloc)
 
    assert(st%wfs_type == type_cmplx)
 
    safe_allocate(stress_tmp(1:st%block_size))
 
    stress_ps_nl = m_zero
 
    do ik = st%d%kpt%start, st%d%kpt%end
 
      if (st%kweights(ik) <= m_epsilon) cycle
 
      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 hamiltonian_elec_copy_and_set_phase(hm, gr, st%d%kpt, st%group%psib(ib, ik), psib)
 
        ! calculate the gradient
        call zderivatives_batch_grad(gr%der, psib, gpsib, set_bc=.false.)
 
 
        ! Get rV_NL |\psi> for all atoms
        do idir = 1, gr%der%dim
          call psib%copy_to(rvnl_psib(idir))
          call batch_set_zero(rvnl_psib(idir))
        end do
        call hm%vnl%zr_vn_local(gr, st%d, gr%der%boundaries%spiral, psib, rvnl_psib)
 
        do idir = 1, ions%space%periodic_dim
          do jdir = idir, ions%space%periodic_dim
            call zmesh_batch_dotp_vector(gr, gpsib(idir), rvnl_psib(jdir), stress_tmp)
 
            do ist = minst, maxst
              stress_ps_nl(idir, jdir) = stress_ps_nl(idir, jdir) &
                + m_two * st%kweights(ik) * st%occ(ist, ik) * real(stress_tmp(ist-minst+1), real64)
            end do
 
          end do
        end do
 
        do idir = 1, gr%der%dim
          call rvnl_psib(idir)%end()
          call gpsib(idir)%end()
        end do
        call psib%end()
      end do
    end do
 
    safe_deallocate_a(stress_tmp)
 
    if (st%parallel_in_states .or. st%d%kpt%parallel) then
      call comm_allreduce(st%st_kpt_mpi_grp, stress_ps_nl)
    end if
 
    ! Symmetrize the kinetic stress tensor
    call upper_triangular_to_hermitian(ions%space%periodic_dim, stress_ps_nl)
 
    ! Symmetrize the stress tensor if we use k-point symmetries
    if (hm%kpoints%use_symmetries) then
      call dsymmetrize_tensor_cart(gr%symm, stress_ps_nl, use_non_symmorphic=.true.)
    end if
 
    ! Add the nonlocal energy
    do idir = 1, ions%space%periodic_dim
      stress_ps_nl(idir, idir) = stress_ps_nl(idir, idir) + hm%energy%extern_non_local
    end do
 
    stress_ps_nl = stress_ps_nl/ions%latt%rcell_volume
 
    call profiling_out("STRESS_FROM_PSEUDO_NL")
    pop_sub(stress_from_pseudo_nonloc)
 
  end subroutine stress_from_pseudo_nonloc
 
 
  ! -------------------------------------------------------
  subroutine stress_from_pseudo_local(gr, st, hm, ions, rho_total, vh, grad_vh, stress_ps_local)
    type(grid_t),      target,           intent(in) :: gr
    type(states_elec_t),              intent(inout) :: st
    type(hamiltonian_elec_t),            intent(in) :: hm
    type(ions_t),                        intent(in) :: ions
    real(real64), contiguous,                intent(inout) :: rho_total(:)
    real(real64),                        intent(in) :: vh(:)
    real(real64),                        intent(in) :: grad_vh(:,:)
    real(real64),                       intent(out) :: stress_ps_local(3, 3)
 
 
    real(real64) :: stress_SR(3, 3), stress_LR(3, 3)
    real(real64) :: energy_ps_SR, charge, zi
    real(real64),  allocatable :: vloc(:), rvloc(:,:), rho_local_lr(:), rho_lr(:)
    real(real64),  allocatable :: grad_rho(:,:), rho_lr_x(:,:), vlr(:), grad_vlr(:,:)
    integer :: idir, jdir, iatom
    type(ps_t), pointer :: spec_ps
 
    call profiling_in("STRESS_FROM_PSEUDO_LOC")
    push_sub(stress_from_pseudo_local)
 
    ! calculate stress from short-range local pseudopotentials
    stress_sr = m_zero
 
    safe_allocate(vloc(1:gr%np))
    vloc = m_zero
    safe_allocate(rvloc(1:gr%np, 1:gr%der%dim))
    rvloc = m_zero
    do iatom = 1, ions%natoms
      call epot_local_pseudopotential_sr(gr, ions, iatom, vloc, rvloc)
    end do
    safe_deallocate_a(vloc)
 
    safe_allocate(grad_rho(1:gr%np,1:gr%der%dim))
    call dderivatives_grad(gr%der, rho_total, grad_rho)
 
    energy_ps_sr = hm%energy%extern_local
    do idir = 1, ions%space%periodic_dim
      do jdir = idir, ions%space%periodic_dim
        stress_sr(idir, jdir) = stress_sr(idir, jdir) &
          +dmf_dotp(gr, rvloc(:, jdir), grad_rho(:, idir))
      end do
      stress_sr(idir,idir) = stress_sr(idir,idir) + energy_ps_sr
    end do
 
    call upper_triangular_to_hermitian(ions%space%periodic_dim, stress_sr)
 
    stress_sr = stress_sr/ions%latt%rcell_volume
 
    safe_deallocate_a(rvloc)
    safe_deallocate_a(grad_rho)
 
 
    ! calculate stress from long-range local pseudopotentials
    stress_lr = m_zero
 
    ! We treat the long-range part of the local potential as the Hartree term
    ! We first sum the long range densities from atoms
    safe_allocate(rho_lr(1:gr%np_part))
    safe_allocate(rho_lr_x(1:gr%np, 1:gr%der%dim))
    rho_lr = m_zero
    rho_lr_x = m_zero
    safe_allocate(rho_local_lr(1:gr%np))
    do iatom = ions%atoms_dist%start, ions%atoms_dist%end
      assert(ions%atom(iatom)%species%is_ps())
      call species_get_long_range_density(ions%atom(iatom)%species, ions%namespace, ions%space, ions%latt, &
        ions%pos(:, iatom), gr, rho_local_lr, nlr_x=rho_lr_x)
 
      call lalg_axpy(gr%np, m_one, rho_local_lr, rho_lr)
    end do
    safe_deallocate_a(rho_local_lr)
 
    if (ions%atoms_dist%parallel) then
      call comm_allreduce(ions%atoms_dist%mpi_grp, rho_lr)
      call comm_allreduce(ions%atoms_dist%mpi_grp, rho_lr_x)
    end if
 
    do idir = 1, ions%space%periodic_dim
      do jdir = idir, ions%space%periodic_dim
        stress_lr(idir, jdir) = stress_lr(idir, jdir) + dmf_dotp(gr, rho_lr_x(:,jdir), grad_vh(:, idir))
      end do
    end do
    safe_deallocate_a(rho_lr_x)
 
    safe_allocate(vlr(1:gr%np_part))
    call dpoisson_solve(hm%psolver, ions%namespace, vlr, rho_lr, all_nodes = .true.)
    safe_deallocate_a(rho_lr)
 
    safe_allocate(grad_vlr(1:gr%np, 1:gr%der%dim))
    call dderivatives_grad(gr%der, vlr, grad_vlr)
    safe_deallocate_a(vlr)
 
    do idir = 1, ions%space%periodic_dim
      do jdir = idir, ions%space%periodic_dim
        stress_lr(idir, jdir) = stress_lr(idir, jdir) - dmf_dotp(gr, grad_vh(:,idir), grad_vlr(:, jdir))/m_two/m_pi
      end do
    end do
 
    call upper_triangular_to_hermitian(ions%space%periodic_dim, stress_lr)
 
    safe_deallocate_a(grad_vlr)
 
    ! Contribution from G=0 component of the long-range part
    !
    if (ions%space%periodic_dim == 3) then
      charge = m_zero
      do iatom = 1, ions%natoms
        charge = charge + ions%atom(iatom)%species%get_zval()
      end do
 
      do iatom = 1, ions%natoms
        select type(spec => ions%atom(iatom)%species)
        type is(pseudopotential_t)
          zi = spec%get_zval()
          spec_ps => spec%ps
 
          do idir = 1, ions%space%periodic_dim
            stress_lr(idir, idir) = stress_lr(idir, idir) &
              + m_two*m_pi*spec_ps%sigma_erf**2*charge*zi /ions%latt%rcell_volume
          end do
        end select
      end do
    end if
 
    stress_lr = stress_lr/ions%latt%rcell_volume
 
    stress_ps_local = stress_sr + stress_lr
 
    call profiling_out("STRESS_FROM_PSEUDO_LOC")
    pop_sub(stress_from_pseudo_local)
 
  end subroutine stress_from_pseudo_local
 
  ! -------------------------------------------------------
  subroutine epot_local_pseudopotential_sr(mesh, ions, iatom, vpsl, rvpsl)
    class(mesh_t),            intent(in)    :: mesh
    type(ions_t),             intent(in)    :: ions
    integer,                  intent(in)    :: iatom
    real(real64),             intent(inout) :: vpsl(:)
    real(real64),             intent(inout) :: rvpsl(:,:)
 
    integer :: ip
    real(real64) :: radius, vl_ip
    type(submesh_t)  :: sphere
    type(ps_t), pointer :: ps
 
    push_sub(epot_local_pseudopotential_sr)
 
    if (.not. ions%atom(iatom)%species%is_ps()) then
      pop_sub(epot_local_pseudopotential_sr)
      return
    endif
 
    call profiling_in("EPOT_LOCAL_PS_SR")
 
    select type(spec=>ions%atom(iatom)%species)
    type is(pseudopotential_t)
 
      ps => spec%ps
 
      radius = ps%vl%x_threshold*1.05_real64
 
      call submesh_init(sphere, ions%space, mesh, ions%latt, ions%pos(:, iatom), radius)
 
      ! Cannot be written (correctly) as a vector expression since for periodic systems,
      ! there can be values ip, jp such that sphere%map(ip) == sphere%map(jp).
      do ip = 1, sphere%np
        vl_ip = spline_eval(ps%vl, sphere%r(ip))
        vpsl(sphere%map(ip)) = vpsl(sphere%map(ip)) + vl_ip
        rvpsl(sphere%map(ip), 1:ions%space%periodic_dim) = rvpsl(sphere%map(ip), 1:ions%space%periodic_dim) &
          + sphere%rel_x(1:ions%space%periodic_dim, ip) * vl_ip
      end do
 
      call submesh_end(sphere)
 
      nullify(ps)
 
    end select
 
    call profiling_out("EPOT_LOCAL_PS_SR")
    pop_sub(epot_local_pseudopotential_sr)
  end subroutine epot_local_pseudopotential_sr
 
 
  ! -------------------------------------------------------
  subroutine stress_from_hubbard(namespace, gr, st, hm, space, rcell_volume, stress_hubbard)
    type(namespace_t),                   intent(in) :: namespace
    type(grid_t),      target,           intent(in) :: gr
    type(states_elec_t),              intent(inout) :: st
    type(hamiltonian_elec_t),            intent(in) :: hm
    type(space_t),                       intent(in) :: space
    real(real64),                        intent(in) :: rcell_volume
    real(real64),                       intent(out) :: stress_hubbard(3, 3)
 
    integer :: ik, ist, idir, jdir
    integer :: ib, minst, maxst
    type(wfs_elec_t) :: psib, rvu_psib(3), gpsib(3)
    complex(real64), allocatable :: stress_tmp(:)
 
    if (hm%lda_u%level == dft_u_none) then
      stress_hubbard = m_zero
      return
    end if
 
    push_sub_with_profile(stress_from_hubbard)
 
    assert(st%wfs_type == type_cmplx)
 
    safe_allocate(stress_tmp(1:st%block_size))
 
    stress_hubbard = m_zero
 
    do ik = st%d%kpt%start, st%d%kpt%end
 
      if (st%kweights(ik) <= m_epsilon) cycle
 
      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 hamiltonian_elec_copy_and_set_phase(hm, gr, st%d%kpt, st%group%psib(ib, ik), psib)
 
        ! calculate the gradient
        call zderivatives_batch_grad(gr%der, psib, gpsib, set_bc=.false.)
 
        ! Get rV_U |\psi> for all atoms
        do idir = 1, gr%der%dim
          call psib%copy_to(rvu_psib(idir))
          call batch_set_zero(rvu_psib(idir))
        end do
 
        call zlda_u_rvu(hm%lda_u, gr, space, hm%d, namespace, psib, rvu_psib)
 
        do idir = 1,3
          do jdir = idir,3
            call zmesh_batch_dotp_vector(gr, gpsib(idir), rvu_psib(jdir), stress_tmp)
 
            do ist = minst, maxst
              stress_hubbard(idir, jdir) = stress_hubbard(idir, jdir) &
                + m_two * st%kweights(ik) * st%occ(ist, ik) * real(stress_tmp(ist-minst+1), real64)
            end do
 
          end do
        end do
 
        do idir = 1, gr%der%dim
          call rvu_psib(idir)%end()
          call gpsib(idir)%end()
        end do
        call psib%end()
      end do
    end do
 
    safe_deallocate_a(stress_tmp)
 
    if (st%parallel_in_states .or. st%d%kpt%parallel) then
      call comm_allreduce(st%st_kpt_mpi_grp, stress_hubbard)
    end if
 
    ! Symmetrize the kinetic stress tensor
    call upper_triangular_to_hermitian(gr%der%dim, stress_hubbard)
 
    ! Symmetrize the stress tensor if we use k-point symmetries
    if (hm%kpoints%use_symmetries) then
      call dsymmetrize_tensor_cart(gr%symm, stress_hubbard)
    end if
 
    ! Add the Hubbard energy
    do idir = 1,3
      stress_hubbard(idir, idir) = stress_hubbard(idir, idir) + hm%energy%int_dft_u
    end do
 
    stress_hubbard = stress_hubbard/rcell_volume
 
    pop_sub_with_profile(stress_from_hubbard)
  end subroutine stress_from_hubbard
 
 
  ! -------------------------------------------------------
  subroutine output_stress(iunit, space_dim, stress_tensors, all_terms)
    integer,           intent(in) :: iunit
    integer,           intent(in) :: space_dim
    type(stress_t),    intent(in) :: stress_tensors
    logical, optional, intent(in) :: all_terms
 
    logical :: write_all_terms
    character(len=16) :: stress_unit
 
    write_all_terms = optional_default(all_terms, .true.)
 
    write(stress_unit, '(4a,i1)') trim(units_abbrev(units_out%energy)), '/', &
      trim(units_abbrev(units_out%length)), '^', space_dim
 
    if (mpi_grp_is_root(mpi_world)) then
 
      if (write_all_terms) then
        write(iunit, '(3a)') 'Kinetic stress tensor [', trim(stress_unit), '] ='
        call print_stress_tensor(iunit, space_dim, stress_tensors%kinetic)
        if (space_dim == 3) then
          write(iunit, '(a, es15.6, 3a)') 'Kinetic pressure sumrule violation: ', &
            units_from_atomic(units_out%energy, stress_tensors%kinetic_sumrule), &
            ' [', trim(units_abbrev(units_out%energy)), ']'
          write(iunit,*)
        end if
 
 
        write(iunit, '(3a)') 'Hartree stress tensor [', trim(stress_unit), '] ='
        call print_stress_tensor(iunit, space_dim, stress_tensors%Hartree)
        if (space_dim == 3) then
          write(iunit, '(a, es15.6, 3a)') 'Hartree pressure sumrule violation: ', &
            units_from_atomic(units_out%energy, stress_tensors%hartree_sumrule), &
            ' [', trim(units_abbrev(units_out%energy)), ']'
          write(iunit,*)
        end if
 
        write(iunit, '(3a)') 'XC stress tensor [', trim(stress_unit), '] ='
        call print_stress_tensor(iunit, space_dim, stress_tensors%xc)
 
        write(iunit, '(3a)') 'Local pseudo. stress tensor [', trim(stress_unit), '] ='
        call print_stress_tensor(iunit, space_dim, stress_tensors%ps_local)
 
        write(iunit, '(3a)') 'Nonlocal pseudo. stress tensor [', trim(stress_unit), '] ='
        call print_stress_tensor(iunit, space_dim, stress_tensors%ps_nl)
 
        write(iunit, '(3a)') 'Ion-ion stress tensor [', trim(stress_unit), '] ='
        call print_stress_tensor(iunit, space_dim, stress_tensors%ion_ion)
 
        write(iunit, '(3a)') 'vdW stress tensor [', trim(stress_unit), '] ='
        call print_stress_tensor(iunit, space_dim, stress_tensors%vdw)
 
        write(iunit, '(3a)') 'Hubbard stress tensor [', trim(stress_unit), '] ='
        call print_stress_tensor(iunit, space_dim, stress_tensors%hubbard)
      end if
 
      write(iunit, '(3a)') 'Total stress tensor [', trim(stress_unit), '] ='
      call print_stress_tensor(iunit, space_dim, stress_tensors%total)
 
    end if
  end subroutine output_stress
 
 
  subroutine output_pressure(iunit, space_dim, total_stress_tensor)
    integer, intent(in) :: iunit
    integer, intent(in) :: space_dim
    real(real64),   intent(in) :: total_stress_tensor(3,3)
 
    integer :: idim
    real(real64) :: pressure
    character(len=16) :: stress_unit
 
    write(stress_unit, '(4a,i1)') trim(units_abbrev(units_out%energy)), '/', &
      trim(units_abbrev(units_out%length)), '^', space_dim
 
    pressure = m_zero
    do idim = 1, space_dim
      pressure = pressure - total_stress_tensor(idim, idim) / real(space_dim, real64)
    end do
 
    write(iunit,'(3a,es16.8)', advance="no") 'Pressure [', trim(stress_unit), '] = ', &
      units_from_atomic(units_out%energy/units_out%length**space_dim, pressure)
    if (space_dim == 3) then
      write(iunit,'(2x,a,f16.8)') 'Pressure [GPa] = ', units_from_atomic(unit_gpa, pressure)
    else
      write(iunit,*)
    end if
 
  end subroutine output_pressure
 
  subroutine print_stress_tensor(ounit, space_dim, tensor)
    integer, intent(in) :: ounit
    integer, intent(in) :: space_dim
    real(real64),   intent(in) :: tensor(3,3)
 
    real(real64)   :: tensor_with_unit(3,3)
    integer :: idim, jdim
 
    tensor_with_unit = units_from_atomic(units_out%energy/units_out%length**space_dim, tensor)
 
    write(ounit,'(a9,2x)', advance="no")"T_{ij}"
    do jdim = 1, space_dim
      write(ounit,'(i18)', advance="no") jdim
    end do
    write(ounit,*)
    do idim = 1, space_dim
      write(ounit,'(i9,2x)', advance="no") idim
      do jdim = 1, space_dim
        write(ounit,'(es18.9)', advance="no") tensor_with_unit(idim, jdim)
      end do
      write(ounit,*)
    end do
    write(ounit,*)
 
  end subroutine print_stress_tensor
 
 
end module stress_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: