!  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/>.

!TODO unify rk and kind=8

subroutine dlf_get_params(nvar,            &
                          nvar2,           &
                          nspec,           &
                          coords,          &
                          coords2,         &
                          spec,            &
                          ierr,            &
                          tolerance,       &
                          printl,          &
                          maxcycle,        &
                          maxene,          &
                          ivar,            &
                          icoord,          &
                          iopt,            &
                          iline,           &
                          maxstep,         &
                          scale_step,      &
                          lbfgs_mem,       &
                          nimages,         &
                          nebk,            &
                          dump,            &
                          restart,         &
                          nz,              &
                          ncons,           &
                          nconn,           &
                          update,          &
                          maxupdate,       &
                          delta,           &
                          soft,            &
                          inithessian,     &
                          carthessian,     &
                          tsrelative,      &
                          maxrot,          &
                          tolrot,          &
                          nframes,         &
                          nmass,           &
                          nweight,         &
                          timestep,        &
                          fric0,           &
                          fricfac,         &
                          fricp,           &
                          imultistate,     &
                          state_i,         &
                          state_j,         &
                          pf_c1,           &
                          pf_c2,           &
                          gp_c3,           &
                          gp_c4,           &
                          ln_t1,           &
                          ln_t2,           &
                          printf,          &
                          tolerance_e,     &
                          distort,         &
                          massweight,      &
                          minstep,         &
                          maxdump,         &
                          task,            &
                          temperature,     &
                          po_pop_size,     &
                          po_radius,       &
                          po_contraction,  &
                          po_tolerance_r,  &
                          po_tolerance_g,  &
                          po_distribution, &
                          po_maxcycle,     &
                          po_init_pop_size,&
                          po_reset,        &
                          po_mutation_rate,&
                          po_death_rate,   &
                          po_scalefac,     &
                          po_nsave,        &
                          ntasks,          &
                          tdlf_farm,       &
                          n_po_scaling,    &
                          neb_climb_test,  &
                          neb_freeze_test, &
                          nzero,           &
                          coupled_states,  &
                          qtsflag,         &
                          imicroiter,      &
                          maxmicrocycle,   &
                          micro_esp_fit)

    use iso_c_binding
    use ChemShellParallelModule, only: chemsh_parallel
    use dlf_parameter_module   , only: rk
    use ChemShellDLFModule     , only: PyFrag, PyFrag2, PyOptions

    implicit none

    ! options (defined in chemsh/tasks/opt/opt.py) are pointers while others are not
    integer , intent(in)    :: nvar
    integer , intent(in)    :: nvar2
    integer , intent(in)    :: nspec
    real(rk), intent(inout) :: coords(nvar)
    real(rk), intent(inout) :: coords2(nvar2)
    integer , intent(inout) :: spec(nspec)
    integer , intent(out)   :: ierr
    real(rk), intent(out)   :: tolerance
    integer , intent(out)   :: printl
    integer , intent(out)   :: maxcycle
    integer , intent(out)   :: maxene
    integer , intent(out)   :: ivar
    integer , intent(inout) :: icoord
    integer , intent(inout) :: iopt
    integer , intent(inout) :: iline
    real(rk), intent(out)   :: maxstep
    real(rk), intent(out)   :: scale_step
    integer , intent(out)   :: lbfgs_mem
    integer , intent(out)   :: nimages
    real(rk), intent(out)   :: nebk
    integer , intent(out)   :: dump
    integer , intent(out)   :: restart
    integer , intent(inout) :: nz
    integer , intent(inout) :: ncons
    integer , intent(inout) :: nconn
    integer , intent(inout) :: update
    integer , intent(out)   :: maxupdate
    real(rk), intent(out)   :: delta
    real(rk), intent(out)   :: soft
    integer , intent(out)   :: inithessian
    integer , intent(inout) :: carthessian
    integer , intent(out)   :: tsrelative
    integer , intent(out)   :: maxrot
    real(rk), intent(out)   :: tolrot
    integer , intent(inout) :: nframes
    integer , intent(inout) :: nmass
    integer , intent(inout) :: nweight
    real(rk), intent(out)   :: timestep
    real(rk), intent(out)   :: fric0
    real(rk), intent(out)   :: fricfac
    real(rk), intent(out)   :: fricp
    integer , intent(inout) :: imultistate
    integer , intent(out)   :: state_i
    integer , intent(out)   :: state_j
    real(rk), intent(out)   :: pf_c1
    real(rk), intent(out)   :: pf_c2
    real(rk), intent(out)   :: gp_c3
    real(rk), intent(out)   :: gp_c4
    real(rk), intent(out)   :: ln_t1
    real(rk), intent(out)   :: ln_t2
    integer , intent(inout) :: printf
    real(rk), intent(out)   :: tolerance_e
    real(rk), intent(out)   :: distort
    integer , intent(inout) :: massweight
    real(rk), intent(out)   :: minstep
    integer , intent(inout) :: maxdump
    integer , intent(inout) :: task
    real(rk), intent(out)   :: temperature
    integer , intent(out)   :: po_pop_size
    real(rk), intent(out)   :: po_radius
    real(rk), intent(out)   :: po_contraction
    real(rk), intent(out)   :: po_tolerance_r
    real(rk), intent(out)   :: po_tolerance_g
    integer , intent(out)   :: po_distribution
    integer , intent(out)   :: po_maxcycle
    integer , intent(out)   :: po_init_pop_size
    integer , intent(out)   :: po_reset
    real(rk), intent(out)   :: po_mutation_rate
    real(rk), intent(out)   :: po_death_rate
    real(rk), intent(out)   :: po_scalefac
    integer , intent(out)   :: po_nsave
    integer , intent(inout) :: ntasks
    integer , intent(inout) :: tdlf_farm
    integer , intent(inout) :: n_po_scaling
    real(rk), intent(out)   :: neb_climb_test
    real(rk), intent(out)   :: neb_freeze_test
    integer , intent(out)   :: nzero
    integer , intent(out)   :: coupled_states
    integer , intent(inout) :: qtsflag
    integer , intent(inout) :: imicroiter
    integer , intent(out)   :: maxmicrocycle
    integer , intent(out)   :: micro_esp_fit

    integer                   :: debug, i, iat, istart, ityp, natoms, natoms2, iparameterise
    integer                   :: nresidues, ncartcons
    integer                   :: save_path

    real(kind=8)    , pointer :: onedimdbl(:), twodimdbl(:,:)
    integer(kind=8) , pointer :: onedimint(:), twodimint(:,:)

    character(len=255)        :: algorithm, coordinate
    character(len=:), pointer :: trust_radius

    ! parameterisation mode
    logical                   :: parameterise

    ! print level
    call PyOptions%get("save_path", save_path)
    if(save_path.eq.1) then
        printf = 4
    else
        printf = 2
    endif
    call PyOptions%get("debug", debug)
    ! TWK 20/12/2021 Set default DL-FIND print level to 2 as in Tcl-ChemShell
    printl = debug + 2
    ! printf should be no smaller than printl
    printf = max(printf, printl)

    ! get natoms
    call PyFrag%get('npts', natoms)

    ierr = 0

    ! let ChemShell set up MPI for DL_FIND 
    ntasks = chemsh_parallel%nworkgroups
    tdlf_farm = 0

! >>> COORDS

    call PyOptions%get("parameterise", iparameterise)
    parameterise = iparameterise.eq.1
    ! atom mode
    if(.not.parameterise) then
        nz     = natoms
        nweight = 0
        ! YL: copying is unavoidable without modifying DLF
        allocate(twodimdbl(3,natoms))
        call PyFrag%get('xyz', twodimdbl)
        coords  = reshape(twodimdbl, (/nvar/))
        deallocate(twodimdbl)
        ivar = 1
    ! parameter mode
    else
        nz     = 0
        allocate(onedimdbl(nvar))
        call PyOptions%get("parameters", onedimdbl)
        coords  = onedimdbl
        deallocate(onedimdbl)
        ivar = 0
    endif

! <<< END OF COORDS

    ! remember: nvarin=natoms*3 and nvarin2=3*natoms*nframes+natoms
    nmass  = natoms
    ntasks = 1
    nconn  = 0
    nresidues = 0

    ! not in use but needed
    task = 0
    po_distribution = 3

    ! TODO
    temperature = 300.0
    qtsflag = 0

    ! number of frames
    call PyOptions%get("nframes", nframes)

    if(.not.parameterise) then
        ! Set up DL-FIND integer specification array
        !    nat     entries of freezing/residue number
        !    nz      entries of nuclear charges (same order as coords)
        !    5*ncons entries of HDLC constraints (typ, atom1,atom2, atom3, atom4)
        !    2*nconn entries of connections (atom1 atom2)
        !    nat     entries of microiterative region specification
        ! i.e. nspec= nat + nz + 5*ncons + 2*nconn + nat
        spec(:) = 0
        allocate(onedimint(natoms))
        ! frozen centres
        call PyFrag%get("frozen", onedimint)
        spec(1:natoms) = onedimint(:)
        ! TWK 17/08/2022 process HDLC residue numbers
        ! Make sure not to overwrite frozen atoms
        ! TODO: delete residues entirely if they become too small (<3 atoms)
        ! as a result of removing frozen atoms from them
        call PyOptions%get("nresidues", nresidues)
        if (nresidues > 0) then
            call PyOptions%get("resIDs", onedimint)
            do i = 1, natoms
                if (spec(i) /= -1) then
                    spec(i) = onedimint(i) + 1
                end if
            end do
        end if
        ! TWK 19/08/2022 process constraints used with Cartesian coordinates
        ! Again make sure not to overwrite frozen atoms, which take precedence.
        call PyOptions%get("ncartcons", ncartcons)
        if (ncartcons > 0) then
            allocate(twodimint(5,ncartcons))
            call PyOptions%get("constraints", twodimint)
            do i = 1, ncartcons
                iat = twodimint(2, i) + 1
                if (spec(iat) /= -1) then
                    spec(iat) = twodimint(1, i)
                end if
            end do
            deallocate(twodimint)
        end if
        ! atomic numbers
        call PyFrag%get("znums", onedimint)
        spec(1+natoms:natoms+natoms) = onedimint(:)
        ! TWK 19/08/2022 process constraints used with HDLC coordinates
        call PyOptions%get("ncons", ncons)
        if (ncons > 0) then
            allocate(twodimint(5,ncons))
            call PyOptions%get("constraints", twodimint)
            istart = 1+natoms+natoms
            do i = 1, ncons
                ityp = twodimint(1,i)
                spec(istart) = ityp 
                ! DL-FIND atom numbers are offset by 1 from Py-ChemShell
                ! BUT we must leave the empty entries at 0...
                spec(istart+1:istart+1+ityp) = twodimint(2:ityp+2,i) + 1
                istart = istart + 5 
            end do
            deallocate(twodimint)
        end if
        ! TODO: connections, microiterative specs
        deallocate(onedimint)
    else
        spec = 0
    endif

! COORDS2

    coords2 = 0.0

    ! auxiliary array ("coords2" in Tcl-ChemShell)
    if(nframes > 0.and.associated(PyFrag2)) then
        call PyFrag2%get("npts", natoms2)
        allocate(twodimdbl(3,natoms2))
        call PyFrag2%get("xyz", twodimdbl)
        if (debug > 0) then
            write(*,'(/" DL_FIND auxiliary coordinates (a.u.):"/)')
            do i = 1, natoms2
                write(*,'(3F16.9)') twodimdbl(:,i)
            enddo
        end if
!        coords2(nvar*nframes+1:) = theory%frag%masses
        coords2(:nvar*nframes) = reshape(twodimdbl, (/nvar*nframes/))
        deallocate(twodimdbl)

    endif

    if(.not.parameterise) then
        ! TODO: to implement nweights
        ! weights
        allocate(onedimdbl(natoms))
        call PyFrag%get("masses", onedimdbl)
        coords2(nvar*nframes+1:) = onedimdbl
        deallocate(onedimdbl)
    endif

! END OF COORDS2
! OPTIONS

    call PyOptions%get("pf_c1", pf_c1)
    call PyOptions%get("tolerance"       , tolerance)
!    call PyOptions%get("printl"          , printl)    ! set by `debug`
    call PyOptions%get("maxcycle"        , maxcycle)
    call PyOptions%get("maxene"          , maxene)
    ! 3D coordinates or sequetial variables mode (1 for atoms mode)
    call PyOptions%get("icoord"          , icoord)
    call PyOptions%get("iopt"            , iopt)
    call PyOptions%get("maxstep"         , maxstep)
    call PyOptions%get("scale_step"      , scale_step)
    call PyOptions%get("lbfgs_mem"       , lbfgs_mem)
    call PyOptions%get("nimages"         , nimages)
    call PyOptions%get("nebk"            , nebk)
    call PyOptions%get("dump"            , dump)
    call PyOptions%get("restart"         , restart)
!    call PyOptions%get("nz"              , nz)
    call PyOptions%get("ncartcons"       , ncartcons)
    call PyOptions%get("ncons"           , ncons)
!    call PyOptions%get("nconn"           , nconn)
    call PyOptions%get("update"          , update)
    call PyOptions%get("maxupdate"       , maxupdate)
    call PyOptions%get("delta"           , delta)
    call PyOptions%get("soft"            , soft)
    call PyOptions%get("init_hessian"    , inithessian)
    call PyOptions%get("cart_hessian"    , carthessian)
    call PyOptions%get("tsrelative"      , tsrelative)
    call PyOptions%get("maxrot"          , maxrot)
    call PyOptions%get("tolrot"          , tolrot)
!    call PyOptions%get("nmass"           , nmass)
!    call PyOptions%get("nweight"         , nweight)
    call PyOptions%get("timestep"        , timestep)
    call PyOptions%get("fric0"           , fric0)
    call PyOptions%get("fricfac"         , fricfac)
    call PyOptions%get("fricp"           , fricp)
!    call PyOptions%get("imultistate"     , imultistate)
    call PyOptions%get("state_i"         , state_i)
    call PyOptions%get("state_j"         , state_j)
    call PyOptions%get("pf_c1"           , pf_c1)
    call PyOptions%get("pf_c2"           , pf_c2)
    call PyOptions%get("gp_c3"           , gp_c3)
    call PyOptions%get("gp_c4"           , gp_c4)
    call PyOptions%get("ln_t1"           , ln_t1)
    call PyOptions%get("ln_t2"           , ln_t2)
    call PyOptions%get("tolerance_e"     , tolerance_e)
    call PyOptions%get("distort"         , distort)
!    call PyOptions%get("massweight"      , massweight)
    call PyOptions%get("minstep"         , minstep)
!    call PyOptions%get("maxdump"         , maxdump)
!    call PyOptions%get("task"            , task)
    call PyOptions%get("temperature"     , temperature)
    call PyOptions%get("po_pop_size"     , po_pop_size)
    call PyOptions%get("po_radius"       , po_radius)
    call PyOptions%get("po_contraction"  , po_contraction)
    call PyOptions%get("po_tolerance_r"  , po_tolerance_r)
    call PyOptions%get("po_tolerance_g"  , po_tolerance_g)
    call PyOptions%get("po_distribution" , po_distribution)
    call PyOptions%get("po_maxcycle"     , po_maxcycle)
    call PyOptions%get("po_init_pop_size", po_init_pop_size)
    call PyOptions%get("po_reset"        , po_reset)
    call PyOptions%get("po_mutation_rate", po_mutation_rate)
    call PyOptions%get("po_death_rate"   , po_death_rate)
    call PyOptions%get("po_scalefac"     , po_scalefac)
    call PyOptions%get("po_nsave"        , po_nsave)
!    call PyOptions%get("ntasks"          , ntasks)
!    call PyOptions%get("tdlf_farm"       , tdlf_farm)
!    call PyOptions%get("n_po_scaling"    , n_po_scaling)
    call PyOptions%get("neb_climb_test"  , neb_climb_test)
    call PyOptions%get("neb_freeze_test" , neb_freeze_test)
    call PyOptions%get("nzero"           , nzero)
    call PyOptions%get("coupled_states"  , coupled_states)
!    call PyOptions%get("qtsflag"         , qtsflag)
!    call PyOptions%get("imicroiter"      , imicroiter)
    call PyOptions%get("maxmicrocycle"   , maxmicrocycle)
    call PyOptions%get("micro_esp_fit"   , micro_esp_fit)
    call PyOptions%get("trust_radius"    , trust_radius)

    if (debug > 0) then
        write(*, '(1X,A)') ">>> DL_FIND specifications:"
        select case(iopt)
            case(0)
                algorithm = "steepest descent"
            case(1)
                algorithm = "conjugate gradient"
            case(2)
                algorithm = "conjugate gradient with CG restart"
            case(3)
                algorithm = "L-BFGS"
            case(9)
                algorithm = "delta test for finite-difference in gradients"
            case(10)
                algorithm = "P-RFO"
            case(11)
                algorithm = "just Hessian and thermal analysis"
            case(12)
                algorithm = "Hessian of all images and qTS rate"
            case(13)
                algorithm = "rate without tunnelling"
            case(20)
                algorithm = "Newton-Raphson (quasi-Newton)"
            case(30)
                algorithm = "damped dynamics"
            case(51)
                algorithm = "random (stochastic) search"
            case(52)
                algorithm = "genetic algorithm"
            case default
        end select
        write(*, '(5X,A,A,1X,A,I2,A)') "algorithm : ", trim(algorithm), "(iopt =", iopt, ")"
        select case(icoord)
            case(0)
                coordinate = "Cartesian"
            case(1)
                coordinate = "HDLC"
            case(2)
                coordinate = "HDLC-TC (total connection)"
            case(3)
                coordinate = "DLC"
            case(4)
                coordinate = "DLC-TC (total connection)"
        end select
        write(*, '(5X,A,A,1X,A,I2,A)') "coordinate: ", trim(coordinate), "(icoord =", icoord, ")"
    end if

    ! trust radius method
    selectcase(trim(trust_radius))
        case("constant")
            iline = 0
        case("scaled")
            iline = 0
        case("energy")
            iline = 1
        case("gradient")
            iline = 2
        case("line_search")
            iline = 3
        case("line")
            iline = 3
    endselect

    if(debug > 3) then 
       print *
       print *, ">>> DL-FIND specifications:"
       print *
       print *, "    nvar     : ", nvar
       print *, "    nvar2    : ", nvar2
       print *, "    nspec    : ", nspec
       print *, "    icoord   : ", icoord
       print *, "    iopt     : ", iopt
       print *, "    iline    : ", iline
       print *, "    nimages  : ", nimages
       print *, "    nframes  : ", nframes
       print *, "    tolerance: ", tolerance
       print *
    endif

    if(debug > 4) then
        print *, " >>> DL_FIND parameters:"
        print *, " coords(nvar)   = ", coords
        print *, " coords2(nvar2) = ", coords2
        print *, " spec(nspec)         = ", spec
        print *, " ierr                = ", ierr
        print *, " printl              = ", printl
        print *, " maxcycle            = ", maxcycle
        print *, " maxene              = ", maxene
        print *, " maxstep             = ", maxstep
        print *, " scale_step          = ", scale_step
        print *, " lbfgs_mem           = ", lbfgs_mem
        print *, " nebk                = ", nebk
        print *, " dump                = ", dump
        print *, " restart             = ", restart
        print *, " nz                  = ", nz
        print *, " ncons               = ", ncons
        print *, " nconn               = ", nconn
        print *, " update              = ", update
        print *, " maxupdate           = ", maxupdate
        print *, " delta               = ", delta
        print *, " soft                = ", soft
        print *, " inithessian         = ", inithessian
        print *, " carthessian         = ", carthessian
        print *, " tsrelative          = ", tsrelative
        print *, " maxrot              = ", maxrot
        print *, " tolrot              = ", tolrot
        print *, " nmass               = ", nmass
        print *, " nweight             = ", nweight
        print *, " timestep            = ", timestep
        print *, " fric0               = ", fric0
        print *, " fricfac             = ", fricfac
        print *, " fricp               = ", fricp
        print *, " imultistate         = ", imultistate
        print *, " state_i             = ", state_i
        print *, " state_j             = ", state_j
        print *, " pf_c1               = ", pf_c1
        print *, " pf_c2               = ", pf_c2
        print *, " gp_c3               = ", gp_c3
        print *, " gp_c4               = ", gp_c4
        print *, " ln_t1               = ", ln_t1
        print *, " ln_t2               = ", ln_t2
        print *, " printf              = ", printf
        print *, " tolerance_e         = ", tolerance_e
        print *, " distort             = ", distort
        print *, " massweight          = ", massweight
        print *, " minstep             = ", minstep
        print *, " maxdump             = ", maxdump
        print *, " task                = ", task
        print *, " temperature         = ", temperature
        print *, " po_pop_size         = ", po_pop_size
        print *, " po_radius           = ", po_radius
        print *, " po_contraction      = ", po_contraction
        print *, " po_tolerance_r      = ", po_tolerance_r
        print *, " po_tolerance_g      = ", po_tolerance_g
        print *, " po_distribution     = ", po_distribution
        print *, " po_maxcycle         = ", po_maxcycle
        print *, " po_init_pop_size    = ", po_init_pop_size
        print *, " po_reset            = ", po_reset
        print *, " po_mutation_rate    = ", po_mutation_rate
        print *, " po_death_rate       = ", po_death_rate
        print *, " po_scalefac         = ", po_scalefac
        print *, " po_nsave            = ", po_nsave
        print *, " ntasks              = ", ntasks
        print *, " tdlf_farm           = ", tdlf_farm
        print *, " n_po_scaling        = ", n_po_scaling
        print *, " neb_climb_test      = ", neb_climb_test
        print *, " neb_freeze_test     = ", neb_freeze_test
        print *, " nzero               = ", nzero
        print *, " coupled_states      = ", coupled_states
        print *, " qtsflag             = ", qtsflag
        print *, " imicroiter          = ", imicroiter
        print *, " maxmicrocycle       = ", maxmicrocycle
        print *, " micro_esp_fit       = ", micro_esp_fit
    endif

endsubroutine dlf_get_params
! YL 07/02/2021: kiter is always -1
subroutine dlf_get_gradient(nvar,      &
                            coords,    &
                            energy,    &
                            gradients, &
                            iimage,    &
                            kiter,     &
                            status)

    use iso_c_binding
    use dlf_parameter_module, only: rk
    use ChemShellDLFModule  , only: PyFrag, PyResult, PyCallback, PyOptions, PySetParams, ierror
    use dlf_convergence     , only: valg, valrmsg, vals, valrmss

    implicit none
    integer          , intent(in)  :: nvar
    real(rk), target , intent(in)  :: coords(nvar)
    real(rk)         , intent(out) :: energy
    real(rk)         , intent(out) :: gradients(nvar)
    integer          , intent(in)  :: iimage
    integer          , intent(in)  :: kiter
    integer          , intent(out) :: status

    integer                        :: i, debug, natoms, iparameterise
    real(rk), pointer              :: onedimdbl(:), twodimdbl(:,:)
    real(rk), pointer              :: temp(:)
    integer(kind=8), pointer       :: onedimint(:)

    real :: start_time, end_time

    logical                        :: parameterise

    status = 0

    call PyOptions%get("debug", debug)

    call PyOptions%get("parameterise", iparameterise)
    parameterise = iparameterise.eq.1
    
    ! atomic mode
    if(.not.parameterise) then

        call PyFrag%get('npts' , natoms)

        allocate(onedimint(natoms))
        call PyFrag%get("znums", onedimint)

        if(debug > 3) then
            allocate(twodimdbl(3,natoms))
            call PyFrag%get('xyz', twodimdbl)
            write(*,'(/" DL_FIND original coords (a.u.):"/)')
            write(*,'(3X,"index",1X,"Z-number",12X,"X",15X,"Y",15X,"Z"/)')
            do i = 1, natoms
                write(*,'(1X,I6,3X,I6,3X,3F16.9)') i-1, onedimint(i), twodimdbl(:,i)
            enddo
            call flush(6)
            deallocate(twodimdbl)
        endif

        ! update coords
        call PyFrag%set('xyz', reshape(coords, (/3,natoms/)))
 
        ! CALLBACK
        if(debug > 0) then
            write(*,'(/" >>> DL_FIND module is calling callback function for gradients evaluation..."/)')
        endif
        
        call cpu_time(start_time)

        ! energy of the previous step
        allocate(onedimdbl(1))
        call PyResult%get('energy', onedimdbl)
        energy = onedimdbl(1)
        ! save the energy, max grad, max step, RMS grad, RMS step
        call PyResult%set('convergences', (/energy,valg,vals,valrmsg,valrmss/))

        ! do calculations on the Python side
        ierror = PyCallback()
        call flush(5)
        call flush(6)
        if(ierror.ne.0) then
    !        print *, ">>> ERROR: DL_FIND failed to run callback function. Exiting..."
            call exit(ierror)
        endif
        call cpu_time(end_time)
        if (debug > 4) then
            print *, ">>> DL_FIND: PyCallback() costs ", end_time - start_time, " seconds"
        endif
        ! END OF CALLBACK 

        ! note the shape!
        allocate(twodimdbl(3,natoms))
    
        call PyResult%get('energy'   , onedimdbl)
        call PyResult%get('gradients', twodimdbl)

        ! update energy and gradients
        ! note: energy in Python is an array!
        energy = onedimdbl(1)

        ! can use c_f_pointer to reshape without copying, too
        ! call c_f_pointer(c_loc(twodimdbl), onedimdbl, [nvar])
        ! gradients = onedimdbl(:)

        if(debug > 3) then 
            write(*,'(/" DL_FIND energy (a.u.) of updated coordinates:"/)')
            write(*,'(4XF16.9)') energy
        endif

        if(debug > 3) then 
            write(*,'(/" DL_FIND gradients (a.u.) of updated coordinates:"/)')
            write(*,'(3X,"index",1X,"Z-number",12X,"X",15X,"Y",15X,"Z"/)')
        endif
    
        do i = 1, natoms
            gradients(i*3-2:i*3) = twodimdbl(:,i)
            if(debug > 3) then
                ! atomic numbers
                write(*,'(1X,I6,3X,I6,3X,3F16.9)') i-1, onedimint(i), gradients(i*3-2:i*3)
            endif
        enddo

        if(debug > 3) then
            write(*,'(//)')
        endif
        call flush(6)

        deallocate(onedimdbl, twodimdbl)
        deallocate(onedimint)

    ! parameterisation mode
    else

        allocate(onedimdbl(nvar))
        call PyOptions%get('parameters', onedimdbl)
        write(*,'(/" DL_FIND original parameters:"/)')
        do i = 1, nvar
            write(*,'(3X,I6,F20.12)') i-1, onedimdbl(i)
        enddo

        ! update parameters
        call PyOptions%set('parameters', coords)
        call flush(6)
        if(PySetParams() /= 0) then
            write(*, '(A)') ">>> ChemShell ERROR: DL_FIND failed to run PySetParams(). Exiting..."
            status = 1
        endif

        ! CALLBACK
        write(*,'(/" >>> DL_FIND is running the callback function for derivatives evaluation..."/)')
        ! do calculations on the Python side

        call cpu_time(start_time)
        call flush(6)

        ! run callback function to calculate the gradients
        if(PyCallback() /= 0) then
            print *, ">>> ERROR: DL_FIND failed to run callback function. Exiting..."
            status = 1
            stop
        endif

        call cpu_time(end_time)
        if(debug > 3) then 
            write (*, '(1X,A,F12.3,1X,A)') ">>> DL_FIND: CPU time for the callback function is", end_time-start_time, "s"
        endif
        call PyOptions%get('derivatives', onedimdbl)
        gradients(:) = onedimdbl(:)

        if(debug > 3) then 
            write(*,'(/" DL_FIND derivatives of updated parameters:"/)')
            write(*,'(3X,"index",7X,"derivative"/)')
            do i = 1, nvar
                write(*,'(1X,I6,3X,F16.9)') i-1, onedimdbl(i)
            enddo
            write(*,'(//)')
            call flush(6)
        endif
        deallocate(onedimdbl)

        ! note: energy in Python is an array!
        allocate(onedimdbl(1))
        call PyResult%get('energy', onedimdbl)
        energy = onedimdbl(1)

        deallocate(onedimdbl)
        ! END OF CALLBACK

    endif

endsubroutine dlf_get_gradient
! TODO: to implement Hessian array to feed DL_FIND
subroutine dlf_get_hessian(nvar, coords, hessian, status)

    use iso_c_binding
    use DLF_Parameter_Module, only: rk
    use ChemShellDLFModule  , only: PyFrag, PyCallback, ierror

    integer              , intent(in)  :: nvar
    real(rk)    , target , intent(in)  :: coords(nvar)
    real(rk)             , intent(out) :: hessian(nvar)
    integer              , intent(out) :: status
    
    integer                            :: natoms
    real(kind=8), pointer              :: twodimdbl(:,:)

    hessian = 0.0

    call PyFrag%get('npts', natoms)

    allocate(twodimdbl(3,natoms))
    call PyFrag%get('xyz', twodimdbl)
    write(*,'(/" DL_FIND dlf_get_hessian: original coords (a.u.):"/)')
    do i = 1, natoms
        write(*,'(3F16.9)') twodimdbl(:,i)
    enddo
    deallocate(twodimdbl)
    
    call PyFrag%set('xyz', reshape(coords, (/3,natoms/)))
    ! CALLBACK 
    write(*,'(/" >>> DL_FIND module is calling callback function for Hessian evaluation..."/)')
    ! do calculations on the Python side
    call cpu_time(start_time)
    ierror = PyCallback()
    if(ierror.ne.0) then
!        print *, ">>> ERROR: DL_FIND failed to run callback function. Exiting..."
        call exit(ierror)
    endif
    call cpu_time(end_time)
    print *, ">>> DL_FIND: PyCallback() costs ", end_time - start_time, " seconds"

! AD HOC
    status = 1
    ! END OF CALLBACK 

endsubroutine dlf_get_hessian
subroutine dlf_get_multistate_gradients()

    print *, "### PLACEHOLDER: dlf_get_multistate_gradients"

endsubroutine dlf_get_multistate_gradients
subroutine dlf_put_coords(nvar, mode, e0, coords0, iam)

    use dlf_parameter_module, only: rk
    use ChemShellDLFModule  , only: PyFrag, PyFrag2, PyResult, PyOptions

    implicit none

    integer                  , intent(in) :: nvar
    integer                  , intent(in) :: mode
    real(kind=8)             , intent(in) :: e0
    real(kind=8)             , intent(in) :: coords0(*)
    integer                  , intent(in) :: iam        ! a dummy argument

    real(rk)                              :: energy(1)
    integer                               :: i, natoms, natoms2, iparamterise
    real(kind=rk)   , pointer             :: coords(:,:)
    real(kind=8)    , pointer             :: onedimdbl(:)
    integer(kind=rk), pointer             :: znums(:)

    logical                               :: paramterise
    integer                               :: debug

    call PyOptions%get("debug", debug)
    call PyOptions%get("parameterise", iparamterise)
    paramterise = iparamterise.eq.1

    ! final energy
    if(mode>0) then
        energy = e0
        call PyResult%set('energy', energy)
    endif

    ! atomic mode
    if(.not.paramterise) then

        call PyFrag%get('npts', natoms)
        selectcase(mode)
            ! TODO
            case(:-1)
            ! update coordinates of frag
            case(1)
                call PyFrag%set('xyz', reshape(coords0(:nvar), (/3,natoms/)))
            ! update coordinates of frag2
            case(2)
                if(associated(PyFrag2)) then
                    ! YL: natoms2 here is the sum of frag2's natoms (TS and product)
                    ! YL 18/09/2018: as of new GNU compilers (7 or later?) setting values with part of the array (see dl_py2f.f90) causes segfault of mismatched array shapes, so use the full coordinates
!                    call PyFrag2%get('npts', natoms2)
                    ! YL 06/01/2022: it was reshape(coords0(:nvar), (/3,natoms2/)) here and might cause an overflow because the calling routine dlf_formstep_set_tsmode() in DL_FIND passes in an array of size nvar; in fact we only need to set the TS coordinates (verified to work with both Intel and GNU)
                    call PyFrag2%set('xyz', reshape(coords0(:nvar), (/3,natoms/)))
                endif
            ! same as mode 1
            case(3:5)
                call PyFrag%set('xyz', reshape(coords0(:nvar), (/3,natoms/)))
        endselect

        allocate(znums(natoms))
        call PyFrag%get('znums', znums)
    
        allocate(coords(3,natoms))
        call PyFrag%get('xyz', coords)

        if(debug > 3) then 
            write(*,'(//" DL_FIND has updated the Cartesian coords (a.u.):"/)')
            write(*,'(3X,"index",1X,"Z-number",12X,"X",15X,"Y",15X,"Z"/)')
            do i = 1, natoms
                write(*,'(1X,I6,3X,I6,3X,3F16.9)') i-1, znums(i), coords(:,i)
            end do
            write(*,'(//)')
            call flush(6)
        endif
    
        deallocate(znums, coords)

    ! parameter mode
    else

        allocate(onedimdbl(nvar))
        call PyOptions%get('parameters', onedimdbl)

        write(*,'(//" DL_FIND has updated the parameters:"/)')
        do i = 1, nvar
            write(*,'(3X,I6,F20.12)') i-1, onedimdbl(i)
        enddo
        write(*,'(//)')

        deallocate(onedimdbl)

    endif

endsubroutine dlf_put_coords

  
subroutine dlf_put_thermal()
endsubroutine dlf_put_thermal


subroutine dlf_error()

    write(*, *) ">>> ChemShell DL_FIND ERROR"
    ! DLF will quit without an error signal so cause this here
    call exit(1)
    ! TODO: return control to python with error message

endsubroutine dlf_error
subroutine dlf_update()
    print *, "### PLACEHOLDER: dlf_update"
endsubroutine dlf_update
subroutine dlf_get_taskfarm(nworkgroups, workgroupsize, workgroupnode, workgroupid, workgroup_comm, counterparts_comm)

    use ChemShellParallelModule, only : chemsh_parallel

    integer, intent(out) :: nworkgroups
    integer, intent(out) :: workgroupsize
    integer, intent(out) :: workgroupnode
    integer, intent(out) :: workgroupid
    integer, intent(out) :: workgroup_comm
    integer, intent(out) :: counterparts_comm
    
    nworkgroups       = chemsh_parallel%nworkgroups
    workgroupsize     =  1
    workgroupnode     = chemsh_parallel%node
    workgroupid       = chemsh_parallel%workgroup
    workgroup_comm    = -1
    counterparts_comm = chemsh_parallel%comm_node

endsubroutine dlf_get_taskfarm
subroutine dlf_put_taskfarm(nworkgroups, workgroupsize, workgroupnode, workgroupid, workgroup_comm, counterparts_comm)

    integer :: nworkgroups
    integer :: workgroupsize
    integer :: workgroupnode
    integer :: workgroupid
    integer :: workgroup_comm
    integer :: counterparts_comm

    print *, "### DL_FIND: dlf_put_taskfarm...."

endsubroutine dlf_put_taskfarm
subroutine dlf_put_procinfo(nprocs, rank, global_comm)

    integer, intent(in) :: nprocs
    integer, intent(in) :: rank
    integer, intent(in) :: global_comm

    print *, "### DL_FIND: dlf_put_procinfo..."

endsubroutine dlf_put_procinfo
subroutine dlf_get_procinfo(nprocs, rank, global_comm)

    use ChemShellParallelModule

    integer, intent(out) :: nprocs
    integer, intent(out) :: rank
    integer, intent(out) :: global_comm

    nprocs      = chemsh_parallel%nworkgroups
    rank        = chemsh_parallel%workgroup
    global_comm = chemsh_parallel%comm_node

endsubroutine dlf_get_procinfo
