module classes_box_volume_change use, intrinsic :: iso_fortran_env, only: DP => REAL64 use data_constants, only: num_dimensions use types_logical_wrapper, only: Logical_Line use procedures_logical_factory, only: logical_create => create use procedures_random_number, only: random_integer use classes_tower_sampler, only: Abstract_Tower_Sampler use procedures_tower_sampler_factory, only: tower_sampler_destroy => destroy use types_environment_wrapper, only: Environment_Wrapper use types_mixture_wrapper, only: Mixture_Wrapper use procedures_mixture_factory, only: mixture_rescale_positions => rescale_positions use types_cells_wrapper, only: Cells_Wrapper use procedures_cells_factory, only: cells_destroy => destroy use procedures_cells_memento, only: cells_memento_save => save, cells_memento_restore => restore use types_short_interactions_wrapper, only: Short_Interactions_Wrapper use procedures_short_interactions_resetter, only: short_interactions_reset => reset use procedures_short_interactions_visitor, only: short_interactions_visit => visit use types_dipolar_interactions_static_wrapper, only: Dipolar_Interactions_Static_Wrapper use procedures_dipolar_interactions_factory, only: dipolar_interactions_destroy => destroy use procedures_dipolar_interactions_visitor, only: dipolar_interactions_visit => visit use classes_dipolar_interactions_facade, only: Abstract_Dipolar_Interactions_Facade use classes_changed_box_size, only: Abstract_Changed_Box_Size use procedures_triangle_observables, only: triangle_observables_diff, triangle_observables_sum use types_observables_energies, only: Concrete_Observables_Energies use procedures_observables_energies_factory, only: observables_energies_create => create, & observables_energies_set => set use types_generating_observables_wrapper, only: Generating_Observables_Wrapper use classes_generating_algorithm, only: Abstract_Generating_Algorithm use procedures_selectors_resetters, only: selectors_reset => reset use procedures_metropolis_algorithm, only: metropolis_algorithm implicit none private type, extends(Abstract_Generating_Algorithm), public :: Box_Volume_Change private type(Environment_Wrapper), pointer :: environment => null() type(Mixture_Wrapper), pointer :: mixture => null() type(Short_Interactions_Wrapper), pointer :: short_interactions => null() class(Abstract_Dipolar_Interactions_Facade), pointer :: dipolar_interactions_facades(:) => & null() class(Abstract_Changed_Box_Size), pointer :: changed_boxes_size(:) => null() logical, allocatable :: have_positions(:, :) class(Abstract_Tower_Sampler), allocatable :: selectors(:) contains procedure :: construct => Concrete_construct procedure :: destroy => Concrete_destroy procedure :: reset_selectors => Concrete_reset_selectors procedure :: get_num_choices => Concrete_get_num_choices procedure :: try => Concrete_try procedure, private :: metropolis_algorithm => Concrete_metropolis_algorithm procedure, private :: acceptation_probability => Concrete_acceptation_probability end type Box_Volume_Change contains subroutine Concrete_construct(this, environment, mixture, short_interactions, & dipolar_interactions_facades, changed_boxes_size, have_positions, selectors) class(Box_Volume_Change), intent(out) :: this type(Environment_Wrapper), target, intent(in) :: environment type(Mixture_Wrapper), target, intent(in) :: mixture type(Short_Interactions_Wrapper), target, intent(in) :: short_interactions class(Abstract_Dipolar_Interactions_Facade), target, intent(in) :: & dipolar_interactions_facades(:) class(Abstract_Changed_Box_Size), target, intent(in) :: changed_boxes_size(:) logical, intent(in) :: have_positions(:, :) class(Abstract_Tower_Sampler), intent(in) :: selectors(:) this%environment => environment this%mixture => mixture this%short_interactions => short_interactions this%dipolar_interactions_facades => dipolar_interactions_facades this%changed_boxes_size => changed_boxes_size allocate(this%have_positions, source=have_positions) allocate(this%selectors, source=selectors) end subroutine Concrete_construct subroutine Concrete_destroy(this) class(Box_Volume_Change), intent(inout) :: this call tower_sampler_destroy(this%selectors) if (allocated(this%have_positions)) deallocate(this%have_positions) this%changed_boxes_size => null() this%dipolar_interactions_facades => null() this%short_interactions => null() this%mixture => null() this%environment => null() end subroutine Concrete_destroy subroutine Concrete_reset_selectors(this) class(Box_Volume_Change), intent(inout) :: this call selectors_reset(this%selectors, this%changed_boxes_size, this%mixture%& average_nums_particles, this%have_positions) end subroutine Concrete_reset_selectors pure integer function Concrete_get_num_choices(this) result(num_choices) class(Box_Volume_Change), intent(in) :: this integer :: i_box num_choices = 0 do i_box = 1, size(this%selectors) num_choices = num_choices + this%selectors(i_box)%get_num_choices() end do end function Concrete_get_num_choices !> @note [[procedures_dipolar_interactions_factory:destroy_static]] is not necessary !> since dipolar_interactions_static is declared inside the subroutine. !> @todo Send a bug report to ifort. subroutine Concrete_try(this, observables) class(Box_Volume_Change), intent(in) :: this type(Generating_Observables_Wrapper), intent(inout) :: observables logical :: success type(Concrete_Observables_Energies) :: new_energies real(DP), dimension(num_dimensions) :: new_box_size, box_size, box_size_ratio integer :: i_box type(Logical_Line), allocatable :: only_resized_triangle(:) type(Cells_Wrapper) :: cells type(Dipolar_Interactions_Static_Wrapper) :: dipolar_interactions_static logical :: reset_real_pair i_box = random_integer(size(this%environment%periodic_boxes)) observables%volumes_change_counter(i_box)%num_hits = observables%& volumes_change_counter(i_box)%num_hits + 1 box_size = this%environment%periodic_boxes(i_box)%get_size() box_size_ratio = this%changed_boxes_size(i_box)%get_ratio() new_box_size = box_size * box_size_ratio call this%dipolar_interactions_facades(i_box)%save(dipolar_interactions_static, & reset_real_pair, new_box_size) call this%environment%periodic_boxes(i_box)%set(new_box_size) call mixture_rescale_positions(this%mixture%components(:, i_box), box_size_ratio) call logical_create(only_resized_triangle, size(this%mixture%components, 1)) call cells_memento_save(cells, only_resized_triangle, this%short_interactions%& visitable_cells_memento, this%short_interactions%cells(i_box)) call short_interactions_reset(this%short_interactions%cells(i_box)%neighbour_cells, & only_resized_triangle, this%short_interactions%cells(i_box)%visitable_cells) call this%dipolar_interactions_facades(i_box)%reset(reset_real_pair) call observables_energies_create(new_energies, size(this%mixture%components, 1)) call this%metropolis_algorithm(success, new_energies, i_box, product(box_size_ratio), & observables%energies(i_box)) if (success) then observables%accessible_domains_size(:, i_box) = this%environment%& accessible_domains(i_box)%get_size() call observables_energies_set(observables%energies(i_box), new_energies) observables%volumes_change_counter(i_box)%num_successes = observables%& volumes_change_counter(i_box)%num_successes + 1 call cells_destroy(cells%visitable_cells) call cells_destroy(cells%neighbour_cells) else call this%environment%periodic_boxes(i_box)%set(box_size) call mixture_rescale_positions(this%mixture%components(:, i_box), 1._DP/box_size_ratio) call cells_memento_restore(this%short_interactions%cells(i_box), only_resized_triangle,& this%short_interactions%visitable_cells_memento, cells) call this%dipolar_interactions_facades(i_box)%restore(dipolar_interactions_static, & reset_real_pair) end if call dipolar_interactions_destroy(dipolar_interactions_static) end subroutine Concrete_try subroutine Concrete_metropolis_algorithm(this, success, new_energies, i_box, box_volume_ratio, & energies) class(Box_Volume_Change), intent(in) :: this logical, intent(out) :: success type(Concrete_Observables_Energies), intent(inout) :: new_energies integer, intent(in) :: i_box real(DP), intent(in) :: box_volume_ratio type(Concrete_Observables_Energies), intent(in) :: energies real(DP) :: delta_energy type(Concrete_Observables_Energies) :: deltas logical :: overlap success = .false. call observables_energies_create(deltas, size(this%mixture%components, 1)) call short_interactions_visit(overlap, new_energies%walls_energies, this%& mixture%components(:, i_box), this%short_interactions%walls_visitors(i_box), this%& short_interactions%wall_pairs) if (overlap) return deltas%walls_energies = new_energies%walls_energies - energies%walls_energies call short_interactions_visit(overlap, new_energies%short_energies, this%& mixture%components(:, i_box), this%short_interactions%cells(i_box)%visitable_cells) if (overlap) return call triangle_observables_diff(deltas%short_energies, new_energies%short_energies, & energies%short_energies) call dipolar_interactions_visit(new_energies%field_energies, this%environment%& external_fields(i_box), this%mixture%components(:, i_box)) deltas%field_energies = new_energies%field_energies - energies%field_energies call this%dipolar_interactions_facades(i_box)%visit(new_energies%dipolar_energies, & new_energies%dipolar_shared_energy, box_volume_ratio, energies%dipolar_energies,& energies%dipolar_shared_energy) call triangle_observables_diff(deltas%dipolar_energies, new_energies%dipolar_energies, & energies%dipolar_energies) deltas%dipolar_shared_energy = new_energies%dipolar_shared_energy - energies%& dipolar_shared_energy delta_energy = sum(deltas%walls_energies + deltas%field_energies) + & triangle_observables_sum(deltas%short_energies) + & triangle_observables_sum(deltas%dipolar_energies) + deltas%dipolar_shared_energy success = metropolis_algorithm(this%acceptation_probability(i_box, box_volume_ratio, & delta_energy)) end subroutine Concrete_metropolis_algorithm !> \[ !> P[V \to V^\prime] = \min \left( 1, e^{-\beta p V \left( \frac{V^\prime}{V} - 1 \right)} !> \left( \frac{V^\prime}{V} \right)^{N+1} !> e^{-\beta [U(\vec{s}^N, V^\prime) - U(\vec{s}^N, V)]} \right) !> \] pure real(DP) function Concrete_acceptation_probability(this, i_box, box_volume_ratio, & delta_energy) result(probability) class(Box_Volume_Change), intent(in) :: this integer, intent(in) :: i_box real(DP), intent(in) :: box_volume_ratio real(DP), intent(in) :: delta_energy integer :: i_component, num_particles num_particles = 0._DP do i_component = 1, size(this%mixture%components, 1) num_particles = num_particles + this%mixture%components(i_component, i_box)%& num_particles%get() end do probability = exp(-this%environment%beta_pressure%get() * & product(this%environment%accessible_domains(i_box)%get_size()) * & (box_volume_ratio - 1._DP)) * box_volume_ratio**(num_particles + 1) * & exp(-delta_energy / this%environment%temperature%get()) probability = min(1._DP, probability) end function Concrete_acceptation_probability end module classes_box_volume_change