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

! YL 01/08/2018: this file should contain ONLY methods called by the (Python) interface
module ChemShellClusterModule

    use iso_c_binding
    use DL_PY2F

    private

    public PyFrag, PyCluster, PyOptions, PyCell, PyField, PyBQs, PyShells
    public get_ndims, cluster_origin
    public MAXCONN
    public debug

    integer(kind=8), parameter :: MAXCONN = 14
    type(dictType) , pointer   :: PyFrag, PyCluster, PyOptions, PyCell, PyField, PyBQs, PyShells
    integer                    :: debug
    integer        , external  :: cluster_cut
    real(kind=8)               :: cluster_origin(3)

! YL: since Intel Compiler 2016, "external" is no longer allowed to be bind(c)
!     "error #8809: An OPTIONAL or EXTERNAL dummy argument to a BIND(C) procedure is not interoperable"
!    integer(c_long), pointer, external :: PyCallback
!    abstract interface
!        integer(c_long) function callback() bind(c)
!            use iso_c_binding
!        endfunction callback
!    endinterface
!
!    procedure(callback), pointer :: PyCallback

    contains

    ! this is a wrapper around the original cluster_cut() of CONSTRUCT
    integer(c_int) function cluster_get_size(natoms_cluster,  &
                                             nshells_cluster, &
                                             chemsh_debug)    &
                   result(ierror)                             &
                   bind(c, name="cluster_get_size")

        implicit none

        integer, intent(out) :: natoms_cluster
        integer, intent(out) :: nshells_cluster
        integer, intent(in)  :: chemsh_debug

        debug = chemsh_debug

        ! set size of new cluster to zero to allow cluster_cut to return
        natoms_cluster  = 0
        nshells_cluster = 0

        ! TODO support fully coordinated mode
        if(cluster_cut(natoms_cluster, nshells_cluster).ne.0) then
            write(*,'(/A)') " >>> ChemShell CLUSTER ERROR: cluster_cut called in cluster_get_size"
        endif
        call flush(6)

        ! DEBUG
        write(*,'(/A)') " The created cluster contains"
        write(*,'(1X,I0,1X,A,I0,1X,A/)') natoms_cluster, "atoms and ", nshells_cluster, "shell centres"

        ierror = 0

    endfunction cluster_get_size
    integer(c_int) function cluster_get_coords(clusterPtr,   &
                                               chemsh_debug) &
                   result(ierror)                            &
                   bind(c, name='cluster_get_coords')

        implicit none
            
        type(PyType)             , intent(inout) :: clusterPtr
        integer                  , intent(in)    :: chemsh_debug

        real(kind=8)                             :: radius_cluster
        real(kind=8)    , pointer                :: coords(:,:), charges(:)
!        real(kind=8)                             :: origin_point(3)
        integer                                  :: natoms_frag, nshells_frag
        integer                                  :: natoms_cluster, nshells_cluster
        integer(kind=8) , pointer                :: znums(:)
        character(len=:), pointer                :: crystal_type

        debug = chemsh_debug

        if(associated(PyCluster)) deallocate(PyCluster)
        allocate(PyCluster, source=ptr2dict(clusterPtr))

        !Get the parameters for the routine from Python
        call PyFrag%get('npts'             , natoms_frag)
        call PyFrag%get('nshells'          , nshells_frag)

        call PyOptions%get('radius_cluster', radius_cluster)
        call PyOptions%get('crystal_type'  , crystal_type)
!        call PyOptions%get('origin'        , origin_point)

        call PyCluster%get('npts'   , natoms_cluster)
        call PyCluster%get('nshells', nshells_cluster)

        ! YL 01/08/2018: what is this?
        ! TODO support fully coordinated mode
        if(cluster_cut(natoms_cluster, nshells_cluster).ne.0) then
            write(*,'(/A)') "Something went wrong getting coordinates of cluster"
        endif
        call flush(6)

        ierror = 0

    endfunction cluster_get_coords
    integer(c_int) function cluster_adjust_charges(clusterPtr,   &
                                                   chemsh_debug) &
                   result(ierror)                                &
                   bind(c, name='cluster_adjust_charges')
      
        implicit none

        type(PyType)               , intent(inout) :: clusterPtr
        integer                    , intent(in)    :: chemsh_debug
        
        integer         , parameter                :: nbqs=0
        integer         , parameter                :: nbqs_cluster=0
        integer                                    :: natoms_frag, nshells_frag
        integer                                    :: natoms_cluster, nshells_cluster
        integer                                    :: dimensionality
        character(len=:), pointer                  :: crystal_type
        character(len=:), pointer                  :: adjust_mode
        
        debug = chemsh_debug

        if(associated(PyCluster)) deallocate(PyCluster)
        allocate(PyCluster, source=ptr2dict(clusterPtr))
        
        ! get the parameters for the routine from Python
        call PyFrag   %get('npts'         , natoms_frag)
        call PyFrag   %get('nshells'      , nshells_frag)
        call PyCell   %get('ndimensions'  , dimensionality)
        call PyCluster%get('npts'         , natoms_cluster)
        call PyCluster%get('nshells'      , nshells_cluster)
        call PyOptions%get('crystal_type' , crystal_type)
        call PyOptions%get('adjust_charge', adjust_mode)

        write(*,'(A)') ' Charges adjustment configurations'
        write(*,'(A,I12)') ' fragment natoms            : ', natoms_frag
        write(*,'(A,I12)') ' fragment nshells           : ', nshells_frag
        write(*,'(A,I12)') ' cluster natoms             : ', natoms_cluster
        write(*,'(A,I12)') ' cluster nshells            : ', nshells_cluster
        write(*,'(A,A)')   ' dimensionality             : ', dimensionality
        write(*,'(A,A)')   ' crystal_type               : ', trim(crystal_type)
        write(*,'(A,A)')   ' adjust mode                : ', trim(adjust_mode)

        ! nbqs_frag=0 and nbqs_cluster = 0
        call cluster_adjust_charge(natoms_cluster, 0, nshells_cluster, &
                                   natoms_frag   , 0, nshells_frag   , &
                                   ierror)

    endfunction cluster_adjust_charges
    integer(c_int) function cluster_calc_electrostatics(fieldPtr,     &
                                                        chemsh_debug) &
                   result(ierror)                                     &
                   bind(c, name='cluster_calc_electrostatics')

        implicit none

        type(PyType), intent(inout) :: fieldPtr
        integer     , intent(in)    :: chemsh_debug

        integer                     :: grid_mode
        integer                     :: natoms_frag, nshells_frag
        integer                     :: ngridpts

        debug = chemsh_debug

        allocate(pyField,source=ptr2dict(fieldPtr))

        call PyFrag%get('npts'   , natoms_frag)
        call PyFrag%get('nshells', nshells_frag)
        call PyOptions%get('grid_mode'  , grid_mode)
        call PyOptions%get('grid_points', ngridpts)

        ! nbqs=0
        call periodic_electrostatics(grid_mode, natoms_frag, 0, nshells_frag, ngridpts, ierror)

    endfunction cluster_calc_electrostatics
    integer(c_int) function cluster_fit_bqs(clusterPtr,   &
                                            fieldPtr  ,   &
                                            chemsh_debug) &
                   result(ierror)                         &
                   bind(c, name='cluster_fit_bqs')

        implicit none

        type(PyType)         , intent(in) :: clusterPtr
        type(PyType)         , intent(in) :: fieldPtr
        integer              , intent(in) :: chemsh_debug

        real(kind=8)                      :: radius_cluster, radius_active
        real(kind=8)                      :: bq_margin, bq_fit_weight, bq_layer
        integer                           :: bq_density
        integer                           :: natoms_frag, nshells_frag
        integer                           :: natoms_cluster, nshells_cluster

        integer                           :: field_npoints
        type(PyType), pointer             :: bqPtr

        debug = chemsh_debug

        ! YL: haven't implemented deallocation but do not forget it!
        if(associated(PyCluster)) deallocate(PyCluster)
        allocate(PyCluster, source=ptr2dict(clusterPtr))
        allocate(PyField  , source=ptr2dict(fieldPtr))
        call PyCluster%get('bqs', bqPtr)
        allocate(PyBQs    , source=ptr2dict(bqPtr))

        call PyField%get('npts', field_npoints)

        !Get the parameters for the routine from Python
        call PyFrag   %get('npts'          , natoms_frag)
        call PyFrag   %get('nshells'       , nshells_frag)
        call PyCluster%get('npts'          , natoms_cluster)
        call PyCluster%get('nshells'       , nshells_cluster)
        call PyOptions%get('radius_cluster', radius_cluster)
        call PyOptions%get('radius_active' , radius_active)
        call PyOptions%get('bq_margin'     , bq_margin)
        call PyOptions%get('bq_density'    , bq_density)
        call PyOptions%get('bq_layer'      , bq_layer)
        ! AJL, 12/06/18: Added so we can constrain the cluster charge
        call PyOptions%get('bq_fit_weight' , bq_fit_weight)

        write(*,'(A)') ' BQ fitting configurations'
        write(*,'(A,i12)')    ' fragment natoms            : ', natoms_frag
        write(*,'(A,i12)')    ' fragment nshells           : ', nshells_frag
        write(*,'(A,i12)')    ' cluster natoms             : ', natoms_cluster
        write(*,'(A,i12)')    ' cluster nshells            : ', nshells_cluster
        write(*,'(A,g20.12)') ' cluster radius             : ', radius_cluster
        write(*,'(A,g20.12)') ' active radius              : ', radius_active
        write(*,'(A,g20.12)') ' BQs margin                 : ', bq_margin
        write(*,'(A,g20.12)') ' BQs layer                  : ', bq_layer
        write(*,'(A,i12)')    ' BQs density                : ', bq_density
        write(*,'(A,g20.12)') ' fitting weight             : ', bq_fit_weight
        
        ! nbqs_frag=0 and nbqs_cluster=0
        call fit_bqs(natoms_cluster, 0, nshells_cluster, &
                     natoms_frag   , 0, nshells_frag   , &
                     radius_cluster, radius_active     , &
                     bq_margin, bq_density, bq_layer, bq_fit_weight, ierror)

    endfunction cluster_fit_bqs
    integer(c_int) function cluster_init(fragPtr   ,   &
                                         optionsPtr,   &
                                         chemsh_debug) &
                   result(ierror)                      &
                   bind(c, name='cluster_init')

        implicit none
        type(PyType)         , intent(in) :: fragPtr
        type(PyType)         , intent(in) :: optionsPtr
        integer              , intent(in) :: chemsh_debug

        type(PyType), pointer             :: cellPtr

        debug = chemsh_debug

        ! frag
        allocate(PyFrag, source=ptr2dict(fragPtr))

        ! options
        allocate(PyOptions, source=ptr2dict(optionsPtr))

        ! cell
        call PyFrag%get('cell', cellPtr)
        allocate(PyCell, source=ptr2dict(cellPtr))

        ierror = 0

    endfunction cluster_init
    integer(c_int) function cluster_finalise(chemsh_debug) &
                   result(ierror)                   &
                   bind(c, name='cluster_finalise')

        implicit none
 
        integer, intent(in) :: chemsh_debug

        debug = chemsh_debug

        deallocate(PyFrag, PyOptions, PyCell)

        ierror = 0

    endfunction cluster_finalise
    integer function get_ndims() result(ndims)

        implicit none

        call PyCell%get('ndimensions', ndims)

    endfunction get_ndims

endmodule ChemShellClusterModule
