ion_dynamics.F90

Go to the documentation of this file.

!! Copyright (C) 2008 X. Andrade
!!
!! 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 ion_dynamics_oct_m
  use accel_oct_m
  use box_parallelepiped_oct_m
  use iso_c_binding
  use debug_oct_m
  use global_oct_m
  use grid_oct_m
  use ions_oct_m
  use, intrinsic :: iso_fortran_env
  use kpoints_oct_m
  use lalg_adv_oct_m
  use lattice_vectors_oct_m
  use loct_math_oct_m
  use math_oct_m
  use messages_oct_m
  use mpi_oct_m
  use multicomm_oct_m
  use namespace_oct_m
  use parser_oct_m
  use poisson_oct_m
  use read_coords_oct_m
  use restart_oct_m
  use space_oct_m
  use species_oct_m
  use symmetrizer_oct_m
  use tdfunction_oct_m
  use profiling_oct_m
  use unit_oct_m
  use unit_system_oct_m
  use varinfo_oct_m
 
  implicit none
 
  private
 
  public ::                                &
    ion_dynamics_t,                        &
    ion_state_t,                           &
    ion_dynamics_init,                     &
    ion_dynamics_end,                      &
    ion_dynamics_propagate,                &
    ion_dynamics_propagate_vel,            &
    ion_dynamics_save_state,               &
    ion_dynamics_restore_state,            &
    ion_dynamics_drive_ions,               &
    ion_dynamics_temperature,              &
    ion_dynamics_freeze,                   &
    ion_dynamics_unfreeze,                 &
    ion_dynamics_verlet_step1,             &
    ion_dynamics_verlet_step2,             &
    ion_dynamics_dump,                     &
    ion_dynamics_load,                     &
    ion_dynamics_box_update,               &
    electrons_lattice_vectors_update
 
 
  integer, parameter ::   &
    THERMO_NONE     = 0,  &
    thermo_scal     = 1,  &
    thermo_nh       = 2
 
  type nose_hoover_t
    private
    real(real64) :: mass
    real(real64) :: pos
    real(real64) :: vel
  end type nose_hoover_t
 
  type ion_td_displacement_t
    private
    logical     :: move
    type(tdf_t) :: fx
    type(tdf_t) :: fy
    type(tdf_t) :: fz
  end type ion_td_displacement_t
 
  type ion_dynamics_t
    private
    logical          :: move_ions
    logical          :: relax_cell
    logical          :: constant_velocity
    integer          :: thermostat
    real(real64)     :: dt
    real(real64)     :: current_temperature
 
    real(real64), allocatable :: oldforce(:, :)
 
    real(real64), allocatable :: old_pos(:, :)
 
    real(real64), allocatable :: cell_force(:)
    real(real64), allocatable :: old_cell_force(:)
    real(real64), allocatable :: cell_vel(:)
    real(real64), allocatable :: initial_rlattice(:,:)
    real(real64), allocatable :: strain(:)
 
    real(real64) :: pressure
 
    logical :: symmetrize = .false.
    type(symmetrizer_t), pointer :: symm
 
    type(nose_hoover_t) :: nh(1:2)
    type(tdf_t) :: temperature_function
 
    logical :: drive_ions
    type(ion_td_displacement_t), allocatable ::  td_displacements(:)
    type(ions_t),     pointer :: ions_t0
 
    real(real64), public :: ionic_scale
  contains
    procedure :: ions_move => ion_dynamics_ions_move 
    procedure :: cell_relax => ion_dynamics_cell_relax 
    procedure :: update_stress => ion_dynamics_update_stress 
    procedure :: is_active => ion_dynamics_is_active 
  end type ion_dynamics_t
 
  type ion_state_t
    private
    real(real64), allocatable :: pos(:, :)
    real(real64), allocatable :: vel(:, :)
    real(real64), allocatable :: old_pos(:, :)
    type(nose_hoover_t) :: nh(1:2)
  end type ion_state_t
 
  type cell_state_t
    private
    real(real64), allocatable :: pos(:, :)
    real(real64), allocatable :: vel(:, :)
    real(real64), allocatable :: old_pos(:, :)
  end type cell_state_t
 
contains
 
  ! ---------------------------------------------------------
  subroutine ion_dynamics_init(this, namespace, ions, symmetrize, symm)
    use iso_c_binding, only: c_ptr
    type(ion_dynamics_t),        intent(out)   :: this
    type(namespace_t),           intent(in)    :: namespace
    type(ions_t),                intent(inout) :: ions
    logical,                     intent(in)    :: symmetrize
    type(symmetrizer_t), optional, target, intent(in)    :: symm
 
    integer :: i, j, iatom, ierr, periodic_dim, ncomp
    real(real64)   :: xx(ions%space%dim), temperature, sigma, kin1, kin2
    type(c_ptr) :: random_gen_pointer
    type(read_coords_info) :: xyz
    character(len=100)  :: temp_function_name
    logical :: have_velocities
 
    type(block_t)      :: blk
    integer            :: ndisp
    character(len=200) :: expression
 
    push_sub(ion_dynamics_init)
 
    have_velocities = .false.
    this%drive_ions = .false.
 
    this%symmetrize = symmetrize
    if (this%symmetrize) then
      assert(present(symm))
      this%symm => symm
    else
      nullify(this%symm)
    end if
 
    !%Variable TDIonicTimeScale
    !%Type float
    !%Default 1.0
    !%Section Time-Dependent::Propagation
    !%Description
    !% This variable defines the factor between the timescale of ionic
    !% and electronic movement. It allows reasonably fast
    !% Born-Oppenheimer molecular-dynamics simulations based on
    !% Ehrenfest dynamics. The value of this variable is equivalent to
    !% the role of <math>\mu</math> in Car-Parrinello. Increasing it
    !% linearly accelerates the time step of the ion
    !% dynamics, but also increases the deviation of the system from the
    !% Born-Oppenheimer surface. The default is 1, which means that both
    !% timescales are the same. Note that a value different than 1
    !% implies that the electrons will not follow physical behaviour.
    !%
    !% According to our tests, values around 10 are reasonable, but it
    !% will depend on your system, mainly on the width of the gap.
    !%
    !% Important: The electronic time step will be the value of
    !% <tt>TDTimeStep</tt> divided by this variable, so if you have determined an
    !% optimal electronic time step (that we can call <i>dte</i>), it is
    !% recommended that you define your time step as:
    !%
    !% <tt>TDTimeStep</tt> = <i>dte</i> * <tt>TDIonicTimeScale</tt>
    !%
    !% so you will always use the optimal electronic time step
    !% (<a href=http://arxiv.org/abs/0710.3321>more details</a>).
    !%End
    call parse_variable(namespace, 'TDIonicTimeScale', m_one, this%ionic_scale)
 
    if (this%ionic_scale <= m_zero) then
      write(message(1),'(a)') 'Input: TDIonicTimeScale must be positive.'
      call messages_fatal(1, namespace=namespace)
    end if
 
    call messages_print_var_value('TDIonicTimeScale', this%ionic_scale, namespace=namespace)
 
 
 
 
    !%Variable IonsConstantVelocity
    !%Type logical
    !%Default no
    !%Section Time-Dependent::Propagation
    !%Description
    !% (Experimental) If this variable is set to yes, the ions will
    !% move with a constant velocity given by the initial
    !% conditions. They will not be affected by any forces.
    !%End
    call parse_variable(namespace, 'IonsConstantVelocity', .false., this%constant_velocity)
    call messages_print_var_value('IonsConstantVelocity', this%constant_velocity, namespace=namespace)
 
    if (this%constant_velocity) then
      call messages_experimental('IonsConstantVelocity', namespace=namespace)
      have_velocities = .true.
      this%drive_ions = .true.
    end if
 
    !%Variable IonsTimeDependentDisplacements
    !%Type block
    !%Section Time-Dependent::Propagation
    !%Description
    !% (Experimental) This variable allows you to specify a
    !% time-dependent function describing the displacement of the ions
    !% from their equilibrium position: <math>r(t) = r_0 + \Delta
    !% r(t)</math>.  Specify the displacements dx(t), dy(t), dz(t) as
    !% follows, for some or all of the atoms:
    !%
    !% <tt>%IonsTimeDependentDisplacements
    !% <br>&nbsp;&nbsp; atom_index | "dx(t)" | "dy(t)" | "dz(t)"
    !% <br>%</tt>
    !%
    !% The displacement functions are time-dependent functions and should match one
    !% of the function names given in the first column of the <tt>TDFunctions</tt> block.
    !% If this block is set, the ions will not be affected by any forces.
    !%End
 
 
    ndisp = 0
    if (parse_block(namespace, 'IonsTimeDependentDisplacements', blk) == 0) then
      call messages_experimental("IonsTimeDependentDisplacements", namespace=namespace)
      ndisp= parse_block_n(blk)
      safe_allocate(this%td_displacements(1:ions%natoms))
      this%td_displacements(1:ions%natoms)%move = .false.
      if (ndisp > 0) this%drive_ions =.true.
 
      do i = 1, ndisp
        call parse_block_integer(blk, i-1, 0, iatom)
        this%td_displacements(iatom)%move = .true.
 
        call parse_block_string(blk, i-1, 1, expression)
        call tdf_read(this%td_displacements(iatom)%fx, namespace, trim(expression), ierr)
        if (ierr /= 0) then
          write(message(1),'(3A)') 'Could not find "', trim(expression), '" in the TDFunctions block:'
          call messages_warning(1, namespace=namespace)
        end if
 
 
        call parse_block_string(blk, i-1, 2, expression)
        call tdf_read(this%td_displacements(iatom)%fy, namespace, trim(expression), ierr)
        if (ierr /= 0) then
          write(message(1),'(3A)') 'Could not find "', trim(expression), '" in the TDFunctions block:'
          call messages_warning(1, namespace=namespace)
        end if
 
        call parse_block_string(blk, i-1, 3, expression)
        call tdf_read(this%td_displacements(iatom)%fz, namespace, trim(expression), ierr)
        if (ierr /= 0) then
          write(message(1),'(3A)') 'Could not find "', trim(expression), '" in the TDFunctions block:'
          call messages_warning(1, namespace=namespace)
        end if
 
      end do
 
      safe_allocate(this%ions_t0)
      this%ions_t0 = ions
 
    end if
 
 
 
 
    !%Variable Thermostat
    !%Type integer
    !%Default none
    !%Section Time-Dependent::Propagation
    !%Description
    !% This variable selects the type of thermostat applied to
    !% control the ionic temperature.
    !%Option none 0
    !% No thermostat is applied. This is the default.
    !%Option velocity_scaling 1
    !% Velocities are scaled to control the temperature.
    !%Option nose_hoover 2
    !% Nose-Hoover thermostat.
    !%End
 
    call parse_variable(namespace, 'Thermostat', thermo_none, this%thermostat)
    if (.not. varinfo_valid_option('Thermostat', this%thermostat)) call messages_input_error(namespace, 'Thermostat')
    call messages_print_var_option('Thermostat', this%thermostat, namespace=namespace)
 
    if (this%thermostat /= thermo_none) then
 
      have_velocities = .true.
 
      if (this%drive_ions) then
        call messages_write('You cannot use a Thermostat and IonsConstantVelocity or IonsTimeDependentDisplacements')
        call messages_write('at the same time.')
        call messages_fatal(namespace=namespace)
      end if
 
      call messages_experimental('Thermostat', namespace=namespace)
 
      !%Variable TemperatureFunction
      !%Type integer
      !%Default "temperature"
      !%Section Time-Dependent::Propagation
      !%Description
      !% If a thermostat is used, this variable indicates the name of the
      !% function in the <tt>TDFunctions</tt> block that will be used to control the
      !% temperature. The values of the temperature are given in
      !% degrees Kelvin.
      !%End
      call parse_variable(namespace, 'TemperatureFunction', 'temperature', temp_function_name)
 
      call tdf_read(this%temperature_function, namespace, temp_function_name, ierr)
 
      if (ierr /= 0) then
        message(1) = "You have enabled a thermostat but Octopus could not find"
        message(2) = "the '"//trim(temp_function_name)//"' function in the TDFunctions block."
        call messages_fatal(2, namespace=namespace)
      end if
 
      if (this%thermostat == thermo_nh) then
        !%Variable ThermostatMass
        !%Type float
        !%Default 1.0
        !%Section Time-Dependent::Propagation
        !%Description
        !% This variable sets the fictitious mass for the Nose-Hoover
        !% thermostat.
        !%End
        call messages_obsolete_variable(namespace, 'NHMass', 'ThermostatMass')
 
        call parse_variable(namespace, 'ThermostatMass', m_one, this%nh(1)%mass)
        this%nh(2)%mass = this%nh(1)%mass
 
        this%nh(1:2)%pos = m_zero
        this%nh(1:2)%vel = m_zero
 
        safe_allocate(this%old_pos(1:ions%space%dim, 1:ions%natoms))
 
        this%old_pos = ions%pos
      end if
 
    end if
 
    !now initialize velocities
 
    !%Variable RandomVelocityTemp
    !%Type float
    !%Default 0.0
    !%Section System::Velocities
    !%Description
    !% If this variable is present, <tt>Octopus</tt> will assign random
    !% velocities to the atoms following a Boltzmann distribution with
    !% temperature given by <tt>RandomVelocityTemp</tt> (in degrees Kelvin).
    !% The seed for the random number generator can be modified by setting
    !% <tt>GSL_RNG_SEED</tt> environment variable.
    !%End
 
    ! we now load the velocities, either from the temperature, from the input, or from a file
    if (parse_is_defined(namespace, 'RandomVelocityTemp')) then
 
      have_velocities = .true.
 
      if (mpi_world%is_root()) then
        call loct_ran_init(random_gen_pointer)
        call parse_variable(namespace, 'RandomVelocityTemp', m_zero, temperature, unit = unit_kelvin)
      end if
 
      do i = 1, ions%natoms
        !generate the velocities in the root node
        if (mpi_world%is_root()) then
          sigma = sqrt(temperature / ions%mass(i))
          do j = 1, 3
            ions%vel(j, i) = loct_ran_gaussian(random_gen_pointer, sigma)
          end do
        end if
        !and send them to the others
        call mpi_world%bcast(ions%vel(:, i), ions%space%dim, mpi_double_precision, 0)
      end do
 
      if (mpi_world%is_root()) then
        call loct_ran_end(random_gen_pointer)
      end if
 
      call ions%update_kinetic_energy()
      kin1 = ions%kinetic_energy
 
      xx = ions%center_of_mass_vel()
      do i = 1, ions%natoms
        ions%vel(:, i) = ions%vel(:, i) - xx
      end do
 
      call ions%update_kinetic_energy()
      kin2 = ions%kinetic_energy
 
      do i = 1, ions%natoms
        ions%vel(:, i) =  sqrt(kin1/kin2)*ions%vel(:, i)
      end do
 
      call ions%update_kinetic_energy()
 
      write(message(1),'(a,f10.4,1x,a)') 'Info: Initial velocities randomly distributed with T =', &
        units_from_atomic(unit_kelvin, temperature), units_abbrev(unit_kelvin)
      write(message(2),'(2x,a,f8.4,1x,a)') '<K>       =', &
        units_from_atomic(units_out%energy, ions%kinetic_energy/ions%natoms), &
        units_abbrev(units_out%energy)
      write(message(3),'(2x,a,f8.4,1x,a)') '3/2 k_B T =', &
        units_from_atomic(units_out%energy, (m_three/m_two)*temperature), &
        units_abbrev(units_out%energy)
      call messages_info(3, namespace=namespace)
 
    else
      !%Variable XYZVelocities
      !%Type string
      !%Section System::Velocities
      !%Description
      !% <tt>Octopus</tt> will try to read the starting velocities of the atoms from the XYZ file
      !% specified by the variable <tt>XYZVelocities</tt>.
      !% Note that you do not need to specify initial velocities if you are not going
      !% to perform ion dynamics; if you are going to allow the ions to move but the velocities
      !% are not specified, they are considered to be null.
      !% Note: It is important for the velocities to maintain the ordering
      !% in which the atoms were defined in the coordinates specifications.
      !%End
 
      !%Variable XSFVelocities
      !%Type string
      !%Section System::Velocities
      !%Description
      !% Like <tt>XYZVelocities</tt> but in XCrySDen format, as in <tt>XSFCoordinates</tt>.
      !%End
 
      !%Variable PDBVelocities
      !%Type string
      !%Section System::Velocities
      !%Description
      !% Like <tt>XYZVelocities</tt> but in PDB format, as in <tt>PDBCoordinates</tt>.
      !%End
 
      !%Variable Velocities
      !%Type block
      !%Section System::Velocities
      !%Description
      !% If <tt>XYZVelocities</tt>, <tt>PDBVelocities</tt>, and <tt>XSFVelocities</tt>
      !% are not present, <tt>Octopus</tt> will try to fetch the initial
      !% atomic velocities from this block. If this block is not present, <tt>Octopus</tt>
      !% will set the initial velocities to zero. The format of this block can be
      !% illustrated by this example:
      !%
      !% <tt>%Velocities
      !% <br>&nbsp;&nbsp;'C'  |      -1.7 | 0.0 | 0.0
      !% <br>&nbsp;&nbsp;'O'  | &nbsp;1.7 | 0.0 | 0.0
      !% <br>%</tt>
      !%
      !% It describes one carbon and one oxygen moving at the relative
      !% velocity of 3.4 velocity units.
      !%
      !% Note: It is important for the velocities to maintain the ordering
      !% in which the atoms were defined in the coordinates specifications.
      !%End
 
      call read_coords_init(xyz)
      call read_coords_read('Velocities', xyz, ions%space, namespace)
      if (xyz%source /= read_coords_err) then
 
        have_velocities = .true.
 
        if (ions%natoms /= xyz%n) then
          write(message(1), '(a,i4,a,i4)') 'I need exactly ', ions%natoms, ' velocities, but I found ', xyz%n
          call messages_fatal(1, namespace=namespace)
        end if
 
        ! copy information
        do i = 1, ions%natoms
          ions%vel(:, i) = xyz%atom(i)%x(1:ions%space%dim)
        end do
 
        call read_coords_end(xyz)
 
      else
        ions%vel = m_zero
      end if
    end if
 
    call ions%update_kinetic_energy()
 
    !%Variable MoveIons
    !%Type logical
    !%Section Time-Dependent::Propagation
    !%Description
    !% This variable controls whether atoms are moved during a time
    !% propagation run. The default is yes when the ion velocity is
    !% set explicitly or implicitly, otherwise is no.
    !%End
    call parse_variable(namespace, 'MoveIons', have_velocities, this%move_ions)
    call messages_print_var_value('MoveIons', this%move_ions, namespace=namespace)
 
    if (this%move_ions .and. ions%space%periodic_dim == 1) then
      call messages_input_error(namespace, 'MoveIons', &
        'Moving ions for a 1D periodic system is not allowed, as forces are incorrect.')
    end if
 
    if (this%ions_move()) then
      safe_allocate(this%oldforce(1:ions%space%dim, 1:ions%natoms))
    end if
 
    if (ions%space%is_periodic()) then
      !%Variable CellDynamics
      !%Type logical
      !%Section Time-Dependent::Propagation
      !%Description
      !% This variable controls whether the cell relaxation is done during a time
      !% propagation run. The default is no.
      !% This is done based on the Parrinello-Rahman equation of motion of the cell,
      !% see Parrinello and Rahman, J. Appl. Pys. 52, 7182 (1981).
      !% Only for periodic systems.
      !%End
      call parse_variable(namespace, 'CellDynamics', .false., this%relax_cell)
      call messages_print_var_value('CellDynamics', this%relax_cell, namespace=namespace)
 
      if (this%cell_relax()) then
        if (accel_is_enabled()) then
          message(1) = "Cell dynamics not supported on GPUs."
          call messages_fatal(1, namespace=namespace)
        end if
 
        periodic_dim = ions%space%periodic_dim
        ncomp = periodic_dim * periodic_dim
        safe_allocate(this%cell_force(1:ncomp))
        this%cell_force = m_zero
        safe_allocate(this%old_cell_force(1:ncomp))
        this%old_cell_force = m_zero
        safe_allocate(this%cell_vel(1:ncomp))
        this%cell_vel = m_zero
        ! We start from identity for the strain
        safe_allocate(this%strain(1:ncomp))
        this%strain = m_zero
        ncomp = 1
        do i = 1, periodic_dim
          do j = i, periodic_dim
            if(i == j) this%strain(ncomp) = m_one
            ncomp = ncomp + 1
          end do
        end do
 
        ! As we work with "reduced" quantities, the initial lattice vectors are used as a reference
        ! and the strain is propagated from this initial reference
        safe_allocate(this%initial_rlattice(1:periodic_dim, 1:periodic_dim))
        this%initial_rlattice(1:periodic_dim, 1:periodic_dim) = ions%latt%rlattice(1:periodic_dim, 1:periodic_dim)
      end if
 
      !%Variable HydrostaticPressure
      !%Type float
      !%Default 0.0
      !%Section Time-Dependent::Propagation
      !%Description
      !% Geometry optimization and molecular dynamics can be performed in presence of an external
      !% hydrostatic pressure. This variable allows to set this value.
      !% Only for periodic systems.
      !%End
      call parse_variable(namespace, 'HydrostaticPressure', m_zero, this%pressure)
    else
      this%relax_cell = .false.
    end if
 
 
    pop_sub(ion_dynamics_init)
  end subroutine ion_dynamics_init
 
 
  ! ---------------------------------------------------------
  subroutine ion_dynamics_end(this)
    type(ion_dynamics_t), intent(inout) :: this
 
    push_sub(ion_dynamics_end)
    safe_deallocate_a(this%oldforce)
 
    if (this%thermostat /= thermo_none) then
      call tdf_end(this%temperature_function)
    end if
 
    if (this%drive_ions .and. allocated(this%td_displacements)) then
      if (any(this%td_displacements(1:this%ions_t0%natoms)%move)) then
        ! ions end cannot be called here, otherwise the species are destroyed twice
        ! call ions_end(this%ions_t0)
        !AFE_DEALLOCATE_P(this%ions_t0)
      end if
      safe_deallocate_a(this%td_displacements)
    end if
 
    safe_deallocate_a(this%cell_force)
    safe_deallocate_a(this%old_cell_force)
    safe_deallocate_a(this%cell_vel)
    safe_deallocate_a(this%initial_rlattice)
 
    pop_sub(ion_dynamics_end)
  end subroutine ion_dynamics_end
 
 
  ! ---------------------------------------------------------
  subroutine ion_dynamics_propagate(this, ions, time, dt, namespace)
    type(ion_dynamics_t), intent(inout) :: this
    type(ions_t),         intent(inout) :: ions
    real(real64),         intent(in)    :: time
    real(real64),         intent(in)    :: dt
    type(namespace_t),    intent(in)    :: namespace
 
    integer      :: iatom
 
    push_sub(ion_dynamics_propagate)
 
    this%dt = dt
 
    if (this%drive_ions) then
 
      call ion_dynamics_propagate_driven_ions(this, ions, time, dt)
 
    else
 
      ! Update the thermostat temperature
      call ion_dynamics_update_temperature(this, time, namespace)
 
      ! Conversion to reduced coordinates
      do iatom = 1, ions%natoms
        ions%pos(:, iatom) = ions%latt%cart_to_red(ions%pos(:, iatom))
        ions%vel(:, iatom) = ions%latt%cart_to_red(ions%vel(:, iatom))
        ions%tot_force(:, iatom) = ions%latt%cart_to_red(ions%tot_force(:, iatom))
      end do
 
      ! Get the new reduced coordinates
      if (this%ions_move()) then
        call ion_dynamics_propagate_ions(this, ions, dt)
      end if
 
      ! Updating the lattice vectors
      if (this%cell_relax()) then
        call ion_dynamics_propagate_cell(this, ions, dt, namespace)
      end if
 
      ! Get the new Cartesian coordinates
      do iatom = 1, ions%natoms
        ions%pos(:, iatom) = ions%latt%red_to_cart(ions%pos(:, iatom))
        ions%vel(:, iatom) = ions%latt%red_to_cart(ions%vel(:, iatom))
        if (allocated(this%oldforce)) then
          this%oldforce(:, iatom) = ions%latt%red_to_cart(this%oldforce(:, iatom))
        end if
        ions%tot_force(:, iatom) = ions%latt%red_to_cart(ions%tot_force(:, iatom))
      end do
 
    end if
 
    call ions%fold_atoms_into_cell()
 
    pop_sub(ion_dynamics_propagate)
  end subroutine ion_dynamics_propagate
 
 
  ! ---------------------------------------------------------
  subroutine ion_dynamics_update_temperature(this, time, namespace)
    type(ion_dynamics_t), intent(inout) :: this
    real(real64),         intent(in)    :: time
    type(namespace_t),    intent(in)    :: namespace
 
    push_sub(ion_dynamics_update_temperature)
 
    ! get the temperature from the tdfunction for the current time
    if (this%thermostat /= thermo_none) then
      this%current_temperature = units_to_atomic(unit_kelvin, tdf(this%temperature_function, time))
 
      if (this%current_temperature < m_zero) then
        write(message(1), '(a, f10.3, 3a, f10.3, 3a)') &
          "Negative temperature (", &
          units_from_atomic(unit_kelvin, this%current_temperature), " ", units_abbrev(unit_kelvin), &
          ") at time ", &
          units_from_atomic(units_out%time, time), " ", trim(units_abbrev(units_out%time)), "."
        call messages_fatal(1, namespace=namespace)
      end if
    else
      this%current_temperature = m_zero
    end if
 
    pop_sub(ion_dynamics_update_temperature)
  end subroutine ion_dynamics_update_temperature
 
  subroutine ion_dynamics_propagate_driven_ions(this, ions, time, dt)
    type(ion_dynamics_t), intent(inout) :: this
    type(ions_t),         intent(inout) :: ions
    real(real64),         intent(in)    :: time
    real(real64),         intent(in)    :: dt
 
    integer :: iatom
    real(real64) :: dr(3)
 
    push_sub(ion_dynamics_propagate_driven_ions)
 
    assert(this%drive_ions)
 
    do iatom = 1, ions%natoms
      if (ions%fixed(iatom)) cycle
 
      if (this%constant_velocity) then
        ions%pos(:, iatom) = ions%pos(:, iatom) + dt*ions%vel(:, iatom)
 
      else if (allocated(this%td_displacements)) then
 
        if (this%td_displacements(iatom)%move) then
          dr(1:3)=(/ real(tdf(this%td_displacements(iatom)%fx, time), real64), &
            real(tdf(this%td_displacements(iatom)%fy, time), real64), &
            real(tdf(this%td_displacements(iatom)%fz, time), real64) /)
 
          ions%pos(:, iatom) = this%ions_t0%pos(:, iatom) + dr(1:ions%space%dim)
        end if
 
      end if
    end do
 
    pop_sub(ion_dynamics_propagate_driven_ions)
  end subroutine ion_dynamics_propagate_driven_ions
 
  ! ---------------------------------------------------------
  subroutine ion_dynamics_propagate_ions(this, ions, dt)
    type(ion_dynamics_t), intent(inout) :: this
    type(ions_t),         intent(inout) :: ions
    real(real64),         intent(in)    :: dt
 
    integer      :: iatom
 
    push_sub(ion_dynamics_propagate_ions)
 
    assert(.not. this%drive_ions)
 
    if (this%thermostat /= thermo_nh) then
      ! integrate using verlet
      do iatom = 1, ions%natoms
        if (ions%fixed(iatom)) cycle
 
        ions%pos(:, iatom) = ions%pos(:, iatom) + dt*ions%vel(:, iatom) + &
          m_half*dt**2 / ions%mass(iatom) * ions%tot_force(:, iatom)
 
        this%oldforce(:, iatom) = ions%tot_force(:, iatom)
      end do
 
    else
      ! for the Nose-Hoover thermostat we use a special integrator
 
      ! The implementation of the Nose-Hoover thermostat is based on
      ! Understanding Molecular Simulations by Frenkel and Smit,
      ! Appendix E, page 540-542.
 
      call nh_chain(this, ions)
 
      do iatom = 1, ions%natoms
        if (ions%fixed(iatom)) cycle
 
        ions%pos(:, iatom) = ions%pos(:, iatom) + m_half*dt*ions%vel(:, iatom)
      end do
 
    end if
 
    pop_sub(ion_dynamics_propagate_ions)
  end subroutine ion_dynamics_propagate_ions
 
  ! ---------------------------------------------------------
  subroutine ion_dynamics_propagate_cell(this, ions, dt, namespace)
    type(ion_dynamics_t), intent(inout) :: this
    type(ions_t),         intent(inout) :: ions
    real(real64),         intent(in)    :: dt
    type(namespace_t),    intent(in)    :: namespace
 
    integer      :: idir, jdir, comp
    real(real64) :: rlattice_change(ions%space%periodic_dim*ions%space%periodic_dim)
 
    push_sub(ion_dynamics_propagate_cell)
 
    rlattice_change = dt * this%cell_vel + m_half*dt**2 * this%cell_force
 
    write(message(1),'(a,3a,a)') '  Cell force [', trim(units_abbrev(units_out%energy/units_out%length**ions%space%dim)), ']'
    do idir = 1, ions%space%periodic_dim
      write(message(1+idir),'(9e18.10)') (units_from_atomic(units_out%energy/units_out%length**ions%space%dim, &
        this%cell_force(jdir + (idir-1)*ions%space%periodic_dim)), &
        jdir = 1, ions%space%periodic_dim)
    end do
    call messages_info(1+ions%space%periodic_dim, namespace=ions%namespace)
 
    write(message(1),'(a,3a,a)') '  Cell vel [', &
      trim(units_abbrev(units_out%energy/units_out%length**ions%space%dim*units_out%time)), ']'
    do idir = 1, ions%space%periodic_dim
      write(message(1+idir),'(9e18.10)') (units_from_atomic(units_out%energy/units_out%length**ions%space%dim*units_out%time, &
        this%cell_vel(jdir+ (idir-1)*ions%space%periodic_dim)), &
        jdir = 1, ions%space%periodic_dim)
    end do
    call messages_info(1+ions%space%periodic_dim, namespace=ions%namespace)
 
 
    comp = 1
    do idir = 1, ions%space%periodic_dim
      do jdir = 1, ions%space%periodic_dim
        ions%latt%rlattice(idir, jdir) = ions%latt%rlattice(idir, jdir) + rlattice_change(comp)
        comp = comp + 1
      end do
    end do
 
    this%old_cell_force = this%cell_force
 
    if (associated(this%symm)) then
      call this%symm%symmetrize_lattice_vectors(ions%space%periodic_dim, &
        this%initial_rlattice, ions%latt%rlattice, this%symmetrize)
    end if
    call ions%update_lattice_vectors(ions%latt, this%symmetrize)
 
    pop_sub(ion_dynamics_propagate_cell)
  end subroutine ion_dynamics_propagate_cell
 
 
  ! ---------------------------------------------------------
  subroutine nh_chain(this, ions)
    type(ion_dynamics_t), intent(inout) :: this
    type(ions_t),         intent(inout) :: ions
 
    real(real64) :: g1, g2, ss, uk, dt, temp
 
    push_sub(nh_chain)
 
    dt = this%dt
 
    call ions%update_kinetic_energy()
    uk = ions%kinetic_energy
 
    temp = this%current_temperature
 
    g2 = (this%nh(1)%mass*this%nh(1)%vel**2 - temp)/this%nh(2)%mass
    this%nh(2)%vel = this%nh(2)%vel + g2*dt/m_four
    this%nh(1)%vel = this%nh(1)%vel*exp(-this%nh(2)%vel*dt/8.0_real64)
 
    g1 = (m_two*uk - m_three*ions%natoms*temp)/this%nh(1)%mass
    this%nh(1)%vel = this%nh(1)%vel + g1*dt/m_four
    this%nh(1)%vel = this%nh(1)%vel*exp(-this%nh(2)%vel*dt/8.0_real64)
    this%nh(1)%pos = this%nh(1)%pos + this%nh(1)%vel*dt/m_two
    this%nh(2)%pos = this%nh(2)%pos + this%nh(2)%vel*dt/m_two
 
    ss = exp(-this%nh(1)%vel*dt/m_two)
 
    ions%vel = ss*ions%vel
 
    uk = uk*ss**2
 
    this%nh(1)%vel = this%nh(1)%vel*exp(-this%nh(2)%vel*dt/8.0_real64)
    g1 = (m_two*uk - m_three*ions%natoms*temp)/this%nh(1)%mass
    this%nh(1)%vel = this%nh(1)%vel + g1*dt/m_four
    this%nh(1)%vel = this%nh(1)%vel*exp(-this%nh(2)%vel*dt/8.0_real64)
 
    g2 = (this%nh(1)%mass*this%nh(1)%vel**2 - temp)/this%nh(2)%mass
    this%nh(2)%vel = this%nh(2)%vel + g2*dt/m_four
 
    pop_sub(nh_chain)
  end subroutine nh_chain
 
 
  ! ---------------------------------------------------------
  subroutine ion_dynamics_propagate_vel(this, ions, atoms_moved)
    type(ion_dynamics_t), intent(inout) :: this
    type(ions_t),         intent(inout) :: ions
    logical, optional,    intent(out)   :: atoms_moved
 
    integer      :: iatom
    real(real64) :: scal
 
    if (.not. ion_dynamics_ions_move(this)) return
    if (this%drive_ions) return
 
    push_sub(ion_dynamics_propagate_vel)
 
    if (present(atoms_moved)) atoms_moved = this%thermostat == thermo_nh
 
    if (this%thermostat /= thermo_nh) then
      ! velocity verlet
 
      do iatom = 1, ions%natoms
        if (ions%fixed(iatom)) cycle
 
        ions%vel(:, iatom) = ions%vel(:, iatom) &
          + this%dt/ions%mass(iatom) * m_half * (this%oldforce(:, iatom) + &
          ions%tot_force(:, iatom))
 
      end do
 
    else
      ! the nose-hoover integration
      do iatom = 1, ions%natoms
        ions%vel(:, iatom) = ions%vel(:, iatom) + this%dt*ions%tot_force(:, iatom) / ions%mass(iatom)
        ions%pos(:, iatom) = ions%pos(:, iatom) + m_half*this%dt*ions%vel(:, iatom)
      end do
 
      call nh_chain(this, ions)
 
    end if
 
    if (this%thermostat == thermo_scal) then
      scal = sqrt(this%current_temperature/ion_dynamics_temperature(ions))
 
      ions%vel = scal*ions%vel
    end if
 
    if (this%cell_relax()) then
      this%cell_vel = this%cell_vel + this%dt * m_half * (this%old_cell_force + this%cell_force)
    end if
 
    pop_sub(ion_dynamics_propagate_vel)
  end subroutine ion_dynamics_propagate_vel
 
 
  ! ---------------------------------------------------------
  subroutine ion_dynamics_verlet_step1(ions, q, v, fold, dt)
    type(ions_t),         intent(in)    :: ions
    real(real64),         intent(inout) :: q(:, :)
    real(real64),         intent(inout) :: v(:, :)
    real(real64),         intent(in)    :: fold(:, :)
    real(real64),         intent(in)    :: dt
 
    integer :: iatom
 
    push_sub(ion_dynamics_verlet_step1)
 
    ! First transform momenta to velocities
    do iatom = 1, ions%natoms
      v(iatom, 1:ions%space%dim) = v(iatom, 1:ions%space%dim) / ions%mass(iatom)
    end do
 
    ! integrate using verlet
    do iatom = 1, ions%natoms
      if (ions%fixed(iatom)) cycle
      q(iatom, 1:ions%space%dim) = q(iatom, 1:ions%space%dim) + dt * v(iatom, 1:ions%space%dim) + &
        m_half*dt**2 / ions%mass(iatom) * fold(iatom, 1:ions%space%dim)
    end do
 
    ! And back to momenta.
    do iatom = 1, ions%natoms
      v(iatom, 1:ions%space%dim) = ions%mass(iatom) * v(iatom, 1:ions%space%dim)
    end do
 
    pop_sub(ion_dynamics_verlet_step1)
  end subroutine ion_dynamics_verlet_step1
 
 
 
  ! ---------------------------------------------------------
  subroutine ion_dynamics_verlet_step2(ions, v, fold, fnew, dt)
    type(ions_t),         intent(in)    :: ions
    real(real64),         intent(inout) :: v(:, :)
    real(real64),         intent(in)    :: fold(:, :)
    real(real64),         intent(in)    :: fnew(:, :)
    real(real64),         intent(in)    :: dt
 
    integer :: iatom
 
    push_sub(ion_dynamics_verlet_step2)
 
    ! First transform momenta to velocities
    do iatom = 1, ions%natoms
      v(iatom, 1:ions%space%dim) = v(iatom, 1:ions%space%dim) / ions%mass(iatom)
    end do
 
    ! velocity verlet
    do iatom = 1, ions%natoms
      if (ions%fixed(iatom)) cycle
      v(iatom, 1:ions%space%dim) = v(iatom, 1:ions%space%dim) &
        + dt / ions%mass(iatom) * m_half * (fold(iatom, 1:ions%space%dim) + fnew(iatom, 1:ions%space%dim))
    end do
 
    ! And back to momenta.
    do iatom = 1, ions%natoms
      v(iatom, 1:ions%space%dim) = ions%mass(iatom) * v(iatom, 1:ions%space%dim)
    end do
 
    pop_sub(ion_dynamics_verlet_step2)
  end subroutine ion_dynamics_verlet_step2
 
 
  ! ---------------------------------------------------------
  subroutine ion_dynamics_save_state(this, ions, state)
    type(ion_dynamics_t), intent(in)    :: this
    type(ions_t),         intent(in)    :: ions
    type(ion_state_t),    intent(out)   :: state
 
    if (.not. this%ions_move()) return
 
    push_sub(ion_dynamics_save_state)
 
    safe_allocate(state%pos(1:ions%space%dim, 1:ions%natoms))
    safe_allocate(state%vel(1:ions%space%dim, 1:ions%natoms))
 
    state%pos = ions%pos
    state%vel = ions%vel
 
    if (this%thermostat == thermo_nh) then
      safe_allocate(state%old_pos(1:ions%space%dim, 1:ions%natoms))
      state%old_pos(1:ions%space%dim, 1:ions%natoms) = this%old_pos(1:ions%space%dim, 1:ions%natoms)
      state%nh(1:2)%pos = this%nh(1:2)%pos
      state%nh(1:2)%vel = this%nh(1:2)%vel
    end if
 
    pop_sub(ion_dynamics_save_state)
  end subroutine ion_dynamics_save_state
 
 
  ! ---------------------------------------------------------
  subroutine ion_dynamics_restore_state(this, ions, state)
    type(ion_dynamics_t), intent(inout) :: this
    type(ions_t),         intent(inout) :: ions
    type(ion_state_t),    intent(inout) :: state
 
    assert(.not. this%cell_relax())
    if (.not. this%ions_move()) return
 
    push_sub(ion_dynamics_restore_state)
 
    ions%pos = state%pos
    ions%vel = state%vel
 
    if (this%thermostat == thermo_nh) then
      this%old_pos(1:ions%space%dim, 1:ions%natoms) = state%old_pos(1:ions%space%dim, 1:ions%natoms)
      this%nh(1:2)%pos = state%nh(1:2)%pos
      this%nh(1:2)%vel = state%nh(1:2)%vel
      safe_deallocate_a(state%old_pos)
    end if
 
    safe_deallocate_a(state%pos)
    safe_deallocate_a(state%vel)
 
    pop_sub(ion_dynamics_restore_state)
  end subroutine ion_dynamics_restore_state
 
 
  ! ---------------------------------------------------------
  logical pure function ion_dynamics_ions_move(this) result(ions_move)
    class(ion_dynamics_t), intent(in)    :: this
 
    ions_move = this%move_ions
 
  end function ion_dynamics_ions_move
 
 
  ! ---------------------------------------------------------
  logical pure function ion_dynamics_drive_ions(this) result(drive_ions)
    type(ion_dynamics_t), intent(in)    :: this
 
    drive_ions = this%drive_ions
 
  end function ion_dynamics_drive_ions
 
  ! ---------------------------------------------------------
  logical pure function ion_dynamics_cell_relax(this) result(cell_dynamics)
    class(ion_dynamics_t), intent(in)    :: this
 
    cell_dynamics = this%relax_cell
 
  end function ion_dynamics_cell_relax
 
  ! ---------------------------------------------------------
  logical pure function ion_dynamics_is_active(this) result(is_active)
    class(ion_dynamics_t), intent(in)    :: this
 
    is_active = this%relax_cell .or. this%move_ions
 
  end function ion_dynamics_is_active
 
 
  ! ---------------------------------------------------------
  real(real64) function ion_dynamics_temperature(ions) result(temperature)
    type(ions_t),          intent(in) :: ions
 
    integer :: iatom
    real(real64) :: kinetic_energy
 
    kinetic_energy = m_zero
    do iatom = 1, ions%natoms
      kinetic_energy = kinetic_energy + &
        m_half * ions%mass(iatom) * sum(ions%vel(:, iatom)**2)
    end do
    temperature = m_two/m_three*kinetic_energy/ions%natoms
 
  end function ion_dynamics_temperature
 
 
  ! ---------------------------------------------------------
  logical function ion_dynamics_freeze(this) result(freeze)
    type(ion_dynamics_t), intent(inout)   :: this
    if (this%move_ions) then
      this%move_ions = .false.
      freeze = .true.
    else
      freeze = .false.
    end if
  end function ion_dynamics_freeze
 
 
  ! ---------------------------------------------------------
  subroutine ion_dynamics_unfreeze(this)
    type(ion_dynamics_t), intent(inout)   :: this
    this%move_ions = .true.
  end subroutine ion_dynamics_unfreeze
 
  ! ---------------------------------------------------------
  subroutine ion_dynamics_dump(this, restart, ierr)
    type(ion_dynamics_t), intent(in)  :: this
    type(restart_t),      intent(in)  :: restart
    integer,              intent(out) :: ierr
 
    push_sub(ion_dynamics_dump)
 
    if (allocated(this%oldforce)) then
      call restart%write_binary("ion_dynamics_oldforce", size(this%oldforce), &
        this%oldforce, ierr)
    end if
 
    if( allocated(this%old_cell_force)) then
      call restart%write_binary("ion_dynamics_old_cell_force", size(this%old_cell_force), &
        this%old_cell_force, ierr)
    end if
 
    pop_sub(ion_dynamics_dump)
  end subroutine ion_dynamics_dump
 
  ! ---------------------------------------------------------
  subroutine ion_dynamics_load(this, restart, ierr)
    type(ion_dynamics_t), intent(inout) :: this
    type(restart_t),      intent(in)    :: restart
    integer,              intent(out)   :: ierr
 
    push_sub(ion_dynamics_load)
 
    if (allocated(this%oldforce)) then
      call restart%read_binary("ion_dynamics_oldforce", size(this%oldforce), &
        this%oldforce, ierr)
    end if
 
    if( allocated(this%old_cell_force)) then
      call restart%read_binary("ion_dynamics_old_cell_force", size(this%old_cell_force), &
        this%old_cell_force, ierr)
    end if
 
    pop_sub(ion_dynamics_load)
  end subroutine ion_dynamics_load
 
  !----------------------------------------------------------
  subroutine ion_dynamics_update_stress(this, space, stress, rlattice, rcell_volume)
    class(ion_dynamics_t), intent(inout) :: this
    class(space_t),        intent(in)    :: space
    real(real64),          intent(in)    :: stress(3,3)
    real(real64),          intent(in)    :: rlattice(:,:)
    real(real64),          intent(in)    :: rcell_volume
 
    integer :: idir, jdir, comp
    real(real64) :: inv_latt(space%periodic_dim, space%periodic_dim), tmp_stress(space%periodic_dim, space%periodic_dim)
    real(real64) :: cell_force(space%periodic_dim, space%periodic_dim)
 
    push_sub(ion_dynamics_update_stress)
 
    ! Get the inverse lattice
    inv_latt = rlattice(1:space%periodic_dim, 1:space%periodic_dim)
    call lalg_inverse(space%periodic_dim, inv_latt, 'dir')
 
    tmp_stress = -stress(1:space%periodic_dim, 1:space%periodic_dim)
    do idir = 1, space%periodic_dim
      tmp_stress(idir, idir) = tmp_stress(idir, idir) - this%pressure/space%periodic_dim
    end do
    cell_force = matmul(tmp_stress, transpose(inv_latt)) * rcell_volume
 
    comp = 1
    do idir = 1, space%periodic_dim
      do jdir = 1, space%periodic_dim
        this%cell_force(comp) = cell_force(idir, jdir)
        comp = comp + 1
      end do
    end do
 
    if (debug%info) then
      write(message(1),'(a,3a,a)') '  Stress tensor [', trim(units_abbrev(units_out%energy/units_out%length**space%dim)), ']'
      do idir = 1, space%periodic_dim
        write(message(1+idir),'(9e18.10)') (units_from_atomic(units_out%energy/units_out%length**space%dim, stress(jdir, idir)), &
          jdir = 1, space%periodic_dim)
      end do
      call messages_info(1+space%periodic_dim, namespace=global_namespace)
    end if
 
    pop_sub(ion_dynamics_update_stress)
  end subroutine ion_dynamics_update_stress
 
  !----------------------------------------------------------
  subroutine ion_dynamics_box_update(namespace, gr, space, new_latt)
    type(namespace_t),        intent(in)    :: namespace
    type(grid_t),             intent(inout) :: gr
    class(space_t),           intent(in)    :: space
    type(lattice_vectors_t),  intent(in)    :: new_latt
 
    integer :: idir
    real(real64) :: length(1:space%dim)
 
    push_sub(ion_dynamics_box_update)
 
    ! Regenerate the box
    select type(box => gr%box)
    type is (box_parallelepiped_t)
      do idir = 1, space%dim
        length(idir) = norm2(new_latt%rlattice(1:space%dim, idir))
      end do
      call box%regenerate(space%dim, new_latt%rlattice, length, namespace)
    class default
      call messages_not_implemented("Grid regeneration for non-parallelepiped boxes", namespace=namespace)
    end select
 
    pop_sub(ion_dynamics_box_update)
  end subroutine ion_dynamics_box_update
 
  !----------------------------------------------------------
  subroutine electrons_lattice_vectors_update(namespace, gr, space, psolver, kpoints, mc, qtot, new_latt)
    type(namespace_t),        intent(in)    :: namespace
    type(grid_t),             intent(inout) :: gr
    class(space_t),           intent(in)    :: space
    type(poisson_t),          intent(inout) :: psolver
    type(kpoints_t),          intent(inout) :: kpoints
    type(multicomm_t),        intent(in)    :: mc
    real(real64),             intent(in)    :: qtot
    type(lattice_vectors_t),  intent(in)    :: new_latt
 
    push_sub(electrons_lattice_vectors_update)
 
    call grid_lattice_vectors_update(gr, space, namespace, mc, new_latt)
 
    !Initialize Poisson solvers
    call poisson_end(psolver)
    call poisson_init(psolver, namespace, space, gr%der, mc, gr%stencil, qtot, verbose=.false.)
 
    call kpoints_lattice_vectors_update(kpoints, new_latt)
 
    pop_sub(electrons_lattice_vectors_update)
 
  end subroutine electrons_lattice_vectors_update
 
 
end module ion_dynamics_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: