In the multisysten framework, the total energy (total_energy) of a calculation consists of various contributions, where we follow the standard definitions of
thermodynamics (see here)
The kinetic energy (kinetic_energy)
The internal energy (internal_energy)
The potential energy (potential_energy)
Kinetic energy
We use the term kinetic energy in the usual sense. The only slight exception is for container systems:
Ideally, for containers, the kinetic energy should be the kinetic energy with respect to the centre of mass of the container.
The kinetic energies of the constituents, relative to the centre of mass frame, should be accounted for as part of the internal energy of the container.
multisystem_update_kinetic_energy()
recursive subroutine multisystem_update_kinetic_energy(this)class(multisystem_t),intent(inout)::thisPUSH_SUB(multisystem_update_kinetic_energy)! We currently do not have the center of mass coordinates implemented for multisystems,
! hence we set the kinetic energy to zero.
! The kinetic energies of the constituents are contributing to the internal energy.
this%kinetic_energy=M_ZEROPOP_SUB(multisystem_update_kinetic_energy)end subroutine multisystem_update_kinetic_energy
Specific systems need their specific routines to calculate the kinetic energy.
For the Maxwell system, this is the so-called electromagnetic energy.
Interaction energies
Everything else is treated as an interaction energy. Here we have to distinguish between the interactions between two different systems (which can be of the same type), and the interaction of a system with itself (intra-interaction).
inter-interactions
The interactions between two different systems do not require any further thought, as by definition no physical self-interaction can occur.
As the method to calculate the interactions are part of the interaction class, it is independent of the actual system and can be implemented
at the level of the classsystem_t.
system_update_potential_energy()
subroutine system_update_potential_energy(this)class(system_t),intent(inout)::thistype(interaction_iterator_t)::iterclass(interaction_t),pointer::interactionPUSH_SUB(system_update_potential_energy)this%potential_energy=M_ZEROcall iter%start(this%interactions)do while(iter%has_next())interaction=>iter%get_next()if(.not.interaction%intra_interaction)then
call interaction%calculate_energy()this%potential_energy=this%potential_energy+interaction%energyend if
end do
POP_SUB(system_update_potential_energy)end subroutine system_update_potential_energy
The exception are containers. Here we need to loop over the constituents. In order to distinguish inter- from intra-interactions, we need to
query each interaction for its interaction partner, and skip the interaction, if the partner is part of the container.
multisystem_update_potential_energy()
subroutine multisystem_update_potential_energy(this)class(multisystem_t),intent(inout)::thistype(system_iterator_t)::system_iterclass(system_t),pointer::system
type(interaction_iterator_t)::interaction_iterclass(interaction_t),pointer::interactiontype(system_list_t)::flat_listPUSH_SUB(multisystem_update_potential_energy)this%potential_energy=M_ZERO!> We need to handle interactions of the container itself:
call system_update_potential_energy(this)!> generate a list of all systems inside the container and its subcontainers:
call this%get_flat_list(flat_list)!> loop over all systems inside the container
call system_iter%start(flat_list)do while(system_iter%has_next())system=>system_iter%get_next()!> Even though we are not using the potential energy of the subsystems here, we need to trigger their calculation
call system%update_potential_energy()!> loop over all interactions and discard those with partners inside the container
call interaction_iter%start(system%interactions)do while(interaction_iter%has_next())interaction=>interaction_iter%get_next()select type(interaction)class is(interaction_with_partner_t)if(.not.flat_list%contains(interaction%partner).and..not.interaction%intra_interaction)then
call interaction%calculate_energy()this%potential_energy=this%potential_energy+interaction%energyend if
class is(interaction_t)call interaction%calculate_energy()this%potential_energy=this%potential_energy+interaction%energyend select
end do
end do
POP_SUB(multisystem_update_potential_energy)end subroutine multisystem_update_potential_energy
intra-interactions
Systems may contain more than one physical particle (e.g. the electrons, a set of ions or container systems). In order to account for the interaction of these particles with other particles of the same system, we decided to treat this case as a system interacting with itself, which we call intra-interaction.
In some cases, such as a single particle, this intra interaction has to be zero, while in other cases with many particles, the interactions have to be calculated, where – of course – the interaction of one particle with itself has to be removed (at least approximatively).
Another important aspect of the implementation is that Octopus deals with one-sided interactions. This has the implication that there is no double counting when calculating
the interaction energies. Both contributions have to be counted:
for each system, we add the the energy of the system in the field of the partner.
system_update_internal_energy()
subroutine system_update_internal_energy(this)class(system_t),intent(inout)::thistype(interaction_iterator_t)::iterclass(interaction_t),pointer::interactionPUSH_SUB(system_update_internal_energy)this%internal_energy=M_ZEROcall iter%start(this%interactions)do while(iter%has_next())interaction=>iter%get_next()if(interaction%intra_interaction)then
call interaction%calculate_energy()this%internal_energy=this%internal_energy+interaction%energyend if
end do
POP_SUB(system_update_internal_energy)end subroutine system_update_internal_energy
For containers, this means that the interaction energy contains the complete interaction energies of all subsystems in the container, plus the energies of the subsystems in the field of all other systems (not part of that container).
Note that the internal_energy of a container is not the sum of the self_energies of the constituents, and also the external energy of a container is not
the sum of the external energies of the constituents!
multisystem_update_internal_energy()
recursive subroutine multisystem_update_internal_energy(this)class(multisystem_t),intent(inout)::thisclass(system_t),pointer::system
class(system_t),pointer::system_2type(system_iterator_t)::system_itertype(system_iterator_t)::system_iter_2PUSH_SUB(multisystem_update_internal_energy)! The internal energy of the multisystem contains the kinetic and internal energies of the consistuents
!TODO: the kinetic energy wrt the center of mass motion should be subtracted.
this%internal_energy=M_ZEROcall system_iter%start(this%list)do while(system_iter%has_next())system=>system_iter%get_next()!> First add the kinetic energies of the subsystems
call system%update_kinetic_energy()this%internal_energy=this%internal_energy+system%kinetic_energy!> First add the internal energies of the subsystems
call system%update_internal_energy()this%internal_energy=this%internal_energy+system%internal_energy!> Now add the (inter-) interactions between the systems in the container.
call system_iter_2%start(this%list)do while(system_iter_2%has_next())system_2=>system_iter_2%get_next()!> exclude self-interactions (intra-interactions) as they are included in the internal energy
!! of the subsystem, which was already added above.
if(.not.associated(system,system_2))then
this%internal_energy=this%internal_energy+multisystem_pair_energy(system,system_2)end if
end do! system_iter_2
end do! system_iter
POP_SUB(multisystem_update_internal_energy)end subroutine multisystem_update_internal_energy
Total energy
The total energy of the whole simulation (top level container) is clearly defined. It contains the sum of all internal and interaction energies.
For a specific system (e.g. electrons), the total energy also is the sum of the internal energy and the interaction energies (including the intra interaction energy).
system_update_total_energy()
subroutine system_update_total_energy(this)class(system_t),intent(inout)::thisPUSH_SUB(system_update_total_energy)call this%update_kinetic_energy()this%total_energy=this%kinetic_energy!> External energy as the sum over interaction energies
call this%update_potential_energy()this%total_energy=this%total_energy+this%potential_energy!> Self energy arises from the interaction with itself
call this%update_internal_energy()this%total_energy=this%total_energy+this%internal_energyPOP_SUB(system_update_total_energy)end subroutine system_update_total_energy
Total energies and interaction energies depend on the state of other systems, and therefore should not be used for convergence criteria for a given system. Only the internal energy is safe to use here.