lennard_jones.F90

Go to the documentation of this file.

!! Copyright (C) 2023 F. Bonafe
!!
!! 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 lennard_jones_oct_m
  use debug_oct_m
  use force_interaction_oct_m
  use global_oct_m
  use interaction_oct_m
  use interaction_partner_oct_m
  use, intrinsic :: iso_fortran_env
  use messages_oct_m
  use namespace_oct_m
  use profiling_oct_m
 
  implicit none
 
  private
  public ::                &
    lennard_jones_t
 
  type, extends(force_interaction_t) :: lennard_jones_t
    private
    real(real64), pointer :: system_pos(:,:) => null() 
    real(real64), public  :: lj_epsilon
    real(real64), public  :: lj_sigma
 
    integer, public :: partner_np = 0
    real(real64), allocatable, public :: partner_pos(:,:)
 
  contains
    procedure :: init => lennard_jones_init
    procedure :: calculate => lennard_jones_calculate
    procedure :: calculate_energy => lennard_jones_calculate_energy
    final :: lennard_jones_finalize
  end type lennard_jones_t
 
  interface lennard_jones_t
    module procedure lennard_jones_constructor
  end interface lennard_jones_t
 
contains
 
  ! ---------------------------------------------------------
  function lennard_jones_constructor(partner) result(this)
    class(interaction_partner_t), target, intent(inout) :: partner
    class(lennard_jones_t),               pointer       :: this
 
    push_sub(lennard_jones_constructor)
 
    allocate(this)
 
    this%label = "lennard_jones"
 
    this%partner => partner
 
    ! The Lennard Jones interaction needs the position of the particles and the parameters
    this%system_quantities = ["position"]
 
    ! The Lennard Jones interaction needs the position of the partner
    this%couplings_from_partner = ["position"]
 
    pop_sub(lennard_jones_constructor)
  end function lennard_jones_constructor
 
  ! ---------------------------------------------------------
  subroutine lennard_jones_init(this, dim, system_np, system_pos, system_eps, system_sigma)
    class(lennard_jones_t),             intent(inout) :: this
    integer,                            intent(in)    :: dim
    integer,                            intent(in)    :: system_np
    real(real64),               target, intent(in)    :: system_pos(:,:)
    real(real64),                       intent(in)    :: system_eps
    real(real64),                       intent(in)    :: system_sigma
 
    push_sub(lennard_jones_init)
 
    this%dim = dim
    this%system_np = system_np
    safe_allocate(this%force(1:dim, 1:system_np))
    this%force = m_zero
 
    this%lj_epsilon = system_eps
    this%lj_sigma = system_sigma
 
    this%system_pos => system_pos
 
    pop_sub(lennard_jones_init)
  end subroutine lennard_jones_init
 
  ! ---------------------------------------------------------
  subroutine lennard_jones_calculate(this)
    class(lennard_jones_t),             intent(inout) :: this
 
    integer :: ip, jp
    real(real64) :: dist, rr(1:this%dim), lj_force
 
    push_sub(lennard_jones_calculate)
 
    assert(allocated(this%partner_pos))
 
    do ip = 1, this%system_np
      do jp = 1, this%partner_np
        if (this%intra_interaction .and. ip == jp ) cycle
 
        ! r_ij = r_i - r_j
        rr(1:this%dim) = this%system_pos(1:this%dim, ip) - this%partner_pos(1:this%dim, jp)
        dist = sqrt(sum(rr(1:this%dim)**2))
 
        ! lj_force = -d U(|r_ij|) / d |r_ij|
        lj_force = 48.0_real64 * this%lj_epsilon * (this%lj_sigma**12 / dist**13 - &
          m_half * this%lj_sigma**6 / dist**7)
 
        ! F_i = - d U(r) / d r_i = - (d U (r) / d |r_ij|) * (d r_ij / d r_i) =
        ! = - d U (r) / d r_ij, because d r_ij / d r_i = 1
        this%force(1:this%dim, ip) = rr(1:this%dim) / dist * lj_force
      end do
    end do
 
    pop_sub(lennard_jones_calculate)
  end subroutine lennard_jones_calculate
 
  ! ---------------------------------------------------------
  subroutine lennard_jones_calculate_energy(this)
    class(lennard_jones_t),             intent(inout) :: this
 
    integer :: ip, jp
    real(real64) :: dist
 
    push_sub(lennard_jones_calculate_energy)
 
    assert(allocated(this%partner_pos))
 
    this%energy = m_zero
    do ip = 1, this%system_np
      do jp = 1, this%partner_np
        if (this%intra_interaction .and. ip == jp ) cycle
 
        dist = sqrt(sum((this%system_pos(1:this%dim, ip) - this%partner_pos(1:this%dim, jp))**2)) + m_epsilon
 
        ! this%energy = this%energy + 0.5 * U(r_ij)
        this%energy = this%energy + m_two * this%lj_epsilon * ( (this%lj_sigma / dist)**12 - &
          (this%lj_sigma / dist)**6 )
      end do
    end do
 
    pop_sub(lennard_jones_calculate_energy)
  end subroutine lennard_jones_calculate_energy
 
 
  ! ---------------------------------------------------------
  subroutine lennard_jones_finalize(this)
    type(lennard_jones_t), intent(inout) :: this
 
    push_sub(lennard_jones_finalize)
 
    this%force = m_zero
    nullify(this%system_pos)
    safe_deallocate_a(this%partner_pos)
    safe_deallocate_a(this%force)
 
    call interaction_end(this)
 
    pop_sub(lennard_jones_finalize)
  end subroutine lennard_jones_finalize
 
end module lennard_jones_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: