dftb.F90

Go to the documentation of this file.

!! Copyright (C) 2020 F. Bonafé, H. Appel
!!
!! This program is free software; you can redistribute it and/or modify
!! it under the terms of the GNU General Public License as published by
!! the Free Software Foundation; either version 2, or (at your option)
!! any later version.
!!
!! This program is distributed in the hope that it will be useful,
!! but WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!! GNU General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with this program; if not, write to the Free Software
!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
!! 02110-1301, USA.
!!
 
#include "global.h"
 
module dftb_oct_m
  use algorithm_oct_m
  use debug_oct_m
#ifdef HAVE_DFTBPLUS
  use dftbplus
#endif
  use global_oct_m
  use interaction_oct_m
  use interaction_surrogate_oct_m
  use interactions_factory_oct_m
  use io_oct_m
  use ion_dynamics_oct_m
  use ions_oct_m
  use, intrinsic :: iso_fortran_env
  use lasers_oct_m
  use messages_oct_m
  use mpi_oct_m
  use multisystem_debug_oct_m
  use namespace_oct_m
  use parser_oct_m
  use profiling_oct_m
  use propagator_oct_m
  use propagator_verlet_oct_m
  use quantity_oct_m
  use species_oct_m
  use system_oct_m
  use tdfunction_oct_m
  use unit_oct_m
  use unit_system_oct_m
  use write_iter_oct_m
 
  implicit none
 
  private
  public ::           &
    dftb_t,    &
    dftb_init
 
  integer, parameter ::     &
    OUTPUT_COORDINATES = 1, &
    output_forces      = 2
 
  type, extends(system_t) :: dftb_t
    integer :: n_atom
    real(real64), allocatable :: coords(:,:), gradients(:,:)
    real(real64), allocatable :: acc(:,:)
    real(real64), allocatable :: tot_force(:,:)
    real(real64), allocatable :: vel(:,:)
    real(real64), allocatable :: prev_tot_force(:,:)
    integer, allocatable :: species(:)
    integer              :: dynamics
    real(real64), allocatable :: mass(:)
    real(real64), allocatable :: atom_charges(:,:)
    character(len=LABEL_LEN), allocatable  :: labels(:)
    real(real64), allocatable :: prev_acc(:,:,:)
    real(real64) :: scc_tolerance
    class(ions_t),     pointer :: ions => null()
    type(c_ptr) :: output_handle(2)
    type(ion_dynamics_t) :: ions_dyn
    class(lasers_t), pointer :: lasers => null()
    logical :: laser_field
    real(real64) :: field(3)
    real(real64) :: energy
#ifdef HAVE_DFTBPLUS
    type(TDftbPlus) :: dftbp
#endif
  contains
    procedure :: init_interaction => dftb_init_interaction
    procedure :: initialize => dftb_initialize
    procedure :: do_algorithmic_operation => dftb_do_algorithmic_operation
    procedure :: output_start => dftb_output_start
    procedure :: output_write => dftb_output_write
    procedure :: output_finish => dftb_output_finish
    procedure :: is_tolerance_reached => dftb_is_tolerance_reached
    procedure :: copy_quantities_to_interaction => dftb_copy_quantities_to_interaction
    procedure :: init_interaction_as_partner => dftb_init_interaction_as_partner
    procedure :: update_interactions_start => dftb_update_interactions_start
    procedure :: restart_write_data => dftb_restart_write_data
    procedure :: restart_read_data => dftb_restart_read_data
    procedure :: update_kinetic_energy => dftb_update_kinetic_energy
    final :: dftb_finalize
  end type dftb_t
 
  interface dftb_t
    procedure dftb_constructor
  end interface dftb_t
 
  integer, parameter :: &
    EHRENFEST = 1,      &
    bo        = 2
 
contains
 
  ! ---------------------------------------------------------
  function dftb_constructor(namespace) result(sys)
    class(dftb_t),       pointer    :: sys
    type(namespace_t),   intent(in) :: namespace
 
    push_sub(dftb_constructor)
 
#ifndef HAVE_DFTBPLUS
    message(1) = "DFTB+ system not available. This feature requires compiling Octopus with the DFTB+ library"
    call messages_fatal(1, namespace=namespace)
#endif
 
    allocate(sys)
 
    call dftb_init(sys, namespace)
 
    pop_sub(dftb_constructor)
  end function dftb_constructor
 
  ! ---------------------------------------------------------
  ! ---------------------------------------------------------
  subroutine dftb_init(this, namespace)
    class(dftb_t),     target, intent(inout) :: this
    type(namespace_t),         intent(in)    :: namespace
 
    integer :: ii, jj, ispec
    character(len=MAX_PATH_LEN) :: slako_dir
    character(len=1), allocatable  :: max_ang_mom(:)
    character(len=LABEL_LEN) :: this_max_ang_mom, this_label
    integer :: n_maxang_block
    type(block_t) :: blk
    real(real64) :: initial_temp
 
#ifdef HAVE_DFTBPLUS
    type(tdftbplusinput) :: input
    type(fnode), pointer :: proot, pgeo, pham, pdftb, pmaxang, pslakos, ptype2files, panalysis
    type(fnode), pointer :: pparseropts
    type(fnode), pointer :: pelecdyn, pperturb, plaser
#endif
 
    push_sub(dftb_init)
 
    this%namespace = namespace
 
    ! Currently this system does not support any interaction
    allocate(this%supported_interactions(0))
    allocate(this%supported_interactions_as_partner(0))
 
    ! Quantities updated by the algorithm
    call this%quantities%add(quantity_t("position", updated_on_demand = .false.))
    call this%quantities%add(quantity_t("velocity", updated_on_demand = .false.))
 
    call messages_print_with_emphasis(msg="DFTB+ System", namespace=namespace)
 
    this%ions => ions_t(namespace)
    this%n_atom = this%ions%natoms
    safe_allocate(this%coords(1:3, 1:this%n_atom))
    safe_allocate(this%acc(1:3, 1:this%n_atom))
    safe_allocate(this%vel(1:3, 1:this%n_atom))
    safe_allocate(this%tot_force(1:3, 1:this%n_atom))
    safe_allocate(this%prev_tot_force(1:3, 1:this%n_atom))
    safe_allocate(this%gradients(1:3, 1:this%n_atom))
    safe_allocate(this%species(1:this%n_atom))
    safe_allocate(this%mass(1:this%n_atom))
    safe_allocate(this%atom_charges(1:this%n_atom, 1))
    safe_allocate(this%labels(1:this%ions%nspecies))
    safe_allocate(max_ang_mom(1:this%ions%nspecies))
    max_ang_mom = "s"
 
    ispec = 1
    this%species(1) = 1
    this%labels(1) = trim(this%ions%atom(1)%label)
 
    this%coords = this%ions%pos
    ! mass is read from the default pseudopotential files
    this%mass = this%ions%mass
    do ii = 1, this%n_atom
      if ((ii > 1) .and. .not. (any(this%labels(1:ispec) == this%ions%atom(ii)%label))) then
        ispec = ispec + 1
        this%labels(ispec) = trim(this%ions%atom(ii)%label)
      end if
      do jj = 1, ispec
        if (trim(this%ions%atom(ii)%label) == trim(this%labels(jj))) then
          this%species(ii) = jj
        end if
      end do
    end do
    this%vel = m_zero
    this%tot_force = m_zero
 
    !%Variable MaxAngularMomentum
    !%Type block
    !%Section DFTBPlusInterface
    !%Description
    !% Specifies the highest angular momentum for each atom type. All orbitals up
    !% to that angular momentum will be included in the calculation.
    !% Possible values for the angular momenta are s, p, d, f.
    !% These are examples:
    !%
    !% <tt>%MaxAngularMomentum
    !% <br>&nbsp;&nbsp;'O'   | 'p'
    !% <br>&nbsp;&nbsp;'H'   | 's'
    !% <br>%</tt>
    !%End
    n_maxang_block = 0
    if (parse_block(namespace, 'MaxAngularMomentum', blk) == 0) then
      n_maxang_block = parse_block_n(blk)
      if (n_maxang_block /= this%ions%nspecies) then
        call messages_input_error(namespace, "MaxAngularMomentum", "Wrong number of species.")
      end if
 
      do ii = 1, n_maxang_block
        call parse_block_string(blk, ii-1, 0, this_label)
        call parse_block_string(blk, ii-1, 1, this_max_ang_mom)
        if (.not. any(["s","p","d","f"] == trim(adjustl(this_max_ang_mom)))) then
          call messages_input_error(namespace, "MaxAngularMomentum", "Wrong maximum angular momentum for element"//trim(this_label))
        end if
        do jj = 1, this%ions%nspecies
          if (trim(adjustl(this_label)) == trim(adjustl(this%labels(jj)))) then
            max_ang_mom(jj) = trim(adjustl(this_max_ang_mom))
          end if
        end do
      end do
    end if
    call parse_block_end(blk)
 
    !%Variable SlakoDir
    !%Type string
    !%Default "./"
    !%Section Execution::IO
    !%Description
    !% Folder containing the Slako files
    !%End
    call parse_variable(namespace, 'SlakoDir', './', slako_dir)
 
 
    ! Dynamics variables
 
    call ion_dynamics_init(this%ions_dyn, namespace, this%ions, .false.)
 
    call parse_variable(namespace, 'TDDynamics', bo, this%dynamics)
    call messages_print_var_option('TDDynamics', this%dynamics, namespace=namespace)
    if (this%dynamics == bo) then
      call ion_dynamics_unfreeze(this%ions_dyn)
    end if
 
    !%Variable InitialIonicTemperature
    !%Type float
    !%Default 0.0
    !%Section DFTBPlusInterface
    !%Description
    !% If this variable is present, the ions will have initial velocities
    !% velocities to the atoms following a Boltzmann distribution with
    !% this temperature (in Kelvin). Used only if <tt>TDDynamics = Ehrenfest</tt>
    !% and  <tt>MoveIons = yes</tt>.
    !%End
    call parse_variable(namespace, 'InitialIonicTemperature', m_zero, initial_temp, unit = unit_kelvin)
 
    this%lasers => lasers_t(namespace)
    call lasers_parse_external_fields(this%lasers)
    if (this%lasers%no_lasers > 0) then
      this%laser_field = .true.
    else
      this%laser_field = .false.
    end if
 
#ifdef HAVE_DFTBPLUS
    call tdftbplus_init(this%dftbp)
 
    call this%dftbp%getEmptyInput(input)
    call input%getRootNode(proot)
    call setchild(proot, "Geometry", pgeo)
    call setchildvalue(pgeo, "Periodic", .false.)
    call setchildvalue(pgeo, "TypeNames", this%labels(1:this%ions%nspecies))
    call setchildvalue(pgeo, "TypesAndCoordinates", reshape(this%species, [1, size(this%species)]), this%coords)
    call setchild(proot, "Hamiltonian", pham)
    call setchild(pham, "Dftb", pdftb)
    call setchildvalue(pdftb, "Scc", .true.)
 
    !%Variable SccTolerance
    !%Type float
    !%Section DFTBPlusInterface
    !%Description
    !% Self-consistent-charges convergence tolerance. Once this
    !% tolerance has been achieved the SCC cycle will stop.
    !%End
    call parse_variable(namespace, 'SccTolerance', 1e-9_real64, this%scc_tolerance)
    call messages_print_var_value('SccTolerance', this%scc_tolerance, namespace=namespace)
    call setchildvalue(pdftb, "SccTolerance", this%scc_tolerance)
 
    ! sub-block inside hamiltonian for the maximum angular momenta
    call setchild(pdftb, "MaxAngularMomentum", pmaxang)
    ! explicitly set the maximum angular momenta for the species
    do ii = 1, this%ions%nspecies
      call setchildvalue(pmaxang, this%labels(ii), max_ang_mom(ii))
    end do
 
    ! get the SK data
    ! You should provide the skfiles as found in the external/slakos/origin/mio-1-1/ folder. These can
    ! be downloaded with the utils/get_opt_externals script
    call setchild(pdftb, "SlaterKosterFiles", pslakos)
    call setchild(pslakos, "Type2FileNames", ptype2files)
    call setchildvalue(ptype2files, "Prefix", slako_dir)
    call setchildvalue(ptype2files, "Separator", "-")
    call setchildvalue(ptype2files, "Suffix", ".skf")
 
    !  set up analysis options
    call setchild(proot, "Analysis", panalysis)
    call setchildvalue(panalysis, "CalculateForces", .true.)
 
    call setchild(proot, "ParserOptions", pparseropts)
    call setchildvalue(pparseropts, "ParserVersion", 5)
 
    if (this%dynamics == ehrenfest) then
      call setchild(proot, "ElectronDynamics", pelecdyn)
      call setchildvalue(pelecdyn, "IonDynamics", this%ions_dyn%is_active())
      if (this%ions_dyn%is_active()) then
        call setchildvalue(pelecdyn, "InitialTemperature", initial_temp)
      end if
 
      ! initialize with wrong arguments for the moment, will be overriden later
      call setchildvalue(pelecdyn, "Steps", 1)
      call setchildvalue(pelecdyn, "TimeStep", m_one)
      call setchild(pelecdyn, "Perturbation", pperturb)
      if (this%laser_field) then
        call setchild(pperturb, "Laser", plaser)
        call setchildvalue(plaser, "PolarizationDirection", [ m_one , m_zero , m_zero ])
        call setchildvalue(plaser, "LaserEnergy", m_one)
        call setchildvalue(pelecdyn, "FieldStrength", m_one)
      else
        call setchild(pperturb, "None", plaser)
      end if
    end if
 
    message(1) = 'Input tree in HSD format:'
    call messages_info(1, namespace=namespace)
    call dumphsd(input%hsdTree, stdout)
 
    ! initialise the DFTB+ calculator
    call this%dftbp%setupCalculator(input)
    call this%dftbp%setGeometry(this%coords)
#endif
 
    pop_sub(dftb_init)
  end subroutine dftb_init
 
  ! ---------------------------------------------------------
  subroutine dftb_init_interaction(this, interaction)
    class(dftb_t),        target, intent(inout) :: this
    class(interaction_t),         intent(inout) :: interaction
 
    push_sub(dftb_init_interaction)
 
    select type (interaction)
    class default
      message(1) = "Trying to initialize an unsupported interaction by DFTB+."
      call messages_fatal(1, namespace=this%namespace)
    end select
 
    pop_sub(dftb_init_interaction)
  end subroutine dftb_init_interaction
 
  ! ---------------------------------------------------------
  subroutine dftb_initialize(this)
    class(dftb_t), intent(inout) :: this
 
    push_sub(dftb_initialize)
 
#ifdef HAVE_DFTBPLUS
    select type (algo => this%algo)
    class is (propagator_t)
      select case (this%dynamics)
      case (bo)
        call this%dftbp%getGradients(this%gradients)
        this%tot_force = -this%gradients
      case (ehrenfest)
        call this%dftbp%getEnergy(this%energy)
        call this%dftbp%initializeTimeProp(algo%dt, this%laser_field, .false.)
      end select
    end select
#endif
 
    pop_sub(dftb_initialize)
  end subroutine dftb_initialize
 
  ! ---------------------------------------------------------
  logical function dftb_do_algorithmic_operation(this, operation, updated_quantities) result(done)
    class(dftb_t),                  intent(inout) :: this
    class(algorithmic_operation_t), intent(in)    :: operation
    character(len=:), allocatable,  intent(out)   :: updated_quantities(:)
 
    integer :: ii, jj, il
    type(tdf_t) :: ff, phi
    complex(real64) :: amp, pol(3)
    real(real64) :: time, omega
 
    push_sub(dftb_do_algorithmic_operation)
    call profiling_in(trim(this%namespace%get())//":"//trim(operation%id))
 
    select type (algo => this%algo)
    class is (propagator_t)
 
      done = .true.
      select case (this%dynamics)
      case (bo)
        ! Born-Oppenheimer dynamics
        select case (operation%id)
        case (store_current_status)
          ! Do nothing
 
        case (verlet_start)
          safe_allocate(this%prev_acc(1:3, 1:this%n_atom, 1))
          do jj = 1, this%n_atom
            this%acc(1:3, jj) = this%tot_force(1:3, jj) / this%mass(jj)
          end do
 
        case (verlet_finish)
          safe_deallocate_a(this%prev_acc)
 
        case (verlet_update_pos)
          do jj = 1, this%n_atom
            this%coords(1:3, jj) = this%coords(1:3, jj) + algo%dt * this%vel(1:3, jj) &
              + m_half * algo%dt**2 * this%acc(1:3, jj)
          end do
          updated_quantities = ["position"]
 
        case (verlet_compute_acc)
          do ii = size(this%prev_acc, dim=3) - 1, 1, -1
            this%prev_acc(1:3, 1:this%n_atom, ii + 1) = this%prev_acc(1:3, 1:this%n_atom, ii)
          end do
          this%prev_acc(1:3, 1:this%n_atom, 1) = this%acc(1:3, 1:this%n_atom)
#ifdef HAVE_DFTBPLUS
          call this%dftbp%setGeometry(this%coords)
          call this%dftbp%getGradients(this%gradients)
          this%tot_force = -this%gradients
#endif
          do jj = 1, this%n_atom
            this%acc(1:3, jj) = this%tot_force(1:3, jj) / this%mass(jj)
          end do
 
        case (verlet_compute_vel)
          this%vel(1:3, 1:this%n_atom) = this%vel(1:3, 1:this%n_atom) &
            + m_half * algo%dt * (this%prev_acc(1:3, 1:this%n_atom, 1) + &
            this%acc(1:3, 1:this%n_atom))
          updated_quantities = ["velocity"]
 
        case default
          done = .false.
        end select
 
      case (ehrenfest)
        ! Ehrenfest dynamics
        select case (operation%id)
        case (store_current_status)
          ! Do nothing
        case (verlet_start)
          !Do nothing
        case (verlet_finish)
          !Do nothing
        case (verlet_update_pos)
          this%field = m_zero
          time = this%iteration%value()
          do il = 1, this%lasers%no_lasers
            ! get properties of laser
            call laser_get_f(this%lasers%lasers(il), ff)
            call laser_get_phi(this%lasers%lasers(il), phi)
            omega = laser_carrier_frequency(this%lasers%lasers(il))
            pol = laser_polarization(this%lasers%lasers(il))
            ! calculate electric field from laser
            amp = tdf(ff, time) * exp(m_zi * (omega*time + tdf(phi, time)))
            this%field(1:3) = this%field(1:3) + real(amp*pol(1:3), real64)
          end do
#ifdef HAVE_DFTBPLUS
          call this%dftbp%setTdElectricField(this%field)
          call this%dftbp%doOneTdStep(this%iteration%counter(), atomnetcharges=this%atom_charges, coord=this%coords,&
            force=this%tot_force, energy=this%energy)
#endif
        case (verlet_compute_acc)
          !Do nothing
        case (verlet_compute_vel)
          !Do nothing
        case default
          done = .false.
        end select
 
      end select
 
    end select
 
    call profiling_out(trim(this%namespace%get())//":"//trim(operation%id))
    pop_sub(dftb_do_algorithmic_operation)
  end function dftb_do_algorithmic_operation
 
  ! ---------------------------------------------------------
  logical function dftb_is_tolerance_reached(this, tol) result(converged)
    class(dftb_t),        intent(in)    :: this
    real(real64),         intent(in)    :: tol
 
    push_sub(dftb_is_tolerance_reached)
 
    ! this routine is never called at present, no reason to be here
    assert(.false.)
    converged = .false.
 
    pop_sub(dftb_is_tolerance_reached)
  end function dftb_is_tolerance_reached
 
  ! ---------------------------------------------------------
  subroutine dftb_output_start(this)
    class(dftb_t), intent(inout) :: this
 
    push_sub(dftb_output_start)
 
    select type (algo => this%algo)
    class is (propagator_t)
      ! Create output handle
      call io_mkdir('td.general', this%namespace)
      if (mpi_grp_is_root(mpi_world)) then
        call write_iter_init(this%output_handle(output_coordinates), 0, algo%dt, &
          trim(io_workpath("td.general/coordinates", this%namespace)))
        call write_iter_init(this%output_handle(output_forces), 0, algo%dt, &
          trim(io_workpath("td.general/forces", this%namespace)))
      end if
 
      ! Output info for first iteration
      call this%output_write()
    end select
 
    pop_sub(dftb_output_start)
  end subroutine dftb_output_start
 
  ! ---------------------------------------------------------
  subroutine dftb_output_finish(this)
    class(dftb_t), intent(inout) :: this
 
    push_sub(dftb_output_finish)
 
    select type (algo => this%algo)
    class is (propagator_t)
      if (mpi_grp_is_root(mpi_world)) then
        call write_iter_end(this%output_handle(output_coordinates))
        call write_iter_end(this%output_handle(output_forces))
      end if
    end select
 
    pop_sub(dftb_output_finish)
  end subroutine dftb_output_finish
 
  ! ---------------------------------------------------------
  subroutine dftb_output_write(this)
    class(dftb_t), intent(inout) :: this
 
    integer :: idir, iat, iout
    character(len=50) :: aux
    character(1) :: out_label(2)
    real(real64) :: tmp(3)
 
    if (.not. mpi_grp_is_root(mpi_world)) return ! only first node outputs
 
    push_sub(dftb_output_write)
    call profiling_in(trim(this%namespace%get())//":"//"OUTPUT_WRITE")
 
    select type (algo => this%algo)
    class is (propagator_t)
      out_label(1) = "x"
      out_label(2) = "f"
 
      if (this%iteration%counter() == 0) then
        ! header
        do iout = 1, 2
          call write_iter_clear(this%output_handle(iout))
          call write_iter_string(this%output_handle(iout),'#####################################################################')
          call write_iter_nl(this%output_handle(iout))
          call write_iter_string(this%output_handle(iout),'# HEADER')
          call write_iter_nl(this%output_handle(iout))
 
          ! first line: column names
          call write_iter_header_start(this%output_handle(iout))
 
          do iat = 1, this%n_atom
            do idir = 1, 3
              write(aux, '(a1,a1,i3,a1,i3,a1)') out_label(iout),'(', iat, ',', idir, ')'
              call write_iter_header(this%output_handle(iout), aux)
            end do
          end do
          call write_iter_nl(this%output_handle(iout))
 
          ! second line: units
          call write_iter_string(this%output_handle(iout), '#[Iter n.]')
          call write_iter_header(this%output_handle(iout), '[' // trim(units_abbrev(units_out%time)) // ']')
        end do
 
        do iat = 1, this%n_atom
          do idir = 1, 3
            call write_iter_header(this%output_handle(output_coordinates), '[' // trim(units_abbrev(units_out%length)) // ']')
            call write_iter_header(this%output_handle(output_forces), '[' // trim(units_abbrev(units_out%force)) // ']')
          end do
        end do
 
        do iout = 1, 2
          call write_iter_nl(this%output_handle(iout))
          call write_iter_string(this%output_handle(iout),'#######################################################################')
          call write_iter_nl(this%output_handle(iout))
        end do
      end if
 
      call write_iter_start(this%output_handle(output_coordinates))
      call write_iter_start(this%output_handle(output_forces))
 
      do iat = 1, this%n_atom
        ! Position
        tmp(1:3) = units_from_atomic(units_out%length, this%coords(1:3, iat))
        call write_iter_double(this%output_handle(output_coordinates), tmp, 3)
        ! Force
        tmp(1:3) = units_from_atomic(units_out%force, this%tot_force(1:3, iat))
        call write_iter_double(this%output_handle(output_forces), tmp, 3)
      end do
 
      call write_iter_nl(this%output_handle(output_coordinates))
      call write_iter_nl(this%output_handle(output_forces))
 
    end select
 
    call profiling_out(trim(this%namespace%get())//":"//"OUTPUT_WRITE")
    pop_sub(dftb_output_write)
  end subroutine dftb_output_write
 
  ! ---------------------------------------------------------
  subroutine dftb_init_interaction_as_partner(partner, interaction)
    class(dftb_t),             intent(in)    :: partner
    class(interaction_surrogate_t), intent(inout) :: interaction
 
    push_sub(dftb_init_interaction_as_partner)
 
    select type (interaction)
    class default
      message(1) = "Unsupported interaction."
      call messages_fatal(1, namespace=partner%namespace)
    end select
 
    pop_sub(dftb_init_interaction_as_partner)
  end subroutine dftb_init_interaction_as_partner
 
  ! ---------------------------------------------------------
  subroutine dftb_copy_quantities_to_interaction(partner, interaction)
    class(dftb_t),          intent(inout) :: partner
    class(interaction_surrogate_t), intent(inout) :: interaction
 
    push_sub(dftb_copy_quantities_to_interaction)
 
    select type (interaction)
    class default
      message(1) = "Unsupported interaction."
      call messages_fatal(1, namespace=partner%namespace)
    end select
 
    pop_sub(dftb_copy_quantities_to_interaction)
  end subroutine dftb_copy_quantities_to_interaction
 
  ! ---------------------------------------------------------
  subroutine dftb_update_interactions_start(this)
    class(dftb_t), intent(inout) :: this
 
    push_sub(dftb_update_interactions_start)
 
    ! Store previous force, as it is used as SCF criterium
    this%prev_tot_force(1:3, 1:this%n_atom) = this%tot_force(1:3, 1:this%n_atom)
 
    pop_sub(dftb_update_interactions_start)
  end subroutine dftb_update_interactions_start
 
  ! ---------------------------------------------------------
  subroutine dftb_update_kinetic_energy(this)
    class(dftb_t), intent(inout) :: this
 
    push_sub(dftb_update_kinetic_energy)
 
    this%kinetic_energy = m_zero
 
    pop_sub(dftb_update_kinetic_energy)
  end subroutine dftb_update_kinetic_energy
 
  ! ---------------------------------------------------------
  subroutine dftb_restart_write_data(this)
    class(dftb_t), intent(inout) :: this
 
    push_sub(dftb_restart_write_data)
 
    message(1) = "DFTB system "//trim(this%namespace%get())//" cannot write restart data."
    call messages_warning(1, namespace=this%namespace)
 
    pop_sub(dftb_restart_write_data)
  end subroutine dftb_restart_write_data
 
  ! ---------------------------------------------------------
  ! this function returns true if restart data could be read
  logical function dftb_restart_read_data(this)
    class(dftb_t), intent(inout) :: this
 
    push_sub(dftb_restart_read_data)
 
    ! restarting not yet supported
    dftb_restart_read_data = .false.
 
    pop_sub(dftb_restart_read_data)
  end function dftb_restart_read_data
 
  ! ---------------------------------------------------------
  subroutine dftb_finalize(this)
    type(dftb_t), intent(inout) :: this
 
    push_sub(dftb_finalize)
 
    safe_deallocate_a(this%coords)
    safe_deallocate_a(this%acc)
    safe_deallocate_a(this%vel)
    safe_deallocate_a(this%tot_force)
    safe_deallocate_a(this%prev_tot_force)
    safe_deallocate_a(this%gradients)
    safe_deallocate_a(this%species)
    safe_deallocate_a(this%mass)
    call ion_dynamics_end(this%ions_dyn)
 
    deallocate(this%ions)
 
    ! No call to safe_deallocate macro here, as it gives an ICE with gfortran
    if (associated(this%lasers)) then
      deallocate(this%lasers)
    end if
 
#ifdef HAVE_DFTBPLUS
    call tdftbplus_destruct(this%dftbp)
#endif
 
    call system_end(this)
 
    pop_sub(dftb_finalize)
  end subroutine dftb_finalize
 
end module dftb_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: