module classes_dlc_structures

use, intrinsic :: iso_fortran_env, only: DP => REAL64
use data_constants, only: num_dimensions, PI
use classes_periodic_box, only: Abstract_Periodic_Box
use classes_reciprocal_lattice, only: Abstract_Reciprocal_Lattice
use types_component_wrapper, only: Component_Wrapper
use types_particle_wrapper, only: Concrete_Particle
use classes_structure_factor, only: Abstract_Structure_Factor
use procedures_dipolar_interactions_micro, only: set_fourier, set_exp_kz, reci_number_1_sym

implicit none

private

    type, extends(Abstract_Structure_Factor), abstract, public :: Abstract_DLC_Structures
    private
        class(Abstract_Periodic_Box), pointer :: periodic_box => null()
        integer :: reci_numbers(2) = 0
        type(Component_Wrapper), pointer :: components(:) => null()
        logical, allocatable :: are_dipolar(:)
        complex(DP), dimension(:, :), allocatable :: structure_p, structure_m
    contains
        procedure :: construct => Abstract_construct
        procedure :: destroy => Abstract_destroy
        procedure :: target => Abstract_target
        procedure :: reset => Abstract_reset
        procedure :: is_dipolar => Abstract_is_dipolar
        procedure :: get_plus => Abstract_get_plus
        procedure :: get_minus => Abstract_get_minus
        procedure :: update_translation => Abstract_update_translation
        procedure :: update_transmutation => Abstract_update_transmutation
        procedure :: update_rotation => Abstract_update_rotation
        procedure :: update_add => Abstract_update_add
        procedure :: update_remove => Abstract_update_remove
        procedure :: update_switch => Abstract_update_switch
        procedure, private :: update_exchange => Abstract_update_exchange
    end type Abstract_DLC_Structures

    type, extends(Abstract_DLC_Structures), public :: Concrete_DLC_Structures

    end type Concrete_DLC_Structures

    type, extends(Abstract_DLC_Structures), public :: Null_DLC_Structures
    contains
        procedure :: construct => Null_construct
        procedure :: destroy => Null_destroy
        procedure :: target => Null_target
        procedure :: reset => Null_reset
        procedure :: is_dipolar => Null_is_dipolar
        procedure :: get_plus => Null_get
        procedure :: get_minus => Null_get
        procedure :: update_translation => Null_update_translation
        procedure :: update_transmutation => Null_update_transmutation
        procedure :: update_switch => Null_update_switch
        procedure, private :: update_exchange => Null_update_exchange
    end type Null_DLC_Structures

contains

!implementation Abstract_DLC_Structures

    subroutine Abstract_construct(this, periodic_box, reciprocal_lattice, components, are_dipolar)
        class(Abstract_DLC_Structures), intent(out) :: this
        class(Abstract_Periodic_Box), intent(in) :: periodic_box
        class(Abstract_Reciprocal_Lattice), intent(in) :: reciprocal_lattice
        type(Component_Wrapper), intent(in) :: components(:)
        logical, intent(in) :: are_dipolar(:)

        call this%target(periodic_box, components)
        this%reci_numbers = reshape(reciprocal_lattice%get_numbers(), [2])
        allocate(this%are_dipolar(size(are_dipolar)))
        this%are_dipolar = are_dipolar
        allocate(this%structure_p(-this%reci_numbers(1):this%reci_numbers(1), &
                                                      0:this%reci_numbers(2)))
        allocate(this%structure_m(-this%reci_numbers(1):this%reci_numbers(1), &
                                                      0:this%reci_numbers(2)))
    end subroutine Abstract_construct

    subroutine Abstract_destroy(this)
        class(Abstract_DLC_Structures), intent(inout) :: this

        if (allocated(this%structure_m)) deallocate(this%structure_m)
        if (allocated(this%structure_p)) deallocate(this%structure_p)
        if (allocated(this%are_dipolar)) deallocate(this%are_dipolar)
        this%components => null()
        this%periodic_box => null()
    end subroutine Abstract_destroy

    subroutine Abstract_target(this, periodic_box, components)
        class(Abstract_DLC_Structures), intent(inout) :: this
        class(Abstract_Periodic_Box), target, intent(in) :: periodic_box
        type(Component_Wrapper), target, intent(in) :: components(:)

        this%periodic_box => periodic_box
        this%components => components
    end subroutine Abstract_target

    subroutine Abstract_reset(this)
        class(Abstract_DLC_Structures), intent(inout) :: this

        type(Concrete_Particle) :: particle
        integer :: i_component, i_particle

        this%structure_p = cmplx(0._DP, 0._DP, DP)
        this%structure_m = cmplx(0._DP, 0._DP, DP)
        do i_component = 1, size(this%components)
            do i_particle = 1, this%components(i_component)%dipole_moments%get_num()
                particle%position = this%components(i_component)%positions%get(i_particle)
                particle%dipole_moment = this%components(i_component)%dipole_moments%&
                    get(i_particle)
                call this%update_add(i_component, particle)
            end do
        end do
    end subroutine Abstract_reset

    pure logical function Abstract_is_dipolar(this, i_component) result(is_dipolar)
        class(Abstract_DLC_Structures), intent(in) :: this
        integer, intent(in) :: i_component

        is_dipolar = this%are_dipolar(i_component)
    end function Abstract_is_dipolar

    !> Structure factors
    !> \[
    !>      S_\pm(\vec{k}_{1:2}) = \sum_{\vec{x}, \vec{\mu}}
    !>          (\pm k_{1:2} \mu_3 + i \vec{k}_{1:2} \cdot \vec{\mu}_{1:2})
    !>          e^{\pm k_{1:2} x_3} e^{i \vec{k}_{1:2} \cdot \vec{x}_{1:2}}
    !> \]
    pure complex(DP) function Abstract_get_plus(this, n_1, n_2) result(structure_p)
        class(Abstract_DLC_Structures), intent(in) :: this
        integer, intent(in) :: n_1, n_2

        structure_p = this%structure_p(n_1, n_2)
    end function Abstract_get_plus

    !> Cf. [[Abstract_get_plus]]
    pure complex(DP) function Abstract_get_minus(this, n_1, n_2) result(structure_m)
        class(Abstract_DLC_Structures), intent(in) :: this
        integer, intent(in) :: n_1, n_2

        structure_m = this%structure_m(n_1, n_2)
    end function Abstract_get_minus

    !> Structure factors update when a particle is translated:
    !> \( (\vec{x}, \vec{\mu}) \to (\vec{x}^\prime, \vec{\mu}) \).
    !> \[
    !>      \Delta S_\pm(\vec{k}_{1:2}) =
    !>          (\pm k_{1:2} \mu_3 + i \vec{k}_{1:2} \cdot \vec{\mu}_{1:2})
    !>          \left(
    !>              e^{\pm k_{1:2} x^\prime_3} e^{i \vec{k}_{1:2} \cdot \vec{x}^\prime_{1:2}} -
    !>              e^{\pm k_{1:2} x_3} e^{i \vec{k}_{1:2} \cdot \vec{x}_{1:2}}
    !>          \right)
    !> \]
    !> Warning: only half wave vectors are updated.
    pure subroutine Abstract_update_translation(this, i_component, new_position, old)
        class(Abstract_DLC_Structures), intent(inout) :: this
        integer, intent(in) :: i_component
        real(DP), intent(in) :: new_position(:)
        type(Concrete_Particle), intent(in) :: old

        real(DP) :: surface_size(2)
        real(DP), dimension(2) :: wave_1_x_position_old, wave_1_x_position_new, wave_vector
        real(DP) :: wave_dot_moment_12, wave_x_moment_3
        integer :: n_1, n_2

        complex(DP) :: fourier_position_new
        complex(DP), dimension(-this%reci_numbers(1):this%reci_numbers(1)) :: fourier_position_new_1
        complex(DP), dimension(-this%reci_numbers(2):this%reci_numbers(2)) :: fourier_position_new_2
        real(DP) :: exp_kz_new
        real(DP), dimension(0:this%reci_numbers(1), 0:this%reci_numbers(2)) :: exp_kz_new_tab
        complex(DP) :: fourier_position_old
        complex(DP), dimension(-this%reci_numbers(1):this%reci_numbers(1)) :: fourier_position_old_1
        complex(DP), dimension(-this%reci_numbers(2):this%reci_numbers(2)) :: fourier_position_old_2
        real(DP) :: exp_kz_old
        real(DP), dimension(0:this%reci_numbers(1), 0:this%reci_numbers(2)) :: exp_kz_old_tab

        if (.not.this%are_dipolar(i_component)) return

        surface_size = reshape(this%periodic_box%get_size(), [2])
        wave_1_x_position_old = 2._DP*PI * old%position(1:2) / surface_size
        call set_fourier(fourier_position_old_1, this%reci_numbers(1), wave_1_x_position_old(1))
        call set_fourier(fourier_position_old_2, this%reci_numbers(2), wave_1_x_position_old(2))
        call set_exp_kz(exp_kz_old_tab, surface_size, old%position(3))
        wave_1_x_position_new = 2._DP*PI * new_position(1:2) / surface_size
        call set_fourier(fourier_position_new_1, this%reci_numbers(1), wave_1_x_position_new(1))
        call set_fourier(fourier_position_new_2, this%reci_numbers(2), wave_1_x_position_new(2))
        call set_exp_kz(exp_kz_new_tab, surface_size, new_position(3))

        do n_2 = 0, this%reci_numbers(2)
            wave_vector(2) = 2._DP*PI * real(n_2, DP) / surface_size(2)
            do n_1 = -reci_number_1_sym(this%reci_numbers, 0, n_2), this%reci_numbers(1)
                wave_vector(1) = 2._DP*PI * real(n_1, DP) / surface_size(1)

                if (n_1**2 + n_2**2 > this%reci_numbers(1)**2) cycle

                fourier_position_old = fourier_position_old_1(n_1) * fourier_position_old_2(n_2)
                exp_kz_old = exp_kz_old_tab(abs(n_1), abs(n_2))
                fourier_position_new = fourier_position_new_1(n_1) * fourier_position_new_2(n_2)
                exp_kz_new = exp_kz_new_tab(abs(n_1), abs(n_2))

                wave_dot_moment_12 = dot_product(wave_vector, old%dipole_moment(1:2))
                wave_x_moment_3 = norm2(wave_vector) * old%dipole_moment(3)

                this%structure_p(n_1, n_2) = this%structure_p(n_1, n_2) + &
                    cmplx(+wave_x_moment_3, wave_dot_moment_12, DP) * &
                    (fourier_position_new * exp_kz_new - fourier_position_old * exp_kz_old)
                this%structure_m(n_1, n_2) = this%structure_m(n_1, n_2) + &
                    cmplx(-wave_x_moment_3, wave_dot_moment_12, DP) * &
                    (fourier_position_new / exp_kz_new - fourier_position_old / exp_kz_old)
            end do
        end do
    end subroutine Abstract_update_translation

    !> Structure factors update when a particle is transmuted:
    !> \( (\vec{x}, \vec{\mu}) \to (\vec{x}, \vec{\mu}^\prime) \).
    !>  \[
    !>      \Delta S_\pm(\vec{k}_{1:2}) =
    !>          \left[
    !>              \pm k_{1:2} (\mu^\prime_3 - \mu_3) +
    !>              i \vec{k}_{1:2} \cdot (\vec{\mu}^\prime_{1:2} - \vec{\mu}_{1:2})
    !>          \right]
    !>          e^{\pm k_{1:2} x_3} e^{i \vec{k}_{1:2} \cdot \vec{x}_{1:2}}
    !>  \]
    !> Warning: only half wave vectors are updated.
    pure subroutine Abstract_update_transmutation(this, ij_components, new_dipole_moment, old)
        class(Abstract_DLC_Structures), intent(inout) :: this
        integer, intent(in) :: ij_components(:)
        real(DP), intent(in) :: new_dipole_moment(:)
        type(Concrete_Particle), intent(in) :: old

        real(DP) :: surface_size(2)
        real(DP), dimension(2) :: wave_1_x_position, wave_vector
        real(DP) :: wave_dot_delta_moment_12, wave_delta_moment_3
        integer :: n_1, n_2

        complex(DP) :: fourier_position
        complex(DP), dimension(-this%reci_numbers(1):this%reci_numbers(1)) :: fourier_position_1
        complex(DP), dimension(-this%reci_numbers(2):this%reci_numbers(2)) :: fourier_position_2
        real(DP) :: exp_kz
        real(DP), dimension(0:this%reci_numbers(1), 0:this%reci_numbers(2)) :: exp_kz_tab

        if (.not.(this%are_dipolar(ij_components(1)) .or. this%are_dipolar(ij_components(2)))) &
            return

        surface_size = reshape(this%periodic_box%get_size(), [2])
        wave_1_x_position = 2._DP*PI * old%position(1:2) / surface_size
        call set_fourier(fourier_position_1, this%reci_numbers(1), wave_1_x_position(1))
        call set_fourier(fourier_position_2, this%reci_numbers(2), wave_1_x_position(2))
        call set_exp_kz(exp_kz_tab, surface_size, old%position(3))

        do n_2 = 0, this%reci_numbers(2)
            wave_vector(2) = 2._DP*PI * real(n_2, DP) / surface_size(2)
            do n_1 = -reci_number_1_sym(this%reci_numbers, 0, n_2), this%reci_numbers(1)
                wave_vector(1) = 2._DP*PI * real(n_1, DP) / surface_size(1)

                if (n_1**2 + n_2**2 > this%reci_numbers(1)**2) cycle

                fourier_position = fourier_position_1(n_1) * fourier_position_2(n_2)
                exp_kz = exp_kz_tab(abs(n_1), abs(n_2))

                wave_dot_delta_moment_12 = dot_product(wave_vector, new_dipole_moment(1:2) - old%&
                    dipole_moment(1:2))
                wave_delta_moment_3 = norm2(wave_vector) * (new_dipole_moment(3) - old%&
                    dipole_moment(3))

                this%structure_p(n_1, n_2) = this%structure_p(n_1, n_2) + &
                    cmplx(+wave_delta_moment_3, wave_dot_delta_moment_12, DP) * &
                    (fourier_position * exp_kz)
                this%structure_m(n_1, n_2) = this%structure_m(n_1, n_2) + &
                    cmplx(-wave_delta_moment_3, wave_dot_delta_moment_12, DP) * &
                    (fourier_position / exp_kz)
            end do
        end do
    end subroutine Abstract_update_transmutation

    pure subroutine Abstract_update_rotation(this, i_component, new_dipole_moment, old)
        class(Abstract_DLC_Structures), intent(inout) :: this
        integer, intent(in) :: i_component
        real(DP), intent(in) :: new_dipole_moment(:)
        type(Concrete_Particle), intent(in) :: old

        call this%update_transmutation([i_component, i_component], new_dipole_moment, old)
    end subroutine Abstract_update_rotation

    !> cf. [[classes_dlc_structures:Abstract_update_exchange]]
    pure subroutine Abstract_update_add(this, i_component, particle)
        class(Abstract_DLC_Structures), intent(inout) :: this
        integer, intent(in) :: i_component
        type(Concrete_Particle), intent(in) :: particle

        call this%update_exchange(i_component, particle, +1._DP)
    end subroutine Abstract_update_add

    !> cf. [[classes_dlc_structures:Abstract_update_exchange]]
    pure subroutine Abstract_update_remove(this, i_component, particle)
        class(Abstract_DLC_Structures), intent(inout) :: this
        integer, intent(in) :: i_component
        type(Concrete_Particle), intent(in) :: particle

        call this%update_exchange(i_component, particle, -1._DP)
    end subroutine Abstract_update_remove

    !> Structure factors update when a particle of coordinates \( (\vec{x}, \vec{\mu}) \) is added
    !> (\( + )\) or removed (\( - \)):
    !> \[
    !>      \Delta S_\pm(\vec{k}_{1:2}) =
    !>          {\bf\pm} (\pm k_{1:2} \mu_3 + i \vec{k}_{1:2} \cdot \vec{\mu}_{1:2})
    !>          e^{\pm k_{1:2} x_3} e^{i \vec{k}_{1:2} \cdot \vec{x}_{1:2}}.
    !> \]
    pure subroutine Abstract_update_exchange(this, i_component, particle, signed)
        class(Abstract_DLC_Structures), intent(inout) :: this
        integer, intent(in) :: i_component
        type(Concrete_Particle), intent(in) :: particle
        real(DP), intent(in) :: signed

        real(DP) :: surface_size(2)
        real(DP), dimension(2) :: wave_1_x_position, wave_vector
        real(DP) :: wave_dot_moment_12, wave_x_moment_3
        integer :: n_1, n_2

        complex(DP) :: fourier_position
        complex(DP), dimension(-this%reci_numbers(1):this%reci_numbers(1)) :: fourier_position_1
        complex(DP), dimension(-this%reci_numbers(2):this%reci_numbers(2)) :: fourier_position_2
        real(DP) :: exp_kz
        real(DP), dimension(0:this%reci_numbers(1), 0:this%reci_numbers(2)) :: exp_kz_tab

        if (.not.this%are_dipolar(i_component)) return

        surface_size = reshape(this%periodic_box%get_size(), [2])
        wave_1_x_position = 2._DP*PI * particle%position(1:2) / surface_size
        call set_fourier(fourier_position_1, this%reci_numbers(1), wave_1_x_position(1))
        call set_fourier(fourier_position_2, this%reci_numbers(2), wave_1_x_position(2))
        call set_exp_kz(exp_kz_tab, surface_size, particle%position(3))

        do n_2 = 0, this%reci_numbers(2)
            wave_vector(2) = 2._DP*PI * real(n_2, DP) / surface_size(2)
            do n_1 = -reci_number_1_sym(this%reci_numbers, 0, n_2), this%reci_numbers(1)
                wave_vector(1) = 2._DP*PI * real(n_1, DP) / surface_size(1)

                if (n_1**2 + n_2**2 > this%reci_numbers(1)**2) cycle

                fourier_position = fourier_position_1(n_1) * fourier_position_2(n_2)
                exp_kz = exp_kz_tab(abs(n_1), abs(n_2))
                wave_dot_moment_12 = dot_product(wave_vector, particle%dipole_moment(1:2))
                wave_x_moment_3 = norm2(wave_vector) * particle%dipole_moment(3)

                this%structure_p(n_1, n_2) = this%structure_p(n_1, n_2) + &
                    signed * cmplx(+wave_x_moment_3, wave_dot_moment_12, DP) * &
                    fourier_position * exp_kz
                this%structure_m(n_1, n_2) = this%structure_m(n_1, n_2) + &
                    signed * cmplx(-wave_x_moment_3, wave_dot_moment_12, DP) * &
                    fourier_position / exp_kz
            end do
        end do
    end subroutine Abstract_update_exchange

    !> Structure factors update when 2 particles of coordinates \( (\vec{x}_1, \vec{\mu}_1) \) and
    !> \( (\vec{x}_2, \vec{\mu}_2) \) are switched.
    !> \[
    !>      \Delta S_\pm(\vec{k}_{1:2}) =
    !>          \left[
    !>              \pm k_{1:2} (\mu_{1, 3} - \mu_{2, 3}) +
    !>              i \vec{k}_{1:2} \cdot (\vec{\mu}_{1, 1:2} - \vec{\mu}_{2, 1:2})
    !>          \right] \\
    !>          \left(
    !>              e^{\pm k_{1:2} x_{2, 3}} e^{i \vec{k}_{1:2} \cdot \vec{x}_{2, 1:2}} -
    !>              e^{\pm k_{1:2} x_{1, 3}} e^{i \vec{k}_{1:2} \cdot \vec{x}_{1, 1:2}}
    !>          \right)
    !> \]
    pure subroutine Abstract_update_switch(this, ij_components, particles)
        class(Abstract_DLC_Structures), intent(inout) :: this
        integer, intent(in) :: ij_components(:)
        type(Concrete_Particle), intent(in) :: particles(:)

        real(DP) :: surface_size(2)
        real(DP), dimension(2) :: wave_1_x_position_1, wave_1_x_position_2, wave_vector
        real(DP) :: wave_dot_delta_moment_12, wave_delta_moment_3
        integer :: n_1, n_2

        complex(DP) :: fourier_position_1
        complex(DP), dimension(-this%reci_numbers(1):this%reci_numbers(1)) :: fourier_position_1_1
        complex(DP), dimension(-this%reci_numbers(2):this%reci_numbers(2)) :: fourier_position_1_2
        real(DP) :: exp_kz_1
        real(DP), dimension(0:this%reci_numbers(1), 0:this%reci_numbers(2)) :: exp_kz_1_tab
        complex(DP) :: fourier_position_2
        complex(DP), dimension(-this%reci_numbers(1):this%reci_numbers(1)) :: fourier_position_2_1
        complex(DP), dimension(-this%reci_numbers(2):this%reci_numbers(2)) :: fourier_position_2_2
        real(DP) :: exp_kz_2
        real(DP), dimension(0:this%reci_numbers(1), 0:this%reci_numbers(2)) :: exp_kz_2_tab

        if (.not.(this%are_dipolar(ij_components(1)) .or. this%are_dipolar(ij_components(2)))) &
            return

        surface_size = reshape(this%periodic_box%get_size(), [2])
        wave_1_x_position_1 = 2._DP*PI * particles(1)%position(1:2) / surface_size
        call set_fourier(fourier_position_1_1, this%reci_numbers(1), wave_1_x_position_1(1))
        call set_fourier(fourier_position_1_2, this%reci_numbers(2), wave_1_x_position_1(2))
        call set_exp_kz(exp_kz_1_tab, surface_size, particles(1)%position(3))
        wave_1_x_position_2 = 2._DP*PI * particles(2)%position(1:2) / surface_size
        call set_fourier(fourier_position_2_1, this%reci_numbers(1), wave_1_x_position_2(1))
        call set_fourier(fourier_position_2_2, this%reci_numbers(2), wave_1_x_position_2(2))
        call set_exp_kz(exp_kz_2_tab, surface_size, particles(2)%position(3))

        do n_2 = 0, this%reci_numbers(2)
            wave_vector(2) = 2._DP*PI * real(n_2, DP) / surface_size(2)
            do n_1 = -reci_number_1_sym(this%reci_numbers, 0, n_2), this%reci_numbers(1)
                wave_vector(1) = 2._DP*PI * real(n_1, DP) / surface_size(1)

                if (n_1**2 + n_2**2 > this%reci_numbers(1)**2) cycle

                fourier_position_1 = fourier_position_1_1(n_1) * fourier_position_1_2(n_2)
                exp_kz_1 = exp_kz_1_tab(abs(n_1), abs(n_2))
                fourier_position_2 = fourier_position_2_1(n_1) * fourier_position_2_2(n_2)
                exp_kz_2 = exp_kz_2_tab(abs(n_1), abs(n_2))

                wave_dot_delta_moment_12 = dot_product(wave_vector, &
                    particles(1)%dipole_moment(1:2) - particles(2)%dipole_moment(1:2))
                wave_delta_moment_3 = norm2(wave_vector) * &
                    (particles(1)%dipole_moment(3) - particles(2)%dipole_moment(3))

                this%structure_p(n_1, n_2) = this%structure_p(n_1, n_2) + &
                    cmplx(+wave_delta_moment_3, wave_dot_delta_moment_12, DP) * &
                    (fourier_position_2 * exp_kz_2 - fourier_position_1 * exp_kz_1)
                this%structure_m(n_1, n_2) = this%structure_m(n_1, n_2) + &
                    cmplx(-wave_delta_moment_3, wave_dot_delta_moment_12, DP) * &
                    (fourier_position_2 / exp_kz_2 - fourier_position_1 / exp_kz_1)
            end do
        end do
    end subroutine Abstract_update_switch

!end implementation Abstract_DLC_Structures

!implementation Null_DLC_Structures

    subroutine Null_construct(this, periodic_box, reciprocal_lattice, components, are_dipolar)
        class(Null_DLC_Structures), intent(out) :: this
        class(Abstract_Periodic_Box), intent(in) :: periodic_box
        class(Abstract_Reciprocal_Lattice), intent(in) :: reciprocal_lattice
        type(Component_Wrapper), intent(in) :: components(:)
        logical, intent(in) :: are_dipolar(:)
    end subroutine Null_construct

    subroutine Null_destroy(this)
        class(Null_DLC_Structures), intent(inout) :: this
    end subroutine Null_destroy

    subroutine Null_target(this, periodic_box, components)
        class(Null_DLC_Structures), intent(inout) :: this
        class(Abstract_Periodic_Box), target, intent(in) :: periodic_box
        type(Component_Wrapper), target, intent(in) :: components(:)
    end subroutine Null_target

    subroutine Null_reset(this)
        class(Null_DLC_Structures), intent(inout) :: this
    end subroutine Null_reset

    pure logical function Null_is_dipolar(this, i_component) result(is_dipolar)
        class(Null_DLC_Structures), intent(in) :: this
        integer, intent(in) :: i_component
        is_dipolar = .false.
    end function Null_is_dipolar

    pure complex(DP) function Null_get(this, n_1, n_2) result(structure_pm)
        class(Null_DLC_Structures), intent(in) :: this
        integer, intent(in) :: n_1, n_2
        structure_pm = cmplx(0._DP, 0._DP, DP)
    end function Null_get

    pure subroutine Null_update_translation(this, i_component, new_position, old)
        class(Null_DLC_Structures), intent(inout) :: this
        integer, intent(in) :: i_component
        real(DP), intent(in) :: new_position(:)
        type(Concrete_Particle), intent(in) :: old
    end subroutine Null_update_translation

    pure subroutine Null_update_transmutation(this, ij_components, new_dipole_moment, old)
        class(Null_DLC_Structures), intent(inout) :: this
        integer, intent(in) :: ij_components(:)
        real(DP), intent(in) :: new_dipole_moment(:)
        type(Concrete_Particle), intent(in) :: old
    end subroutine Null_update_transmutation

    pure subroutine Null_update_exchange(this, i_component, particle, signed)
        class(Null_DLC_Structures), intent(inout) :: this
        integer, intent(in) :: i_component
        type(Concrete_Particle), intent(in) :: particle
        real(DP), intent(in) :: signed
    end subroutine Null_update_exchange

    pure subroutine Null_update_switch(this, ij_components, particles)
        class(Null_DLC_Structures), intent(inout) :: this
        integer, intent(in) :: ij_components(:)
        type(Concrete_Particle), intent(in) :: particles(:)
    end subroutine Null_update_switch

!end implementation Null_DLC_Structures

end module classes_dlc_structures