f77_interface.F

Go to the documentation of this file.

!--------------------------------------------------------------------------------------------------!
!   CP2K: A general program to perform molecular dynamics simulations                              !
!   Copyright 2000-2026 CP2K developers group <https://cp2k.org>                                   !
!                                                                                                  !
!   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
!--------------------------------------------------------------------------------------------------!
 
! **************************************************************************************************
!> \brief interface to use cp2k as library
!> \note
!>      useful additions for the future would be:
!>      - string(path) based set/get of simple values (to change the new
!>        input during the run and extract more data (energy types for example).
!>      - set/get of a subset of atoms
!> \par History
!>      07.2004 created [fawzi]
!>      11.2004 parallel version [fawzi]
!> \author fawzi & Johanna
! **************************************************************************************************
MODULE f77_interface
   USE base_hooks,                      ONLY: cp_abort_hook,&
                                              cp_warn_hook,&
                                              timeset_hook,&
                                              timestop_hook
   USE bibliography,                    ONLY: add_all_references
   USE cell_methods,                    ONLY: init_cell
   USE cell_types,                      ONLY: cell_type
   USE cp2k_info,                       ONLY: get_runtime_info
   USE cp_dbcsr_api,                    ONLY: dbcsr_finalize_lib,&
                                              dbcsr_init_lib
   USE cp_dlaf_utils_api,               ONLY: cp_dlaf_finalize,&
                                              cp_dlaf_free_all_grids
   USE cp_error_handling,               ONLY: cp_error_handling_setup
   USE cp_files,                        ONLY: open_file
   USE cp_log_handling,                 ONLY: &
        cp_add_default_logger, cp_default_logger_stack_size, cp_failure_level, &
        cp_get_default_logger, cp_logger_create, cp_logger_release, cp_logger_retain, &
        cp_logger_type, cp_rm_default_logger, cp_to_string
   USE cp_output_handling,              ONLY: cp_iterate
   USE cp_output_handling_openpmd,      ONLY: cp_openpmd_output_finalize
   USE cp_result_methods,               ONLY: get_results,&
                                              test_for_result
   USE cp_result_types,                 ONLY: cp_result_type
   USE cp_subsys_types,                 ONLY: cp_subsys_get,&
                                              cp_subsys_set,&
                                              cp_subsys_type,&
                                              unpack_subsys_particles
   USE dbm_api,                         ONLY: dbm_library_finalize,&
                                              dbm_library_init
   USE dbt_api,                         ONLY: dbt_library_finalize,&
                                              dbt_library_init
   USE eip_environment,                 ONLY: eip_init
   USE eip_environment_types,           ONLY: eip_env_create,&
                                              eip_environment_type
   USE embed_main,                      ONLY: embed_create_force_env
   USE embed_types,                     ONLY: embed_env_type
   USE environment,                     ONLY: cp2k_finalize,&
                                              cp2k_init,&
                                              cp2k_read,&
                                              cp2k_setup
   USE fist_main,                       ONLY: fist_create_force_env
   USE force_env_methods,               ONLY: force_env_calc_energy_force,&
                                              force_env_create
   USE force_env_types,                 ONLY: &
        force_env_get, force_env_get_frc, force_env_get_natom, force_env_get_nparticle, &
        force_env_get_pos, force_env_get_vel, force_env_release, force_env_retain, force_env_set, &
        force_env_type, multiple_fe_list
   USE fp_types,                        ONLY: fp_env_create,&
                                              fp_env_read,&
                                              fp_env_write,&
                                              fp_type
   USE global_types,                    ONLY: global_environment_type,&
                                              globenv_create,&
                                              globenv_release
   USE grid_api,                        ONLY: grid_library_finalize,&
                                              grid_library_init
   USE input_constants,                 ONLY: &
        do_eip, do_embed, do_fist, do_ipi, do_mixed, do_nnp, do_qmmm, do_qmmmx, do_qs, do_sirius
   USE input_cp2k_check,                ONLY: check_cp2k_input
   USE input_cp2k_force_eval,           ONLY: create_force_eval_section
   USE input_cp2k_read,                 ONLY: empty_initial_variables,&
                                              read_input
   USE input_enumeration_types,         ONLY: enum_i2c,&
                                              enumeration_type
   USE input_keyword_types,             ONLY: keyword_get,&
                                              keyword_type
   USE input_section_types,             ONLY: &
        section_get_keyword, section_release, section_type, section_vals_duplicate, &
        section_vals_get, section_vals_get_subs_vals, section_vals_release, &
        section_vals_remove_values, section_vals_retain, section_vals_type, section_vals_val_get, &
        section_vals_write
   USE ipi_environment,                 ONLY: ipi_init
   USE ipi_environment_types,           ONLY: ipi_environment_type
   USE kinds,                           ONLY: default_path_length,&
                                              default_string_length,&
                                              dp
   USE libint_wrapper,                  ONLY: cp_libint_static_cleanup,&
                                              cp_libint_static_init
   USE machine,                         ONLY: default_output_unit,&
                                              m_chdir,&
                                              m_getcwd,&
                                              m_memory
   USE message_passing,                 ONLY: mp_comm_type,&
                                              mp_comm_world,&
                                              mp_para_env_release,&
                                              mp_para_env_type,&
                                              mp_query_thread_level,&
                                              mp_world_finalize,&
                                              mp_world_init
   USE metadynamics_types,              ONLY: meta_env_type
   USE metadynamics_utils,              ONLY: metadyn_read
   USE mixed_environment_types,         ONLY: mixed_environment_type
   USE mixed_main,                      ONLY: mixed_create_force_env
   USE mp_perf_env,                     ONLY: add_mp_perf_env,&
                                              get_mp_perf_env,&
                                              mp_perf_env_release,&
                                              mp_perf_env_retain,&
                                              mp_perf_env_type,&
                                              rm_mp_perf_env
   USE nnp_environment,                 ONLY: nnp_init
   USE nnp_environment_types,           ONLY: nnp_type
   USE offload_api,                     ONLY: offload_get_chosen_device,&
                                              offload_get_device_count,&
                                              offload_init,&
                                              offload_set_chosen_device
   USE periodic_table,                  ONLY: init_periodic_table
   USE pw_fpga,                         ONLY: pw_fpga_finalize,&
                                              pw_fpga_init
   USE pw_gpu,                          ONLY: pw_gpu_finalize,&
                                              pw_gpu_init
   USE pwdft_environment,               ONLY: pwdft_init
   USE pwdft_environment_types,         ONLY: pwdft_env_create,&
                                              pwdft_environment_type
   USE qmmm_create,                     ONLY: qmmm_env_create
   USE qmmm_types,                      ONLY: qmmm_env_type
   USE qmmmx_create,                    ONLY: qmmmx_env_create
   USE qmmmx_types,                     ONLY: qmmmx_env_type
   USE qs_environment,                  ONLY: qs_init
   USE qs_environment_types,            ONLY: get_qs_env,&
                                              qs_env_create,&
                                              qs_environment_type
   USE reference_manager,               ONLY: remove_all_references
   USE sirius_interface,                ONLY: cp_sirius_finalize,&
                                              cp_sirius_init,&
                                              cp_sirius_is_initialized
   USE string_table,                    ONLY: string_table_allocate,&
                                              string_table_deallocate
   USE timings,                         ONLY: add_timer_env,&
                                              get_timer_env,&
                                              rm_timer_env,&
                                              timer_env_release,&
                                              timer_env_retain,&
                                              timings_register_hooks
   USE timings_types,                   ONLY: timer_env_type
   USE virial_types,                    ONLY: virial_type
   USE xc_gauxc_interface,              ONLY: cp_gauxc_status_type,&
                                              gauxc_check_status,&
                                              gauxc_finalize,&
                                              gauxc_init
   USE xc_libxc,                        ONLY: libxc_release_workers
#include "./base/base_uses.f90"
 
   IMPLICIT NONE
   PRIVATE
 
   LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .true.
   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'f77_interface'
 
! **************************************************************************************************
   TYPE f_env_p_type
      TYPE(f_env_type), POINTER :: f_env => null()
   END TYPE f_env_p_type
 
! **************************************************************************************************
   TYPE f_env_type
      INTEGER :: id_nr = 0
      LOGICAL :: is_dirty = .true.
      TYPE(force_env_type), POINTER      :: force_env => null()
      TYPE(cp_logger_type), POINTER      :: logger => null()
      TYPE(timer_env_type), POINTER      :: timer_env => null()
      TYPE(mp_perf_env_type), POINTER    :: mp_perf_env => null()
      CHARACTER(len=default_path_length) :: my_path = "", old_path = ""
   END TYPE f_env_type
 
   TYPE(f_env_p_type), DIMENSION(:), POINTER, SAVE :: f_envs
   TYPE(mp_para_env_type), POINTER, SAVE :: default_para_env
   LOGICAL, SAVE :: module_initialized = .false.
   INTEGER, SAVE :: last_f_env_id = 0, n_f_envs = 0
 
   PUBLIC :: default_para_env
   PUBLIC :: init_cp2k, finalize_cp2k, configure_dbcsr
   PUBLIC :: create_force_env, destroy_force_env, set_pos, get_pos, &
             get_force, calc_energy_force, get_energy, get_stress_tensor, &
             calc_energy, calc_force, check_input, get_natom, get_nparticle, &
             f_env_add_defaults, f_env_rm_defaults, f_env_type, &
             f_env_get_from_id, &
             set_vel, set_cell, get_cell, get_qmmm_cell, get_result_r1
CONTAINS
 
! **************************************************************************************************
!> \brief returns the position of the force env corresponding to the given id
!> \param env_id the id of the requested environment
!> \return ...
!> \author fawzi
!> \note
!>      private utility function
! **************************************************************************************************
   FUNCTION get_pos_of_env(env_id) RESULT(res)
      INTEGER, INTENT(in)                                :: env_id
      INTEGER                                            :: res
 
      INTEGER                                            :: env_pos, isub
 
      env_pos = -1
      DO isub = 1, n_f_envs
         IF (f_envs(isub)%f_env%id_nr == env_id) THEN
            env_pos = isub
         END IF
      END DO
      res = env_pos
   END FUNCTION get_pos_of_env
 
! **************************************************************************************************
!> \brief initializes cp2k, needs to be called once before using any of the
!>      other functions when using cp2k as library
!> \param init_mpi if the mpi environment should be initialized
!> \param ierr returns a number different from 0 if there was an error
!> \param mpi_comm an existing mpi communicator (if not given mp_comm_world
!>      will be used)
!> \author fawzi
! **************************************************************************************************
   SUBROUTINE init_cp2k(init_mpi, ierr, mpi_comm)
      LOGICAL, INTENT(in)                                :: init_mpi
      INTEGER, INTENT(out)                               :: ierr
      TYPE(mp_comm_type), INTENT(in), OPTIONAL           :: mpi_comm
 
      INTEGER                                            :: offload_chosen_device, &
                                                            offload_device_count, unit_nr
      TYPE(cp_gauxc_status_type)                         :: gauxc_status
      TYPE(cp_logger_type), POINTER                      :: logger
 
      IF (.NOT. module_initialized) THEN
         ! install error handler hooks
         CALL cp_error_handling_setup()
 
         ! install timming handler hooks
         CALL timings_register_hooks()
 
         ! get runtime information
         CALL get_runtime_info()
 
         ! Intialize CUDA/HIP before MPI
         ! Needed for HIP on ALPS & LUMI
         CALL offload_init()
 
         ! re-create the para_env and log with correct (reordered) ranks
         ALLOCATE (default_para_env)
         IF (init_mpi) THEN
            ! get the default system wide communicator
            CALL mp_world_init(default_para_env)
         ELSE
            CALL mp_query_thread_level()
            IF (PRESENT(mpi_comm)) THEN
               default_para_env = mpi_comm
            ELSE
               default_para_env = mp_comm_world
            END IF
         END IF
 
         CALL string_table_allocate()
         CALL add_mp_perf_env()
         CALL add_timer_env()
 
         IF (default_para_env%is_source()) THEN
            unit_nr = default_output_unit
         ELSE
            unit_nr = -1
         END IF
         NULLIFY (logger)
 
         CALL cp_logger_create(logger, para_env=default_para_env, &
                               default_global_unit_nr=unit_nr, &
                               close_global_unit_on_dealloc=.false.)
         CALL cp_add_default_logger(logger)
         CALL cp_logger_release(logger)
 
         ALLOCATE (f_envs(0))
         module_initialized = .true.
         ierr = 0
 
         ! Initialize mathematical constants
         CALL init_periodic_table()
 
         ! Init the bibliography
         CALL add_all_references()
 
         offload_device_count = offload_get_device_count()
 
         ! Select active offload device when available.
         IF (offload_device_count > 0) THEN
            offload_chosen_device = mod(default_para_env%mepos, offload_device_count)
            CALL offload_set_chosen_device(offload_chosen_device)
         END IF
 
         CALL configure_dbcsr(default_para_env, unit_nr)
#if defined(__parallel)
         CALL gauxc_init(default_para_env%get_handle(), gauxc_status)
#else
         CALL gauxc_init(status=gauxc_status)
#endif
         CALL gauxc_check_status(gauxc_status)
         CALL pw_fpga_init()
         CALL pw_gpu_init()
         CALL grid_library_init()
         CALL dbm_library_init()
         CALL dbt_library_init()
#if defined(__LIBINT)
         CALL cp_libint_static_init()
#endif
      ELSE
         ierr = cp_failure_level
      END IF
 
      !sample peak memory
      CALL m_memory()
 
   END SUBROUTINE init_cp2k
 
! **************************************************************************************************
!> \brief Configure DBCSR with CP2K's hooks and selected accelerator device.
!>        For an active library the DBCSR API only updates the output unit.
!>        Initialization belongs to init_cp2k, except for FARMING which must
!>        temporarily initialize DBCSR on each job's own MPI communicator.
!> \param mpi_comm communicator for initialization (ignored if already initialized)
!> \param output_unit diagnostic output unit, negative on non-I/O ranks
! **************************************************************************************************
   SUBROUTINE configure_dbcsr(mpi_comm, output_unit)
      CLASS(mp_comm_type), INTENT(IN)                    :: mpi_comm
      INTEGER, INTENT(IN)                                :: output_unit
 
      INTEGER, POINTER                                   :: active_device_id
      INTEGER, TARGET                                    :: chosen_device
 
      NULLIFY (active_device_id)
      IF (offload_get_device_count() > 0) THEN
         chosen_device = offload_get_chosen_device()
         active_device_id => chosen_device
      END IF
      CALL dbcsr_init_lib(mpi_comm%get_handle(), timeset_hook, timestop_hook, &
                          cp_abort_hook, cp_warn_hook, io_unit=output_unit, &
                          accdrv_active_device_id=active_device_id)
   END SUBROUTINE configure_dbcsr
 
! **************************************************************************************************
!> \brief cleanup after you have finished using this interface
!> \param finalize_mpi if the mpi environment should be finalized
!> \param ierr returns a number different from 0 if there was an error
!> \author fawzi
! **************************************************************************************************
   SUBROUTINE finalize_cp2k(finalize_mpi, ierr)
      LOGICAL, INTENT(in)                                :: finalize_mpi
      INTEGER, INTENT(out)                               :: ierr
 
      INTEGER                                            :: ienv
      TYPE(cp_gauxc_status_type)                         :: gauxc_status
 
!sample peak memory
 
      CALL m_memory()
 
      IF (.NOT. module_initialized) THEN
         ierr = cp_failure_level
      ELSE
         DO ienv = n_f_envs, 1, -1
            CALL destroy_force_env(f_envs(ienv)%f_env%id_nr, ierr=ierr)
            cpassert(ierr == 0)
         END DO
         DEALLOCATE (f_envs)
 
#if defined(__LIBINT)
         CALL cp_libint_static_cleanup()
#endif
 
         ! Release the cached LibXC workers, because their device buffers are
         ! owned by the offload memory pool which is cleared below
         CALL libxc_release_workers()
 
         ! Finalize libraries (Offload)
         CALL dbt_library_finalize()
         CALL dbm_library_finalize()
         CALL grid_library_finalize()
         CALL pw_gpu_finalize()
         CALL pw_fpga_finalize()
         IF (cp_sirius_is_initialized()) CALL cp_sirius_finalize()
         CALL gauxc_finalize(gauxc_status)
         CALL gauxc_check_status(gauxc_status)
         ! Finalize the DBCSR library
         CALL dbcsr_finalize_lib()
 
         ! Finalize DLA-Future and pika runtime; if already finalized does nothing
         CALL cp_dlaf_free_all_grids()
         CALL cp_dlaf_finalize()
 
         CALL mp_para_env_release(default_para_env)
         CALL cp_rm_default_logger()
 
         ! Deallocate the bibliography
         CALL remove_all_references()
         CALL rm_timer_env()
         CALL rm_mp_perf_env()
         CALL string_table_deallocate(0)
         CALL cp_openpmd_output_finalize()
         IF (finalize_mpi) THEN
            CALL mp_world_finalize()
         END IF
 
         ierr = 0
      END IF
   END SUBROUTINE finalize_cp2k
 
! **************************************************************************************************
!> \brief deallocates a f_env
!> \param f_env the f_env to deallocate
!> \author fawzi
! **************************************************************************************************
   RECURSIVE SUBROUTINE f_env_dealloc(f_env)
      TYPE(f_env_type), POINTER                          :: f_env
 
      INTEGER                                            :: ierr
 
      cpassert(ASSOCIATED(f_env))
      CALL force_env_release(f_env%force_env)
      CALL cp_logger_release(f_env%logger)
      CALL timer_env_release(f_env%timer_env)
      CALL mp_perf_env_release(f_env%mp_perf_env)
      IF (f_env%old_path /= f_env%my_path) THEN
         CALL m_chdir(f_env%old_path, ierr)
         cpassert(ierr == 0)
      END IF
   END SUBROUTINE f_env_dealloc
 
! **************************************************************************************************
!> \brief createates a f_env
!> \param f_env the f_env to createate
!> \param force_env the force_environment to be stored
!> \param timer_env the timer env to be stored
!> \param mp_perf_env the mp performance environment to be stored
!> \param id_nr ...
!> \param logger ...
!> \param old_dir ...
!> \author fawzi
! **************************************************************************************************
   SUBROUTINE f_env_create(f_env, force_env, timer_env, mp_perf_env, id_nr, logger, old_dir)
      TYPE(f_env_type), POINTER                          :: f_env
      TYPE(force_env_type), POINTER                      :: force_env
      TYPE(timer_env_type), POINTER                      :: timer_env
      TYPE(mp_perf_env_type), POINTER                    :: mp_perf_env
      INTEGER, INTENT(in)                                :: id_nr
      TYPE(cp_logger_type), POINTER                      :: logger
      CHARACTER(len=*), INTENT(in)                       :: old_dir
 
      ALLOCATE (f_env)
      f_env%force_env => force_env
      CALL force_env_retain(f_env%force_env)
      f_env%logger => logger
      CALL cp_logger_retain(logger)
      f_env%timer_env => timer_env
      CALL timer_env_retain(f_env%timer_env)
      f_env%mp_perf_env => mp_perf_env
      CALL mp_perf_env_retain(f_env%mp_perf_env)
      f_env%id_nr = id_nr
      CALL m_getcwd(f_env%my_path)
      f_env%old_path = old_dir
   END SUBROUTINE f_env_create
 
! **************************************************************************************************
!> \brief ...
!> \param f_env_id ...
!> \param f_env ...
! **************************************************************************************************
   SUBROUTINE f_env_get_from_id(f_env_id, f_env)
      INTEGER, INTENT(in)                                :: f_env_id
      TYPE(f_env_type), POINTER                          :: f_env
 
      INTEGER                                            :: f_env_pos
 
      NULLIFY (f_env)
      f_env_pos = get_pos_of_env(f_env_id)
      IF (f_env_pos < 1) THEN
         cpabort("invalid env_id "//cp_to_string(f_env_id))
      ELSE
         f_env => f_envs(f_env_pos)%f_env
      END IF
 
   END SUBROUTINE f_env_get_from_id
 
! **************************************************************************************************
!> \brief adds the default environments of the f_env to the stack of the
!>      defaults, and returns a new error and sets failure to true if
!>      something went wrong
!> \param f_env_id the f_env from where to take the defaults
!> \param f_env will contain the f_env corresponding to f_env_id
!> \param handle ...
!> \author fawzi
!> \note
!>      The following routines need to be synchronized wrt. adding/removing
!>      of the default environments (logging, performance,error):
!>      environment:cp2k_init, environment:cp2k_finalize,
!>      f77_interface:f_env_add_defaults, f77_interface:f_env_rm_defaults,
!>      f77_interface:create_force_env, f77_interface:destroy_force_env
! **************************************************************************************************
   SUBROUTINE f_env_add_defaults(f_env_id, f_env, handle)
      INTEGER, INTENT(in)                                :: f_env_id
      TYPE(f_env_type), POINTER                          :: f_env
      INTEGER, INTENT(out), OPTIONAL                     :: handle
 
      INTEGER                                            :: f_env_pos, ierr
      TYPE(cp_logger_type), POINTER                      :: logger
 
      NULLIFY (f_env)
      f_env_pos = get_pos_of_env(f_env_id)
      IF (f_env_pos < 1) THEN
         cpabort("invalid env_id "//cp_to_string(f_env_id))
      ELSE
         f_env => f_envs(f_env_pos)%f_env
         logger => f_env%logger
         cpassert(ASSOCIATED(logger))
         CALL m_getcwd(f_env%old_path)
         IF (f_env%old_path /= f_env%my_path) THEN
            CALL m_chdir(trim(f_env%my_path), ierr)
            cpassert(ierr == 0)
         END IF
         CALL add_mp_perf_env(f_env%mp_perf_env)
         CALL add_timer_env(f_env%timer_env)
         CALL cp_add_default_logger(logger)
         IF (PRESENT(handle)) handle = cp_default_logger_stack_size()
      END IF
   END SUBROUTINE f_env_add_defaults
 
! **************************************************************************************************
!> \brief removes the default environments of the f_env to the stack of the
!>      defaults, and sets ierr accordingly to the failuers stored in error
!>      It also releases the error
!> \param f_env the f_env from where to take the defaults
!> \param ierr variable that will be set to a number different from 0 if
!>        error contains an error (otherwise it will be set to 0)
!> \param handle ...
!> \author fawzi
!> \note
!>      The following routines need to be synchronized wrt. adding/removing
!>      of the default environments (logging, performance,error):
!>      environment:cp2k_init, environment:cp2k_finalize,
!>      f77_interface:f_env_add_defaults, f77_interface:f_env_rm_defaults,
!>      f77_interface:create_force_env, f77_interface:destroy_force_env
! **************************************************************************************************
   SUBROUTINE f_env_rm_defaults(f_env, ierr, handle)
      TYPE(f_env_type), POINTER                          :: f_env
      INTEGER, INTENT(out), OPTIONAL                     :: ierr
      INTEGER, INTENT(in), OPTIONAL                      :: handle
 
      INTEGER                                            :: ierr2
      TYPE(cp_logger_type), POINTER                      :: d_logger, logger
      TYPE(mp_perf_env_type), POINTER                    :: d_mp_perf_env
      TYPE(timer_env_type), POINTER                      :: d_timer_env
 
      IF (ASSOCIATED(f_env)) THEN
         IF (PRESENT(handle)) THEN
            cpassert(handle == cp_default_logger_stack_size())
         END IF
 
         logger => f_env%logger
         d_logger => cp_get_default_logger()
         d_timer_env => get_timer_env()
         d_mp_perf_env => get_mp_perf_env()
         cpassert(ASSOCIATED(logger))
         cpassert(ASSOCIATED(d_logger))
         cpassert(ASSOCIATED(d_timer_env))
         cpassert(ASSOCIATED(d_mp_perf_env))
         cpassert(ASSOCIATED(logger, d_logger))
         ! CPASSERT(ASSOCIATED(d_timer_env, f_env%timer_env))
         cpassert(ASSOCIATED(d_mp_perf_env, f_env%mp_perf_env))
         IF (f_env%old_path /= f_env%my_path) THEN
            CALL m_chdir(trim(f_env%old_path), ierr2)
            cpassert(ierr2 == 0)
         END IF
         IF (PRESENT(ierr)) THEN
            ierr = 0
         END IF
         CALL cp_rm_default_logger()
         CALL rm_timer_env()
         CALL rm_mp_perf_env()
      ELSE
         IF (PRESENT(ierr)) THEN
            ierr = 0
         END IF
      END IF
   END SUBROUTINE f_env_rm_defaults
 
! **************************************************************************************************
!> \brief creates a new force environment using the given input, and writing
!>        the output to the given output unit
!> \param new_env_id will contain the id of the newly created environment
!> \param input_declaration ...
!> \param input_path where to read the input (if the input is given it can
!>        a virtual path)
!> \param output_path filename (or name of the unit) for the output
!> \param mpi_comm the mpi communicator to be used for this environment
!>        it will not be freed when you get rid of the force_env
!> \param output_unit if given it should be the unit for the output
!>        and no file is open (should be valid on the processor with rank 0)
!> \param owns_out_unit if the output unit should be closed upon destroing
!>        of the force_env (defaults to true if not default_output_unit)
!> \param input the parsed input, if given and valid it is used
!>        instead of parsing from file
!> \param ierr will return a number different from 0 if there was an error
!> \param work_dir ...
!> \param initial_variables key-value list of initial preprocessor variables
!> \author fawzi
!> \note
!>      The following routines need to be synchronized wrt. adding/removing
!>      of the default environments (logging, performance,error):
!>      environment:cp2k_init, environment:cp2k_finalize,
!>      f77_interface:f_env_add_defaults, f77_interface:f_env_rm_defaults,
!>      f77_interface:create_force_env, f77_interface:destroy_force_env
! **************************************************************************************************
   RECURSIVE SUBROUTINE create_force_env(new_env_id, input_declaration, input_path, &
                                         output_path, mpi_comm, output_unit, owns_out_unit, &
                                         input, ierr, work_dir, initial_variables)
      INTEGER, INTENT(out)                               :: new_env_id
      TYPE(section_type), POINTER                        :: input_declaration
      CHARACTER(len=*), INTENT(in)                       :: input_path
      CHARACTER(len=*), INTENT(in), OPTIONAL             :: output_path
 
      CLASS(mp_comm_type), INTENT(IN), OPTIONAL          :: mpi_comm
      INTEGER, INTENT(in), OPTIONAL                      :: output_unit
      LOGICAL, INTENT(in), OPTIONAL                      :: owns_out_unit
      TYPE(section_vals_type), OPTIONAL, POINTER         :: input
      INTEGER, INTENT(out), OPTIONAL                     :: ierr
      CHARACTER(len=*), INTENT(in), OPTIONAL             :: work_dir
      CHARACTER(len=*), DIMENSION(:, :), OPTIONAL        :: initial_variables
 
      CHARACTER(len=*), PARAMETER                        :: routinen = 'create_force_env'
 
      CHARACTER(len=default_path_length)                 :: old_dir, wdir
      INTEGER :: handle, i, ierr2, iforce_eval, isubforce_eval, k, method_name_id, my_group, &
                 nforce_eval, ngroups, nsubforce_size, unit_nr
      INTEGER, DIMENSION(:), POINTER                     :: group_distribution, i_force_eval, &
                                                            lgroup_distribution
      LOGICAL :: check, do_qmmm_force_mixing, multiple_subsys, my_owns_out_unit, &
                 use_motion_section, use_multiple_para_env
      TYPE(cp_logger_type), POINTER                      :: logger, my_logger
      TYPE(mp_para_env_type), POINTER                    :: my_para_env, para_env
      TYPE(eip_environment_type), POINTER                :: eip_env
      TYPE(embed_env_type), POINTER                      :: embed_env
      TYPE(enumeration_type), POINTER                    :: enum
      TYPE(f_env_p_type), DIMENSION(:), POINTER          :: f_envs_old
      TYPE(force_env_type), POINTER                      :: force_env, my_force_env
      TYPE(fp_type), POINTER                             :: fp_env
      TYPE(global_environment_type), POINTER             :: globenv
      TYPE(ipi_environment_type), POINTER                :: ipi_env
      TYPE(keyword_type), POINTER                        :: keyword
      TYPE(meta_env_type), POINTER                       :: meta_env
      TYPE(mixed_environment_type), POINTER              :: mixed_env
      TYPE(mp_perf_env_type), POINTER                    :: mp_perf_env
      TYPE(nnp_type), POINTER                            :: nnp_env
      TYPE(pwdft_environment_type), POINTER              :: pwdft_env
      TYPE(qmmm_env_type), POINTER                       :: qmmm_env
      TYPE(qmmmx_env_type), POINTER                      :: qmmmx_env
      TYPE(qs_environment_type), POINTER                 :: qs_env
      TYPE(section_type), POINTER                        :: section
      TYPE(section_vals_type), POINTER :: fe_section, force_env_section, force_env_sections, &
                                          fp_section, input_file, qmmm_section, qmmmx_section, root_section, subsys_section, &
                                          wrk_section
      TYPE(timer_env_type), POINTER                      :: timer_env
 
      cpassert(ASSOCIATED(input_declaration))
      NULLIFY (para_env, force_env, timer_env, mp_perf_env, globenv, meta_env, &
               fp_env, eip_env, pwdft_env, mixed_env, qs_env, qmmm_env, embed_env)
      new_env_id = -1
      IF (PRESENT(mpi_comm)) THEN
         ALLOCATE (para_env)
         para_env = mpi_comm
      ELSE
         para_env => default_para_env
         CALL para_env%retain()
      END IF
 
      CALL timeset(routinen, handle)
 
      CALL m_getcwd(old_dir)
      wdir = old_dir
      IF (PRESENT(work_dir)) THEN
         IF (work_dir /= " ") THEN
            CALL m_chdir(work_dir, ierr2)
            IF (ierr2 /= 0) THEN
               IF (PRESENT(ierr)) ierr = ierr2
               RETURN
            END IF
            wdir = work_dir
         END IF
      END IF
 
      IF (PRESENT(output_unit)) THEN
         unit_nr = output_unit
      ELSE
         IF (para_env%is_source()) THEN
            IF (output_path == "__STD_OUT__") THEN
               unit_nr = default_output_unit
            ELSE
               CALL open_file(file_name=output_path, &
                              file_status="UNKNOWN", &
                              file_action="WRITE", &
                              file_position="APPEND", &
                              unit_number=unit_nr)
            END IF
         ELSE
            unit_nr = -1
         END IF
      END IF
 
      my_owns_out_unit = unit_nr /= default_output_unit
      IF (PRESENT(owns_out_unit)) my_owns_out_unit = owns_out_unit
      CALL globenv_create(globenv)
      CALL cp2k_init(para_env, output_unit=unit_nr, globenv=globenv, input_file_name=input_path, &
                     wdir=wdir)
      logger => cp_get_default_logger()
      ! warning this is dangerous, I did not check that all the subfunctions
      ! support it, the program might crash upon error
 
      NULLIFY (input_file)
      IF (PRESENT(input)) input_file => input
      IF (.NOT. ASSOCIATED(input_file)) THEN
         IF (PRESENT(initial_variables)) THEN
            input_file => read_input(input_declaration, input_path, initial_variables, para_env=para_env)
         ELSE
            input_file => read_input(input_declaration, input_path, empty_initial_variables, para_env=para_env)
         END IF
      ELSE
         CALL section_vals_retain(input_file)
      END IF
 
      CALL check_cp2k_input(input_declaration, input_file, para_env=para_env, output_unit=unit_nr)
 
      root_section => input_file
      CALL section_vals_retain(root_section)
 
      IF (n_f_envs + 1 > SIZE(f_envs)) THEN
         f_envs_old => f_envs
         ALLOCATE (f_envs(n_f_envs + 10))
         DO i = 1, n_f_envs
            f_envs(i)%f_env => f_envs_old(i)%f_env
         END DO
         DO i = n_f_envs + 1, SIZE(f_envs)
            NULLIFY (f_envs(i)%f_env)
         END DO
         DEALLOCATE (f_envs_old)
      END IF
 
      CALL cp2k_read(root_section, para_env, globenv)
 
      CALL cp2k_setup(root_section, para_env, globenv)
      ! Group Distribution
      ALLOCATE (group_distribution(0:para_env%num_pe - 1))
      group_distribution = 0
      lgroup_distribution => group_distribution
      ! Setup all possible force_env
      force_env_sections => section_vals_get_subs_vals(root_section, "FORCE_EVAL")
      CALL section_vals_val_get(root_section, "MULTIPLE_FORCE_EVALS%MULTIPLE_SUBSYS", &
                                l_val=multiple_subsys)
      CALL multiple_fe_list(force_env_sections, root_section, i_force_eval, nforce_eval)
      ! Enforce the deletion of the subsys (unless not explicitly required)
      IF (.NOT. multiple_subsys) THEN
         DO iforce_eval = 2, nforce_eval
            wrk_section => section_vals_get_subs_vals(force_env_sections, "SUBSYS", &
                                                      i_rep_section=i_force_eval(iforce_eval))
            CALL section_vals_remove_values(wrk_section)
         END DO
      END IF
      nsubforce_size = nforce_eval - 1
      use_multiple_para_env = .false.
      use_motion_section = .true.
      DO iforce_eval = 1, nforce_eval
         NULLIFY (force_env_section, my_force_env, subsys_section)
         ! Reference subsys from the first ordered force_eval
         IF (.NOT. multiple_subsys) THEN
            subsys_section => section_vals_get_subs_vals(force_env_sections, "SUBSYS", &
                                                         i_rep_section=i_force_eval(1))
         END IF
         ! Handling para_env in case of multiple force_eval
         IF (use_multiple_para_env) THEN
            ! Check that the order of the force_eval is the correct one
            CALL section_vals_val_get(force_env_sections, "METHOD", i_val=method_name_id, &
                                      i_rep_section=i_force_eval(1))
            IF ((method_name_id /= do_mixed) .AND. (method_name_id /= do_embed)) THEN
               CALL cp_abort(__location__, &
                             "In case of multiple force_eval the MAIN force_eval (the first in the list of FORCE_EVAL_ORDER or "// &
                             "the one omitted from that order list) must be a MIXED_ENV type calculation. Please check your "// &
                             "input file and possibly correct the MULTIPLE_FORCE_EVAL%FORCE_EVAL_ORDER. ")
            END IF
 
            IF (method_name_id == do_mixed) THEN
               check = ASSOCIATED(force_env%mixed_env%sub_para_env)
               cpassert(check)
               ngroups = force_env%mixed_env%ngroups
               my_group = lgroup_distribution(para_env%mepos)
               isubforce_eval = iforce_eval - 1
               ! If task not allocated on this procs skip setup..
               IF (modulo(isubforce_eval - 1, ngroups) /= my_group) cycle
               my_para_env => force_env%mixed_env%sub_para_env(my_group + 1)%para_env
               my_logger => force_env%mixed_env%sub_logger(my_group + 1)%p
               CALL cp_rm_default_logger()
               CALL cp_add_default_logger(my_logger)
            END IF
            IF (method_name_id == do_embed) THEN
               check = ASSOCIATED(force_env%embed_env%sub_para_env)
               cpassert(check)
               ngroups = force_env%embed_env%ngroups
               my_group = lgroup_distribution(para_env%mepos)
               isubforce_eval = iforce_eval - 1
               ! If task not allocated on this procs skip setup..
               IF (modulo(isubforce_eval - 1, ngroups) /= my_group) cycle
               my_para_env => force_env%embed_env%sub_para_env(my_group + 1)%para_env
               my_logger => force_env%embed_env%sub_logger(my_group + 1)%p
               CALL cp_rm_default_logger()
               CALL cp_add_default_logger(my_logger)
            END IF
         ELSE
            my_para_env => para_env
         END IF
 
         ! Initialize force_env_section
         ! No need to allocate one more force_env_section if only 1 force_eval
         ! is provided.. this is in order to save memory..
         IF (nforce_eval > 1) THEN
            CALL section_vals_duplicate(force_env_sections, force_env_section, &
                                        i_force_eval(iforce_eval), i_force_eval(iforce_eval))
            IF (iforce_eval /= 1) use_motion_section = .false.
         ELSE
            force_env_section => force_env_sections
            use_motion_section = .true.
         END IF
         CALL section_vals_val_get(force_env_section, "METHOD", i_val=method_name_id)
 
         IF (method_name_id == do_qmmm) THEN
            qmmmx_section => section_vals_get_subs_vals(force_env_section, "QMMM%FORCE_MIXING")
            CALL section_vals_get(qmmmx_section, explicit=do_qmmm_force_mixing)
            IF (do_qmmm_force_mixing) THEN
               method_name_id = do_qmmmx
            END IF ! QMMM Force-Mixing has its own (hidden) method_id
         END IF
 
         SELECT CASE (method_name_id)
         CASE (do_fist)
            CALL fist_create_force_env(my_force_env, root_section, my_para_env, globenv, &
                                       force_env_section=force_env_section, subsys_section=subsys_section, &
                                       use_motion_section=use_motion_section)
 
         CASE (do_qs)
            ALLOCATE (qs_env)
            CALL qs_env_create(qs_env, globenv)
            CALL qs_init(qs_env, my_para_env, root_section, globenv=globenv, force_env_section=force_env_section, &
                         subsys_section=subsys_section, use_motion_section=use_motion_section)
            CALL force_env_create(my_force_env, root_section, qs_env=qs_env, para_env=my_para_env, globenv=globenv, &
                                  force_env_section=force_env_section)
 
         CASE (do_qmmm)
            qmmm_section => section_vals_get_subs_vals(force_env_section, "QMMM")
            ALLOCATE (qmmm_env)
            CALL qmmm_env_create(qmmm_env, root_section, my_para_env, globenv, &
                                 force_env_section, qmmm_section, subsys_section, use_motion_section)
            CALL force_env_create(my_force_env, root_section, qmmm_env=qmmm_env, para_env=my_para_env, &
                                  globenv=globenv, force_env_section=force_env_section)
 
         CASE (do_qmmmx)
            ALLOCATE (qmmmx_env)
            CALL qmmmx_env_create(qmmmx_env, root_section, my_para_env, globenv, &
                                  force_env_section, subsys_section, use_motion_section)
            CALL force_env_create(my_force_env, root_section, qmmmx_env=qmmmx_env, para_env=my_para_env, &
                                  globenv=globenv, force_env_section=force_env_section)
 
         CASE (do_eip)
            ALLOCATE (eip_env)
            CALL eip_env_create(eip_env)
            CALL eip_init(eip_env, root_section, my_para_env, force_env_section=force_env_section, &
                          subsys_section=subsys_section)
            CALL force_env_create(my_force_env, root_section, eip_env=eip_env, para_env=my_para_env, &
                                  globenv=globenv, force_env_section=force_env_section)
 
         CASE (do_sirius)
            IF (.NOT. cp_sirius_is_initialized()) THEN
               IF (unit_nr > 0) WRITE (unit=unit_nr, fmt="(T2,A)", advance="NO") "SIRIUS| "
               CALL cp_sirius_init()
            END IF
            ALLOCATE (pwdft_env)
            CALL pwdft_env_create(pwdft_env)
            CALL pwdft_init(pwdft_env, root_section, my_para_env, force_env_section=force_env_section, &
                            subsys_section=subsys_section, use_motion_section=use_motion_section)
            CALL force_env_create(my_force_env, root_section, pwdft_env=pwdft_env, para_env=my_para_env, &
                                  globenv=globenv, force_env_section=force_env_section)
 
         CASE (do_mixed)
            ALLOCATE (mixed_env)
            CALL mixed_create_force_env(mixed_env, root_section, my_para_env, &
                                        force_env_section=force_env_section, n_subforce_eval=nsubforce_size, &
                                        use_motion_section=use_motion_section)
            CALL force_env_create(my_force_env, root_section, mixed_env=mixed_env, para_env=my_para_env, &
                                  globenv=globenv, force_env_section=force_env_section)
            !TODO: the sub_force_envs should really be created via recursion
            use_multiple_para_env = .true.
            CALL cp_add_default_logger(logger) ! just to get the logger swapping started
            lgroup_distribution => my_force_env%mixed_env%group_distribution
 
         CASE (do_embed)
            ALLOCATE (embed_env)
            CALL embed_create_force_env(embed_env, root_section, my_para_env, &
                                        force_env_section=force_env_section, n_subforce_eval=nsubforce_size, &
                                        use_motion_section=use_motion_section)
            CALL force_env_create(my_force_env, root_section, embed_env=embed_env, para_env=my_para_env, &
                                  globenv=globenv, force_env_section=force_env_section)
            !TODO: the sub_force_envs should really be created via recursion
            use_multiple_para_env = .true.
            CALL cp_add_default_logger(logger) ! just to get the logger swapping started
            lgroup_distribution => my_force_env%embed_env%group_distribution
 
         CASE (do_nnp)
            ALLOCATE (nnp_env)
            CALL nnp_init(nnp_env, root_section, my_para_env, force_env_section=force_env_section, &
                          subsys_section=subsys_section, use_motion_section=use_motion_section)
            CALL force_env_create(my_force_env, root_section, nnp_env=nnp_env, para_env=my_para_env, &
                                  globenv=globenv, force_env_section=force_env_section)
 
         CASE (do_ipi)
            ALLOCATE (ipi_env)
            CALL ipi_init(ipi_env, root_section, my_para_env, force_env_section=force_env_section, &
                          subsys_section=subsys_section)
            CALL force_env_create(my_force_env, root_section, ipi_env=ipi_env, para_env=my_para_env, &
                                  globenv=globenv, force_env_section=force_env_section)
 
         CASE DEFAULT
            CALL create_force_eval_section(section)
            keyword => section_get_keyword(section, "METHOD")
            CALL keyword_get(keyword, enum=enum)
            CALL cp_abort(__location__, &
                          "Invalid METHOD <"//trim(enum_i2c(enum, method_name_id))// &
                          "> was specified, ")
            CALL section_release(section)
         END SELECT
 
         NULLIFY (meta_env, fp_env)
         IF (use_motion_section) THEN
            ! Metadynamics Setup
            fe_section => section_vals_get_subs_vals(root_section, "MOTION%FREE_ENERGY")
            CALL metadyn_read(meta_env, my_force_env, root_section, my_para_env, fe_section)
            CALL force_env_set(my_force_env, meta_env=meta_env)
            ! Flexible Partition Setup
            fp_section => section_vals_get_subs_vals(root_section, "MOTION%FLEXIBLE_PARTITIONING")
            ALLOCATE (fp_env)
            CALL fp_env_create(fp_env)
            CALL fp_env_read(fp_env, fp_section)
            CALL fp_env_write(fp_env, fp_section)
            CALL force_env_set(my_force_env, fp_env=fp_env)
         END IF
         ! Handle multiple force_eval
         IF (nforce_eval > 1 .AND. iforce_eval == 1) THEN
            ALLOCATE (my_force_env%sub_force_env(nsubforce_size))
            ! Nullify subforce_env
            DO k = 1, nsubforce_size
               NULLIFY (my_force_env%sub_force_env(k)%force_env)
            END DO
         END IF
         ! Reference the right force_env
         IF (iforce_eval == 1) THEN
            force_env => my_force_env
         ELSE
            force_env%sub_force_env(iforce_eval - 1)%force_env => my_force_env
         END IF
         ! Multiple para env for sub_force_eval
         IF (.NOT. use_multiple_para_env) THEN
            lgroup_distribution = iforce_eval
         END IF
         ! Release force_env_section
         IF (nforce_eval > 1) CALL section_vals_release(force_env_section)
      END DO
      IF (use_multiple_para_env) THEN
         CALL cp_rm_default_logger()
      END IF
      DEALLOCATE (group_distribution)
      DEALLOCATE (i_force_eval)
      timer_env => get_timer_env()
      mp_perf_env => get_mp_perf_env()
      CALL para_env%max(last_f_env_id)
      last_f_env_id = last_f_env_id + 1
      new_env_id = last_f_env_id
      n_f_envs = n_f_envs + 1
      CALL f_env_create(f_envs(n_f_envs)%f_env, logger=logger, &
                        timer_env=timer_env, mp_perf_env=mp_perf_env, force_env=force_env, &
                        id_nr=last_f_env_id, old_dir=old_dir)
      CALL force_env_release(force_env)
      CALL globenv_release(globenv)
      CALL section_vals_release(root_section)
      CALL mp_para_env_release(para_env)
      CALL f_env_rm_defaults(f_envs(n_f_envs)%f_env, ierr=ierr)
      CALL timestop(handle)
 
   END SUBROUTINE create_force_env
 
! **************************************************************************************************
!> \brief deallocates the force_env with the given id
!> \param env_id the id of the force_env to remove
!> \param ierr will contain a number different from 0 if
!> \param q_finalize ...
!> \author fawzi
!> \note
!>      The following routines need to be synchronized wrt. adding/removing
!>      of the default environments (logging, performance,error):
!>      environment:cp2k_init, environment:cp2k_finalize,
!>      f77_interface:f_env_add_defaults, f77_interface:f_env_rm_defaults,
!>      f77_interface:create_force_env, f77_interface:destroy_force_env
! **************************************************************************************************
   RECURSIVE SUBROUTINE destroy_force_env(env_id, ierr, q_finalize)
      INTEGER, INTENT(in)                                :: env_id
      INTEGER, INTENT(out)                               :: ierr
      LOGICAL, INTENT(IN), OPTIONAL                      :: q_finalize
 
      INTEGER                                            :: env_pos, i
      TYPE(f_env_type), POINTER                          :: f_env
      TYPE(global_environment_type), POINTER             :: globenv
      TYPE(mp_para_env_type), POINTER                    :: para_env
      TYPE(section_vals_type), POINTER                   :: root_section
 
      NULLIFY (f_env)
      CALL f_env_add_defaults(env_id, f_env)
      env_pos = get_pos_of_env(env_id)
      n_f_envs = n_f_envs - 1
      DO i = env_pos, n_f_envs
         f_envs(i)%f_env => f_envs(i + 1)%f_env
      END DO
      NULLIFY (f_envs(n_f_envs + 1)%f_env)
 
      CALL force_env_get(f_env%force_env, globenv=globenv, &
                         root_section=root_section, para_env=para_env)
 
      cpassert(ASSOCIATED(globenv))
      NULLIFY (f_env%force_env%globenv)
      CALL f_env_dealloc(f_env)
      IF (PRESENT(q_finalize)) THEN
         CALL cp2k_finalize(root_section, para_env, globenv, f_env%old_path, q_finalize)
      ELSE
         CALL cp2k_finalize(root_section, para_env, globenv, f_env%old_path)
      END IF
      CALL section_vals_release(root_section)
      CALL globenv_release(globenv)
      DEALLOCATE (f_env)
      ierr = 0
   END SUBROUTINE destroy_force_env
 
! **************************************************************************************************
!> \brief returns the number of atoms in the given force env
!> \param env_id id of the force_env
!> \param n_atom ...
!> \param ierr will return a number different from 0 if there was an error
!> \date   22.11.2010 (MK)
!> \author fawzi
! **************************************************************************************************
   SUBROUTINE get_natom(env_id, n_atom, ierr)
 
      INTEGER, INTENT(IN)                                :: env_id
      INTEGER, INTENT(OUT)                               :: n_atom, ierr
 
      TYPE(f_env_type), POINTER                          :: f_env
 
      n_atom = 0
      NULLIFY (f_env)
      CALL f_env_add_defaults(env_id, f_env)
      n_atom = force_env_get_natom(f_env%force_env)
      CALL f_env_rm_defaults(f_env, ierr)
 
   END SUBROUTINE get_natom
 
! **************************************************************************************************
!> \brief returns the number of particles in the given force env
!> \param env_id id of the force_env
!> \param n_particle ...
!> \param ierr will return a number different from 0 if there was an error
!> \author Matthias Krack
!>
! **************************************************************************************************
   SUBROUTINE get_nparticle(env_id, n_particle, ierr)
 
      INTEGER, INTENT(IN)                                :: env_id
      INTEGER, INTENT(OUT)                               :: n_particle, ierr
 
      TYPE(f_env_type), POINTER                          :: f_env
 
      n_particle = 0
      NULLIFY (f_env)
      CALL f_env_add_defaults(env_id, f_env)
      n_particle = force_env_get_nparticle(f_env%force_env)
      CALL f_env_rm_defaults(f_env, ierr)
 
   END SUBROUTINE get_nparticle
 
! **************************************************************************************************
!> \brief gets a cell
!> \param env_id id of the force_env
!> \param cell the array with the cell matrix
!> \param per periodicity
!> \param ierr will return a number different from 0 if there was an error
!> \author Joost VandeVondele
! **************************************************************************************************
   SUBROUTINE get_cell(env_id, cell, per, ierr)
 
      INTEGER, INTENT(IN)                                :: env_id
      REAL(kind=dp), DIMENSION(3, 3)                     :: cell
      INTEGER, DIMENSION(3), OPTIONAL                    :: per
      INTEGER, INTENT(OUT)                               :: ierr
 
      TYPE(cell_type), POINTER                           :: cell_full
      TYPE(f_env_type), POINTER                          :: f_env
 
      NULLIFY (f_env)
      CALL f_env_add_defaults(env_id, f_env)
      NULLIFY (cell_full)
      CALL force_env_get(f_env%force_env, cell=cell_full)
      cpassert(ASSOCIATED(cell_full))
      cell = cell_full%hmat
      IF (PRESENT(per)) per(:) = cell_full%perd(:)
      CALL f_env_rm_defaults(f_env, ierr)
 
   END SUBROUTINE get_cell
 
! **************************************************************************************************
!> \brief gets the qmmm cell
!> \param env_id id of the force_env
!> \param cell the array with the cell matrix
!> \param ierr will return a number different from 0 if there was an error
!> \author Holly Judge
! **************************************************************************************************
   SUBROUTINE get_qmmm_cell(env_id, cell, ierr)
 
      INTEGER, INTENT(IN)                                :: env_id
      REAL(kind=dp), DIMENSION(3, 3)                     :: cell
      INTEGER, INTENT(OUT)                               :: ierr
 
      TYPE(cell_type), POINTER                           :: cell_qmmm
      TYPE(f_env_type), POINTER                          :: f_env
      TYPE(qmmm_env_type), POINTER                       :: qmmm_env
 
      NULLIFY (f_env)
      CALL f_env_add_defaults(env_id, f_env)
      NULLIFY (cell_qmmm)
      CALL force_env_get(f_env%force_env, qmmm_env=qmmm_env)
      CALL get_qs_env(qmmm_env%qs_env, cell=cell_qmmm)
      cpassert(ASSOCIATED(cell_qmmm))
      cell = cell_qmmm%hmat
      CALL f_env_rm_defaults(f_env, ierr)
 
   END SUBROUTINE get_qmmm_cell
 
! **************************************************************************************************
!> \brief gets a result from CP2K that is a real 1D array
!> \param env_id id of the force_env
!> \param description the tag of the result
!> \param N ...
!> \param RESULT ...
!> \param res_exist ...
!> \param ierr will return a number different from 0 if there was an error
!> \author Joost VandeVondele
! **************************************************************************************************
   SUBROUTINE get_result_r1(env_id, description, N, RESULT, res_exist, ierr)
      INTEGER                                            :: env_id
      CHARACTER(LEN=default_string_length)               :: description
      INTEGER                                            :: n
      REAL(kind=dp), DIMENSION(1:N)                      :: result
      LOGICAL, OPTIONAL                                  :: res_exist
      INTEGER                                            :: ierr
 
      INTEGER                                            :: nres
      LOGICAL                                            :: exist_res
      TYPE(cp_result_type), POINTER                      :: results
      TYPE(cp_subsys_type), POINTER                      :: subsys
      TYPE(f_env_type), POINTER                          :: f_env
 
      NULLIFY (f_env, subsys, results)
      CALL f_env_add_defaults(env_id, f_env)
 
      CALL force_env_get(f_env%force_env, subsys=subsys)
      CALL cp_subsys_get(subsys, results=results)
      ! first test for the result
      IF (PRESENT(res_exist)) THEN
         res_exist = test_for_result(results, description=description)
         exist_res = res_exist
      ELSE
         exist_res = .true.
      END IF
      ! if existing (or assuming the existence) read the results
      IF (exist_res) THEN
         CALL get_results(results, description=description, n_rep=nres)
         CALL get_results(results, description=description, values=result, nval=nres)
      END IF
 
      CALL f_env_rm_defaults(f_env, ierr)
 
   END SUBROUTINE get_result_r1
 
! **************************************************************************************************
!> \brief gets the forces of the particles
!> \param env_id id of the force_env
!> \param frc the array where to write the forces
!> \param n_el number of positions (3*nparticle) just to check
!> \param ierr will return a number different from 0 if there was an error
!> \date   22.11.2010 (MK)
!> \author fawzi
! **************************************************************************************************
   SUBROUTINE get_force(env_id, frc, n_el, ierr)
 
      INTEGER, INTENT(IN)                                :: env_id, n_el
      REAL(kind=dp), DIMENSION(1:n_el)                   :: frc
      INTEGER, INTENT(OUT)                               :: ierr
 
      TYPE(f_env_type), POINTER                          :: f_env
 
      NULLIFY (f_env)
      CALL f_env_add_defaults(env_id, f_env)
      CALL force_env_get_frc(f_env%force_env, frc, n_el)
      CALL f_env_rm_defaults(f_env, ierr)
 
   END SUBROUTINE get_force
 
! **************************************************************************************************
!> \brief gets the stress tensor
!> \param env_id id of the force_env
!> \param stress_tensor the array where to write the stress tensor
!> \param ierr will return a number different from 0 if there was an error
!> \param available whether the potential stress was calculated
!> \author Ole Schuett
! **************************************************************************************************
   SUBROUTINE get_stress_tensor(env_id, stress_tensor, ierr, available)
 
      INTEGER, INTENT(IN)                                :: env_id
      REAL(kind=dp), DIMENSION(3, 3), INTENT(OUT)        :: stress_tensor
      INTEGER, INTENT(OUT)                               :: ierr
      LOGICAL, INTENT(OUT), OPTIONAL                     :: available
 
      TYPE(cell_type), POINTER                           :: cell
      TYPE(cp_subsys_type), POINTER                      :: subsys
      TYPE(f_env_type), POINTER                          :: f_env
      TYPE(virial_type), POINTER                         :: virial
 
      NULLIFY (f_env, subsys, virial, cell)
      stress_tensor(:, :) = 0.0_dp
 
      CALL f_env_add_defaults(env_id, f_env)
      CALL force_env_get(f_env%force_env, subsys=subsys, cell=cell)
      CALL cp_subsys_get(subsys, virial=virial)
      IF (PRESENT(available)) available = virial%pv_availability
      IF (virial%pv_availability) THEN
         stress_tensor(:, :) = virial%pv_virial(:, :)/cell%deth
      END IF
      CALL f_env_rm_defaults(f_env, ierr)
 
   END SUBROUTINE get_stress_tensor
 
! **************************************************************************************************
!> \brief gets the positions of the particles
!> \param env_id id of the force_env
!> \param pos the array where to write the positions
!> \param n_el number of positions (3*nparticle) just to check
!> \param ierr will return a number different from 0 if there was an error
!> \date   22.11.2010 (MK)
!> \author fawzi
! **************************************************************************************************
   SUBROUTINE get_pos(env_id, pos, n_el, ierr)
 
      INTEGER, INTENT(IN)                                :: env_id, n_el
      REAL(kind=dp), DIMENSION(1:n_el)                   :: pos
      INTEGER, INTENT(OUT)                               :: ierr
 
      TYPE(f_env_type), POINTER                          :: f_env
 
      NULLIFY (f_env)
      CALL f_env_add_defaults(env_id, f_env)
      CALL force_env_get_pos(f_env%force_env, pos, n_el)
      CALL f_env_rm_defaults(f_env, ierr)
 
   END SUBROUTINE get_pos
 
! **************************************************************************************************
!> \brief gets the velocities of the particles
!> \param env_id id of the force_env
!> \param vel the array where to write the velocities
!> \param n_el number of velocities (3*nparticle) just to check
!> \param ierr will return a number different from 0 if there was an error
!> \author fawzi
!> date    22.11.2010 (MK)
! **************************************************************************************************
   SUBROUTINE get_vel(env_id, vel, n_el, ierr)
 
      INTEGER, INTENT(IN)                                :: env_id, n_el
      REAL(kind=dp), DIMENSION(1:n_el)                   :: vel
      INTEGER, INTENT(OUT)                               :: ierr
 
      TYPE(f_env_type), POINTER                          :: f_env
 
      NULLIFY (f_env)
      CALL f_env_add_defaults(env_id, f_env)
      CALL force_env_get_vel(f_env%force_env, vel, n_el)
      CALL f_env_rm_defaults(f_env, ierr)
 
   END SUBROUTINE get_vel
 
! **************************************************************************************************
!> \brief sets a new cell
!> \param env_id id of the force_env
!> \param new_cell the array with the cell matrix
!> \param ierr will return a number different from 0 if there was an error
!> \author Joost VandeVondele
! **************************************************************************************************
   SUBROUTINE set_cell(env_id, new_cell, ierr)
 
      INTEGER, INTENT(IN)                                :: env_id
      REAL(kind=dp), DIMENSION(3, 3)                     :: new_cell
      INTEGER, INTENT(OUT)                               :: ierr
 
      TYPE(cell_type), POINTER                           :: cell
      TYPE(cp_subsys_type), POINTER                      :: subsys
      TYPE(f_env_type), POINTER                          :: f_env
 
      NULLIFY (f_env, cell, subsys)
      CALL f_env_add_defaults(env_id, f_env)
      f_env%is_dirty = .true.
      NULLIFY (cell)
      CALL force_env_get(f_env%force_env, cell=cell)
      cpassert(ASSOCIATED(cell))
      cell%hmat = new_cell
      CALL init_cell(cell)
      CALL force_env_get(f_env%force_env, subsys=subsys)
      CALL cp_subsys_set(subsys, cell=cell)
      CALL f_env_rm_defaults(f_env, ierr)
 
   END SUBROUTINE set_cell
 
! **************************************************************************************************
!> \brief sets the positions of the particles
!> \param env_id id of the force_env
!> \param new_pos the array with the new positions
!> \param n_el number of positions (3*nparticle) just to check
!> \param ierr will return a number different from 0 if there was an error
!> \date   22.11.2010 updated (MK)
!> \author fawzi
! **************************************************************************************************
   SUBROUTINE set_pos(env_id, new_pos, n_el, ierr)
 
      INTEGER, INTENT(IN)                                :: env_id, n_el
      REAL(kind=dp), DIMENSION(1:n_el)                   :: new_pos
      INTEGER, INTENT(OUT)                               :: ierr
 
      TYPE(cp_subsys_type), POINTER                      :: subsys
      TYPE(f_env_type), POINTER                          :: f_env
 
      NULLIFY (f_env)
      CALL f_env_add_defaults(env_id, f_env)
      NULLIFY (subsys)
      CALL force_env_get(f_env%force_env, subsys=subsys)
      CALL unpack_subsys_particles(subsys=subsys, r=new_pos)
      f_env%is_dirty = .true.
      CALL f_env_rm_defaults(f_env, ierr)
 
   END SUBROUTINE set_pos
 
! **************************************************************************************************
!> \brief sets the velocities of the particles
!> \param env_id id of the force_env
!> \param new_vel the array with the new velocities
!> \param n_el number of velocities (3*nparticle) just to check
!> \param ierr will return a number different from 0 if there was an error
!> \date   22.11.2010 updated (MK)
!> \author fawzi
! **************************************************************************************************
   SUBROUTINE set_vel(env_id, new_vel, n_el, ierr)
 
      INTEGER, INTENT(IN)                                :: env_id, n_el
      REAL(kind=dp), DIMENSION(1:n_el)                   :: new_vel
      INTEGER, INTENT(OUT)                               :: ierr
 
      TYPE(cp_subsys_type), POINTER                      :: subsys
      TYPE(f_env_type), POINTER                          :: f_env
 
      NULLIFY (f_env)
      CALL f_env_add_defaults(env_id, f_env)
      NULLIFY (subsys)
      CALL force_env_get(f_env%force_env, subsys=subsys)
      CALL unpack_subsys_particles(subsys=subsys, v=new_vel)
      f_env%is_dirty = .true.
      CALL f_env_rm_defaults(f_env, ierr)
 
   END SUBROUTINE set_vel
 
! **************************************************************************************************
!> \brief updates the energy and the forces of given force_env
!> \param env_id id of the force_env that you want to update
!> \param calc_force if the forces should be updated, if false the forces
!>        might be wrong.
!> \param ierr will return a number different from 0 if there was an error
!> \author fawzi
! **************************************************************************************************
   RECURSIVE SUBROUTINE calc_energy_force(env_id, calc_force, ierr)
 
      INTEGER, INTENT(in)                                :: env_id
      LOGICAL, INTENT(in)                                :: calc_force
      INTEGER, INTENT(out)                               :: ierr
 
      TYPE(cp_logger_type), POINTER                      :: logger
      TYPE(f_env_type), POINTER                          :: f_env
 
      NULLIFY (f_env)
      CALL f_env_add_defaults(env_id, f_env)
      logger => cp_get_default_logger()
      CALL cp_iterate(logger%iter_info) ! add one to the iteration count
      f_env%is_dirty = .true.
      CALL force_env_calc_energy_force(f_env%force_env, calc_force=calc_force)
      f_env%is_dirty = .false.
      CALL f_env_rm_defaults(f_env, ierr)
 
   END SUBROUTINE calc_energy_force
 
! **************************************************************************************************
!> \brief returns the energy of the last configuration calculated
!> \param env_id id of the force_env that you want to update
!> \param e_pot the potential energy of the system
!> \param ierr will return a number different from 0 if there was an error
!> \author fawzi
! **************************************************************************************************
   SUBROUTINE get_energy(env_id, e_pot, ierr)
 
      INTEGER, INTENT(in)                                :: env_id
      REAL(kind=dp), INTENT(out)                         :: e_pot
      INTEGER, INTENT(out)                               :: ierr
 
      TYPE(f_env_type), POINTER                          :: f_env
 
      NULLIFY (f_env)
      CALL f_env_add_defaults(env_id, f_env)
      CALL force_env_get(f_env%force_env, potential_energy=e_pot)
      CALL f_env_rm_defaults(f_env, ierr)
 
   END SUBROUTINE get_energy
 
! **************************************************************************************************
!> \brief returns the energy of the configuration given by the positions
!>      passed as argument
!> \param env_id id of the force_env that you want to update
!> \param pos array with the positions
!> \param n_el number of elements in pos (3*natom)
!> \param e_pot the potential energy of the system
!> \param ierr will return a number different from 0 if there was an error
!> \author fawzi
!> \note
!>      utility call
! **************************************************************************************************
   RECURSIVE SUBROUTINE calc_energy(env_id, pos, n_el, e_pot, ierr)
 
      INTEGER, INTENT(IN)                                :: env_id, n_el
      REAL(kind=dp), DIMENSION(1:n_el), INTENT(IN)       :: pos
      REAL(kind=dp), INTENT(OUT)                         :: e_pot
      INTEGER, INTENT(OUT)                               :: ierr
 
      REAL(kind=dp), DIMENSION(1)                        :: dummy_f
 
      CALL calc_force(env_id, pos, n_el, e_pot, dummy_f, 0, ierr)
 
   END SUBROUTINE calc_energy
 
! **************************************************************************************************
!> \brief returns the energy of the configuration given by the positions
!>      passed as argument
!> \param env_id id of the force_env that you want to update
!> \param pos array with the positions
!> \param n_el_pos number of elements in pos (3*natom)
!> \param e_pot the potential energy of the system
!> \param force array that will contain the forces
!> \param n_el_force number of elements in force (3*natom). If 0 the
!>        forces are not calculated
!> \param ierr will return a number different from 0 if there was an error
!> \author fawzi
!> \note
!>      utility call, but actually it could be a better and more efficient
!>      interface to connect to other codes if cp2k would be deeply
!>      refactored
! **************************************************************************************************
   RECURSIVE SUBROUTINE calc_force(env_id, pos, n_el_pos, e_pot, force, n_el_force, ierr)
 
      INTEGER, INTENT(in)                                :: env_id, n_el_pos
      REAL(kind=dp), DIMENSION(1:n_el_pos), INTENT(in)   :: pos
      REAL(kind=dp), INTENT(out)                         :: e_pot
      INTEGER, INTENT(in)                                :: n_el_force
      REAL(kind=dp), DIMENSION(1:n_el_force), &
         INTENT(inout)                                   :: force
      INTEGER, INTENT(out)                               :: ierr
 
      LOGICAL                                            :: calc_f
 
      calc_f = (n_el_force /= 0)
      CALL set_pos(env_id, pos, n_el_pos, ierr)
      IF (ierr == 0) CALL calc_energy_force(env_id, calc_f, ierr)
      IF (ierr == 0) CALL get_energy(env_id, e_pot, ierr)
      IF (calc_f .AND. (ierr == 0)) CALL get_force(env_id, force, n_el_force, ierr)
 
   END SUBROUTINE calc_force
 
! **************************************************************************************************
!> \brief performs a check of the input
!> \param input_declaration ...
!> \param input_file_path the path of the input file to check
!> \param output_file_path path of the output file (to which it is appended)
!>        if it is "__STD_OUT__" the default_output_unit is used
!> \param echo_input if the parsed input should be written out with all the
!>        defaults made explicit
!> \param mpi_comm the mpi communicator (if not given it uses the default
!>        one)
!> \param initial_variables key-value list of initial preprocessor variables
!> \param ierr error control, if different from 0 there was an error
!> \author fawzi
! **************************************************************************************************
   SUBROUTINE check_input(input_declaration, input_file_path, output_file_path, &
                          echo_input, mpi_comm, initial_variables, ierr)
      TYPE(section_type), POINTER                        :: input_declaration
      CHARACTER(len=*), INTENT(in)                       :: input_file_path, output_file_path
      LOGICAL, INTENT(in), OPTIONAL                      :: echo_input
      TYPE(mp_comm_type), INTENT(in), OPTIONAL           :: mpi_comm
      CHARACTER(len=default_path_length), &
         DIMENSION(:, :), INTENT(IN)                     :: initial_variables
      INTEGER, INTENT(out)                               :: ierr
 
      INTEGER                                            :: unit_nr
      LOGICAL                                            :: my_echo_input
      TYPE(cp_logger_type), POINTER                      :: logger
      TYPE(mp_para_env_type), POINTER                    :: para_env
      TYPE(section_vals_type), POINTER                   :: input_file
 
      my_echo_input = .false.
      IF (PRESENT(echo_input)) my_echo_input = echo_input
 
      IF (PRESENT(mpi_comm)) THEN
         ALLOCATE (para_env)
         para_env = mpi_comm
      ELSE
         para_env => default_para_env
         CALL para_env%retain()
      END IF
      IF (para_env%is_source()) THEN
         IF (output_file_path == "__STD_OUT__") THEN
            unit_nr = default_output_unit
         ELSE
            CALL open_file(file_name=output_file_path, file_status="UNKNOWN", &
                           file_action="WRITE", file_position="APPEND", &
                           unit_number=unit_nr)
         END IF
      ELSE
         unit_nr = -1
      END IF
 
      NULLIFY (logger)
      CALL cp_logger_create(logger, para_env=para_env, &
                            default_global_unit_nr=unit_nr, &
                            close_global_unit_on_dealloc=.false.)
      CALL cp_add_default_logger(logger)
      CALL cp_logger_release(logger)
 
      input_file => read_input(input_declaration, input_file_path, initial_variables=initial_variables, &
                               para_env=para_env)
      CALL check_cp2k_input(input_declaration, input_file, para_env=para_env, output_unit=unit_nr)
      IF (my_echo_input .AND. (unit_nr > 0)) THEN
         CALL section_vals_write(input_file, &
                                 unit_nr=unit_nr, &
                                 hide_root=.true., &
                                 hide_defaults=.false.)
      END IF
      CALL section_vals_release(input_file)
 
      CALL cp_logger_release(logger)
      CALL mp_para_env_release(para_env)
      ierr = 0
      CALL cp_rm_default_logger()
 
   END SUBROUTINE check_input
 
END MODULE f77_interface