!  Copyright (C) 2017 The authors of Py-ChemShell
!
!  This file is part of Py-ChemShell.
!
!  Py-ChemShell is free software: you can redistribute it and/or modify
!  it under the terms of the GNU Lesser General Public License as
!  published by the Free Software Foundation, either version 3 of the
!  License, or (at your option) any later version.
!
!  Py-ChemShell 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 Lesser General Public License for more details.
!
!  You should have received a copy of the GNU Lesser General Public
!  License along with Py-ChemShell.  If not, see
!  <http://www.gnu.org/licenses/>.

! Replace CP2K's main program (in CP2K/src/start/cp2k.F)
! with a ChemShell-specific startup routine
! YL 21/02/2021: upgraded to support v8.1
! TWK 10/10/2022: modified file handling to support v2022.1

subroutine cp2k(input, output, initialise, is_root, communicator)

!    use ChemShellParallelModule, only : chemsh_parallel
!    use omp_lib            , only : omp_set_num_threads
!    use cp2k_info          , only : cp2k_flags
    use cp2k_runs          , only : run_input
    use f77_interface      , only : finalize_cp2k, init_cp2k
    use input_cp2k         , only : create_cp2k_root_section
    use input_cp2k_read    , only : empty_initial_variables
    use input_section_types, only : section_release, section_type
    use kinds              , only : default_path_length
    use cp_files           , only : open_file, close_file, get_unit_number
!    use machine            , only : default_output_unit
!    use message_passing    , only : mp_comm_world
!    use timings            , only : add_timer_env
! YL 09/01/2020: a local copy of the original source file
#include "cp2k_base_uses.f90"

    implicit none
 
    character(len=*), intent(in)         :: input
    character(len=*), intent(in)         :: output
    logical         , intent(in)         :: initialise
    logical         , intent(in)         :: is_root
    integer(kind=4) , intent(in)         :: communicator
 
    character(len=default_path_length)   :: input_file_name, output_file_name
    integer                              :: output_unit
    integer                              :: ierror, linput, loutput, ifunit
    type(section_type), pointer          :: input_declaration

    ! set default behaviour for the command line switches
    input_file_name             = repeat(' ', default_path_length)
    output_file_name            = repeat(' ', default_path_length)
    linput                      = len(trim(input))
    loutput                     = len(trim(output))
    input_file_name (1:linput)  = input (1:linput)
    output_file_name(1:loutput) = output(1:loutput)
    call flush(6)
 
    if(initialise) then
        if(is_root) then
            write(*, '(/A/)') " ChemShell/CP2K: initialising..."
            call flush(6)
        endif
        call init_cp2k(.false., ierror)
    endif
 
    ! create a section instance
    nullify(input_declaration)
    call create_cp2k_root_section(input_declaration)
    call flush(6)

    ! run via input file
    if(is_root) then
        ifunit = get_unit_number(output_file_name)
        ! TWK 10/10/2022 CP2K no longer opens the output file in run_input,
        ! so we do it here
        call open_file(file_name=output_file_name, file_status="UNKNOWN", &
            file_action="WRITE", file_position="APPEND", skip_get_unit_number=.TRUE., &
            unit_number=ifunit)
    end if
    if(is_root) then
        write(*, '(/A,A,A/)') " ChemShell/CP2K: running input file ", trim(input_file_name), "..."
        call flush(6)
    endif
    call run_input(input_declaration, input_file_name, output_file_name, empty_initial_variables, communicator)
    call flush(6)

    ! YL 03/03/2021: because we archive the output file at every cycle of geometry
    !                geometry optimisation or MD simulation to keep each single-
    !                point calculation separate, each time we have to close the file unit
    !                which is not closed by the run_input() function (CP2K isn't designed
    !                to work this way); it looks safe to use the root process to close it
    !                (and get the unit number) but we can't get the local MPI environment
    !                instance para_env and its I/O node para_env%ionode; in fact it works
    !                OK without closing the output file if compiled with gfortran, but
    !                with ifort it won't work
    if(is_root) then
        call close_file(ifunit)
        call flush(6)
    endif

    ! clear the section instance
    call section_release(input_declaration)
    call flush(6)
 
    ! YL 01/03/2021: we don't finailse the CP2K instance for now because knowing
    !                when to do it needs an extra argument to be passed in
!    call finalize_cp2k(finalize_mpi=.false., ierror=ierror)
!    cpassert(ierror == 0)

endsubroutine
