classical_particles.F90

Go to the documentation of this file.

!! Copyright (C) 2002-2006 M. Marques, A. Castro, A. Rubio, G. Bertsch
!! Copyright (C) 2019 N. Tancogne-Dejean
!! Copyright (C) 2020 M. Oliveira
!! Copyright (C) 2021 S. Ohlmann, I. Albar, A. Obzhirov, A. Geondzhian, M. Lawan
!! Copyright (C) 2022 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 classical_particles_oct_m
  use algorithm_oct_m
  use debug_oct_m
  use force_interaction_oct_m
  use global_oct_m
  use interaction_surrogate_oct_m
  use interaction_oct_m
  use io_oct_m
  use, intrinsic :: iso_fortran_env
  use lalg_adv_oct_m
  use messages_oct_m
  use mpi_oct_m
  use multisystem_debug_oct_m
  use namespace_oct_m
  use parser_oct_m
  use profiling_oct_m
  use propagator_data_classical_particles_oct_m
  use propagator_beeman_oct_m
  use propagator_exp_mid_2step_oct_m
  use propagator_oct_m
  use propagator_verlet_oct_m
  use quantity_oct_m
  use space_oct_m
  use system_oct_m
  use unit_oct_m
  use unit_system_oct_m
  use utils_oct_m
  use write_iter_oct_m
 
  implicit none
 
  private
  public ::                                             &
    classical_particles_t,                              &
    classical_particles_init,                           &
    classical_particles_copy,                           &
    classical_particles_end,                            &
    classical_particles_init_interaction,               &
    classical_particles_restart_write_data,             &
    classical_particles_update_quantity,                &
    classical_particles_init_interaction_as_partner,    &
    classical_particles_copy_quantities_to_interaction
 
  integer, parameter ::     &
    OUTPUT_COORDINATES = 1, &
    output_energy      = 2
 
  type, extends(system_t), abstract :: classical_particles_t
    private
    type(c_ptr), public :: output_handle(2)
    type(space_t), public :: space
    integer, public :: np
    real(real64), allocatable, public :: mass(:)
    real(real64), allocatable, public :: pos(:,:)
    real(real64), allocatable, public :: vel(:,:)
    real(real64), allocatable, public :: tot_force(:,:)
    real(real64), allocatable, public :: lj_epsilon(:)
    real(real64), allocatable, public :: lj_sigma(:)
    logical, allocatable, public :: fixed(:)
    !                                            !! The default is to let the particles move.
    type(propagator_data_t),public :: prop_data
  contains
    procedure :: do_algorithmic_operation => classical_particles_do_algorithmic_operation
    procedure :: is_tolerance_reached => classical_particles_is_tolerance_reached
    procedure :: copy_quantities_to_interaction => classical_particles_copy_quantities_to_interaction
    procedure :: update_interactions_start => classical_particles_update_interactions_start
    procedure :: update_interactions_finish => classical_particles_update_interactions_finish
    procedure :: output_start => classical_particles_output_start
    procedure :: output_write => classical_particles_output_write
    procedure :: output_finish => classical_particles_output_finish
    procedure :: restart_write_data => classical_particles_restart_write_data
    procedure :: restart_read_data => classical_particles_restart_read_data
    procedure :: update_kinetic_energy => classical_particles_update_kinetic_energy
    procedure :: center_of_mass => classical_particles_center_of_mass
    procedure :: center_of_mass_vel => classical_particles_center_of_mass_vel
    procedure :: center => classical_particles_center
    procedure :: axis_large => classical_particles_axis_large
    procedure :: axis_inertia => classical_particles_axis_inertia
  end type classical_particles_t
 
 
contains
 
  ! ---------------------------------------------------------
  ! ---------------------------------------------------------
  subroutine classical_particles_init(this, np)
    class(classical_particles_t), intent(inout) :: this
    integer,                      intent(in)    :: np
 
    push_sub(classical_particles_init)
 
    this%np = np
    safe_allocate(this%mass(1:np))
    this%mass = m_zero
    safe_allocate(this%pos(1:this%space%dim, 1:np))
    this%pos = m_zero
    safe_allocate(this%vel(1:this%space%dim, 1:np))
    this%vel = m_zero
    safe_allocate(this%tot_force(1:this%space%dim, 1:np))
    this%tot_force = m_zero
    safe_allocate(this%fixed(1:np))
    this%fixed = .false. ! By default we let the particles move.
    safe_allocate(this%lj_epsilon(1:np))
    this%lj_epsilon = m_zero
    safe_allocate(this%lj_sigma(1:np))
    this%lj_sigma = m_zero
 
    ! No interaction directly supported by the classical particles (but classes
    ! that extend this one can add their own)
    allocate(this%supported_interactions(0))
    allocate(this%supported_interactions_as_partner(0))
 
    ! Quantities updated by the algorithm
    call this%quantities%add(quantity_t("position", updated_on_demand = .false.))
    call this%quantities%add(quantity_t("velocity", updated_on_demand = .false.))
 
    ! Other quantities
    call this%quantities%add(quantity_t("mass", updated_on_demand = .false., always_available=.true.))
 
    pop_sub(classical_particles_init)
  end subroutine classical_particles_init
 
  ! ---------------------------------------------------------
  subroutine classical_particles_copy(this, cp_in)
    class(classical_particles_t), intent(out) :: this
    class(classical_particles_t), intent(in)  :: cp_in
 
    push_sub(classical_particles_copy)
 
    this%np = cp_in%np
    safe_allocate_source_a(this%mass,      cp_in%mass)
    safe_allocate_source_a(this%pos,       cp_in%pos)
    safe_allocate_source_a(this%vel,       cp_in%vel)
    safe_allocate_source_a(this%tot_force, cp_in%tot_force)
    safe_allocate_source_a(this%fixed,     cp_in%fixed)
 
    this%kinetic_energy    = cp_in%kinetic_energy
 
    this%quantities = cp_in%quantities
    this%supported_interactions = cp_in%supported_interactions
 
    this%prop_data = cp_in%prop_data
 
    pop_sub(classical_particles_copy)
  end subroutine classical_particles_copy
 
  ! ---------------------------------------------------------
  subroutine classical_particles_init_interaction(this, interaction)
    class(classical_particles_t), target, intent(inout) :: this
    class(interaction_surrogate_t),       intent(inout) :: interaction
 
    push_sub(classical_particles_init_interaction)
 
    select type (interaction)
    class default
      message(1) = "Trying to initialize an unsupported interaction by classical particles."
      call messages_fatal(1, namespace=this%namespace)
    end select
 
    pop_sub(classical_particles_init_interaction)
  end subroutine classical_particles_init_interaction
 
  ! ---------------------------------------------------------
  logical function classical_particles_do_algorithmic_operation(this, operation, updated_quantities) result(done)
    class(classical_particles_t),    intent(inout) :: this
    class(algorithmic_operation_t),  intent(in)    :: operation
    character(len=:), allocatable,   intent(out)   :: updated_quantities(:)
 
    integer :: ii, ip
    real(real64), allocatable :: tmp_pos(:,:,:), tmp_vel(:,:,:)
    real(real64) :: factor
 
    push_sub(classical_particles_do_algorithmic_operation)
    call profiling_in(trim(this%namespace%get())//":"//trim(operation%id))
 
    done = .true.
    select type (prop => this%algo)
    class is (propagator_t)
 
      select case (operation%id)
      case (store_current_status)
        this%prop_data%save_pos(:, 1:this%np) = this%pos(:, 1:this%np)
        this%prop_data%save_vel(:, 1:this%np) = this%vel(:, 1:this%np)
 
      case (verlet_start)
        if (.not. this%prop_data%initialized) then
          call this%prop_data%initialize(prop, this%space%dim, this%np)
          do ip = 1, this%np
            if (this%fixed(ip)) then
              this%prop_data%acc(:, ip) = m_zero
            else
              this%prop_data%acc(:, ip) = this%tot_force(:, ip) / this%mass(ip)
            end if
          end do
        end if
 
      case (verlet_finish)
        call this%prop_data%end()
 
      case (beeman_finish)
        call this%prop_data%end()
 
      case (verlet_update_pos)
        this%pos(:, 1:this%np) = this%pos(:, 1:this%np) + prop%dt * this%vel(:, 1:this%np) &
          + m_half * prop%dt**2 * this%prop_data%acc(:, 1:this%np)
        updated_quantities = ["position"]
 
      case (verlet_compute_acc, beeman_compute_acc)
        do ii = size(this%prop_data%prev_acc, dim=3) - 1, 1, -1
          this%prop_data%prev_acc(:, 1:this%np, ii + 1) = this%prop_data%prev_acc(:, 1:this%np, ii)
        end do
        do ip = 1, this%np
          this%prop_data%prev_acc(:, ip, 1) = this%prop_data%acc(:, ip)
          if (this%fixed(ip)) then
            this%prop_data%acc(:, ip) = m_zero
          else
            this%prop_data%acc(:, ip) = this%tot_force(:, ip) / this%mass(ip)
          end if
        end do
 
      case (verlet_compute_vel)
        this%vel(:, 1:this%np) = this%vel(:, 1:this%np) &
          + m_half * prop%dt * (this%prop_data%prev_acc(:, 1:this%np, 1) + this%prop_data%acc(:, 1:this%np))
        updated_quantities = ["velocity"]
 
      case (beeman_start)
        if (.not. this%prop_data%initialized) then
          call this%prop_data%initialize(prop, this%space%dim, this%np)
          do ip = 1, this%np
            if (this%fixed(ip)) then
              this%prop_data%acc(:, ip) = m_zero
            else
              this%prop_data%acc(:, ip) = this%tot_force(:, ip) / this%mass(ip)
            end if
            this%prop_data%prev_acc(:, ip, 1) = this%prop_data%acc(:, ip)
          end do
        end if
 
      case (beeman_predict_pos)
        this%pos(:, 1:this%np) = this%pos(:, 1:this%np) + prop%dt * this%vel(:, 1:this%np) + &
          m_one/6.0_real64 * prop%dt**2 * (m_four*this%prop_data%acc(:, 1:this%np) - this%prop_data%prev_acc(:, 1:this%np, 1))
 
        if (.not. prop%predictor_corrector) then
          updated_quantities = ["position"]
        end if
 
      case (beeman_predict_vel)
        this%vel(:, 1:this%np) = this%vel(:, 1:this%np) + &
          m_one/6.0_real64 * prop%dt * ( m_two * this%prop_data%acc(:, 1:this%np) + &
          5.0_real64 * this%prop_data%prev_acc(:, 1:this%np, 1) - this%prop_data%prev_acc(:, 1:this%np, 2))
        updated_quantities = ["velocity"]
 
      case (beeman_correct_pos)
        this%pos(:, 1:this%np) = this%prop_data%save_pos(:, 1:this%np) &
          + prop%dt * this%prop_data%save_vel(:, 1:this%np) &
          + m_one/6.0_real64 * prop%dt**2 * (this%prop_data%acc(:, 1:this%np) &
          + m_two * this%prop_data%prev_acc(:, 1:this%np, 1))
        updated_quantities = ["position"]
 
      case (beeman_correct_vel)
        this%vel(:, 1:this%np) = this%prop_data%save_vel(:, 1:this%np) &
          + m_half * prop%dt * (this%prop_data%acc(:, 1:this%np) + this%prop_data%prev_acc(:, 1:this%np, 1))
        updated_quantities = ["velocity"]
 
      case (expmid_2step_start)
        if (.not. this%prop_data%initialized) then
          call this%prop_data%initialize(prop, this%space%dim, this%np)
          this%prop_data%prev_pos(:, 1:this%np, 1) = this%pos(:, 1:this%np)
          this%prop_data%prev_vel(:, 1:this%np, 1) = this%vel(:, 1:this%np)
        end if
 
      case (expmid_2step_finish)
        call this%prop_data%end()
 
      case (expmid_2step_predict_dt_2)
        this%pos(:, 1:this%np) = 1.5_real64*this%prop_data%save_pos(:, 1:this%np) &
          - m_half*this%prop_data%prev_pos(:, 1:this%np, 1)
        this%vel(:, 1:this%np) = 1.5_real64*this%prop_data%save_vel(:, 1:this%np) &
          - m_half*this%prop_data%prev_vel(:, 1:this%np, 1)
        this%prop_data%prev_pos(:, 1:this%np, 1) = this%prop_data%save_pos(:, 1:this%np)
        this%prop_data%prev_vel(:, 1:this%np, 1) = this%prop_data%save_vel(:, 1:this%np)
        updated_quantities = ["position", "velocity"]
 
      case (update_hamiltonian)
        do ip = 1, this%np
          if (this%fixed(ip)) then
            this%prop_data%hamiltonian_elements(:, ip) = m_zero
          else
            this%prop_data%hamiltonian_elements(:, ip) = this%tot_force(:, ip) / (this%mass(ip) * this%pos(:, ip) + m_epsilon)
          end if
        end do
 
      case (expmid_2step_predict_dt)
        safe_allocate(tmp_pos(1:this%space%dim, 1:this%np, 1:2))
        safe_allocate(tmp_vel(1:this%space%dim, 1:this%np, 1:2))
        ! apply exponential - at some point this could use the machinery of
        !   exponential_apply (but this would require a lot of boilerplate code
        !   like a Hamiltonian class etc)
        ! prop_data%save_pos/vel contain the state at t - this is the state we want to
        !   apply the Hamiltonian to
        tmp_pos(:, 1:this%np, 1) = this%prop_data%save_pos(:, 1:this%np)
        tmp_vel(:, 1:this%np, 1) = this%prop_data%save_vel(:, 1:this%np)
        this%pos(:, 1:this%np) = this%prop_data%save_pos(:, 1:this%np)
        this%vel(:, 1:this%np) = this%prop_data%save_vel(:, 1:this%np)
        ! compute exponential with Taylor expansion
        factor = m_one
        do ii = 1, 4
          factor = factor * prop%dt / ii
          do ip = 1, this%np
            ! apply hamiltonian
            tmp_pos(:, ip, 2) = tmp_vel(:, ip, 1)
            tmp_vel(:, ip, 2) = this%prop_data%hamiltonian_elements(:, ip) * tmp_pos(:, ip, 1)
            ! swap temporary variables
            tmp_pos(:, ip, 1) = tmp_pos(:, ip, 2)
            tmp_vel(:, ip, 1) = tmp_vel(:, ip, 2)
            ! accumulate components of Taylor expansion
            this%pos(:, ip) = this%pos(:, ip) + factor * tmp_pos(:, ip, 1)
            this%vel(:, ip) = this%vel(:, ip) + factor * tmp_vel(:, ip, 1)
          end do
        end do
        safe_deallocate_a(tmp_pos)
        safe_deallocate_a(tmp_vel)
        updated_quantities = ["position", "velocity"]
 
      case (expmid_2step_correct_dt_2)
        ! only correct for dt/2 if not converged yet
        if (.not. this%is_tolerance_reached(prop%scf_tol)) then
          this%pos(:, 1:this%np) = m_half*(this%pos(:, 1:this%np) + this%prop_data%save_pos(:, 1:this%np))
          this%vel(:, 1:this%np) = m_half*(this%vel(:, 1:this%np) + this%prop_data%save_vel(:, 1:this%np))
          updated_quantities = ["position", "velocity"]
        end if
 
      case default
        done = .false.
      end select
 
    class default
      done = .false.
    end select
 
    call profiling_out(trim(this%namespace%get())//":"//trim(operation%id))
    pop_sub(classical_particles_do_algorithmic_operation)
  end function classical_particles_do_algorithmic_operation
 
  ! ---------------------------------------------------------
  logical function classical_particles_is_tolerance_reached(this, tol) result(converged)
    class(classical_particles_t), intent(in)    :: this
    real(real64),                 intent(in)    :: tol
 
    integer :: ip
    real(real64) :: change, max_change
 
    push_sub(classical_particles_is_tolerance_reached)
 
    ! Here we put the criterion that maximum acceleration change is below the tolerance
    max_change = m_zero
    do ip = 1, this%np
      change = sum((this%prop_data%prev_tot_force(1:this%space%dim, ip) - this%tot_force(1:this%space%dim, ip))**2) / &
        this%mass(ip)
      if (change > max_change) then
        max_change = change
      end if
    end do
    converged = max_change < tol**2
 
    write(message(1), '(a, e13.6, a, e13.6)') "Debug: -- Maximum change in acceleration  ", &
      sqrt(max_change), " and tolerance ", tol
    call messages_info(1, namespace=this%namespace, debug_only=.true.)
 
    pop_sub(classical_particles_is_tolerance_reached)
  end function classical_particles_is_tolerance_reached
 
  ! ---------------------------------------------------------
  subroutine classical_particles_update_quantity(this, label)
    class(classical_particles_t), intent(inout) :: this
    character(len=*),             intent(in)    :: label
 
    push_sub(classical_particles_update_quantity)
 
    select case (label)
    case default
      message(1) = "Incompatible quantity."
      call messages_fatal(1, namespace=this%namespace)
    end select
 
    pop_sub(classical_particles_update_quantity)
  end subroutine classical_particles_update_quantity
 
  ! ---------------------------------------------------------
  subroutine classical_particles_init_interaction_as_partner(partner, interaction)
    class(classical_particles_t), intent(in)    :: partner
    class(interaction_surrogate_t), intent(inout) :: interaction
 
    push_sub(classical_particles_init_interaction_as_partner)
 
    select type (interaction)
    class default
      message(1) = "Unsupported interaction."
      call messages_fatal(1, namespace=partner%namespace)
    end select
 
    pop_sub(classical_particles_init_interaction_as_partner)
  end subroutine classical_particles_init_interaction_as_partner
 
  ! ---------------------------------------------------------
  subroutine classical_particles_copy_quantities_to_interaction(partner, interaction)
    class(classical_particles_t),         intent(inout) :: partner
    class(interaction_surrogate_t),       intent(inout) :: interaction
 
    push_sub(classical_particles_copy_quantities_to_interaction)
 
    select type (interaction)
    class default
      message(1) = "Unsupported interaction."
      call messages_fatal(1, namespace=partner%namespace)
    end select
 
    pop_sub(classical_particles_copy_quantities_to_interaction)
  end subroutine classical_particles_copy_quantities_to_interaction
 
  ! ---------------------------------------------------------
  subroutine classical_particles_update_interactions_start(this)
    class(classical_particles_t), intent(inout) :: this
 
    push_sub(classical_particles_update_interactions_start)
 
    ! Store previous force, as it is used as SCF criterium
    select type (prop => this%algo)
    class is (propagator_t)
      if (prop%predictor_corrector) then
        this%prop_data%prev_tot_force(1:this%space%dim, 1:this%np) = this%tot_force(1:this%space%dim, 1:this%np)
      end if
    end select
 
    pop_sub(classical_particles_update_interactions_start)
  end subroutine classical_particles_update_interactions_start
 
  ! ---------------------------------------------------------
  subroutine classical_particles_update_interactions_finish(this)
    class(classical_particles_t), intent(inout) :: this
 
    type(interaction_iterator_t) :: iter
 
    push_sub(classical_particles_update_interactions_finish)
 
    ! Compute the total force acting on the classical particles
    this%tot_force(1:this%space%dim, 1:this%np) = m_zero
    call iter%start(this%interactions)
    do while (iter%has_next())
      select type (interaction => iter%get_next())
      class is (force_interaction_t)
        this%tot_force(1:this%space%dim, 1:this%np) = this%tot_force(1:this%space%dim, 1:this%np) + &
          interaction%force(1:this%space%dim, 1:this%np)
      end select
    end do
 
    pop_sub(classical_particles_update_interactions_finish)
  end subroutine classical_particles_update_interactions_finish
 
  ! ---------------------------------------------------------
  subroutine classical_particles_output_start(this)
    class(classical_particles_t), intent(inout) :: this
 
    integer :: iteration
    real(real64) :: dt
 
    push_sub(classical_particles_output_start)
 
    select type (algo => this%algo)
    class is (propagator_t)
      dt = algo%dt
    end select
 
    iteration = this%iteration%counter()
    if (iteration > 0) then
      ! add one if restarting because the output routine is only called at the end of the timestep
      iteration = iteration + 1
    end if
    ! Create output handle
    call io_mkdir('td.general', this%namespace)
    if (mpi_grp_is_root(this%grp)) then
      call write_iter_init(this%output_handle(output_coordinates), iteration, dt, &
        trim(io_workpath("td.general/coordinates", this%namespace)))
      call write_iter_init(this%output_handle(output_energy), iteration, dt, &
        trim(io_workpath("td.general/energy", this%namespace)))
    end if
 
    ! Output info for first iteration
    if (iteration == 0) then
      call this%output_write()
    end if
 
    pop_sub(classical_particles_output_start)
  end subroutine classical_particles_output_start
 
  ! ---------------------------------------------------------
  subroutine classical_particles_output_finish(this)
    class(classical_particles_t), intent(inout) :: this
 
    push_sub(classical_particles_output_finish)
 
    if (mpi_grp_is_root(this%grp)) then
      call write_iter_end(this%output_handle(output_coordinates))
      call write_iter_end(this%output_handle(output_energy))
    end if
 
    pop_sub(classical_particles_output_finish)
  end subroutine classical_particles_output_finish
 
  ! ---------------------------------------------------------
  subroutine classical_particles_output_write(this)
    class(classical_particles_t), intent(inout) :: this
 
    integer :: idir, ii, ip
    character(len=50) :: aux
    real(real64) :: tmp(this%space%dim)
 
    if (.not. mpi_grp_is_root(this%grp)) return ! only first node outputs
 
    push_sub(classical_particles_output_write)
 
    if (this%iteration%counter() == 0) then
      ! header
      call write_iter_clear(this%output_handle(output_coordinates))
      call write_iter_string(this%output_handle(output_coordinates), &
        '############################################################################')
      call write_iter_nl(this%output_handle(output_coordinates))
      call write_iter_string(this%output_handle(output_coordinates),'# HEADER')
      call write_iter_nl(this%output_handle(output_coordinates))
 
      ! first line: column names
      call write_iter_header_start(this%output_handle(output_coordinates))
      do ip = 1, this%np
        do idir = 1, this%space%dim
          write(aux, '(a2,i3,a1,i3,a1)') 'x(', ip, ',', idir, ')'
          call write_iter_header(this%output_handle(output_coordinates), aux)
        end do
      end do
      do ip = 1, this%np
        do idir = 1, this%space%dim
          write(aux, '(a2,i3,a1,i3,a1)') 'v(', ip, ',', idir, ')'
          call write_iter_header(this%output_handle(output_coordinates), aux)
        end do
      end do
      do ip = 1, this%np
        do idir = 1, this%space%dim
          write(aux, '(a2,i3,a1,i3,a1)') 'f(', ip, ',', idir, ')'
          call write_iter_header(this%output_handle(output_coordinates), aux)
        end do
      end do
      call write_iter_nl(this%output_handle(output_coordinates))
 
      ! second line: units
      call write_iter_string(this%output_handle(output_coordinates), '#[Iter n.]')
      call write_iter_header(this%output_handle(output_coordinates), '[' // trim(units_abbrev(units_out%time)) // ']')
      do ip = 1, this%np
        do idir = 1, this%space%dim
          call write_iter_header(this%output_handle(output_coordinates), '[' // trim(units_abbrev(units_out%length)) // ']')
        end do
      end do
      do ip = 1, this%np
        do idir = 1, this%space%dim
          call write_iter_header(this%output_handle(output_coordinates), '[' // trim(units_abbrev(units_out%velocity)) // ']')
        end do
      end do
      do ip = 1, this%np
        do idir = 1, this%space%dim
          call write_iter_header(this%output_handle(output_coordinates), '[' // trim(units_abbrev(units_out%force)) // ']')
        end do
      end do
      call write_iter_nl(this%output_handle(output_coordinates))
 
      call write_iter_string(this%output_handle(output_coordinates), &
        '############################################################################')
      call write_iter_nl(this%output_handle(output_coordinates))
    end if
 
    call write_iter_start(this%output_handle(output_coordinates))
 
    ! Position
    do ip = 1, this%np
      tmp(1:this%space%dim) = units_from_atomic(units_out%length, this%pos(1:this%space%dim, ip))
      call write_iter_double(this%output_handle(output_coordinates), tmp, this%space%dim)
    end do
    ! Velocity
    do ip = 1, this%np
      tmp(1:this%space%dim) = units_from_atomic(units_out%velocity, this%vel(1:this%space%dim, ip))
      call write_iter_double(this%output_handle(output_coordinates), tmp, this%space%dim)
    end do
    ! Force
    do ip = 1, this%np
      tmp(1:this%space%dim) = units_from_atomic(units_out%force, this%tot_force(1:this%space%dim, ip))
      call write_iter_double(this%output_handle(output_coordinates), tmp, this%space%dim)
    end do
 
    call write_iter_nl(this%output_handle(output_coordinates))
 
    ! Energies
    if (this%iteration%counter() == 0) then
      ! header
      call write_iter_clear(this%output_handle(output_energy))
      call write_iter_string(this%output_handle(output_energy), &
        '############################################################################')
      call write_iter_nl(this%output_handle(output_energy))
      call write_iter_string(this%output_handle(output_energy),'# HEADER')
      call write_iter_nl(this%output_handle(output_energy))
 
      ! first line: column names
      call write_iter_header_start(this%output_handle(output_energy))
      call write_iter_header(this%output_handle(output_energy), 'Total')
      call write_iter_header(this%output_handle(output_energy), 'Kinetic')
      call write_iter_header(this%output_handle(output_energy), 'Potential')
      call write_iter_header(this%output_handle(output_energy), 'Internal')
      call write_iter_nl(this%output_handle(output_energy))
 
      ! second line: units
      call write_iter_string(this%output_handle(output_energy), '#[Iter n.]')
      call write_iter_header(this%output_handle(output_energy), '[' // trim(units_abbrev(units_out%time)) // ']')
      do ii = 1, 4
        call write_iter_header(this%output_handle(output_energy), '[' // trim(units_abbrev(units_out%energy)) // ']')
      end do
      call write_iter_nl(this%output_handle(output_energy))
 
      call write_iter_string(this%output_handle(output_energy), &
        '############################################################################')
      call write_iter_nl(this%output_handle(output_energy))
    end if
    call write_iter_start(this%output_handle(output_energy))
 
    call write_iter_double(this%output_handle(output_energy), units_from_atomic(units_out%energy, this%total_energy), 1)
    call write_iter_double(this%output_handle(output_energy), units_from_atomic(units_out%energy, this%kinetic_energy), 1)
    call write_iter_double(this%output_handle(output_energy), units_from_atomic(units_out%energy, this%potential_energy), 1)
    call write_iter_double(this%output_handle(output_energy), units_from_atomic(units_out%energy, this%internal_energy), 1)
    call write_iter_nl(this%output_handle(output_energy))
 
    pop_sub(classical_particles_output_write)
  end subroutine classical_particles_output_write
 
  ! ---------------------------------------------------------
  subroutine classical_particles_restart_write_data(this)
    class(classical_particles_t), intent(inout) :: this
 
    integer :: restart_file_unit
 
    push_sub(classical_particles_restart_write_data)
 
    if (mpi_grp_is_root(this%grp)) then
      call write_iter_flush(this%output_handle(output_coordinates))
      call write_iter_flush(this%output_handle(output_energy))
    end if
 
    call profiling_in(trim(this%namespace%get())//":"//"RESTART_WRITE")
 
    call io_mkdir('restart/'//td_dir, this%namespace, parents=.true.)
    restart_file_unit = io_open('restart/'//td_dir// 'restart_classical_particles', this%namespace, action='write')
    write(restart_file_unit, *) this%np
    write(restart_file_unit, *) this%mass(:)
    write(restart_file_unit, *) this%pos(:,:)
    write(restart_file_unit, *) this%vel(:,:)
    write(restart_file_unit, *) this%tot_force(:,:)
    call io_close(restart_file_unit)
 
    if (this%iteration%counter() > 0) then
      ! only initialized after the first time step
      call this%prop_data%restart_write(this%namespace, this%algo)
    end if
 
    message(1) = "Successfully wrote restart data for system "//trim(this%namespace%get())
    call messages_info(1, namespace=this%namespace)
 
    call profiling_out(trim(this%namespace%get())//":"//"RESTART_WRITE")
    pop_sub(classical_particles_restart_write_data)
  end subroutine classical_particles_restart_write_data
 
  ! ---------------------------------------------------------
  logical function classical_particles_restart_read_data(this)
    class(classical_particles_t), intent(inout) :: this
 
    integer :: restart_file_unit
 
    push_sub(classical_particles_restart_read_data)
    call profiling_in(trim(this%namespace%get())//":"//"RESTART_READ")
 
    call io_mkdir('restart/'//td_dir, this%namespace, parents=.true.)
    restart_file_unit = io_open('restart/'//td_dir// 'restart_classical_particles', this%namespace, action='read', die=.false.)
    if (restart_file_unit /= -1) then
      read(restart_file_unit, *) this%np
      read(restart_file_unit, *) this%mass(:)
      read(restart_file_unit, *) this%pos(:,:)
      read(restart_file_unit, *) this%vel(:,:)
      read(restart_file_unit, *) this%tot_force(:,:)
      call io_close(restart_file_unit)
      call this%prop_data%initialize(this%algo, this%space%dim, this%np)
      if (this%iteration%counter() > 0) then
        ! only initialized after the first time step
        classical_particles_restart_read_data = this%prop_data%restart_read(this%namespace, this%algo)
      else
        classical_particles_restart_read_data = .true.
      end if
    else
      ! could not open file
      classical_particles_restart_read_data = .false.
    end if
 
    if (classical_particles_restart_read_data) then
      message(1) = "Successfully read restart data for system "//trim(this%namespace%get())
      ! namespace should be added here at some point
      call messages_info(1)
    else
      message(1) = "Failed to read restart data for system "//trim(this%namespace%get())
      ! namespace should be added here at some point
      call messages_info(1)
    end if
 
    call profiling_out(trim(this%namespace%get())//":"//"RESTART_READ")
    pop_sub(classical_particles_restart_read_data)
  end function classical_particles_restart_read_data
 
  ! ---------------------------------------------------------
  subroutine classical_particles_update_kinetic_energy(this)
    class(classical_particles_t),     intent(inout) :: this
 
    integer :: ip
 
    push_sub(classical_particles_update_kinetic_energy)
 
    this%kinetic_energy = m_zero
    do ip = 1, this%np
      this%kinetic_energy = this%kinetic_energy + m_half * this%mass(ip) * sum(this%vel(:, ip)**2)
    end do
 
    pop_sub(classical_particles_update_kinetic_energy)
  end subroutine classical_particles_update_kinetic_energy
 
  ! ---------------------------------------------------------
  function classical_particles_center_of_mass(this, mask, pseudo) result(pos)
    class(classical_particles_t), intent(in) :: this
    logical,            optional, intent(in) :: mask(:)
    logical,            optional, intent(in) :: pseudo
    real(real64) :: pos(this%space%dim)
 
    real(real64) :: mass, total_mass
    integer :: ip
 
    push_sub(classical_particles_center_of_mass)
 
    assert(.not. this%space%is_periodic())
 
    pos = m_zero
    total_mass = m_zero
    mass = m_one
    do ip = 1, this%np
      if (present(mask)) then
        if (.not. mask(ip)) cycle
      end if
      if (.not. optional_default(pseudo, .false.)) then
        mass = this%mass(ip)
      end if
      pos = pos + mass*this%pos(:, ip)
      total_mass = total_mass + mass
    end do
    pos = pos/total_mass
 
    pop_sub(classical_particles_center_of_mass)
  end function classical_particles_center_of_mass
 
  ! ---------------------------------------------------------
  function classical_particles_center_of_mass_vel(this) result(vel)
    class(classical_particles_t), intent(in) :: this
    real(real64) :: vel(this%space%dim)
 
    real(real64) :: mass, total_mass
    integer :: ip
 
    push_sub(classical_particles_center_of_mass_vel)
 
    vel = m_zero
    total_mass = m_zero
    do ip = 1, this%np
      mass = this%mass(ip)
      total_mass = total_mass + mass
      vel = vel + mass*this%vel(:, ip)
    end do
    vel = vel/total_mass
 
    pop_sub(classical_particles_center_of_mass_vel)
  end function classical_particles_center_of_mass_vel
 
  ! ---------------------------------------------------------
  function classical_particles_center(this) result(pos)
    class(classical_particles_t), intent(in) :: this
    real(real64) :: pos(this%space%dim)
 
    real(real64) :: xmin(this%space%dim), xmax(this%space%dim)
    integer  :: ip, idir
 
    push_sub(classical_particles_center)
 
    xmin =  m_huge
    xmax = -m_huge
    do ip = 1, this%np
      do idir = 1, this%space%dim
        if (this%pos(idir, ip) > xmax(idir)) xmax(idir) = this%pos(idir, ip)
        if (this%pos(idir, ip) < xmin(idir)) xmin(idir) = this%pos(idir, ip)
      end do
    end do
 
    pos = (xmax + xmin)/m_two
 
    pop_sub(classical_particles_center)
  end function classical_particles_center
 
  ! ---------------------------------------------------------
  subroutine classical_particles_axis_large(this, x, x2)
    class(classical_particles_t), intent(in)  :: this
    real(real64),                 intent(out) :: x(this%space%dim)
    real(real64),                 intent(out) :: x2(this%space%dim)
 
    integer  :: ip, jp
    real(real64) :: rmax, r, r2
 
    push_sub(classical_particles_axis_large)
 
    assert(.not. this%space%is_periodic())
 
    ! first get the further apart atoms
    rmax = -m_huge
    do ip = 1, this%np
      do jp = 1, this%np/2 + 1
        r = norm2(this%pos(:, ip) - this%pos(:, jp))
        if (r > rmax) then
          rmax = r
          x = this%pos(:, ip) - this%pos(:, jp)
        end if
      end do
    end do
    x  = x /norm2(x)
 
    ! now let us find out what is the second most important axis
    rmax = -m_huge
    do ip = 1, this%np
      r2 = sum(x * this%pos(:, ip))
      r = norm2(this%pos(:, ip) - r2*x)
      if (r > rmax) then
        rmax = r
        x2 = this%pos(:, ip) - r2*x
      end if
    end do
 
    pop_sub(classical_particles_axis_large)
  end subroutine classical_particles_axis_large
 
  ! ---------------------------------------------------------
  subroutine classical_particles_axis_inertia(this, x, x2, pseudo)
    class(classical_particles_t), intent(in)  :: this
    real(real64),                 intent(out) :: x(this%space%dim)
    real(real64),                 intent(out) :: x2(this%space%dim)
    logical,                      intent(in)  :: pseudo
 
    real(real64) :: mass, tinertia(this%space%dim, this%space%dim), eigenvalues(this%space%dim)
    integer :: ii, jj, ip
    type(unit_t) :: unit
 
    push_sub(classical_particles_axis_inertia)
 
    assert(.not. this%space%is_periodic())
 
    ! first calculate the inertia tensor
    tinertia = m_zero
    mass = m_one
    do ip = 1, this%np
      if (.not. pseudo) mass = this%mass(ip)
      do ii = 1, this%space%dim
        tinertia(ii, :) = tinertia(ii, :) - mass*this%pos(ii, ip)*this%pos(:, ip)
        tinertia(ii, ii) = tinertia(ii, ii) + mass*sum(this%pos(:, ip)**2)
      end do
    end do
 
    unit = units_out%length**2
    ! note: we always use amu for atomic masses, so no unit conversion to/from atomic is needed.
    if (pseudo) then
      write(message(1),'(a)') 'Moment of pseudo-inertia tensor [' // trim(units_abbrev(unit)) // ']'
    else
      write(message(1),'(a)') 'Moment of inertia tensor [amu*' // trim(units_abbrev(unit)) // ']'
    end if
    call messages_info(1, namespace=this%namespace)
    call output_tensor(tinertia, this%space%dim, unit, write_average = .true., namespace=this%namespace)
 
    call lalg_eigensolve(this%space%dim, tinertia, eigenvalues)
 
    write(message(1),'(a,6f25.6)') 'Eigenvalues: ', units_from_atomic(unit, eigenvalues)
    call messages_info(1, namespace=this%namespace)
 
    ! make a choice to fix the sign of the axis.
    do ii = 1, 2
      jj = maxloc(abs(tinertia(:,ii)), dim = 1)
      if (tinertia(jj,ii) < m_zero) tinertia(:,ii) = -tinertia(:,ii)
    end do
    x  = tinertia(:,1)
    x2 = tinertia(:,2)
 
    pop_sub(classical_particles_axis_inertia)
  end subroutine classical_particles_axis_inertia
 
  ! ---------------------------------------------------------
  subroutine classical_particles_end(this)
    class(classical_particles_t), intent(inout) :: this
 
    push_sub(classical_particles_end)
 
    safe_deallocate_a(this%mass)
    safe_deallocate_a(this%pos)
    safe_deallocate_a(this%vel)
    safe_deallocate_a(this%tot_force)
    safe_deallocate_a(this%fixed)
    safe_deallocate_a(this%lj_epsilon)
    safe_deallocate_a(this%lj_sigma)
 
    call system_end(this)
 
    pop_sub(classical_particles_end)
  end subroutine classical_particles_end
 
end module classical_particles_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: