System: Electrons

  type, extends(system_t) :: electrons_t
    ! Components are public by default
    type(ions_t),     pointer    :: ions => NULL() !< the ion component of the system
    type(grid_t)                 :: gr    !< the mesh
    type(states_elec_t)          :: st    !< the states
    type(v_ks_t)                 :: ks    !< the Kohn-Sham potentials
    type(output_t)               :: outp  !< the output
    type(multicomm_t)            :: mc    !< index and domain communicators
    type(hamiltonian_elec_t)     :: hm    !< the Hamiltonian
    type(td_t)                   :: td    !< everything related to time propagation
    type(current_t)              :: current_calculator
    type(scf_t)                  :: scf   !< SCF for BOMD

    type(kpoints_t) :: kpoints                   !< the k-points

    logical :: generate_epot

    type(states_elec_t)          :: st_copy  !< copy of the states

    ! At the moment this is not treated as an external potential
    class(lasers_t), pointer :: lasers => null()      !< lasers
    class(gauge_field_t), pointer :: gfield => null()      !< gauge field

    ! List with all the external partners
    ! This will become a list of interactions in the future
    type(partner_list_t) :: ext_partners

    !TODO: have a list of self interactions
    type(xc_interaction_t), pointer   :: xc_interaction => null()

    logical :: ions_propagated = .false.

    logical :: needs_current = .false.   !< If current is needed by an interaction
  contains
    procedure :: init_interaction => electrons_init_interaction
    procedure :: init_parallelization => electrons_init_parallelization
    procedure :: init_algorithm => electrons_init_algorithm
    procedure :: initial_conditions => electrons_initial_conditions
    procedure :: do_algorithmic_operation => electrons_do_algorithmic_operation
    procedure :: is_tolerance_reached => electrons_is_tolerance_reached
    procedure :: update_quantity => electrons_update_quantity
    procedure :: update_exposed_quantity => electrons_update_exposed_quantity
    procedure :: init_interaction_as_partner => electrons_init_interaction_as_partner
    procedure :: copy_quantities_to_interaction => electrons_copy_quantities_to_interaction
    procedure :: output_start => electrons_output_start
    procedure :: output_write => electrons_output_write
    procedure :: output_finish => electrons_output_finish
    procedure :: process_is_slave  => electrons_process_is_slave
    procedure :: restart_write_data => electrons_restart_write_data
    procedure :: restart_read_data => electrons_restart_read_data
    procedure :: update_kinetic_energy => electrons_update_kinetic_energy
    procedure :: propagation_start => electrons_propagation_start
    final :: electrons_finalize
  end type electrons_t

  function electrons_constructor(namespace, generate_epot) result(sys)
    class(electrons_t), pointer    :: sys
    type(namespace_t),  intent(in) :: namespace
    logical,  optional, intent(in) :: generate_epot

    integer :: iatom
    type(lattice_vectors_t) :: latt_inp
    type(profile_t), save :: prof

    PUSH_SUB(electrons_constructor)
    call profiling_in(prof,"ELECTRONS_CONSTRUCTOR")

    SAFE_ALLOCATE(sys)

    sys%namespace = namespace

    call messages_obsolete_variable(sys%namespace, 'SystemName')

    call space_init(sys%space, sys%namespace)
    call sys%space%write_info(sys%namespace)
    if (sys%space%has_mixed_periodicity()) then
      call messages_experimental('Support for mixed periodicity systems')
    end if

    sys%ions => ions_t(sys%namespace, latt_inp=latt_inp)

    call grid_init_stage_1(sys%gr, sys%namespace, sys%space, sys%ions%symm, latt_inp, sys%ions%natoms, sys%ions%pos)

    if (sys%space%is_periodic()) then
      call sys%ions%latt%write_info(sys%namespace)
    end if

    ! Sanity check for atomic coordinates
    do iatom = 1, sys%ions%natoms
      if (.not. sys%gr%box%contains_point(sys%ions%pos(:, iatom))) then
        if (sys%space%periodic_dim /= sys%space%dim) then
          ! FIXME: This could fail for partial periodicity systems
          ! because contains_point is too strict with atoms close to
          ! the upper boundary to the cell.
          write(message(1), '(a,i5,a)') "Atom ", iatom, " is outside the box."
          call messages_warning(1, namespace=sys%namespace)
        end if
      end if
    end do

    ! we need k-points for periodic systems
    call kpoints_init(sys%kpoints, sys%namespace, sys%gr%symm, sys%space%dim, sys%space%periodic_dim, sys%ions%latt)

    call states_elec_init(sys%st, sys%namespace, sys%space, sys%ions%val_charge(), sys%kpoints)
    call sys%st%write_info(sys%namespace)
    ! if independent particles in N dimensions are being used, need to initialize them
    !  after masses are set to 1 in grid_init_stage_1 -> derivatives_init
    call modelmb_copy_masses (sys%st%modelmbparticles, sys%gr%der%masses)
    call elf_init(sys%namespace)

    sys%generate_epot = optional_default(generate_epot, .true.)

    call sys%supported_interactions_as_partner%add(CURRENT_TO_MXLL_FIELD)
    call sys%supported_interactions%add(MXLL_FIELD_TO_MEDIUM)
    sys%quantities(CURRENT)%required = .true.
    sys%quantities(CURRENT)%updated_on_demand = .false.
    call current_init(sys%current_calculator, sys%namespace)

    call profiling_out(prof)
    POP_SUB(electrons_constructor)
  end function electrons_constructor