Factories

  recursive function system_factory_create(this, namespace, name, type) result(system)
    class(system_factory_t), intent(in) :: this       !< the system factory
    type(namespace_t),       intent(in) :: namespace  !< namespace of the parent
    character(len=*),        intent(in) :: name       !< name of the system (will define namespace)
    integer,                 intent(in) :: type       !< type of the system to create
    class(system_t),         pointer    :: system     !< pointer to newly created system

    PUSH_SUB(system_factory_create)

    !%Variable Systems
    !%Type block
    !%Section System
    !%Description
    !% List of systems that will be treated in the calculation.
    !% The first column should be a string containing the system name.
    !% The second column should be the system type. See below for a list of
    !% available system types.
    !%Option electronic 1
    !% An electronic system. (not fully implemented yet)
    !%Option maxwell 2
    !% A maxwell system.
    !%Option classical_particle 3
    !% A classical particle. Used for testing purposes only.
    !%Option charged_particle 4
    !% A charged classical particle.
    !%Option dftbplus 5
    !% A DFTB+ system
    !%Option linear_medium 6
    !% A linear medium for classical electrodynamics.
    !%Option matter 7
    !% A matter system containing electrons and classical ions.
    !%Option dispersive_medium 8
    !% (Experimental) A dispersive medium for classical electrodynamics.
    !%Option multisystem 9
    !% A system containing other systems.
    !%End
    select case (type)
    case (SYSTEM_MULTISYSTEM)
      system => multisystem_basic_t(namespace_t(name, parent=namespace), this)
    case (SYSTEM_ELECTRONIC)
      system => electrons_t(namespace_t(name, parent=namespace))
    case (SYSTEM_MAXWELL)
      system => maxwell_t(namespace_t(name, parent=namespace))
    case (SYSTEM_CLASSICAL_PARTICLE)
      system => classical_particle_t(namespace_t(name, parent=namespace))
    case (SYSTEM_CHARGED_PARTICLE)
      system => charged_particle_t(namespace_t(name, parent=namespace))
    case (SYSTEM_DFTBPLUS)
      system => dftb_t(namespace_t(name, parent=namespace))
    case (SYSTEM_LINEAR_MEDIUM)
      system => linear_medium_t(namespace_t(name, parent=namespace))
    case (SYSTEM_MATTER)
      system => matter_t(namespace_t(name, parent=namespace))
    case (SYSTEM_DISPERSIVE_MEDIUM)
      system => dispersive_medium_t(namespace_t(name, parent=namespace))
      call messages_experimental('dispersive_medium', namespace=namespace)
    case default
      system => null()
    end select

    POP_SUB(system_factory_create)
  end function system_factory_create

  recursive subroutine interactions_factory_abst_create_interactions(this, system, available_partners)
    class(interactions_factory_abst_t),    intent(in)    :: this
    class(system_t),                       intent(inout) :: system
    class(partner_list_t),         target, intent(in)    :: available_partners

    type(integer_list_t) :: interactions_to_create
    type(integer_iterator_t) :: interaction_iter
    integer :: interaction_type
    type(partner_list_t) :: partners, partners_flat_list
    type(partner_iterator_t) :: partner_iter
    class(interaction_partner_t), pointer :: partner
    type(system_iterator_t) :: iter
    class(system_t), pointer :: subsystem

    integer :: il, ic, mode
    type(block_t) :: blk
    character(len=MAX_NAMESPACE_LEN) :: input_name

    PUSH_SUB(interactions_factory_abst_create_interactions)

    ! Make a copy of the interactions list so that we can modify it
    interactions_to_create = system%supported_interactions

    ! Get the list of partners as a flat list
    call flatten_partner_list(available_partners, partners_flat_list)

    ! Parse input. The variable name and description should be given by the
    ! factory, as different factories might have different options.
    if (parse_block(system%namespace, this%block_name(), blk) == 0) then

      ! Loop over all interactions specified in the input file
      do il = 0, parse_block_n(blk) - 1
        ! Read the interaction type (first column)
        call parse_block_integer(blk, il, 0, interaction_type)

        ! Sanity check: the interaction type must be known and must not be mistaken for an interaction mode
        if (.not. varinfo_valid_option(this%block_name(), interaction_type) .or. &
          any(interaction_type == (/ALL_PARTNERS, ONLY_PARTNERS, NO_PARTNERS, ALL_EXCEPT/))) then
          call messages_input_error(system%namespace, this%block_name(), details="Unknown interaction type", row=il, column=0)
        end if

        ! Ignore interactions that are not supported by this system
        if (.not. interactions_to_create%has(interaction_type)) cycle

        ! Read how this interaction should be treated (second column)
        call parse_block_integer(blk, il, 1, mode)

        ! Create list of partners for this interaction taking into account the selected mode
        select case (mode)
        case (ALL_PARTNERS)
          ! Use all available partners
          partners = partners_flat_list
        case (NO_PARTNERS)
          ! No partners for this interaction
          call partners%empty()
        case (ONLY_PARTNERS)
          ! Start with an empty list. We will add only the select partners bellow
          call partners%empty()
        case (ALL_EXCEPT)
          ! Start with full list. We will remove the select partners bellow
          partners = partners_flat_list
        case default
          call messages_input_error(system%namespace, this%block_name(), "Unknown interaction mode", row=il, column=1)
        end select

        if (mode == ONLY_PARTNERS .or. mode == ALL_EXCEPT) then
          ! In these two cases we need to read the names of the selected
          ! partners (remaining columns) and handled them appropriately
          do ic = 2, parse_block_cols(blk, il) - 1
            call parse_block_string(blk, il, ic, input_name)

            ! Loop over available partners and either add them or remove them
            ! from the list depending on the selected mode
            call partner_iter%start(partners_flat_list)
            do while (partner_iter%has_next())
              partner => partner_iter%get_next()
              if (partner%namespace%is_contained_in(input_name)) then
                select case (mode)
                case (ONLY_PARTNERS)
                  call partners%add(partner)
                case (ALL_EXCEPT)
                  call partners%delete(partner)
                end select
              end if
            end do
          end do

        end if

        ! Now actually create the interactions for the selected partners
        call create_interactions(this, system, partners, interaction_type)

        ! Remove this interaction type from the list, as it has just been handled
        call interactions_to_create%delete(interaction_type)
      end do
      call parse_block_end(blk)
    end if

    ! Loop over all the remaining interactions supported by the system
    call interaction_iter%start(interactions_to_create)
    do while (interaction_iter%has_next())
      interaction_type = interaction_iter%get_next()

      ! Check what is the default mode for this interaction type (all or none)
      select case (this%default_mode(system%namespace, interaction_type))
      case (ALL_PARTNERS)
        partners = partners_flat_list
      case (NO_PARTNERS)
        call partners%empty()
      case default
        message(1) = "Default interaction mode can only be all_partners or no_partners."
        call messages_fatal(1, namespace=system%namespace)
      end select

      call create_interactions(this, system, partners, interaction_type)
    end do

    ! All systems need to be connected to make sure they remain synchronized.
    ! We enforce that be adding a ghost interaction between all systems
    call create_interactions(this, system, partners_flat_list)

    ! If the system is a multisystem, then we also need to create the interactions for the subsystems
    select type (system)
    class is (multisystem_t)
      call iter%start(system%list)
      do while (iter%has_next())
        subsystem => iter%get_next()
        call this%create_interactions(subsystem, available_partners)
      end do
    end select

    POP_SUB(interactions_factory_abst_create_interactions)
  contains

    recursive subroutine flatten_partner_list(partners, flat_list)
      class(partner_list_t),  intent(in)    :: partners
      class(partner_list_t),  intent(inout) :: flat_list

      class(interaction_partner_t), pointer :: partner
      type(partner_iterator_t) :: iterator

      PUSH_SUB(interactions_factory_abst_create_interactions.flatten_partner_list)

      call iterator%start(partners)
      do while (iterator%has_next())
        partner => iterator%get_next()

        call flat_list%add(partner)

        select type (partner)
        class is (multisystem_t)
          ! Also include the subsystems of a multisystem
          call flatten_partner_list(partner%list, flat_list)
        end select

      end do

      POP_SUB(interactions_factory_abst_create_interactions.flatten_partner_list)
    end subroutine flatten_partner_list

  function interactions_factory_create(this, type, partner) result(interaction)
    class(interactions_factory_t),         intent(in)    :: this
    integer,                               intent(in)    :: type
    class(interaction_partner_t),  target, intent(inout) :: partner
    class(interaction_t),                  pointer       :: interaction

    PUSH_SUB(interactions_factory_create)

    !%Variable Interactions
    !%Type block
    !%Section System
    !%Description
    !% This input option controls the interactions between systems. It basically
    !% allows to select which systems will interact with another system through
    !% a given interaction type. The format of the block is the following:
    !%
    !%  <br>%<tt>Namespace.Interactions
    !%   <br>&nbsp;&nbsp;interaction_type | interaction_mode | ...
    !%  <br>%</tt>
    !%
    !% Here is an example to better understand how this works:
    !%
    !%  <br>%<tt>SystemA.Interactions
    !%   <br>&nbsp;&nbsp;gravity | all_except | "SystemB"
    !%  <br>%</tt>
    !%
    !% This means that SystemA and all the systems that belong to the same
    !% namespace (i.e., all its subsystems) will interact through gravity with
    !% all interaction partners that are also able to interact through gravity,
    !% except with SystemB. Note that the opposite is not true so, although
    !% clearly unphysical, this will not prevent SystemB from feeling the
    !% gravity from SystemA (in <tt>Octopus</tt> the interactions are always
    !% one-sided).
    !%
    !% NB: Each interaction type should only appear once in the block. Any
    !% further instances beyond the first will be ignored.
    !%
    !% Available modes and interaction types:
    !%Option no_partners -1
    !%  (interaction mode)
    !% Do not interact with any partner.
    !%Option all_partners -2
    !%  (interaction mode)
    !% Interact with all available partners.
    !%Option only_partners -3
    !%  (interaction mode)
    !% Interact only with some specified partners. A list of partner names must
    !% be given.
    !%Option all_except -4
    !%  (interaction mode)
    !% Interact with all available partners except with some specified
    !% partners. A list of partner names to exclude must be given.
    !%Option gravity 1
    !%  (interaction type)
    !% Gravity interaction between two masses.
    !%Option lorentz_force 2
    !%  (interaction type)
    !% Lorentz force resulting from an EM field acting on a moving charge.
    !%Option coulomb_force 3
    !%  (interaction type)
    !% Coulomb force between two charged particles.
    !%Option linear_medium_to_em_field 4
    !%  (interaction type)
    !% Linear medium for propagation of EM fields.
    !%Option current_to_mxll_field 5
    !%  (interaction type)
    !% Drude dispersive linear medium for propagation of EM fields.
    !%Option maxwell_e_field 6
    !%  (interaction type)
    !% Electric field resulting from the Maxwell solver.
    !%Option maxwell_b_field 7
    !%  (interaction type)
    !% Magnetic field resulting from the Maxwell solver.
    !%Option maxwell_vector_potential 8
    !%  (interaction type)
    !% Vector potential resulting from the Maxwell solver.
    !%Option lennard_jones 9
    !%  (interaction type)
    !% Force resulting from a Lennard Jones potential between classical particles.
    !%End
    select case (type)
    case (GRAVITY)
      interaction => gravity_t(partner, .false.)
    case (COULOMB_FORCE)
      interaction => coulomb_force_t(partner, .false.)
    case (LORENTZ_FORCE)
      interaction => lorentz_force_t(partner)
    case (LINEAR_MEDIUM_TO_EM_FIELD)
      interaction => linear_medium_to_em_field_t(partner)
    case (CURRENT_TO_MXLL_FIELD)
      interaction => current_to_mxll_field_t(partner)
    case (MXLL_E_FIELD_TO_MATTER)
      interaction => mxll_e_field_to_matter_t(partner)
    case (MXLL_B_FIELD_TO_MATTER)
      interaction => mxll_b_field_to_matter_t(partner)
    case (MXLL_VEC_POT_TO_MATTER)
      interaction => mxll_vec_pot_to_matter_t(partner)
    case (LENNARD_JONES)
      interaction => lennard_jones_t(partner, .false.)
    case default
      ! This should never happen, as this is handled in
      ! interactions_factory_abst_create_interactions
    end select

    POP_SUB(interactions_factory_create)
  end function interactions_factory_create