ions.F90

Go to the documentation of this file.

!! Copyright (C) 2002-2006 M. Marques, A. Castro, A. Rubio, G. Bertsch
!! Copyright (C) 2021 M. Oliveira
!!
!! 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 ions_oct_m
  use atom_oct_m
  use blas_oct_m
  use calc_mode_par_oct_m
  use charged_particles_oct_m
  use comm_oct_m
  use iso_c_binding
  use debug_oct_m
  use distributed_oct_m
  use global_oct_m
  use interaction_oct_m
  use interaction_surrogate_oct_m
  use io_binary_oct_m
  use io_oct_m
  use ion_interaction_oct_m
  use, intrinsic :: iso_fortran_env
  use jellium_oct_m
  use lalg_adv_oct_m
  use lattice_vectors_oct_m
  use math_oct_m
  use messages_oct_m
  use multicomm_oct_m
  use mpi_oct_m
  use namespace_oct_m
  use parser_oct_m
  use phonon_modes_oct_m
  use profiling_oct_m
  use pseudopotential_oct_m
  use read_coords_oct_m
  use space_oct_m
  use species_oct_m
  use species_factory_oct_m
  use spglib_f08
  use symm_op_oct_m
  use symmetries_oct_m
  use tdfunction_oct_m
  use unit_oct_m
  use unit_system_oct_m
  use varinfo_oct_m
  use wigner_distribution_oct_m
 
  implicit none
 
  private
  public :: ions_t
 
  type, extends(charged_particles_t) :: ions_t
    ! Components are public by default
 
    type(lattice_vectors_t) :: latt
 
    integer                   :: natoms
    type(atom_t), allocatable :: atom(:)
 
    type(symmetries_t) :: symm
 
    type(distributed_t) :: atoms_dist
 
    integer, allocatable :: map_symm_atoms(:,:)
    integer, allocatable :: inv_map_symm_atoms(:,:)
 
    real(real64), allocatable :: equilibrium_pos(:,:)
    !                                                 !! for multitrajectory runs.
    !                                                 !!
    !                                                 !! This is necessary, as the displacements are added before we
    !                                                 !! generate the grid, but the grid should always be constructed
    !                                                 !! for the equilibrium positions.
    real(real64), allocatable :: pos_displacements(:,:)
    real(real64), allocatable :: vel_displacements(:,:)
 
    ! Information about the species
    integer                                :: nspecies
    class(species_wrapper_t), allocatable  :: species(:)
    logical                                :: only_user_def
    logical,          private              :: species_time_dependent
 
    logical                 :: force_total_enforce
    type(ion_interaction_t) :: ion_interaction
 
    logical,     private :: apply_global_force
    type(tdf_t), private :: global_force_function
  contains
    procedure :: copy => ions_copy
    generic   :: assignment(=) => copy
    procedure :: partition => ions_partition
    procedure :: init_interaction => ions_init_interaction
    procedure :: initialize => ions_initialize
    procedure :: update_quantity => ions_update_quantity
    procedure :: init_interaction_as_partner => ions_init_interaction_as_partner
    procedure :: copy_quantities_to_interaction => ions_copy_quantities_to_interaction
    procedure :: fold_atoms_into_cell => ions_fold_atoms_into_cell
    procedure :: min_distance => ions_min_distance
    procedure :: has_time_dependent_species => ions_has_time_dependent_species
    procedure :: val_charge => ions_val_charge
    procedure :: dipole => ions_dipole
    procedure :: translate => ions_translate
    procedure :: rotate => ions_rotate
    procedure :: global_force => ions_global_force
    procedure :: write_xyz => ions_write_xyz
    procedure :: read_xyz => ions_read_xyz
    procedure :: write_poscar => ions_write_poscar
    procedure :: write_crystal => ions_write_crystal
    procedure :: write_bild_forces_file => ions_write_bild_forces_file
    procedure :: write_vtk_geometry => ions_write_vtk_geometry
    procedure :: update_lattice_vectors => ions_update_lattice_vectors
    procedure :: symmetrize_atomic_coord => ions_symmetrize_atomic_coord
    procedure :: generate_mapping_symmetric_atoms => ions_generate_mapping_symmetric_atoms
    procedure :: print_spacegroup => ions_print_spacegroup
    procedure :: current => ions_current
    procedure :: abs_current => ions_abs_current
    procedure :: init_random_displacements => ions_init_random_displacements
    procedure :: single_mode_displacements => ions_single_mode_displacements
    final :: ions_finalize
  end type ions_t
 
  interface ions_t
    procedure ions_constructor
  end interface ions_t
 
contains
 
  ! ---------------------------------------------------------
  function ions_constructor(namespace, print_info, latt_inp) result(ions)
    type(namespace_t),                 intent(in)    :: namespace
    logical,                 optional, intent(in)    :: print_info
    type(lattice_vectors_t), optional, intent(out)   :: latt_inp
    class(ions_t), pointer :: ions
 
    type(read_coords_info) :: xyz
    integer :: ia, ierr, idir
    character(len=100)  :: function_name
    real(real64) :: mindist
    real(real64), allocatable :: factor(:)
    integer, allocatable :: site_type(:)
    logical, allocatable :: spherical_site(:)
    real(real64), parameter :: threshold = 1e-5_real64
    type(species_factory_t) :: factory
    real(real64) :: T
    integer :: displacement_mode
    real(real64) :: amp_pos, amp_vel
 
    logical :: in_ensemble
 
    push_sub_with_profile(ions_constructor)
 
    allocate(ions)
 
    ions%namespace = namespace
 
    ions%space = space_t(namespace)
    in_ensemble = parse_is_defined(namespace, "NumberOfReplicas")
 
 
    call species_factory_init(factory, namespace)
 
    ! initialize geometry
    call read_coords_init(xyz)
 
    ! load positions of the atoms
    call read_coords_read('Coordinates', xyz, ions%space, namespace)
 
    if (xyz%n < 1) then
      message(1) = "Coordinates have not been defined."
      call messages_fatal(1, namespace=namespace)
    end if
 
    ! Initialize parent class
    call charged_particles_init(ions, xyz%n)
 
    ! copy information from xyz to ions
    ions%natoms = xyz%n
    safe_allocate(ions%atom(1:ions%natoms))
    do ia = 1, ions%natoms
      call atom_init(ions%atom(ia), ions%space%dim, xyz%atom(ia)%label)
      ions%pos(:,ia) = xyz%atom(ia)%x(1:ions%space%dim)
      if (bitand(xyz%flags, xyz_flags_move) /= 0) then
        ions%fixed(ia) = .not. xyz%atom(ia)%move
      end if
    end do
 
    if (allocated(xyz%latvec)) then
      ! Build lattice vectors from the XSF input
      ions%latt = lattice_vectors_t(namespace, ions%space, xyz%latvec)
    else
      ! Build lattice vectors from input file
      ions%latt = lattice_vectors_t(namespace, ions%space)
    end if
 
    ! Convert coordinates to Cartesian in case we have reduced coordinates
    if (xyz%source == read_coords_reduced) then
      do ia = 1, ions%natoms
        ions%pos(:, ia) = ions%latt%red_to_cart(ions%pos(:, ia))
      end do
    end if
 
    call read_coords_end(xyz)
 
    ! Save the positions read from the input before applying deviations from the initial_input
    safe_allocate_source_a(ions%equilibrium_pos, ions%pos)
 
    call ions_init_species(ions, factory, print_info=print_info)
 
    ! Set the masses and charges. This needs to be done after initializing the species.
    do ia = 1, ions%natoms
      ions%mass(ia) = ions%atom(ia)%species%get_mass()
      ions%charge(ia) = ions%atom(ia)%species%get_zval()
    end do
 
    !%Variable InitialDisplacementMode
    !%Type integer
    !%Default 0
    !%Section System
    !%Description
    !% When this variable is set and non-zero, the phonon modes file will be read, and ionic
    !% positions and velocities will be displaced according to a single phonon mode.
    !% The amplitudes can be defined with the variables:
    !% - InitialDisplacementsAmplitudePos
    !% - InitialDisplacementsAmplitudeVel
    !%End
    call parse_variable(namespace, 'InitialDisplacementMode', 0, displacement_mode)
 
    !%Variable InitialDisplacementAmplitudePos
    !%Type float
    !%Default 1.0
    !%Section System
    !%Description
    !% Amplitude for initial position displacement according to the phonon mode, defined by InitialDisplacementMode.
    !% The value corresponds to the normal distributed canonical coordinates.
    !%End
    call parse_variable(namespace, 'InitialDisplacementAmplitudePos', 1.0_real64, amp_pos)
 
    !%Variable InitialDisplacementAmplitudeVel
    !%Type float
    !%Default 1.0
    !%Section System
    !%Description
    !% Amplitude for initial velocity displacement according to the phonon mode, defined by InitialDisplacementMode.
    !% The value corresponds to the normal distributed canonical coordinates.
    !%End
    call parse_variable(namespace, 'InitialDisplacementAmplitudeVel', 1.0_real64, amp_vel)
 
    if (displacement_mode>0) then
 
      message(1) = "Random displacements are disabled when InitialDisplacementMode > 0."
      call messages_warning(1, namespace=namespace)
 
      safe_allocate(ions%pos_displacements(1:ions%space%dim, 1:ions%natoms))
      safe_allocate(ions%vel_displacements(1:ions%space%dim, 1:ions%natoms))
 
      call ions%single_mode_displacements(displacement_mode, amp_pos, amp_vel)
 
      ! apply initial displacements
      ions%pos = ions%pos + ions%pos_displacements
 
      ! note that the velocity displacement needs to be applied in ions%initialize()
 
    end if
 
    ! Now we can apply displacements to the positions, if requested.
    ! The random displacements are disabled if a single mode is specified,
    if(in_ensemble .and. displacement_mode==0) then
 
      !%Variable EnsembleTemperature
      !%Type float
      !%Default 0
      !%Section Multi-Trajectory
      !%Description
      !% The temperature for an ensemble calculation in Kelvin
      !%End
      call parse_variable(namespace, 'EnsembleTemperature', 0.0_real64, t)
 
      safe_allocate(ions%pos_displacements(1:ions%space%dim, 1:ions%natoms))
      safe_allocate(ions%vel_displacements(1:ions%space%dim, 1:ions%natoms))
      call ions%init_random_displacements(t)
 
      ! apply initial displacements
      ions%pos = ions%pos + ions%pos_displacements
 
      ! note that the velocity displacement needs to be applied in ions%initialize()
 
    end if
 
    call ions_fold_atoms_into_cell(ions)
    call distributed_nullify(ions%atoms_dist, ions%natoms)
 
    call messages_obsolete_variable(namespace, 'PDBClassical')
 
    if (present(latt_inp)) then
      ! The lattice as read from the input might be needed by some other part of the code, so we save it
      latt_inp = ions%latt
    end if
 
    ! Now that we have processed the atomic coordinates, we renormalize the
    ! lattice parameters along the non-periodic dimensions
    if (ions%space%has_mixed_periodicity()) then
      safe_allocate(factor(ions%space%dim))
      do idir = 1, ions%space%periodic_dim
        factor(idir) = m_one
      end do
      do idir = ions%space%periodic_dim + 1, ions%space%dim
        factor(idir) = m_one/norm2(ions%latt%rlattice(1:ions%space%dim, idir))
      end do
      call ions%latt%scale(factor)
      safe_deallocate_a(factor)
    end if
 
    ! Check that atoms are not too close
    if (ions%natoms > 1) then
      mindist = ions_min_distance(ions, real_atoms_only = .false.)
      if (mindist < threshold) then
        write(message(1), '(a)') "Some of the atoms seem to sit too close to each other."
        write(message(2), '(a)') "Please review your input files and the output geometry (in 'static/')."
        write(message(3), '(a, f12.6, 1x, a)') "Minimum distance = ", &
          units_from_atomic(units_out%length, mindist), trim(units_abbrev(units_out%length))
        call messages_warning(3, namespace=namespace)
 
        ! then write out the geometry, whether asked for or not in Output variable
        call io_mkdir(static_dir, namespace)
        call ions%write_xyz(trim(static_dir)//'/geometry')
      end if
 
      if (ions_min_distance(ions, real_atoms_only = .true.) < threshold) then
        message(1) = "It cannot be correct to run with physical atoms so close."
        call messages_fatal(1, namespace=namespace)
      end if
    end if
 
    !Initialize symmetries
    safe_allocate(spherical_site(1:ions%natoms))
    safe_allocate(site_type(1:ions%natoms))
    do ia = 1, ions%natoms
      select type(spec => ions%atom(ia)%species)
      type is(jellium_slab_t)
        spherical_site(ia) = .false.
      type is(species_user_defined_t)
        spherical_site(ia) = .false.
      type is(species_from_file_t)
        spherical_site(ia) = .false.
      type is(species_charge_density_t)
        spherical_site(ia) = .false.
      class default
        spherical_site(ia) = .true.
      end select
 
      site_type(ia) = ions%atom(ia)%species%get_index()
    end do
 
    ions%symm = symmetries_t(ions%namespace, ions%space, ions%latt, ions%natoms, ions%pos, site_type, spherical_site)
 
    safe_deallocate_a(spherical_site)
    safe_deallocate_a(site_type)
 
    ! Generate the mapping of symmetric atoms
    call ions%generate_mapping_symmetric_atoms()
 
    call ion_interaction_init(ions%ion_interaction, namespace, ions%space, ions%natoms)
 
    !%Variable ForceTotalEnforce
    !%Type logical
    !%Default no
    !%Section Hamiltonian
    !%Description
    !% (Experimental) If this variable is set to "yes", then the sum
    !% of the total forces will be enforced to be zero.
    !%End
    call parse_variable(namespace, 'ForceTotalEnforce', .false., ions%force_total_enforce)
    if (ions%force_total_enforce) call messages_experimental('ForceTotalEnforce', namespace=namespace)
 
    !%Variable TDGlobalForce
    !%Type string
    !%Section Time-Dependent
    !%Description
    !% If this variable is set, a global time-dependent force will be
    !% applied to the ions in the x direction during a time-dependent
    !% run. This variable defines the base name of the force, that
    !% should be defined in the <tt>TDFunctions</tt> block. This force
    !% does not affect the electrons directly.
    !%End
 
    if (parse_is_defined(namespace, 'TDGlobalForce')) then
 
      ions%apply_global_force = .true.
 
      call parse_variable(namespace, 'TDGlobalForce', 'none', function_name)
      call tdf_read(ions%global_force_function, namespace, trim(function_name), ierr)
 
      if (ierr /= 0) then
        call messages_write("You have enabled the GlobalForce option but Octopus could not find")
        call messages_write("the '"//trim(function_name)//"' function in the TDFunctions block.")
        call messages_fatal(namespace=namespace)
      end if
 
    else
 
      ions%apply_global_force = .false.
 
    end if
 
    call species_factory_end(factory)
 
    pop_sub_with_profile(ions_constructor)
  end function ions_constructor
 
  ! ---------------------------------------------------------
  subroutine ions_init_species(ions, factory, print_info)
    type(ions_t),            intent(inout) :: ions
    type(species_factory_t), intent(in)    :: factory
    logical, optional,       intent(in)    :: print_info
 
    logical :: print_info_, spec_user_defined
    integer :: i, j, k, ispin
    class(species_t), pointer :: spec
 
    push_sub_with_profile(ions_init_species)
 
    print_info_ = .true.
    if (present(print_info)) then
      print_info_ = print_info
    end if
    ! First, count the species
    ions%nspecies = 0
    atoms1:  do i = 1, ions%natoms
      do j = 1, i - 1
        if (atom_same_species(ions%atom(j), ions%atom(i))) cycle atoms1
      end do
      ions%nspecies = ions%nspecies + 1
    end do atoms1
 
    ! Allocate the species structure.
    allocate(ions%species(1:ions%nspecies))
 
    ! Now, read the data.
    k = 0
    ions%only_user_def = .true.
    atoms2: do i = 1, ions%natoms
      do j = 1, i - 1
        if (atom_same_species(ions%atom(j), ions%atom(i))) cycle atoms2
      end do
      k = k + 1
      ions%species(k)%s => factory%create_from_input(ions%namespace, ions%atom(j)%get_label(), k)
      ! these are the species which do not represent atoms
      select type(spec => ions%species(k)%s)
      class is(jellium_t)
 
      class default
        ions%only_user_def = .false.
      end select
 
      select type(spec => ions%species(k)%s)
      type is(pseudopotential_t)
        if (ions%space%dim /= 3) then
          message(1) = "Pseudopotentials may only be used with Dimensions = 3."
          call messages_fatal(1, namespace=ions%namespace)
        end if
 
      type is(jellium_slab_t)
        if (ions%space%is_periodic() .and. ions%space%periodic_dim /= 2) then
          message(1) = "Periodic jelium slab can only be used if PeriodicDim = 2"
          call messages_fatal(1, namespace=ions%namespace)
        end if
      end select
 
    end do atoms2
 
    ! Reads the spin components. This is read here, as well as in states_init,
    ! to be able to pass it to the pseudopotential initializations subroutine.
    call parse_variable(ions%namespace, 'SpinComponents', 1, ispin)
    if (.not. varinfo_valid_option('SpinComponents', ispin)) call messages_input_error(ions%namespace, 'SpinComponents')
    ispin = min(2, ispin)
 
    if (print_info_) then
      call messages_print_with_emphasis(msg="Species", namespace=ions%namespace)
    end if
    do i = 1, ions%nspecies
      spec => ions%species(i)%s
      call spec%build(ions%namespace, ispin, ions%space%dim, print_info=print_info_)
    end do
    if (print_info_) then
      call messages_print_with_emphasis(namespace=ions%namespace)
    end if
 
    !%Variable SpeciesTimeDependent
    !%Type logical
    !%Default no
    !%Section System::Species
    !%Description
    !% When this variable is set, the potential defined in the block <tt>Species</tt> is calculated
    !% and applied to the Hamiltonian at each time step. You must have at least one <tt>species_user_defined</tt>
    !% type of species to use this.
    !%End
    call parse_variable(ions%namespace, 'SpeciesTimeDependent', .false., ions%species_time_dependent)
    ! we must have at least one user defined species in order to have time dependency
    spec_user_defined = .false.
    do i = 1,ions%nspecies
      select type(spec=>ions%species(i)%s)
      type is(species_user_defined_t)
        spec_user_defined = .true.
      end select
    end do
    if (ions%species_time_dependent .and. .not. spec_user_defined) then
      call messages_input_error(ions%namespace, 'SpeciesTimeDependent')
    end if
 
    !  assign species
    do i = 1, ions%natoms
      do j = 1, ions%nspecies
        if (atom_same_species(ions%atom(i), ions%species(j)%s)) then
          call atom_set_species(ions%atom(i), ions%species(j)%s)
          exit
        end if
      end do
    end do
 
    pop_sub_with_profile(ions_init_species)
  end subroutine ions_init_species
 
  !--------------------------------------------------------------
  subroutine ions_copy(ions_out, ions_in)
    class(ions_t),     intent(out) :: ions_out
    class(ions_t),     intent(in)  :: ions_in
 
    push_sub(ions_copy)
 
    call charged_particles_copy(ions_out, ions_in)
 
    ions_out%latt = ions_in%latt
 
    ions_out%natoms = ions_in%natoms
    safe_allocate(ions_out%atom(1:ions_out%natoms))
    ions_out%atom = ions_in%atom
 
    ions_out%nspecies = ions_in%nspecies
    allocate(ions_out%species(1:ions_out%nspecies))
    ions_out%species = ions_in%species
 
    ions_out%only_user_def     = ions_in%only_user_def
 
    call distributed_copy(ions_in%atoms_dist, ions_out%atoms_dist)
 
    safe_allocate(ions_out%map_symm_atoms(1:ions_in%natoms, 1:ions_in%symm%nops + ions_in%symm%nops_nonsymmorphic))
    ions_out%map_symm_atoms = ions_in%map_symm_atoms
    safe_allocate(ions_out%inv_map_symm_atoms(1:ions_in%natoms, 1:ions_in%symm%nops + ions_in%symm%nops_nonsymmorphic))
    ions_out%inv_map_symm_atoms = ions_in%inv_map_symm_atoms
 
 
    pop_sub(ions_copy)
  end subroutine ions_copy
 
  ! ---------------------------------------------------------
  subroutine ions_partition(this, mc)
    class(ions_t),               intent(inout) :: this
    type(multicomm_t),           intent(in)    :: mc
 
    push_sub(ions_partition)
 
    call distributed_init(this%atoms_dist, this%natoms, mc%group_comm(p_strategy_states), "atoms")
 
    call ion_interaction_init_parallelization(this%ion_interaction, this%natoms, mc)
 
    call mpi_grp_init(this%grp, mc%master_comm)
 
    pop_sub(ions_partition)
  end subroutine ions_partition
 
  ! ---------------------------------------------------------
  subroutine ions_init_interaction(this, interaction)
    class(ions_t),        target, intent(inout) :: this
    class(interaction_t),         intent(inout) :: interaction
 
    push_sub(ions_init_interaction)
 
    select type (interaction)
    class default
      call charged_particles_init_interaction(this, interaction)
    end select
 
    pop_sub(ions_init_interaction)
  end subroutine ions_init_interaction
 
  ! ---------------------------------------------------------
  !
  subroutine ions_init_random_displacements(ions, T)
    class(ions_t), intent(inout) :: ions
    real(real64),  intent(in)    :: T
 
 
    integer(int64)    :: seed
 
    type(phonon_modes_t)         :: phonons
    type(wigner_distribution_t)  :: wigner
 
    real(real64), allocatable :: sigmaP(:)
    real(real64), allocatable :: muP(:)
    real(real64), allocatable :: sigmaQ(:)
    real(real64), allocatable :: muQ(:)
    real(real64), allocatable :: genP(:)
    real(real64), allocatable :: genQ(:)
 
 
    real(real64) :: beta_half
    real(real64), allocatable :: pos_displacements(:), vel_displacements(:)
    real(real64), allocatable :: l_m(:), prefactor(:)
 
    integer :: imode, idim, iatom, i, iline
    integer :: num_real_modes
 
    real(real64), parameter :: low_temperature_tolerance = 1.0e-6_real64
 
 
    push_sub(ions_init_random_displacements)
 
    ! create a unique seed for each replica, based on a hash of the full (unique) namespace string.
    seed = ions%namespace%get_hash32()
 
    message(1) = "Create initial random displacements for ions."
    call messages_info(1, namespace = ions%namespace)
 
    write(message(1), '("namespace = ",A)') ions%namespace%get()
    write(message(2), '("seed = ",I0)') seed
    call messages_info(2, namespace = ions%namespace, debug_only=.true.)
 
 
    ! initialize the phonons (read info from file)
    call phonons%init(ions%namespace, ions%space%dim, ions%natoms, ions%space%periodic_dim > 0)
 
    num_real_modes = phonons%num_modes
 
    if (num_real_modes == 0) then
 
      ions%pos_displacements = m_zero
      ions%vel_displacements = m_zero
 
      pop_sub(ions_init_random_displacements)
      return
    end if
 
 
    ! initialize wigner distribution for phonons%num_modes modes and with given seed
    call wigner%init(num_real_modes, seed)
 
    ! set up widths and shifts for the Wigner function
    safe_allocate(sigmap(1:num_real_modes))
    safe_allocate(sigmaq(1:num_real_modes))
    safe_allocate(mup(1:num_real_modes))
    safe_allocate(muq(1:num_real_modes))
    safe_allocate(genp(1:num_real_modes))
    safe_allocate(genq(1:num_real_modes))
 
    safe_allocate(pos_displacements(1:(ions%space%dim*ions%natoms)))
    safe_allocate(vel_displacements(1:(ions%space%dim*ions%natoms)))
    safe_allocate(prefactor(1:(ions%space%dim*ions%natoms)))
 
    safe_allocate(l_m(1:num_real_modes))
 
 
    if (t < low_temperature_tolerance) then
      ! TODO: Implement proper low T limit
      sigmap = 1.0_real64/2.0_real64
      sigmaq = 1.0_real64/2.0_real64
    else
      beta_half = m_one / (2 * p_kb * t)
      sigmap = 1.0_real64/sqrt((2.0_real64 * tanh(beta_half * phonons%frequencies)))
      sigmaq = 1.0_real64/sqrt((2.0_real64 * tanh(beta_half * phonons%frequencies)))
    end if
    mup = m_zero
    muq = m_zero
 
    l_m = sqrt(2.0_real64/(unit_amu%factor * phonons%frequencies)) * phonons%alpha
 
    i = 1
    do iatom = 1, ions%natoms
      do idim = 1, ions%space%dim
        prefactor(i) = 1.0_real64/sqrt(ions%mass(iatom)/unit_amu%factor * phonons%num_super)
        i=i+1
      end do
    end do
 
    ! get generalized mode coords from Wigner function
    genq = wigner%get(sigmaq, muq, wigner_q) * l_m
    genp = wigner%get(sigmap, mup, wigner_p) / (l_m * unit_amu%factor) ! TODO: Check!!
 
    ! construct initial displacements and velocities for ions (velocities must be applied later)
 
    ions%pos_displacements = m_zero
    ions%vel_displacements = m_zero
 
    ! implement eq.(9) from Kevins paper for periodic systems
 
    do imode = 1, num_real_modes
 
      pos_displacements = prefactor(:) * phonons%eigenvectors(:, imode) * genq(imode)
      vel_displacements = prefactor(:) * phonons%eigenvectors(:, imode) * genp(imode)
 
      iline = 1
      write(message(iline), '("Displacements for mode ",I4)') imode
      do iatom=1, ions%natoms
        iline = iline+1
        write(message(iline), '(3E15.5)') pos_displacements((iatom-1)*ions%space%dim+1:(iatom-1)*ions%space%dim+ions%space%dim)
      end do
      iline = iline+1
      write(message(iline), '("Velocities for mode ",I4)') imode
      do iatom=1, ions%natoms
        iline = iline+1
        write(message(iline), '(3E15.5)') vel_displacements((iatom-1)*ions%space%dim+1:(iatom-1)*ions%space%dim+ions%space%dim)
      end do
 
      call messages_info(iline, namespace=ions%namespace)
 
      call blas_axpy(phonons%dim, 1.0_real64, pos_displacements(1), 1, ions%pos_displacements(1,1), 1)
      call blas_axpy(phonons%dim, 1.0_real64, vel_displacements(1), 1, ions%vel_displacements(1,1), 1)
 
    end do
 
    safe_deallocate_a(sigmap)
    safe_deallocate_a(sigmaq)
    safe_deallocate_a(mup)
    safe_deallocate_a(muq)
    safe_deallocate_a(genp)
    safe_deallocate_a(genq)
    safe_deallocate_a(l_m)
 
    safe_deallocate_a(prefactor)
    safe_deallocate_a(pos_displacements)
    safe_deallocate_a(vel_displacements)
 
    pop_sub(ions_init_random_displacements)
 
  end subroutine ions_init_random_displacements
 
 
  !
  subroutine ions_single_mode_displacements(ions, mode, amplitude_pos, amplitude_vel)
    class(ions_t), intent(inout) :: ions
    integer,       intent(in)    :: mode
    real(real64), optional, intent(in) :: amplitude_pos
    real(real64), optional, intent(in) :: amplitude_vel
 
    type(phonon_modes_t)  :: phonons
    integer               :: num_real_modes, i, iatom, idim
    real(real64)          :: l_m, genP, genQ
 
    real(real64), allocatable :: pos_displacements(:), vel_displacements(:), prefactor(:)
 
 
    push_sub(ions_single_mode_displacements)
 
    write(message(1), '("Create displacements for mode ", I4)') mode
    call messages_info(1, namespace=ions%namespace)
 
    ! initialize the phonons (read info from file)
    call phonons%init(ions%namespace, ions%space%dim, ions%natoms, ions%space%periodic_dim > 0)
 
    num_real_modes = phonons%num_modes
 
    if (mode > num_real_modes) then
      write(message(1), '("Requested mode ",I0," exceeds number of available modes (",I0,")")') mode, num_real_modes
      call messages_fatal(1, namespace=ions%namespace)
 
      pop_sub(ions_single_mode_displacements)
      return
    end if
 
    ions%pos_displacements = m_zero
    ions%vel_displacements = m_zero
 
    l_m = sqrt(2.0_real64/(unit_amu%factor * phonons%frequencies(mode))) * phonons%alpha(mode)
 
    safe_allocate(pos_displacements(1:(ions%space%dim*ions%natoms)))
    safe_allocate(vel_displacements(1:(ions%space%dim*ions%natoms)))
    safe_allocate(prefactor(1:(ions%space%dim*ions%natoms)))
 
    i = 1
    do iatom = 1, ions%natoms
      do idim=1, ions%space%dim
        prefactor(i) = 1.0_real64/sqrt(ions%mass(iatom)/unit_amu%factor * phonons%num_super)
        i=i+1
      end do
    end do
 
    ! get generalized mode coords from Wigner function
    genq = amplitude_pos * l_m
    genp = amplitude_vel / (l_m * unit_amu%factor)! TODO: Check!!
 
    pos_displacements = prefactor(:) * phonons%eigenvectors(:, mode) * genq
    vel_displacements = prefactor(:) * phonons%eigenvectors(:, mode) * genp
 
    write(message(1), '("Displacements for mode ",I4)') mode
    call messages_info(1, namespace=ions%namespace)
    do iatom=1, ions%natoms
      write(message(1), '(3E15.5)') pos_displacements((iatom-1)*ions%space%dim+1:(iatom-1)*ions%space%dim+ions%space%dim)
      call messages_info(1, namespace=ions%namespace)
    end do
    write(message(1), '("Velocities for mode ",I4)') mode
    call messages_info(1, namespace=ions%namespace)
    do iatom=1, ions%natoms
      write(message(1), '(3E15.5)') vel_displacements((iatom-1)*ions%space%dim+1:(iatom-1)*ions%space%dim+ions%space%dim)
      call messages_info(1, namespace=ions%namespace)
    end do
 
 
    call blas_axpy(phonons%dim, 1.0_real64, pos_displacements(1), 1, ions%pos_displacements(1,1), 1)
    call blas_axpy(phonons%dim, 1.0_real64, vel_displacements(1), 1, ions%vel_displacements(1,1), 1)
 
    safe_deallocate_a(pos_displacements)
    safe_deallocate_a(vel_displacements)
    safe_deallocate_a(prefactor)
 
    pop_sub(ions_single_mode_displacements)
 
  end subroutine ions_single_mode_displacements
 
  ! ---------------------------------------------------------
  subroutine ions_initialize(this)
    class(ions_t), intent(inout) :: this
 
    push_sub(ions_initialize)
 
    ! At this point, we need to apply initial velocities from the random displacements.
    if(allocated(this%vel_displacements)) then
      this%vel = this%vel + this%vel_displacements
    end if
 
    pop_sub(ions_initialize)
  end subroutine ions_initialize
 
  ! ---------------------------------------------------------
  subroutine ions_update_quantity(this, label)
    class(ions_t),    intent(inout) :: this
    character(len=*), intent(in)    :: label
 
    push_sub(ions_update_quantity)
 
    select case (label)
    case default
      ! Other quantities should be handled by the parent class
      call charged_particles_update_quantity(this, label)
    end select
 
    pop_sub(ions_update_quantity)
  end subroutine ions_update_quantity
 
  ! ---------------------------------------------------------
  subroutine ions_init_interaction_as_partner(partner, interaction)
    class(ions_t),                intent(in)    :: partner
    class(interaction_surrogate_t), intent(inout) :: interaction
 
    push_sub(ions_init_interaction_as_partner)
 
    select type (interaction)
    class default
      call charged_particles_init_interaction_as_partner(partner, interaction)
    end select
 
    pop_sub(ions_init_interaction_as_partner)
  end subroutine ions_init_interaction_as_partner
 
  ! ---------------------------------------------------------
  subroutine ions_copy_quantities_to_interaction(partner, interaction)
    class(ions_t),          intent(inout) :: partner
    class(interaction_surrogate_t), intent(inout) :: interaction
 
    push_sub(ions_copy_quantities_to_interaction)
 
    select type (interaction)
    class default
      ! Other interactions should be handled by the parent class
      call charged_particles_copy_quantities_to_interaction(partner, interaction)
    end select
 
    pop_sub(ions_copy_quantities_to_interaction)
  end subroutine ions_copy_quantities_to_interaction
 
  ! ---------------------------------------------------------
  subroutine ions_fold_atoms_into_cell(this)
    class(ions_t),       intent(inout) :: this
 
    integer :: iatom
 
    push_sub(ions_fold_atoms_into_cell)
 
    do iatom = 1, this%natoms
      this%pos(:, iatom) = this%latt%fold_into_cell(this%pos(:, iatom))
    end do
 
    pop_sub(ions_fold_atoms_into_cell)
  end subroutine ions_fold_atoms_into_cell
 
  ! ---------------------------------------------------------
  real(real64) function ions_min_distance(this, real_atoms_only) result(rmin)
    class(ions_t),      intent(in) :: this
    logical,  optional, intent(in) :: real_atoms_only
 
    integer :: iatom, jatom, idir
    real(real64)   :: xx(this%space%dim)
    logical :: real_atoms_only_
    class(species_t), pointer :: species
 
    if (this%natoms == 1 .and. .not. this%space%is_periodic()) then
      rmin = m_huge
      return
    end if
 
    push_sub(ions_min_distance)
 
    real_atoms_only_ = optional_default(real_atoms_only, .false.)
 
    ! Without this line, valgrind complains about a conditional jump on uninitialized variable,
    ! as atom_get_species has an intent(out) that causes a call to the finalizer (with Ifort)
    nullify(species)
 
    rmin = huge(rmin)
    do iatom = 1, this%natoms
      call atom_get_species(this%atom(iatom), species)
      select type(species)
      class is(jellium_t)
        if (real_atoms_only_) cycle
      end select
      do jatom = iatom + 1, this%natoms
        call atom_get_species(this%atom(jatom), species)
        select type(species)
        class is(jellium_t)
          if (real_atoms_only_) cycle
        end select
        xx = abs(this%pos(:, iatom) - this%pos(:, jatom))
        if (this%space%is_periodic()) then
          xx = this%latt%cart_to_red(xx)
          do idir = 1, this%space%periodic_dim
            xx(idir) = xx(idir) - floor(xx(idir) + m_half)
          end do
          xx = this%latt%red_to_cart(xx)
        end if
        rmin = min(norm2(xx), rmin)
      end do
    end do
 
    if (.not. (this%only_user_def .and. real_atoms_only_)) then
      ! what if the nearest neighbors are periodic images?
      do idir = 1, this%space%periodic_dim
        rmin = min(rmin, norm2(this%latt%rlattice(:,idir)))
      end do
    end if
 
    ! To avoid numerical instabilities, we round this to 6 digits only
    if(rmin < huge(rmin)/1e6_real64) then
      rmin = anint(rmin*1e6_real64)*1.0e-6_real64
    end if
 
    pop_sub(ions_min_distance)
  end function ions_min_distance
 
  ! ---------------------------------------------------------
  logical function ions_has_time_dependent_species(this) result(time_dependent)
    class(ions_t),     intent(in) :: this
 
    push_sub(ions_has_time_dependent_species)
 
    time_dependent = this%species_time_dependent
 
    pop_sub(ions_has_time_dependent_species)
  end function ions_has_time_dependent_species
 
  ! ---------------------------------------------------------
  real(real64) function ions_val_charge(this, mask) result(val_charge)
    class(ions_t),              intent(in) :: this
    logical,          optional, intent(in) :: mask(:)
 
    push_sub(ions_val_charge)
 
    if (present(mask)) then
      val_charge = -sum(this%charge, mask=mask)
    else
      val_charge = -sum(this%charge)
    end if
 
    pop_sub(ions_val_charge)
  end function ions_val_charge
 
  ! ---------------------------------------------------------
  function ions_dipole(this, mask) result(dipole)
    class(ions_t),               intent(in) :: this
    logical,           optional, intent(in) :: mask(:)
    real(real64) :: dipole(this%space%dim)
 
    integer :: ia
 
    push_sub(ions_dipole)
 
    dipole = m_zero
    do ia = 1, this%natoms
      if (present(mask)) then
        if (.not. mask(ia)) cycle
      end if
      dipole = dipole + this%charge(ia)*this%pos(:, ia)
    end do
    dipole = p_proton_charge*dipole
 
    pop_sub(ions_dipole)
  end function ions_dipole
 
  ! ---------------------------------------------------------
  subroutine ions_translate(this, xx)
    class(ions_t),     intent(inout) :: this
    real(real64),      intent(in)    :: xx(this%space%dim)
 
    integer  :: iatom
 
    push_sub(ions_translate)
 
    do iatom = 1, this%natoms
      this%pos(:, iatom) = this%pos(:, iatom) - xx
    end do
 
    pop_sub(ions_translate)
  end subroutine ions_translate
 
  ! ---------------------------------------------------------
  subroutine ions_rotate(this, from, from2, to)
    class(ions_t),     intent(inout) :: this
    real(real64),      intent(in)    :: from(this%space%dim)
    real(real64),      intent(in)    :: from2(this%space%dim)
    real(real64),      intent(in)    :: to(this%space%dim)
 
    integer :: iatom
    real(real64) :: m1(3, 3), m2(3, 3)
    real(real64) :: m3(3, 3), f2(3), per(3)
    real(real64) :: alpha, r
 
    push_sub(ions_rotate)
 
    if (this%space%dim /= 3) then
      call messages_not_implemented("ions_rotate in other than 3 dimensions", namespace=this%namespace)
    end if
 
    ! initialize matrices
    m1 = diagonal_matrix(3, m_one)
 
    ! rotate the to-axis to the z-axis
    if (abs(to(2)) >  1d-150) then
      alpha = atan2(to(2), to(1))
      call rotate(m1, alpha, 3)
    end if
    alpha = atan2(norm2(to(1:2)), to(3))
    call rotate(m1, -alpha, 2)
 
    ! get perpendicular to z and from
    f2 = matmul(m1, from)
    per(1) = -f2(2)
    per(2) =  f2(1)
    per(3) = m_zero
    r = norm2(per)
    if (r > m_zero) then
      per = per/r
    else
      per(2) = m_one
    end if
 
    ! rotate perpendicular axis to the y-axis
    m2 = diagonal_matrix(3, m_one)
    alpha = atan2(per(1), per(2))
    call rotate(m2, -alpha, 3)
 
    ! rotate from => to (around the y-axis)
    m3 = diagonal_matrix(3, m_one)
    alpha = acos(min(sum(from*to), m_one))
    call rotate(m3, -alpha, 2)
 
    ! join matrices
    m2 = matmul(transpose(m2), matmul(m3, m2))
 
    ! rotate around the z-axis to get the second axis
    per = matmul(m2, matmul(m1, from2))
    alpha = atan2(per(1), per(2))
    call rotate(m2, -alpha, 3) ! second axis is now y
 
    ! get combined transformation
    m1 = matmul(transpose(m1), matmul(m2, m1))
 
    ! now transform the coordinates
    ! it is written in this way to avoid what I consider a bug in the Intel compiler
    do iatom = 1, this%natoms
      f2 = this%pos(:, iatom)
      this%pos(:, iatom) = matmul(m1, f2)
    end do
 
    pop_sub(ions_rotate)
  contains
 
    ! ---------------------------------------------------------
    subroutine rotate(m, angle, dir)
      real(real64),   intent(inout) :: m(3, 3)
      real(real64),   intent(in)    :: angle
      integer, intent(in)    :: dir
 
      real(real64) :: aux(3, 3), ca, sa
 
      push_sub(ions_rotate.rotate)
 
      ca = cos(angle)
      sa = sin(angle)
 
      aux = m_zero
      select case (dir)
      case (1)
        aux(1, 1) = m_one
        aux(2, 2) = ca
        aux(3, 3) = ca
        aux(2, 3) = sa
        aux(3, 2) = -sa
      case (2)
        aux(2, 2) = m_one
        aux(1, 1) = ca
        aux(3, 3) = ca
        aux(1, 3) = sa
        aux(3, 1) = -sa
      case (3)
        aux(3, 3) = m_one
        aux(1, 1) = ca
        aux(2, 2) = ca
        aux(1, 2) = sa
        aux(2, 1) = -sa
      end select
 
      m = matmul(aux, m)
 
      pop_sub(ions_rotate.rotate)
    end subroutine rotate
 
  end subroutine ions_rotate
 
  ! ---------------------------------------------------------
  function ions_global_force(this, time) result(force)
    class(ions_t),        intent(in)    :: this
    real(real64),         intent(in)    :: time
    real(real64) :: force(this%space%dim)
 
    push_sub(ions_global_force)
 
    force = m_zero
 
    if (this%apply_global_force) then
      force(1) = tdf(this%global_force_function, time)
    end if
 
    pop_sub(ions_global_force)
  end function ions_global_force
 
  ! ---------------------------------------------------------
  subroutine ions_write_xyz(this, fname, append, comment, reduce_coordinates)
    class(ions_t),              intent(in) :: this
    character(len=*),           intent(in) :: fname
    logical,          optional, intent(in) :: append
    character(len=*), optional, intent(in) :: comment
    logical,          optional, intent(in) :: reduce_coordinates
 
    integer :: iatom, idim, iunit
    character(len=6) position
    character(len=19) :: frmt
    real(real64) :: red(this%space%dim)
 
    if (.not. this%grp%is_root()) return
 
    push_sub(ions_write_xyz)
 
    position = 'asis'
    if (present(append)) then
      if (append) position = 'append'
    end if
    if(.not.optional_default(reduce_coordinates, .false.)) then
      iunit = io_open(trim(fname)//'.xyz', this%namespace, action='write', position=position)
    else
      iunit = io_open(trim(fname)//'.xyz_red', this%namespace, action='write', position=position)
    end if
 
    write(iunit, '(i4)') this%natoms
    if (present(comment)) then
      write(iunit, '(1x,a)') comment
    else
      write(iunit, '(1x,a,a)') 'units: ', trim(units_abbrev(units_out%length_xyz_file))
    end if
 
    write(unit=frmt, fmt="(a5,i2.2,a4,i2.2,a6)") "(6x,a", label_len, ",2x,", this%space%dim,"f15.9)"
    do iatom = 1, this%natoms
      if(.not.optional_default(reduce_coordinates, .false.)) then
        write(unit=iunit, fmt=frmt) this%atom(iatom)%label, &
          (units_from_atomic(units_out%length_xyz_file, this%pos(idim, iatom)), idim=1, this%space%dim)
      else
        red = this%latt%cart_to_red(this%pos(:, iatom))
        write(unit=iunit, fmt=frmt) this%atom(iatom)%label, (red(idim), idim=1, this%space%dim)
      end if
    end do
    call io_close(iunit)
 
    pop_sub(ions_write_xyz)
  end subroutine ions_write_xyz
 
  !
  subroutine ions_write_poscar(this, fname, comment)
    class(ions_t),              intent(in) :: this
    character(len=*), optional, intent(in) :: fname
    character(len=*), optional, intent(in) :: comment
 
    integer                        :: iatom, idim, iunit
    character(len=:), allocatable  :: fname_, comment_
    character(len=10) :: format_string
 
    push_sub(ions_write_poscar)
 
    if (.not. this%grp%is_root()) then
      pop_sub(ions_write_poscar)
      return
    end if
 
    comment_ = optional_default(comment, "")
    fname_ = optional_default(fname, "POSCAR")
 
    iunit = io_open(trim(fname_), this%namespace, action='write')
 
    write(iunit, '(A)') comment_ ! mandatory comment
    write(iunit, '("1.0")')      ! scaling factor: we use 1 as lattice vectors are absolute
 
    write(format_string, '("(",I1,A,")")') this%space%dim, 'F15.10'
    do idim=1, this%space%dim
      write(iunit, format_string) units_from_atomic(unit_angstrom, this%latt%rlattice(:,idim))
    end do
 
    do iatom = 1, this%natoms
      write(iunit, '(A)', advance='NO') trim(this%atom(iatom)%label)//" "
    end do
    write(iunit, '("")')
    write(iunit, '(A)') repeat("1 ", this%natoms)
    write(iunit, '("direct")')
    write(format_string, '("(",I1,A,")")') this%space%dim, 'F15.10'
    do iatom=1, this%natoms
      write(iunit, format_string) this%latt%cart_to_red(this%pos(:, iatom))
    end do
 
    call io_close(iunit)
 
    pop_sub(ions_write_poscar)
  end subroutine ions_write_poscar
 
  ! ---------------------------------------------------------
  subroutine ions_read_xyz(this, fname, comment)
    use iso_fortran_env, only: real64
    class(ions_t),              intent(inout) :: this
    character(len=*),           intent(in)    :: fname
    character(len=*), optional, intent(in)    :: comment
 
    integer :: iatom, idir, iunit, ios
    character(len=256) :: line
    character(len=LABEL_LEN) :: label
    character(len=:), allocatable :: coordstr
    real(real64) :: tmp(this%space%dim)
 
    push_sub(ions_read_xyz)
 
    iunit = io_open(trim(fname)//'.xyz', this%namespace, action='read', position='rewind')
 
    ! --- First two lines: natoms + optional comment ------------------------
    read(iunit, '(i4)') this%natoms
    if (present(comment)) then
      read(iunit, *)
    else
      read(iunit, *)
    end if
 
    ! --- Read atom lines ---------------------------------------------------
    do iatom = 1, this%natoms
 
      ! Read entire line (safe even if long)
      read(iunit,'(A)', iostat=ios) line
      if (ios /= 0) then
        call io_close(iunit)
        message(1) = "Error reading XYZ atom line"
        call messages_fatal(1, namespace=this%namespace)
      end if
 
      ! Extract fixed-length label (preserves spaces inside)
      if (len_trim(line) < label_len) then
        call io_close(iunit)
        message(1) = "XYZ file: atom line too short for label"
        call messages_fatal(1, namespace=this%namespace)
      end if
      label = line(1:label_len)
 
      ! Extract remaining part for coordinates
      if (len(line) > label_len) then
        coordstr = adjustl(line(label_len+1:))
      else
        call io_close(iunit)
        message(1) = "XYZ file: no coordinate data after label"
        call messages_fatal(1, namespace=this%namespace)
      end if
 
      ! Parse coordinates in a flexible (list-directed) way
      read(coordstr, *, iostat=ios) (tmp(idir), idir=1, this%space%dim)
      if (ios /= 0) then
        call io_close(iunit)
        message(1) = "XYZ file: malformed coordinate fields"
        call messages_fatal(1, namespace=this%namespace)
      end if
 
      ! Convert and store
      this%pos(:, iatom) = units_to_atomic(units_out%length_xyz_file, tmp)
 
    end do
 
    call io_close(iunit)
 
    pop_sub(ions_read_xyz)
  end subroutine ions_read_xyz
 
  ! ----------------------------------------------------------------
  subroutine ions_write_crystal(this, dir)
    class(ions_t),           intent(in) :: this
    character(len=*),        intent(in) :: dir
 
    type(lattice_iterator_t) :: latt_iter
    real(real64) :: radius, pos(this%space%dim)
    integer :: iatom, icopy, iunit
 
    push_sub(ions_write_crystal)
 
    radius = maxval(m_half*norm2(this%latt%rlattice(:,1:this%space%periodic_dim), dim=1))*(m_one + m_epsilon)
    latt_iter = lattice_iterator_t(this%latt, radius)
 
    if (this%grp%is_root()) then
 
      iunit = io_open(trim(dir)//'/crystal.xyz', this%namespace, action='write')
 
      write(iunit, '(i9)') this%natoms*latt_iter%n_cells
      write(iunit, '(a)') '#generated by Octopus'
 
      do iatom = 1, this%natoms
        do icopy = 1, latt_iter%n_cells
          pos = units_from_atomic(units_out%length, this%pos(:, iatom) + latt_iter%get(icopy))
          write(iunit, '(a, 99f12.6)') this%atom(iatom)%label, pos
        end do
      end do
 
      call io_close(iunit)
    end if
 
    pop_sub(ions_write_crystal)
  end subroutine ions_write_crystal
 
  ! ---------------------------------------------------------
  subroutine ions_write_bild_forces_file(this, dir, fname)
    class(ions_t),      intent(in) :: this
    character(len=*),   intent(in) :: dir, fname
 
    integer :: iunit, iatom, idir
    real(real64) :: force(this%space%dim), center(this%space%dim)
    character(len=20) frmt
 
    if (.not. this%grp%is_root()) return
 
    push_sub(ions_write_bild_forces_file)
 
    call io_mkdir(dir, this%namespace)
    iunit = io_open(trim(dir)//'/'//trim(fname)//'.bild', this%namespace, action='write', &
      position='asis')
 
    write(frmt,'(a,i0,a)')'(a,2(', this%space%dim,'f16.6,1x))'
 
    write(iunit, '(a)')'.comment : force vectors in ['//trim(units_abbrev(units_out%force))//']'
    write(iunit, *)
    write(iunit, '(a)')'.color red'
    write(iunit, *)
    do iatom = 1, this%natoms
      center = units_from_atomic(units_out%length, this%pos(:, iatom))
      force = units_from_atomic(units_out%force, this%tot_force(:, iatom))
      write(iunit, '(a,1x,i4,1x,a2,1x,a6,1x,f10.6,a)') '.comment :', iatom, trim(this%atom(iatom)%label), &
        'force:', norm2(force), '['//trim(units_abbrev(units_out%force))//']'
      write(iunit,fmt=trim(frmt)) '.arrow', (center(idir), idir = 1, this%space%dim), &
        (center(idir) + force(idir), idir = 1, this%space%dim)
      write(iunit,*)
    end do
 
    call io_close(iunit)
 
    pop_sub(ions_write_bild_forces_file)
  end subroutine ions_write_bild_forces_file
 
  ! -----------------------------------------------------
  subroutine ions_write_vtk_geometry(this, filename, ascii)
    class(ions_t),              intent(in) :: this
    character(len=*),           intent(in) :: filename
    logical,          optional, intent(in) :: ascii
 
    integer :: iunit, iatom, ierr
    logical :: ascii_
    real(real64), allocatable :: data(:, :)
    integer, allocatable :: idata(:, :)
    character(len=MAX_PATH_LEN) :: fullname
 
    push_sub(ions_write_vtk_geometry)
 
    assert(this%space%dim == 3)
 
    ascii_ = optional_default(ascii, .true.)
 
    fullname = trim(filename)//".vtk"
 
    iunit = io_open(trim(fullname), this%namespace, action='write')
 
    write(iunit, '(1a)') '# vtk DataFile Version 3.0 '
    write(iunit, '(6a)') 'Generated by octopus ', trim(conf%version), ' -  git: ', &
      trim(conf%git_commit), " configuration: ",  trim(conf%config_time)
 
    if (ascii_) then
      write(iunit, '(1a)') 'ASCII'
    else
      write(iunit, '(1a)') 'BINARY'
    end if
 
    write(iunit, '(1a)') 'DATASET POLYDATA'
 
    write(iunit, '(a,i9,a)') 'POINTS ', this%natoms, ' double'
 
    if (ascii_) then
      do iatom = 1, this%natoms
        write(iunit, '(3f12.6)') this%pos(1:3, iatom)
      end do
    else
      call io_close(iunit)
      safe_allocate(data(1:3, 1:this%natoms))
      do iatom = 1, this%natoms
        data(1:3, iatom) = this%pos(1:3, iatom)
      end do
      call io_binary_write(io_workpath(fullname, this%namespace), i4_to_i8(3*this%natoms), data, &
        ierr, nohead = .true., fendian = io_binary_is_little_endian())
      safe_deallocate_a(data)
      iunit = io_open(trim(fullname), this%namespace, action='write', position = 'append')
      write(iunit, '(1a)') ''
    end if
 
    write(iunit, '(a,2i9)') 'VERTICES ', this%natoms, 2*this%natoms
 
    if (ascii_) then
      do iatom = 1, this%natoms
        write(iunit, '(2i9)') 1, iatom - 1
      end do
    else
      call io_close(iunit)
      safe_allocate(idata(1:2, 1:this%natoms))
      do iatom = 1, this%natoms
        idata(1, iatom) = 1
        idata(2, iatom) = iatom - 1
      end do
      call io_binary_write(io_workpath(fullname, this%namespace), i4_to_i8(2*this%natoms), idata, &
        ierr, nohead = .true., fendian = io_binary_is_little_endian())
      safe_deallocate_a(idata)
      iunit = io_open(trim(fullname), this%namespace, action='write', position = 'append')
      write(iunit, '(1a)') ''
    end if
 
    write(iunit, '(a,i9)') 'POINT_DATA', this%natoms
    write(iunit, '(a)') 'SCALARS element integer'
    write(iunit, '(a)') 'LOOKUP_TABLE default'
 
    if (ascii_) then
      do iatom = 1, this%natoms
        write(iunit, '(i9)') nint(this%atom(iatom)%species%get_z())
      end do
    else
      call io_close(iunit)
 
      safe_allocate(idata(1:this%natoms, 1))
 
      do iatom = 1, this%natoms
        idata(iatom, 1) = nint(this%atom(iatom)%species%get_z())
      end do
 
      call io_binary_write(io_workpath(fullname, this%namespace), i4_to_i8(this%natoms), idata, &
        ierr, nohead = .true., fendian = io_binary_is_little_endian())
 
      safe_deallocate_a(idata)
 
      iunit = io_open(trim(fullname), this%namespace, action='write', position = 'append')
      write(iunit, '(1a)') ''
    end if
 
    call io_close(iunit)
 
    pop_sub(ions_write_vtk_geometry)
  end subroutine ions_write_vtk_geometry
 
  ! ---------------------------------------------------------
  function ions_current(this) result(current)
    class(ions_t),   intent(in) :: this
    real(real64) :: current(this%space%dim)
 
    integer :: iatom
 
    push_sub(ions_current)
 
    current = m_zero
    do iatom = this%atoms_dist%start, this%atoms_dist%end
      current = current + p_proton_charge*this%atom(iatom)%species%get_zval()*this%vel(:,iatom)
    end do
 
    if (this%atoms_dist%parallel) then
      call comm_allreduce(this%atoms_dist%mpi_grp, current, dim=this%space%dim)
    end if
 
    pop_sub(ions_current)
  end function ions_current
 
  ! ---------------------------------------------------------
  function ions_abs_current(this) result(abs_current)
    class(ions_t),   intent(in) :: this
    real(real64) :: abs_current(this%space%dim)
 
    integer :: iatom
 
    push_sub(ions_abs_current)
 
    abs_current = m_zero
    do iatom = this%atoms_dist%start, this%atoms_dist%end
      abs_current = abs_current + abs(p_proton_charge*this%atom(iatom)%species%get_zval()*this%vel(:,iatom))
    end do
 
    if (this%atoms_dist%parallel) then
      call comm_allreduce(this%atoms_dist%mpi_grp, abs_current, dim=this%space%dim)
    end if
 
    pop_sub(ions_abs_current)
  end function ions_abs_current
 
  ! ---------------------------------------------------------
  subroutine ions_finalize(ions)
    type(ions_t), intent(inout) :: ions
 
    integer :: i
 
    push_sub(ions_finalize)
 
    call distributed_end(ions%atoms_dist)
 
    call ion_interaction_end(ions%ion_interaction)
 
    safe_deallocate_a(ions%atom)
    ions%natoms=0
 
    if(allocated(ions%species)) then
      do i = 1, ions%nspecies
        !SAFE_ DEALLOCATE_P(ions%species(i)%s)
        if(associated(ions%species(i)%s)) deallocate(ions%species(i)%s)
      end do
      deallocate(ions%species)
    end if
    ions%nspecies=0
 
    safe_deallocate_a(ions%pos_displacements)
    safe_deallocate_a(ions%vel_displacements)
 
    call charged_particles_end(ions)
 
    safe_deallocate_a(ions%map_symm_atoms)
    safe_deallocate_a(ions%inv_map_symm_atoms)
 
 
    pop_sub(ions_finalize)
  end subroutine ions_finalize
 
  !-------------------------------------------------------------------
  subroutine ions_update_lattice_vectors(ions, latt, symmetrize)
    class(ions_t),           intent(inout) :: ions
    type(lattice_vectors_t), intent(in)    :: latt
    logical,                 intent(in)    :: symmetrize
 
    push_sub(ions_update_lattice_vectors)
 
    ! Update the lattice vectors
    call ions%latt%update(latt%rlattice)
 
    ! Regenerate symmetries in Cartesian space
    if (symmetrize) then
      call symmetries_update_lattice_vectors(ions%symm, latt, ions%space%dim)
    end if
 
    pop_sub(ions_update_lattice_vectors)
  end subroutine ions_update_lattice_vectors
 
  !-------------------------------------------------------------------
  subroutine ions_generate_mapping_symmetric_atoms(ions)
    class(ions_t),    intent(inout) :: ions
 
    integer :: iatom, iop, iatom_symm, dim4symms
    real(real64) :: ratom(ions%space%dim)
 
    push_sub(ions_generate_mapping_symmetric_atoms)
 
    safe_allocate(ions%map_symm_atoms(1:ions%natoms, 1:ions%symm%nops + ions%symm%nops_nonsymmorphic))
    safe_allocate(ions%inv_map_symm_atoms(1:ions%natoms, 1:ions%symm%nops + ions%symm%nops_nonsymmorphic))
 
    ! In the 4D case, symmetries are only defined in 3D, see symmetries.F90
    dim4symms = min(3, ions%space%dim)
    ratom = m_zero
 
    do iop = 1, ions%symm%nops
      do iatom = 1, ions%natoms
        !We find the atom that correspond to this one, once symmetry is applied
        ratom(1:dim4symms) = symm_op_apply_inv_cart(ions%symm%ops(iop), ions%pos(:, iatom))
 
        ratom(:) = ions%latt%fold_into_cell(ratom(:))
 
        ! find iatom_symm
        do iatom_symm = 1, ions%natoms
          if (all(abs(ratom(:) - ions%pos(:, iatom_symm)) < symprec)) exit
        end do
 
        if (iatom_symm > ions%natoms) then
          write(message(1),'(a,i6)') 'Internal error: could not find symmetric partner for atom number', iatom
          write(message(2),'(a,i3,a)') 'with symmetry operation number ', iop, '.'
          call messages_fatal(2, namespace=ions%namespace)
        end if
 
        ions%map_symm_atoms(iatom, iop) = iatom_symm
        ions%inv_map_symm_atoms(iatom_symm, iop) = iatom
      end do
    end do
 
    ! Also build a map for non-symmorphic operations
    do iop = 1, ions%symm%nops_nonsymmorphic
      do iatom = 1, ions%natoms
        !We find the atom that correspond to this one, once symmetry is applied
        ratom(1:dim4symms) = symm_op_apply_inv_cart(ions%symm%non_symmorphic_ops(iop), ions%pos(:, iatom))
 
        ratom(:) = ions%latt%fold_into_cell(ratom(:))
 
        ! find iatom_symm
        do iatom_symm = 1, ions%natoms
          if (all(abs(ratom(:) - ions%pos(:, iatom_symm)) < symprec)) exit
        end do
 
        if (iatom_symm > ions%natoms) then
          write(message(1),'(a,i6)') 'Internal error: could not find symmetric partner for atom number', iatom
          write(message(2),'(a,i3,a)') 'with symmetry operation number ', iop, '.'
          call messages_fatal(2, namespace=ions%namespace)
        end if
 
        ions%map_symm_atoms(iatom, ions%symm%nops+iop) = iatom_symm
        ions%inv_map_symm_atoms(iatom_symm, ions%symm%nops+iop) = iatom
      end do
    end do
 
 
    pop_sub(ions_generate_mapping_symmetric_atoms)
  end subroutine ions_generate_mapping_symmetric_atoms
 
  !-------------------------------------------------------------------
  subroutine ions_symmetrize_atomic_coord(ions)
    class(ions_t),    intent(inout) :: ions
 
    integer :: iatom, iop, iatom_sym
    real(real64) :: ratom(ions%space%dim)
    real(real64), allocatable :: new_pos(:,:)
 
    push_sub(ions_symmetrize_atomic_coord)
 
    safe_allocate(new_pos(1:ions%space%dim, 1:ions%natoms))
 
    do iatom = 1, ions%natoms
      new_pos(:, iatom) = m_zero
 
      ! Symmorphic operations
      do iop = 1, ions%symm%nops
        iatom_sym = ions%inv_map_symm_atoms(iatom, iop)
        ratom = symm_op_apply_inv_cart(ions%symm%ops(iop), ions%pos(:, iatom_sym))
        ratom = ions%latt%fold_into_cell(ratom)
        new_pos(:, iatom) = new_pos(:, iatom) + ratom
      end do
 
      ! Non-symmorphic operations
      do iop = 1, ions%symm%nops_nonsymmorphic
        iatom_sym = ions%inv_map_symm_atoms(iatom, iop + ions%symm%nops)
        ratom = symm_op_apply_inv_cart(ions%symm%non_symmorphic_ops(iop), ions%pos(:, iatom_sym))
        ratom = ions%latt%fold_into_cell(ratom)
        new_pos(:, iatom) = new_pos(:, iatom) + ratom
      end do
 
      new_pos(:, iatom) = new_pos(:, iatom) / (ions%symm%nops + ions%symm%nops_nonsymmorphic)
    end do
 
    ions%pos = new_pos
 
    pop_sub(ions_symmetrize_atomic_coord)
  end subroutine ions_symmetrize_atomic_coord
 
  subroutine ions_print_spacegroup(ions)
    class(ions_t),    intent(inout) :: ions
 
    type(spglibdataset) :: spg_dataset
    character(len=11) :: symbol
    integer, allocatable :: site_type(:)
    integer :: space_group, ia
 
    if(.not. ions%space%is_periodic()) return
 
    push_sub(ions_print_spacegroup)
 
    safe_allocate(site_type(1:ions%natoms))
    do ia = 1, ions%natoms
      site_type(ia) = ions%atom(ia)%species%get_index()
    end do
 
    spg_dataset = symmetries_get_spg_dataset(ions%namespace, ions%space, ions%latt, ions%natoms, ions%pos, site_type)
 
    safe_deallocate_a(site_type)
 
    if (spg_dataset%spglib_error /= 0) then
      pop_sub(ions_print_spacegroup)
      return
    end if
 
    space_group = spg_dataset%spacegroup_number
    symbol = spg_dataset%international_symbol
 
    write(message(1),'(a, i4)') 'Info: Space group No. ', space_group
    write(message(2),'(2a)') 'Info: International: ', trim(symbol)
    call messages_info(2, namespace=ions%namespace)
 
    pop_sub(ions_print_spacegroup)
  end subroutine ions_print_spacegroup
 
end module ions_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: