classical_particle.F90

Go to the documentation of this file.

!! Copyright (C) 2019 N. Tancogne-Dejean
!! Copyright (C) 2020 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_particle_oct_m
  use algorithm_oct_m
  use classical_particles_oct_m
  use debug_oct_m
  use distribution_gaussian_oct_m
  use global_oct_m
  use interaction_oct_m
  use interaction_surrogate_oct_m
  use interaction_enum_oct_m
  use gravity_oct_m
  use io_oct_m
  use lennard_jones_oct_m
  use messages_oct_m
  use namespace_oct_m
  use parser_oct_m
  use profiling_oct_m
  use propagator_oct_m
  use quantity_oct_m
  use quickrnd_oct_m
  use space_oct_m
  use system_oct_m
  use unit_oct_m
  use unit_system_oct_m
 
  implicit none
 
  private
  public ::                  &
    classical_particle_t,    &
    classical_particle_init
 
  integer, parameter ::     &
    OUTPUT_COORDINATES = 1, &
    output_energy      = 2
 
  type, extends(classical_particles_t) :: classical_particle_t
  contains
    procedure :: init_interaction => classical_particle_init_interaction
    procedure :: initialize => classical_particle_initialize
    procedure :: update_quantity => classical_particle_update_quantity
    procedure :: init_interaction_as_partner => classical_particle_init_interaction_as_partner
    procedure :: copy_quantities_to_interaction => classical_particle_copy_quantities_to_interaction
    final :: classical_particle_finalize
  end type classical_particle_t
 
  interface classical_particle_t
    procedure classical_particle_constructor
  end interface classical_particle_t
 
contains
 
  ! ---------------------------------------------------------
  function classical_particle_constructor(namespace) result(sys)
    class(classical_particle_t), pointer    :: sys
    type(namespace_t),           intent(in) :: namespace
 
    push_sub(classical_particle_constructor)
 
    allocate(sys)
 
    call classical_particle_init(sys, namespace)
 
    pop_sub(classical_particle_constructor)
  end function classical_particle_constructor
 
  ! ---------------------------------------------------------
  ! ---------------------------------------------------------
  subroutine classical_particle_init(this, namespace)
    class(classical_particle_t), intent(inout) :: this
    type(namespace_t),           intent(in)    :: namespace
 
    push_sub(classical_particle_init)
 
    this%namespace = namespace
 
    this%space = space_t(namespace)
    if (this%space%is_periodic()) then
      call messages_not_implemented('Classical particle for periodic systems')
    end if
 
    call messages_print_with_emphasis(msg="Classical Particle", namespace=namespace)
 
    call classical_particles_init(this, 1)
 
    !%Variable ParticleMass
    !%Type float
    !%Section ClassicalParticles
    !%Description
    !% Mass of classical particle in Kg.
    !%End
    call parse_variable(namespace, 'ParticleMass', m_one, this%mass(1))
    call messages_print_var_value('ParticleMass', this%mass(1), namespace=namespace)
 
    this%supported_interactions = [this%supported_interactions, gravity, lennard_jones]
    this%supported_interactions_as_partner = [this%supported_interactions_as_partner, gravity, lennard_jones]
 
    !%Variable LennardJonesEpsilon
    !%Type float
    !%Section ClassicalParticles
    !%Description
    !% Epsilon parameter (dispersion energy) of Lennard-Jones interaction for this species.
    !% In case two particles have a different epsilon, the combination rule will be computed
    !% <math>\epsilon_{12} = \sqrt{\epsilon_1 + \epsilon_2}</math>.
    !%End
    call parse_variable(namespace, 'LennardJonesEpsilon', m_one, this%lj_epsilon(1))
    call messages_print_var_value('LennardJonesEpsilon', this%lj_epsilon(1), namespace=namespace)
 
    !%Variable LennardJonesSigma
    !%Type float
    !%Section ClassicalParticles
    !%Description
    !% Sigma parameter (particle size) of Lennard-Jones interaction for this species.
    !% In case two particles have a different sigma, the combination rule will be computed
    !% <math>\sigma_{12} = (\sigma_1 + \sigma_2) / 2</math>.
    !%End
    call parse_variable(namespace, 'LennardJonesSigma', m_one, this%lj_sigma(1))
    call messages_print_var_value('LennardJonesSigma', this%lj_sigma(1), namespace=namespace)
 
    call messages_print_with_emphasis(namespace=namespace)
 
    pop_sub(classical_particle_init)
  end subroutine classical_particle_init
 
  ! ---------------------------------------------------------
  subroutine classical_particle_init_interaction(this, interaction)
    class(classical_particle_t), target, intent(inout) :: this
    class(interaction_t),                intent(inout) :: interaction
 
    push_sub(classical_particle_init_interaction)
 
    select type (interaction)
    type is (gravity_t)
      call interaction%init(this%space%dim, 1, this%mass, this%pos)
    type is (lennard_jones_t)
      if (.not. (parse_is_defined(this%namespace, 'LennardJonesSigma') .and. &
        parse_is_defined(this%namespace, 'LennardJonesEpsilon') )) then
        write(message(1),'(a,es9.2)') 'Using default value for Lennard-Jones parameter.'
        call messages_warning(1, namespace=this%namespace)
      end if
 
      call interaction%init(this%space%dim, 1, this%pos, this%lj_epsilon(1), this%lj_sigma(1))
    class default
      call classical_particles_init_interaction(this, interaction)
    end select
 
    pop_sub(classical_particle_init_interaction)
  end subroutine classical_particle_init_interaction
 
  ! ---------------------------------------------------------
  subroutine classical_particle_initialize(this)
    class(classical_particle_t), intent(inout) :: this
 
    integer :: n_rows, idir
    type(block_t) :: blk
 
    real(real64) :: random_pos(1:this%space%dim), width_pos(1:this%space%dim)
    real(real64) :: random_vel(1:this%space%dim), width_vel(1:this%space%dim)
    integer(int64) :: seed
    logical :: in_ensemble
 
    push_sub(classical_particle_initialize)
 
    in_ensemble = parse_is_defined(this%namespace, "NumberOfReplicas")
 
    !%Variable ParticleInitialPosition
    !%Type block
    !%Section ClassicalParticles
    !%Description
    !% Initial position of classical particle, in Km.
    !%End
    this%pos = m_zero
    if (parse_block(this%namespace, 'ParticleInitialPosition', blk) == 0) then
      n_rows = parse_block_n(blk)
      if (n_rows > 1) call messages_input_error(this%namespace, 'ParticleInitialPosition')
 
      do idir = 1, this%space%dim
        call parse_block_float(blk, 0, idir - 1, this%pos(idir, 1))
      end do
      call parse_block_end(blk)
    end if
    call messages_print_var_value('ParticleInitialPosition', this%pos(1:this%space%dim, 1), namespace=this%namespace)
 
    !%Variable ParticleInitialVelocity
    !%Type block
    !%Section ClassicalParticles
    !%Description
    !% Initial velocity of classical particle in Km/s.
    !%End
    this%vel = m_zero
    if (parse_block(this%namespace, 'ParticleInitialVelocity', blk) == 0) then
      n_rows = parse_block_n(blk)
      if (n_rows > 1) call messages_input_error(this%namespace, 'ParticleInitialVelocity')
      do idir = 1, this%space%dim
        call parse_block_float(blk, 0, idir - 1, this%vel(idir, 1))
      end do
      call parse_block_end(blk)
    end if
    call messages_print_var_value('ParticleInitialVelocity', this%vel(1:this%space%dim, 1), namespace=this%namespace)
 
    if (in_ensemble) then
 
      seed = this%namespace%get_hash32()
      width_pos = m_zero
      width_vel = m_zero
 
      !%Variable ParticlePositionDistributionWidth
      !%Type block
      !%Section ClassicalParticles
      !%Description
      !% Width of the position random displacements of classical particle in Km.
      !%End
      if (parse_block(this%namespace, 'ParticlePositionDistributionWidth', blk) == 0) then
        n_rows = parse_block_n(blk)
        if (n_rows > 1) call messages_input_error(this%namespace, 'ParticlePositionDistributionWidth')
        do idir = 1, this%space%dim
          call parse_block_float(blk, 0, idir - 1, width_pos(idir))
        end do
        call parse_block_end(blk)
      end if
      call messages_print_var_value('ParticlePositionDistributionWidth', width_pos(1:this%space%dim), namespace=this%namespace)
 
      !%Variable ParticleVelocityDistributionWidth
      !%Type block
      !%Section ClassicalParticles
      !%Description
      !% Width of the velocity random distribution of classical particle in Km/s.
      !%End
      if (parse_block(this%namespace, 'ParticleVelocityDistributionWidth', blk) == 0) then
        n_rows = parse_block_n(blk)
        if (n_rows > 1) call messages_input_error(this%namespace, 'ParticleVelocityDistributionWidth')
        do idir = 1, this%space%dim
          call parse_block_float(blk, 0, idir - 1, width_vel(idir))
        end do
        call parse_block_end(blk)
      end if
      call messages_print_var_value('ParticleVelocityDistributionWidth', width_vel(1:this%space%dim), namespace=this%namespace)
 
      call shiftseed(seed, 1000_int64)
      call quickrnd(seed, this%space%dim, random_pos)
      call quickrnd(seed, this%space%dim, random_vel)
 
 
      this%pos(:,1) = this%pos(:,1) + (random_pos - 0.5_real64) * width_pos
      this%vel(:,1) = this%vel(:,1) + (random_vel - 0.5_real64) * width_vel
 
    end if
 
    pop_sub(classical_particle_initialize)
  end subroutine classical_particle_initialize
 
  ! ---------------------------------------------------------
  subroutine classical_particle_update_quantity(this, label)
    class(classical_particle_t), intent(inout) :: this
    character(len=*),            intent(in)    :: label
 
    push_sub(classical_particle_update_quantity)
 
    select case (label)
    case default
      ! Other quantities should be handled by the parent class
      call classical_particles_update_quantity(this, label)
    end select
 
    pop_sub(classical_particle_update_quantity)
  end subroutine classical_particle_update_quantity
 
  ! ---------------------------------------------------------
  subroutine classical_particle_init_interaction_as_partner(partner, interaction)
    class(classical_particle_t),    intent(in)    :: partner
    class(interaction_surrogate_t), intent(inout) :: interaction
 
    push_sub(classical_particle_init_interaction_as_partner)
 
    select type (interaction)
    type is (gravity_t)
      interaction%partner_np = 1
      safe_allocate(interaction%partner_mass(1))
      safe_allocate(interaction%partner_pos(1:partner%space%dim, 1))
      interaction%partner_pos = m_zero
 
    type is (lennard_jones_t)
      interaction%partner_np = 1
      safe_allocate(interaction%partner_pos(1:partner%space%dim, 1))
      ! in case the LennardJones epsilon and sigma of system and partner are different,
      ! we compute the combination rules with arithmetic and geometric means:
      ! (they give back the original parameters if they happen to be equal):
      ! sigma_12 = (sigma_1 + sigma_2)/2,    epsilon_12 = sqrt(epsilon_1 * epsilon_2)
      interaction%lj_sigma = m_half * (partner%lj_sigma(1) + interaction%lj_sigma)
      interaction%lj_epsilon = sqrt(partner%lj_epsilon(1) * interaction%lj_epsilon)
 
    class default
      ! Other interactions should be handled by the parent class
      call classical_particles_init_interaction_as_partner(partner, interaction)
    end select
 
    pop_sub(classical_particle_init_interaction_as_partner)
  end subroutine classical_particle_init_interaction_as_partner
 
  ! ---------------------------------------------------------
  subroutine classical_particle_copy_quantities_to_interaction(partner, interaction)
    class(classical_particle_t),          intent(inout) :: partner
    class(interaction_surrogate_t),       intent(inout) :: interaction
 
    push_sub(classical_particle_copy_quantities_to_interaction)
 
    select type (interaction)
    type is (gravity_t)
      interaction%partner_mass(1) = partner%mass(1)
      interaction%partner_pos(:,1) = partner%pos(:,1)
 
    type is (lennard_jones_t)
      interaction%partner_pos(:,1) = partner%pos(:,1)
 
    class default
      message(1) = "Unsupported interaction."
      call messages_fatal(1, namespace=partner%namespace)
    end select
 
    pop_sub(classical_particle_copy_quantities_to_interaction)
  end subroutine classical_particle_copy_quantities_to_interaction
 
  ! ---------------------------------------------------------
  subroutine classical_particle_finalize(this)
    type(classical_particle_t), intent(inout) :: this
 
    push_sub(classical_particle_finalize)
 
    call classical_particles_end(this)
 
    pop_sub(classical_particle_finalize)
  end subroutine classical_particle_finalize
 
end module classical_particle_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: