forces.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 forces_oct_m
  use accel_oct_m
  use batch_oct_m
  use batch_ops_oct_m
  use born_charges_oct_m
  use boundaries_oct_m
  use comm_oct_m
  use debug_oct_m
  use density_oct_m
  use derivatives_oct_m
  use epot_oct_m
  use ext_partner_list_oct_m
  use gauge_field_oct_m
  use global_oct_m
  use grid_oct_m
  use hamiltonian_elec_oct_m
  use interaction_partner_oct_m
  use io_oct_m
  use ions_oct_m
  use, intrinsic :: iso_fortran_env
  use kpoints_oct_m
  use lalg_basic_oct_m
  use lasers_oct_m
  use lda_u_oct_m
  use linear_response_oct_m
  use math_oct_m
  use mesh_oct_m
  use mesh_function_oct_m
  use messages_oct_m
  use mpi_oct_m
  use namespace_oct_m
  use nonlocal_pseudopotential_oct_m
  use parser_oct_m
  use profiling_oct_m
  use projector_oct_m
  use ps_oct_m
  use pseudopotential_oct_m
  use space_oct_m
  use species_oct_m
  use species_pot_oct_m
  use states_abst_oct_m
  use states_elec_oct_m
  use states_elec_dim_oct_m
  use symm_op_oct_m
  use symmetrizer_oct_m
  use unit_oct_m
  use unit_system_oct_m
  use utils_oct_m
  use v_ks_oct_m
  use vdw_ts_oct_m
  use wfs_elec_oct_m
 
  implicit none
 
  private
  public ::                    &
    forces_calculate,          &
    dforces_from_potential,    &
    zforces_from_potential,    &
    dforces_derivative,        &
    zforces_derivative,        &
    dforces_born_charges,      &
    zforces_born_charges,      &
    total_force_calculate,     &
    forces_costate_calculate,  &
    forces_write_info
 
 
contains
 
  ! ---------------------------------------------------------
  subroutine total_force_calculate(space, gr, ions, ep, st, kpoints, x, lda_u)
    class(space_t),      intent(in)    :: space
    type(grid_t),        intent(in)    :: gr
    type(ions_t),        intent(in)    :: ions
    type(epot_t),        intent(in)    :: ep
    type(states_elec_t), intent(in)    :: st
    type(kpoints_t),     intent(in)    :: kpoints
    real(real64),        intent(inout) :: x(:)
    integer,             intent(in)    :: lda_u
 
 
    call profiling_in("FORCES")
    push_sub(total_force_calculate)
 
    x = m_zero
    if (states_are_real(st)) then
      call dtotal_force_from_potential(space, gr, ions, ep, st, kpoints, x, lda_u)
    else
      call ztotal_force_from_potential(space, gr, ions, ep, st, kpoints, x, lda_u)
    end if
 
    pop_sub(total_force_calculate)
    call profiling_out("FORCES")
  end subroutine total_force_calculate
 
  ! -------------------------------------------------------
 
  subroutine forces_costate_calculate(gr, namespace, ions, hm, psi, chi, ff, qq)
    type(grid_t),             intent(in)    :: gr
    type(namespace_t),        intent(in)    :: namespace
    type(ions_t),             intent(inout) :: ions
    type(hamiltonian_elec_t), intent(in)    :: hm
    type(states_elec_t),      intent(in)    :: psi
    type(states_elec_t),      intent(in)    :: chi
    real(real64),             intent(inout) :: ff(:, :)
    real(real64),             intent(in)    :: qq(:, :)
 
    integer :: jatom, idim, jdim
    integer, target :: j, ist, ik, iatom
    real(real64) :: rr, w2r, w1r, xx(ions%space%dim), dq, pdot3p, pdot3m, pdot3p2, pdot3m2, dforce1, dforce2
    complex(real64), allocatable :: zpsi(:, :), derpsi(:, :, :)
    real(real64), allocatable :: forceks1p(:), forceks1m(:), forceks1p2(:), forceks1m2(:), dforceks1(:)
 
    call profiling_in("FORCES")
    push_sub(forces_costate_calculate)
 
    ! FIXME: is the next section not basically the same as the routine ion_internamespace, action_calculate?
 
    ff = m_zero
    do iatom = 1, ions%natoms
      do jatom = 1, ions%natoms
        if (jatom == iatom) cycle
        xx(1:ions%space%dim) = ions%pos(:, jatom) - ions%pos(:, iatom)
        rr = norm2(xx(1:ions%space%dim))
        w1r = - ions%charge(iatom) * ions%charge(jatom) / rr**2
        w2r = m_two * ions%charge(iatom) * ions%charge(jatom) / rr**3
        do idim = 1, ions%space%dim
          do jdim = 1, ions%space%dim
            ff(iatom, idim) = ff(iatom, idim) + (qq(jatom, jdim) - qq(iatom, jdim)) * w2r * (m_one/rr**2) * xx(idim) * xx(jdim)
            ff(iatom, idim) = ff(iatom, idim) - (qq(jatom, jdim) - qq(iatom, jdim)) * w1r * (m_one/rr**3) * xx(idim) * xx(jdim)
            if (jdim == idim) then
              ff(iatom, idim) = ff(iatom, idim) + (qq(jatom, jdim) - qq(iatom, jdim)) * w1r * (m_one/rr)
            end if
          end do
        end do
      end do
    end do
 
    safe_allocate(derpsi(1:gr%np_part, 1:ions%space%dim, 1:psi%d%dim))
 
    dq = 0.001_real64
 
    safe_allocate(forceks1p(1:ions%space%dim))
    safe_allocate(forceks1m(1:ions%space%dim))
    safe_allocate(forceks1p2(1:ions%space%dim))
    safe_allocate(forceks1m2(1:ions%space%dim))
    safe_allocate(dforceks1(1:ions%space%dim))
    safe_allocate(zpsi(1:gr%np_part, 1:psi%d%dim))
 
    do ist = 1, psi%nst
      do ik = 1, psi%nik
        derpsi = m_z0
        call states_elec_get_state(psi, gr, ist, ik, zpsi)
        call zderivatives_grad(gr%der, zpsi(:, 1), derpsi(:, :, 1))
        do iatom = 1, ions%natoms
          do j = 1, ions%space%dim
            call force1(ions%pos(j, iatom) + dq, forceks1p, pdot3p)
            call force1(ions%pos(j, iatom) - dq, forceks1m, pdot3m)
            call force1(ions%pos(j, iatom) + dq/m_two, forceks1p2, pdot3p2)
            call force1(ions%pos(j, iatom) - dq/m_two, forceks1m2, pdot3m2)
            dforceks1 = ((m_four/m_three) * (forceks1p2 - forceks1m2) - (m_one / 6.0_real64) * (forceks1p - forceks1m)) / dq
            dforce1 = sum(qq(iatom, :) * dforceks1(:))
            dforce2 = ((m_four/m_three) * (pdot3p2 - pdot3m2) - (m_one / 6.0_real64) * (pdot3p - pdot3m)) / dq
            ff(iatom, j) = ff(iatom, j) - m_two * psi%occ(ist, ik) * dforce1 + m_two * dforce2
          end do
        end do
      end do
    end do
 
    safe_deallocate_a(zpsi)
    safe_deallocate_a(forceks1p)
    safe_deallocate_a(forceks1m)
    safe_deallocate_a(forceks1p2)
    safe_deallocate_a(forceks1m2)
    safe_deallocate_a(dforceks1)
    safe_deallocate_a(derpsi)
 
    pop_sub(forces_costate_calculate)
    call profiling_out("FORCES")
 
  contains
 
    subroutine force1(qq, res, pdot3)
      real(real64),      intent(in) :: qq
      real(real64), contiguous, intent(inout) :: res(:)
      real(real64),      intent(inout) :: pdot3
 
      integer :: m
      real(real64) :: qold
      complex(real64), allocatable :: viapsi(:, :), zpsi(:, :)
 
      qold = ions%pos(j, iatom)
      ions%pos(j, iatom) = qq
      safe_allocate(viapsi(1:gr%np_part, 1:psi%d%dim))
      safe_allocate(zpsi(1:gr%np_part, 1:psi%d%dim))
      viapsi = m_z0
      call states_elec_get_state(psi, gr, ist, ik, zpsi)
      call zhamiltonian_elec_apply_atom(hm, namespace, ions%space, ions%latt, ions%atom(iatom)%species, &
        ions%pos(:, iatom), iatom, gr, zpsi, viapsi)
 
      res(:) = m_zero
      do m = 1, ubound(res, 1)
        res(m) = real( zmf_dotp(gr, viapsi(:, 1), derpsi(:, m, 1), reduce = .false.), real64)
      end do
      if (gr%parallel_in_domains) call gr%allreduce(res)
 
      call states_elec_get_state(chi, gr, ist, ik, zpsi)
      pdot3 = real(m_zi * zmf_dotp(gr, zpsi(:, 1), viapsi(:, 1)), real64)
      ions%pos(j, iatom) = qold
 
      safe_deallocate_a(viapsi)
    end subroutine force1
 
  end subroutine forces_costate_calculate
  ! ---------------------------------------------------------
 
 
  ! ---------------------------------------------------------
  subroutine forces_calculate(gr, namespace, ions, hm, ext_partners, st, ks, vhxc_old, t, dt)
    type(grid_t),        intent(in)    :: gr
    type(namespace_t),   intent(in)    :: namespace
    type(ions_t),        intent(inout) :: ions
    type(hamiltonian_elec_t), intent(inout) :: hm
    type(partner_list_t),intent(in)    :: ext_partners
    type(states_elec_t), intent(inout) :: st
    type(v_ks_t),        intent(in)    :: ks
    real(real64),     optional, intent(in)    :: vhxc_old(:,:)
    real(real64),     optional, intent(in)    :: t
    real(real64),     optional, intent(in)    :: dt
 
    integer :: j, iatom, idir
    real(real64) :: xx(ions%space%dim), time, global_force(ions%space%dim)
    real(real64), allocatable :: force(:, :), force_loc(:, :), force_nl(:, :), force_u(:, :)
    real(real64), allocatable :: force_nlcc(: ,:)
    real(real64), allocatable :: force_scf(:, :)
    type(lasers_t), pointer :: lasers
 
    call profiling_in("FORCES")
    push_sub(forces_calculate)
 
    time = m_zero
    if (present(t)) time = t
 
    !We initialize the different components of the force to zero
    do iatom = 1, ions%natoms
      ions%atom(iatom)%f_ii(1:ions%space%dim) = m_zero
      ions%atom(iatom)%f_vdw(1:ions%space%dim) = m_zero
      ions%atom(iatom)%f_loc(1:ions%space%dim) = m_zero
      ions%atom(iatom)%f_nl(1:ions%space%dim) = m_zero
      ions%atom(iatom)%f_u(1:ions%space%dim) = m_zero
      ions%atom(iatom)%f_fields(1:ions%space%dim) = m_zero
      ions%atom(iatom)%f_nlcc(1:ions%space%dim) = m_zero
      ions%atom(iatom)%f_scf(1:ions%space%dim) = m_zero
      ions%atom(iatom)%f_photons(1:ions%space%dim) = m_zero
    end do
 
    ! the ion-ion and vdw terms are already calculated
    ! if we use vdw TS, we need to compute it now
    if (ks%vdw%vdw_correction == option__vdwcorrection__vdw_ts) then
      call vdw_ts_force_calculate(ks%vdw%vdw_ts, hm%ep%vdw_forces, ions, gr, st%d%nspin, st%rho)
    end if
 
    do iatom = 1, ions%natoms
      ions%tot_force(:, iatom) = hm%ep%fii(1:ions%space%dim, iatom) + hm%ep%vdw_forces(1:ions%space%dim, iatom)
      if (ks%has_photons) then
        ions%tot_force(:, iatom) = ions%tot_force(:, iatom) &
          - p_proton_charge*ions%charge(iatom)*hm%ep%photon_forces(1:ions%space%dim)
      end if
      ions%atom(iatom)%f_ii(1:ions%space%dim) = hm%ep%fii(1:ions%space%dim, iatom)
      ions%atom(iatom)%f_vdw(1:ions%space%dim) = hm%ep%vdw_forces(1:ions%space%dim, iatom)
      ions%atom(iatom)%f_photons(1:ions%space%dim) = - p_proton_charge*ions%charge(iatom)*hm%ep%photon_forces(1:ions%space%dim)
    end do
 
    if (present(t)) then
      global_force = ions%global_force(time)
 
      ! the ion-ion term is already calculated
      do iatom = 1, ions%natoms
        ions%tot_force(:, iatom) = ions%tot_force(:, iatom) + global_force
        ions%atom(iatom)%f_ii(1:ions%space%dim) = ions%atom(iatom)%f_ii(1:ions%space%dim) + global_force
      end do
    end if
 
    safe_allocate(force(1:ions%space%dim, 1:ions%natoms))
    safe_allocate(force_loc(1:ions%space%dim, 1:ions%natoms))
    safe_allocate(force_nl(1:ions%space%dim, 1:ions%natoms))
    safe_allocate(force_u(1:ions%space%dim, 1:ions%natoms))
    safe_allocate(force_scf(1:ions%space%dim, 1:ions%natoms))
    safe_allocate(force_nlcc(1:ions%space%dim, 1:ions%natoms))
 
    if (states_are_real(st)) then
      call dforces_from_potential(gr, namespace, ions%space, ions, hm, st, force, force_loc, force_nl, force_u)
    else
      call zforces_from_potential(gr, namespace, ions%space, ions, hm, st, force, force_loc, force_nl, force_u)
    end if
 
    if (allocated(st%rho_core)) then
      call forces_from_nlcc(gr, ions, st%d%spin_channels, hm%vxc, force_nlcc)
    else
      force_nlcc(:, :) = m_zero
    end if
    if (present(vhxc_old)) then
      call forces_from_scf(gr, ions, st%d%spin_channels, hm%vhxc, vhxc_old, force_scf)
    else
      force_scf = m_zero
    end if
 
    if (ions%force_total_enforce) then
      call forces_set_total_to_zero(ions, force)
      call forces_set_total_to_zero(ions, force_loc)
      call forces_set_total_to_zero(ions, force_nl)
      call forces_set_total_to_zero(ions, force_u)
      call forces_set_total_to_zero(ions, force_scf)
      call forces_set_total_to_zero(ions, force_nlcc)
    end if
 
    do iatom = 1, ions%natoms
      do idir = 1, ions%space%dim
        ions%tot_force(idir, iatom) = ions%tot_force(idir, iatom) + force(idir, iatom) &
          + force_scf(idir, iatom) + force_nlcc(idir, iatom)
        ions%atom(iatom)%f_loc(idir) = force_loc(idir, iatom)
        ions%atom(iatom)%f_nl(idir) = force_nl(idir, iatom)
        ions%atom(iatom)%f_u(idir) = force_u(idir, iatom)
        ions%atom(iatom)%f_nlcc(idir) = force_nlcc(idir, iatom)
        ions%atom(iatom)%f_scf(idir) = force_scf(idir, iatom)
      end do
    end do
 
    safe_deallocate_a(force)
    safe_deallocate_a(force_loc)
    safe_deallocate_a(force_nl)
    safe_deallocate_a(force_u)
    safe_deallocate_a(force_nlcc)
    safe_deallocate_a(force_scf)
 
    !\todo forces due to the magnetic fields (static and time-dependent)
    lasers => list_get_lasers(ext_partners)
    if (present(t) .and. associated(lasers)) then
      do j = 1, lasers%no_lasers
        select case (laser_kind(lasers%lasers(j)))
        case (e_field_electric)
          xx(:) = m_zero
          call laser_field(lasers%lasers(j), xx, t)
          do iatom = 1, ions%natoms
            ! Here the proton charge is +1, since the electric field has the usual sign.
            ions%tot_force(:, iatom) = ions%tot_force(:, iatom) + ions%charge(iatom)*xx(:)
            ions%atom(iatom)%f_fields(1:ions%space%dim) = ions%atom(iatom)%f_fields(1:ions%space%dim) &
              + ions%charge(iatom)*xx(:)
          end do
 
        case (e_field_vector_potential)
          ! Forces are correctly calculated only if the time-dependent
          ! vector potential has no spatial dependence.
          ! The full force taking account of the spatial dependence of A should be:
          ! F = q [- dA/dt + v x \nabla x A]
 
          !TODO: Add the gauge-field here
          xx(:) = m_zero
          call laser_electric_field(lasers%lasers(j), xx(:), t, dt) !convert in E field (E = -dA/ c dt)
          do iatom = 1, ions%natoms
            ! Also here the proton charge is +1
            ions%tot_force(:, iatom) = ions%tot_force(:, iatom) + ions%charge(iatom)*xx(:)
            ions%atom(iatom)%f_fields(1:ions%space%dim) = ions%atom(iatom)%f_fields(1:ions%space%dim) &
              + ions%charge(iatom)*xx(:)
          end do
 
          if (lasers_with_nondipole_field(lasers)) then
            write(message(1),'(a)') 'The forces are currently not supported for nondipole '
            write(message(2),'(a)') 'strong-field approximation Hamiltonian approach.'
            call messages_fatal(2, namespace=namespace)
          end if
 
        case (e_field_magnetic, e_field_scalar_potential)
          write(message(1),'(a)') 'The forces are currently not properly calculated if time-dependent'
          write(message(2),'(a)') 'magnetic fields are present.'
          call messages_fatal(2, namespace=namespace)
        end select
      end do
    end if
 
    if (allocated(hm%ep%e_field)) then
      do iatom = 1, ions%natoms
        ! Here the proton charge is +1, since the electric field has the usual sign.
        ions%tot_force(:, iatom) = ions%tot_force(:, iatom) + ions%charge(iatom)*hm%ep%e_field(1:ions%space%dim)
        ions%atom(iatom)%f_fields(1:ions%space%dim) = ions%atom(iatom)%f_fields(1:ions%space%dim) &
          + ions%charge(iatom)*hm%ep%e_field(1:ions%space%dim)
      end do
    end if
 
    if (allocated(hm%ep%b_field) .or. allocated(hm%ep%a_static)) then
      write(message(1),'(a)') 'The forces are currently not properly calculated if static'
      write(message(2),'(a)') 'magnetic fields or static vector potentials are present.'
      call messages_fatal(2, namespace=namespace)
    end if
 
    if (list_has_gauge_field(ext_partners)) then
      write(message(1),'(a)') 'The forces are currently not properly calculated if gauge-field'
      write(message(2),'(a)') 'is applied.'
      call messages_fatal(2, namespace=namespace)
    end if
 
    ! For periodic systems, we substract the averaged force
    if(ions%space%is_periodic()) then
      call forces_set_total_to_zero(ions, ions%tot_force)
    end if
 
    pop_sub(forces_calculate)
    call profiling_out("FORCES")
 
  end subroutine forces_calculate
 
  ! ----------------------------------------------------------------------
 
  subroutine forces_set_total_to_zero(ions, force)
    type(ions_t),        intent(in)    :: ions
    real(real64),        intent(inout) :: force(:, :)
 
    real(real64), allocatable :: total_force(:)
    integer :: iatom
 
    push_sub(forces_set_total_to_zero)
 
    safe_allocate(total_force(1:ions%space%dim))
 
    total_force(1:ions%space%dim) = sum(force, dim=2) / ions%natoms
 
    do iatom = 1, ions%natoms
      force(1:ions%space%dim, iatom) = force(1:ions%space%dim, iatom) - total_force(1:ions%space%dim)
    end do
 
    safe_deallocate_a(total_force)
    pop_sub(forces_set_total_to_zero)
  end subroutine forces_set_total_to_zero
 
  ! ----------------------------------------------------------------------
  subroutine forces_compute_total_torque(ions, total_torque)
    type(ions_t),        intent(in)    :: ions
    real(real64),        intent(inout) :: total_torque(:)
 
    real(real64) :: center_of_mass(ions%space%dim), rr(3), ff(3)
    integer :: iatom
 
    push_sub(forces_compute_total_torque)
 
    center_of_mass = ions%center_of_mass()
 
    total_torque = m_zero
    rr = m_zero
    ff = m_zero
    do iatom = 1, ions%natoms
      rr(1:ions%space%dim) = ions%pos(1:ions%space%dim, iatom) - center_of_mass
      ff(1:ions%space%dim) = ions%tot_force(1:ions%space%dim, iatom)
      total_torque(1:3) = total_torque(1:3) + dcross_product(rr, ff)
    end do
 
    pop_sub(forces_compute_total_torque)
  end subroutine forces_compute_total_torque
 
 
  ! ----------------------------------------------------------------------
 
  subroutine forces_write_info(iunit, ions, dir, namespace)
    integer,             intent(in)    :: iunit
    type(ions_t),        intent(in)    :: ions
    character(len=*),    intent(in)    :: dir
    type(namespace_t),   intent(in)    :: namespace
 
    integer :: iatom, idir, ii, iunit2
    real(real64) :: torque(1:3)
 
    if (.not. mpi_grp_is_root(mpi_world)) return
 
    push_sub(forces_write_info)
 
    write(iunit,'(3a)') 'Forces on the ions [', trim(units_abbrev(units_out%force)), "]"
    write(iunit,'(a,10x,99(14x,a))') ' Ion', (index2axis(idir), idir = 1, ions%space%dim)
    do iatom = 1, ions%natoms
      write(iunit,'(i4,a10,10es17.8)') iatom, trim(ions%atom(iatom)%species%get_label()), &
        units_from_atomic(units_out%force, ions%tot_force(:, iatom))
    end do
    write(iunit,'(1x,100a1)') ("-", ii = 1, 13 + ions%space%dim * 15)
    write(iunit,'(a14, 10es17.8)') " Max abs force", &
      (units_from_atomic(units_out%force, maxval(abs(ions%tot_force(idir, 1:ions%natoms)))), idir=1, ions%space%dim)
    write(iunit,'(a14, 10es17.8)') " Total force", &
      (units_from_atomic(units_out%force, sum(ions%tot_force(idir, 1:ions%natoms))), idir=1, ions%space%dim)
 
    if (ions%space%dim == 2 .or. ions%space%dim == 3 .and. .not. ions%space%is_periodic()) then
      call forces_compute_total_torque(ions, torque)
      write(iunit,'(a14, 10es17.8)') ' Total torque', &
        (units_from_atomic(units_out%force*units_out%length, torque(idir)), idir = 1, 3)
    end if
    write(iunit,'(1x)')
 
 
    iunit2 = io_open(trim(dir)//'/forces', namespace, action='write', position='asis')
    write(iunit2,'(a)') &
      ' # Total force (x,y,z) Ion-Ion (x,y,z) VdW (x,y,z) Local (x,y,z) NL (x,y,z)' // &
      ' Fields (x,y,z) Hubbard(x,y,z) SCF(x,y,z) NLCC(x,y,z) Phot (x,y,z)'
    do iatom = 1, ions%natoms
      write(iunit2,'(i4,a10,30es17.8)') iatom, trim(ions%atom(iatom)%species%get_label()), &
        (units_from_atomic(units_out%force, ions%tot_force(idir, iatom)), idir=1, ions%space%dim), &
        (units_from_atomic(units_out%force, ions%atom(iatom)%f_ii(idir)), idir=1, ions%space%dim), &
        (units_from_atomic(units_out%force, ions%atom(iatom)%f_vdw(idir)), idir=1, ions%space%dim), &
        (units_from_atomic(units_out%force, ions%atom(iatom)%f_loc(idir)), idir=1, ions%space%dim), &
        (units_from_atomic(units_out%force, ions%atom(iatom)%f_nl(idir)), idir=1, ions%space%dim), &
        (units_from_atomic(units_out%force, ions%atom(iatom)%f_fields(idir)), idir=1, ions%space%dim), &
        (units_from_atomic(units_out%force, ions%atom(iatom)%f_u(idir)), idir=1, ions%space%dim), &
        (units_from_atomic(units_out%force, ions%atom(iatom)%f_scf(idir)), idir=1, ions%space%dim), &
        (units_from_atomic(units_out%force, ions%atom(iatom)%f_nlcc(idir)), idir=1, ions%space%dim), &
        (units_from_atomic(units_out%force, ions%atom(iatom)%f_photons(idir)), idir=1, ions%space%dim)
    end do
    call io_close(iunit2)
 
    pop_sub(forces_write_info)
 
  end subroutine forces_write_info
 
  ! ----------------------------------------------------------------------
  ! This routine add the contribution to the forces from the nonlinear core correction
  ! see Eq. 9 of Kronik et al., J. Chem. Phys. 115, 4322 (2001)
  subroutine forces_from_nlcc(mesh, ions, spin_channels, vxc, force_nlcc)
    class(mesh_t),                  intent(in)    :: mesh
    type(ions_t),                   intent(inout) :: ions
    integer,                        intent(in)    :: spin_channels
    real(real64),                   intent(in)    :: vxc(:,:)
    real(real64), contiguous,              intent(out)   :: force_nlcc(:, :)
 
    integer :: is, iatom, idir
    real(real64), allocatable :: drho(:,:)
 
 
    push_sub(forces_from_nlcc)
 
    call profiling_in("FORCES_NLCC")
 
    safe_allocate(drho(1:mesh%np, 1:ions%space%dim))
 
    force_nlcc = m_zero
 
    do iatom = ions%atoms_dist%start, ions%atoms_dist%end
      call species_get_nlcc_grad(ions%atom(iatom)%species, ions%space, ions%latt, ions%pos(:, iatom), mesh, drho)
 
      do idir = 1, ions%space%dim
        do is = 1, spin_channels
          force_nlcc(idir, iatom) = force_nlcc(idir, iatom) &
            - dmf_dotp(mesh, drho(:,idir), vxc(1:mesh%np, is), reduce = .false.)/spin_channels
        end do
      end do
    end do
 
    safe_deallocate_a(drho)
 
    if (ions%atoms_dist%parallel) call dforces_gather(ions, force_nlcc)
 
    if (mesh%parallel_in_domains) then
      call profiling_in("FORCES_COMM")
      call mesh%allreduce(force_nlcc)
      call profiling_out("FORCES_COMM")
    end if
 
    call profiling_out("FORCES_NLCC")
 
    pop_sub(forces_from_nlcc)
  end subroutine forces_from_nlcc
 
  ! Implementation of the term from Chan et al.,  Phys. Rev. B 47, 4771 (1993).
  ! Here we make the approximation that the "atomic densities" are just the one
  ! from the pseudopotential.
  ! NTD : No idea if this is good or bad, but this is easy to implement
  !       and works well in practice
  subroutine forces_from_scf(mesh, ions, spin_channels, vhxc, vhxc_old, force_scf)
    class(mesh_t),                  intent(in)    :: mesh
    type(ions_t),                   intent(inout) :: ions
    integer,                        intent(in)    :: spin_channels
    real(real64),                   intent(in)    :: vhxc(:,:)
    real(real64),                   intent(in)    :: vhxc_old(:,:)
    real(real64), contiguous,              intent(out)   :: force_scf(:, :)
 
    integer :: is, iatom, idir
    real(real64), allocatable :: dvhxc(:,:), drho(:,:,:)
 
    push_sub(forces_from_scf)
 
    call profiling_in("FORCES_SCF")
 
    safe_allocate(dvhxc(1:mesh%np, 1:spin_channels))
    safe_allocate(drho(1:mesh%np, 1:spin_channels, 1:ions%space%dim))
 
    !We average over spin channels
    do is = 1, spin_channels
      dvhxc(1:mesh%np, is) = vhxc(1:mesh%np, is) - vhxc_old(1:mesh%np, is)
    end do
 
    force_scf = m_zero
 
    do iatom = ions%atoms_dist%start, ions%atoms_dist%end
      select type(spec=>ions%atom(iatom)%species)
      type is(pseudopotential_t)
 
        if (ps_has_density(spec%ps)) then
 
          call species_atom_density_grad(spec, ions%namespace, ions%space, ions%latt, &
            ions%pos(:, iatom), mesh, spin_channels, drho)
 
          do idir = 1, ions%space%dim
            do is = 1, spin_channels
              force_scf(idir, iatom) = force_scf(idir, iatom) &
                - dmf_dotp(mesh, drho(:,is,idir), dvhxc(:,is), reduce = .false.)
            end do
          end do
        end if
      end select
    end do
 
    safe_deallocate_a(dvhxc)
    safe_deallocate_a(drho)
 
    if (ions%atoms_dist%parallel) call dforces_gather(ions, force_scf)
 
    if (mesh%parallel_in_domains) then
      call profiling_in("FORCES_COMM")
      call mesh%allreduce(force_scf)
      call profiling_out("FORCES_COMM")
    end if
 
    call profiling_out("FORCES_SCF")
 
    pop_sub(forces_from_scf)
  end subroutine forces_from_scf
 
!---------------------------------------------------------------------------
  subroutine total_force_from_local_potential(mesh, space, vpsl, gdensity, force)
    class(mesh_t),                  intent(in)    :: mesh
    class(space_t),                 intent(in)    :: space
    real(real64),                   intent(in)    :: vpsl(:)
    real(real64),                   intent(in)    :: gdensity(:, :)
    real(real64),                   intent(inout) :: force(:)
 
    integer            :: idir
    real(real64)       :: force_tmp(1:space%dim)
 
    push_sub(total_force_from_local_potential)
 
    do idir = 1, space%dim
      force_tmp(idir) = dmf_dotp(mesh, vpsl(1:mesh%np), gdensity(:, idir), reduce = .false.)
    end do
 
    if (mesh%parallel_in_domains) call mesh%allreduce(force_tmp)
 
    force(1:space%dim) = force(1:space%dim) + force_tmp(1:space%dim)
 
    pop_sub(total_force_from_local_potential)
  end subroutine total_force_from_local_potential
 
 
#include "undef.F90"
#include "real.F90"
#include "forces_inc.F90"
 
#include "undef.F90"
#include "complex.F90"
#include "forces_inc.F90"
 
end module forces_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: