!  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 CONSTRUCT code

subroutine get_frac_coords(periodic_frag, frac_coords, ierror)

    use ChemShellClusterModule
    
    integer, intent(inout) :: periodic_frag, ierror
    real   , intent(inout) :: frac_coords(:,:)

    print *, "### PLACEHOLDER: get_frac_coords"

endsubroutine get_frac_coords
integer function set_grid_electrostatics(ngridpts, cart_grid, pot, efield) result(ierror)

    integer              , intent(inout) :: ngridpts
    real(kind=8)         , intent(inout) :: cart_grid(3,ngridpts)
    real(kind=8), pointer, intent(inout) :: pot(:)
    real(kind=8), pointer, intent(inout) :: efield(:)

    print *, "### PLACEHOLDER: set_grid_electrostatics"

    ierror = 0

endfunction set_grid_electrostatics
! pass the values back to Python
integer function set_cluster_origin(origin) result(ierror)

    use ChemShellClusterModule

    real(kind=8), intent(in) :: origin(3)

    call PyOptions%set('origin', origin)

    ierror = 0

endfunction set_cluster_origin
! YL FIXME: not in use
subroutine get_cluster_origin(origin, ierror)

    use ChemShellClusterModule

    implicit none
    
    real(kind=8), intent(out) :: origin(3)
    integer     , intent(out) :: ierror

    ! TODO Store and recover the origin already chosen in cluster_cut
    ! TODO It should not be possible to change it at fitting stage!

    ! YL 15/11/2019: merged origin_atom and origin_point
    ! cluster cutting needs fractional origin, so this is wrong:
    call PyOptions%get('origin', origin(1:3))

    ierror = 0

endsubroutine get_cluster_origin
subroutine get_unit_cell(periodic_frag, a, b, c, alpha, beta, gamma, ierror)

    integer, intent(inout) :: periodic_frag
    real, intent(inout) :: a, b, c, alpha, beta, gamma
    integer, intent(inout) :: ierror

    print *, "### FORTRAN: get_unit_cell"

endsubroutine get_unit_cell
subroutine get_valence_matrix(no_species, valence, ierror)
    integer, intent(inout) :: no_species
    real   , intent(inout) :: valence(:,:)
    integer, intent(inout) :: ierror

    print *, "### FORTRAN: get_valence_matrix"

endsubroutine get_valence_matrix
subroutine set_cart_coords(periodic_frag, some_switch, cart_coords, ierror)
    integer, intent(inout) :: periodic_frag, some_switch
    real(kind=8), intent(inout) :: cart_coords(:,:)
    integer, intent(inout) :: ierror

    print *, "### FORTRAN: set_cart_coords"

endsubroutine set_cart_coords
subroutine get_cart_coords(switch, coords, ncores, nshells, ierror)

    use ChemShellClusterModule, only: PyFrag, PyCluster, PyShells
    use DL_PY2F
   
    implicit none

    integer     , intent(in)  :: switch
    integer     , intent(in)  :: ncores,nshells
    integer     , intent(out) :: ierror
    real(kind=8), intent(out) :: coords(3,ncores+nshells)

    integer                   :: i
    real(kind=8)              :: shellcoords(3,nshells)
    type(PyType), pointer     :: shellsPtr

!real(kind=8), pointer              :: twodimdbl(:,:)
! YL TODO: use twodimdbl to receive coordinates and shell positions, respectively

    selectcase(switch)
        case(1)
            !get coord matrix (FRAG_get_coord_matrix)
            !print *,'Extracting atom coordinates...'
            call PyFrag%get('xyz', coords(1:3,1:ncores))
            if(nshells > 0) then
                !print *,'Extracting shells..., switch =', switch
                call PyFrag%get('shells', shellsPtr)
                if(associated(PyShells)) deallocate(PyShells)
                allocate(PyShells, source=ptr2dict(shellsPtr))
                call PyShells%get('coords', shellcoords)
                coords(:,ncores+1:ncores+nshells) = shellcoords(:,1:nshells)
            endif 
        case(2)
            call PyCluster%get('xyz', coords(1:3,1:ncores))
            if(nshells > 0) then
                !print *,'Extracting shells..., switch =', switch
                call PyCluster%get('shells', shellsPtr)
                if(associated(PyShells)) deallocate(PyShells)
                allocate(PyShells, source=ptr2dict(shellsPtr))
                call PyShells%get('coords', shellcoords)
                coords(:,ncores+1:ncores+nshells) = shellcoords(:,1:nshells)
            endif
        case default
            print *,"Invalid mode switch in get_cart_coords"
            STOP
    endselect
!    deallocate(PyShells)
    ierror = 0

endsubroutine get_cart_coords
subroutine get_grid_coords(ngridpts, cart_grid, ierror)

    use ChemShellClusterModule

    implicit none

    integer, intent(inout) :: ngridpts
    real   , intent(inout) :: cart_grid(3,ngridpts)
    integer, intent(inout) :: ierror


    ierror =-1

    print *, "### FORTRAN: get_grid_coords"

    ierror = 0

endsubroutine get_grid_coords
subroutine set_cart_charges(dummy, cart_charges, ncent, ierror)

    use ChemShellClusterModule, only: PyCluster, PyShells
    use DL_PY2F

    implicit none  
  
    ! FIXME fragment mode missing
    integer, intent(in) :: dummy
    integer, intent(in) :: ncent
    real(kind=8), intent(in) :: cart_charges(ncent)
    integer, intent(inout) :: ierror
    integer  :: ncores, nshells
    type(PyType), pointer  :: shellsPtr

    ! core charges
    call PyCluster%get('npts', ncores)
    call PyCluster%set('charges', cart_charges(1:ncores))
    ! shell charges
    call PyCluster%get('nshells', nshells)
    if(nshells > 0) then
        call PyCluster%get("shells", shellsPtr)
        if(associated(PyShells)) deallocate(PyShells)
        allocate(PyShells, source=ptr2dict(shellsPtr))
        call PyShells%set("charges", cart_charges(ncores+1:ncent))
    endif
    ierror = 0

endsubroutine set_cart_charges
subroutine get_cart_charges(switch, charges, ncent, ierror)

    use ChemShellClusterModule, only: PyFrag, PyCluster, PyShells
    use DL_PY2F
    implicit none

    integer     , intent(in)  :: switch
    integer     , intent(in)  :: ncent
    real(kind=8), intent(out) :: charges(ncent)
    integer     , intent(out) :: ierror
    integer                   :: ncores, nshells
    type(PyType), pointer  :: shellsPtr

    selectcase(switch)
    case(1)
        ! cores
        call PyFrag%get('npts', ncores)
        call PyFrag%get('charges', charges(1:ncores))
        ! shells
        call PyFrag%get('nshells', nshells)
        if(nshells > 0) then
            call PyFrag%get("shells", shellsPtr)
            if(associated(PyShells)) deallocate(PyShells)
            allocate(PyShells, source=ptr2dict(shellsPtr))
            call PyShells%get("charges", charges(ncores+1:ncent))
        endif
    case(2)
        ! cores
        call PyCluster%get('npts', ncores)
        call PyCluster%get('charges', charges(1:ncores))
        ! shells
        call PyCluster%get('nshells', nshells)
        if(nshells > 0) then
            call PyCluster%get("shells", shellsPtr)
            if(associated(PyShells)) deallocate(PyShells)
            allocate(PyShells, source=ptr2dict(shellsPtr))
            call PyShells%get("charges", charges(ncores+1:ncent))
        endif
    case default
        print *,"Invalid mode switch in get_cart_coords"
        STOP
    endselect

    ierror = 0

endsubroutine get_cart_charges
subroutine get_xtal_map(cluster_frag, xtal_map, ncores,nshells,ierror)
    use ChemShellClusterModule
    use DL_PY2F

    implicit none
    
    integer, intent(in)    :: cluster_frag
    integer,intent(in)     :: nshells
    integer,intent(in)     :: ncores
    integer, intent(inout) :: xtal_map(ncores+nshells)
    integer, intent(inout) :: ierror

    type(PyType), pointer  :: shellsPtr
    
!   From interface.c: 
!  /* Notes:
!      - map was used in xtal to indicate which asymmetric cell sites
!        correspond to the resulting primitive cell sites. As there is
!        no symmetry, here it simply counts from 1 upwards 
!      - charges are not handled by frac_to_cart. Copy straight from 
!        xatoms, etc. instead
!  */
    
    call pyCluster%get('map',xtal_map(1:ncores))
    if(nshells > 0) then
       call PyCluster%get("shells", shellsPtr)
       if(associated(PyShells)) deallocate(PyShells)
       allocate(PyShells, source=ptr2dict(shellsPtr))
       call PyShells%get('map', xtal_map(ncores+1:ncores+nshells))
    endif
 
    ierror = 0

endsubroutine get_xtal_map
subroutine get_cell_vectors(periodic_frag, vect, ierror)

    use ChemShellClusterModule
    use DL_PY2F
 
    implicit none
   
    integer                , intent(inout) :: periodic_frag
    real(kind=8)           , intent(out)   :: vect(3,3)
    integer                , intent(out)   :: ierror
   
!    type(PyType)  , pointer                :: cellPtr 

! YL: it is odd that allocating PyCell fails here but succeeds in ClusterCreateModule%init()
!    call PyFrag%get('cell', cellPtr)
!    allocate(PyCell, source=ptr2dict(cellPtr))

    call PyCell%get('vectors', vect)

    ierror = 0

endsubroutine get_cell_vectors
subroutine set_cart_unit_cell(periodic_frag, vect, cart_coords, ierror)

    integer, intent(inout) :: periodic_frag
    real   , intent(inout) :: vect(3,3)
    real   , intent(inout) :: cart_coords(:,:)
    integer, intent(inout) :: ierror
    
    print *, "### FORTRAN: set_cart_unit_cell"

endsubroutine set_cart_unit_cell
subroutine get_point_core(ifrag, point_core, nshells, ierror)

    use ChemShellClusterModule, only : PyFrag, PyCluster, PyShells
    use DL_PY2F

    implicit none

    integer, intent(in)    :: ifrag
    integer, intent(in)    :: nshells
    integer, intent(inout) :: point_core(nshells)
    integer, intent(inout) :: ierror

    integer                :: i
    integer                :: shellarray(nshells)
    real(kind=8)       :: charges(nshells)
    type(PyType), pointer  :: shellsPtr

    ierror = 1
    
    selectcase(ifrag)
        case(1)
           ! periodic fragment
           call PyFrag%get("shells", shellsPtr)
        case(2)
           ! cluster
           call PyCluster%get("shells", shellsPtr)
        case default
            print *,"Invalid fragment number in get_point_core"
            STOP       
    endselect
     
    if(associated(PyShells)) deallocate(PyShells)
    allocate(PyShells, source=ptr2dict(shellsPtr))

    call PyShells%get("charges", charges)

    call PyShells%get("coreatoms", shellarray) 

    ! shift by 1 to Fortran convention 
    point_core(:) = shellarray(:) + 1

    ierror = 0

endsubroutine get_point_core
subroutine get_znum(periodic_frag, znum_core, ncores,ierror)

    use ChemShellClusterModule, only : PyFrag

    implicit none
    integer, intent(in)    :: periodic_frag
    integer, intent(inout) :: ierror
    integer,intent(in)     :: ncores
    integer(kind=8)        :: znums(ncores)
    integer,intent(inout)  :: znum_core(ncores)
    
    call PyFrag%get('znums',znums) 
    znum_core(:) = int(znums(:))
    ierror = 0 

endsubroutine get_znum
subroutine set_cluster_bqs(n_polar_bq,      &
                           n_lsq_bq,        &
                           r_polar_bq,      &
                           r_lsq_bq,        &
                           charge_polar_bq, &
                           charge_lsq_bq,   &
                           ierror)

    use ChemShellClusterModule
    use DL_PY2F

    implicit none

    integer, intent(in) :: n_polar_bq
    integer, intent(in) :: n_lsq_bq
    real(kind=8), intent(in) :: r_polar_bq(3,n_polar_bq)
    real(kind=8), intent(in) :: r_lsq_bq(3,n_lsq_bq)
    real(kind=8), intent(in) :: charge_polar_bq(n_polar_bq)
    real(kind=8) , intent(in) :: charge_lsq_bq(n_lsq_bq)
    real(kind=8) :: bq_array(3,n_polar_bq+n_lsq_bq)
    integer, intent(inout) :: ierror
    type(PyType), pointer     :: bqPtr
    integer :: nbqs

    ierror = 1

    call PyCluster%get('bqs', bqPtr)
    allocate(PyBQs, source=ptr2dict(bqPtr))
    call PyBQs%get('npts', nbqs)
    write(*, '(/1X,A,I12)') 'Number of point charges    : ', nbqs

    !Set the coordinates
    bq_array(:,1:n_polar_bq) = r_polar_bq(:,:)
    bq_array(:,n_polar_bq+1:nbqs) = r_lsq_bq(:,:)
    call PyBQs%set('xyz', bq_array)

    !set the charges
    bq_array(1,1:n_polar_bq) = charge_polar_bq(:)
    bq_array(1,n_polar_bq+1:nbqs) = charge_lsq_bq(:)    
    call PyBQs%set('charges', bq_array(1,:))

    deallocate(PyBQs)

    ierror = 0

endsubroutine set_cluster_bqs
subroutine set_electrostatics(periodic_frag, &
                              cart_coords,   &
                              pot_core,      &
                              pot_shel,      &
                              efield_core,   &
                              efield_shel,   &
                              ncores,nshells,ierror)

    use ChemShellClusterModule
    implicit none
    integer     , intent(in)    :: ncores,nshells
    integer     , intent(in)    :: periodic_frag
    real(kind=8), intent(in)    :: cart_coords(3,ncores+nshells)
    real(kind=8), intent(in)    :: pot_core(ncores)
    real(kind=8), intent(in)    :: pot_shel(nshells)
    real(kind=8), intent(in)    :: efield_core(3,ncores)
    real(kind=8), intent(in)    :: efield_shel(3,nshells)
    integer     , intent(inout) :: ierror
    integer                            :: npoints
    real(kind=8), allocatable      :: efield_tmp(:,:)
    real(kind=8), allocatable      :: pot_tmp(:)

    call pyField%get('npts', npoints)
    call pyField%set('xyz' , cart_coords)
    
    allocate(efield_tmp(1:3,1:npoints)) 
    allocate(pot_tmp(1:npoints))

    efield_tmp(1:3,1:ncores) = efield_core(1:3,1:ncores)
    pot_tmp(1:ncores)        = pot_core(1:ncores)

    if(nshells > 0) then
        efield_tmp(1:3,1+ncores:npoints) = efield_shel(1:3,1:nshells)
        pot_tmp(1+ncores:npoints)        = pot_shel(1:nshells)
    endif 

    call pyField%set('efield',efield_tmp)
    call pyField%set('pot',pot_tmp)

    deallocate(efield_tmp)
    deallocate(pot_tmp)
    ierror = 0

endsubroutine set_electrostatics
subroutine get_electrostatics(periodic_frag, &
                              pot_core,      &
                              pot_shel,      &
                              efield_core,   &
                              efield_shel,   &
                              ncores,nshells,ierror)

    use ChemShellClusterModule

    implicit none

    integer, intent(in)             :: periodic_frag
    integer, intent(in)             :: ncores,nshells
    real(kind=8)  , intent(out) :: pot_core(ncores)
    real(kind=8)  , intent(out) :: pot_shel(nshells)
    real(kind=8)  , intent(out) :: efield_core(3,ncores)
    real(kind=8)  , intent(out) :: efield_shel(3,nshells)
    integer       , intent(inout)          :: ierror

    integer                       :: i, j, npoints
    real(kind=8), allocatable   :: efield_tmp(:,:)
    real(kind=8), allocatable   :: pot_tmp(:)

    call pyField%get('npts', npoints)
    write(*, '(/A)')    ' Electrostatics'
    write(*, '(A,I12)') ' Number of points           : ', npoints
    if(nshells == 0) then
        call pyField%get('efield', efield_core(1:3,1:ncores))
        call pyField%get('pot'   , pot_core)
    else
        allocate(efield_tmp(1:3,1:npoints))
        allocate(pot_tmp(1:npoints))

        efield_tmp = 0.0
        pot_tmp    = 0.0

        call pyField%get('efield', efield_tmp)
        efield_core(1:3,1:ncores)  = efield_tmp(1:3,1:ncores)
        efield_shel(1:3,1:nshells) = efield_tmp(1:3,1+ncores:npoints)
        call pyField%get('pot'   , pot_tmp)
        pot_core(1:ncores)  = pot_tmp(1:ncores)
        pot_shel(1:nshells) = pot_tmp(1+ncores:npoints)

        deallocate(efield_tmp)
        deallocate(pot_tmp)

    endif

    ierror = 0 

endsubroutine get_electrostatics
subroutine periodic_neighbour_list(natoms_frag,          &
                                   coordination,    &
                                   point_neighbour, &
                                   v_neighbour,     &
                                   ierror)

    use ChemShellClusterModule

    implicit none
    
    integer, intent(in)    :: natoms_frag
    integer, intent(out)   :: coordination(natoms_frag)
    integer, intent(out)   :: point_neighbour(14,natoms_frag)
    real(kind=8), intent(inout) :: v_neighbour(3,14,natoms_frag)
    integer, intent(out)   :: ierror

    integer(kind=8), pointer   :: nneighbours(:)
    integer(kind=8), pointer   :: neighbours(:,:)
    real(kind=8)   , pointer   :: neighb_vectors(:,:)

    ierror = 1

    ! neighbours list
    allocate(nneighbours(natoms_frag))
    allocate(neighbours(MAXCONN, natoms_frag))
    allocate(neighb_vectors(3, MAXCONN*natoms_frag))

    call PyOptions%get('nneighbours', nneighbours)
    call PyOptions%get('neighbours', neighbours)
    call PyOptions%get('neighb_vectors', neighb_vectors)

    ! Neighbours are set at initialisation so simply read them in
    coordination    = nneighbours
    point_neighbour = neighbours
    v_neighbour     = reshape(neighb_vectors, (/3, 14, natoms_frag/))

    deallocate(nneighbours, neighbours, neighb_vectors)

    ierror = 0

endsubroutine periodic_neighbour_list
subroutine get_cluster_connectivity(ncores, connectivity, ierror)

    use ChemShellClusterModule

    implicit none

    integer, intent(in)  :: ncores    
    integer, intent(out) :: connectivity(14, ncores)
    integer, intent(out) :: ierror

    integer                :: i
 
    ! YL 05/05/2018: FIXME
    connectivity = 0
    call PyCluster%get('conn', connectivity) 

    ! Convert to construct Fortran convention
    connectivity = connectivity + 1
    
 !   do i = 1,ncores 
 !       print *,connectivity(:, i) 
 !   enddo

    ierror = 0

endsubroutine get_cluster_connectivity
subroutine setCluster(cart_cluster_cores, cart_cluster_shells,       &
                      charges_cluster_cores, charges_cluster_shells, &
                      map_cluster_cores, map_cluster_shells,         &
                      coreatoms, ierror)

    use ChemShellClusterModule
    use DL_PY2F

    implicit none

    real(kind=8)         , intent(in)  :: cart_cluster_cores(3,*)
    real(kind=8)         , intent(in)  :: cart_cluster_shells(3,*)
    real(kind=8)         , intent(in)  :: charges_cluster_cores(*)
    real(kind=8)         , intent(in)  :: charges_cluster_shells(*)
    integer              , intent(in)  :: coreatoms(*)
    integer              , intent(in)  :: map_cluster_cores(*)
    integer              , intent(in)  :: map_cluster_shells(*)
    integer              , intent(out) :: ierror

    integer     , pointer              :: znums_frag(:), znums_cluster(:)
!    real(kind=8) :: displace(3,max_shel)

    type(PyType), pointer              :: shellsPtr
    integer                            :: i
    integer                            :: natoms_frag, nshells_frag, natoms_cluster, nshells_cluster

    ierror = -1

    call pyFrag%get('npts'   , natoms_frag)
    call pyFrag%get('nshells', nshells_frag)
    call pyCluster%get('npts'   , natoms_cluster)
    call pyCluster%get('nshells', nshells_cluster)

    if(nshells_frag > 0) then
        call PyCluster%get('shells', shellsPtr)
        if(associated(PyShells)) deallocate(PyShells)
        allocate(PyShells, source=ptr2dict(shellsPtr))
    endif

    ! map fragment's znums to cluster's znums
    allocate(znums_frag(natoms_frag))
    allocate(znums_cluster(natoms_cluster))
    call pyFrag%get('znums', znums_frag) 
    znums_cluster(1:natoms_cluster) = znums_frag(map_cluster_cores(1:natoms_cluster))

    ! put the coordinates into the cluster fragment 
    call PyCluster%set('xyz'    , cart_cluster_cores(1:3,1:natoms_cluster))
    call PyCluster%set('charges', charges_cluster_cores(1:natoms_cluster))
    call PyCluster%set('znums'  , znums_cluster)
    call PyCluster%set('map'    , map_cluster_cores(1:natoms_cluster))

    if(nshells_frag > 0) then
        call PyShells%set('coords'   , cart_cluster_shells(1:3,1:nshells_cluster))
        ! YL 05/05/2018: do NOT set displacements for they will be calculated in Python
!        call PyShells%set('displace' , displace(1:3,1:nshells_cluster))
        call PyShells%set('charges'  , charges_cluster_shells(1:nshells_cluster))
        call PyShells%set('coreatoms', coreatoms(1:nshells_cluster)-1)
        call PyShells%set('map'      , map_cluster_shells(1:nshells_cluster))
    endif

    deallocate(znums_frag, znums_cluster)

    ierror = 0

endsubroutine setCluster
!  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: better call this module ClusterOptionsModule to store fragment and other information

!    private
!    public init, finalise
!    public get_cluster_size
!    public PyFrag, PyCluster, PyOptions, PyCell, PyField, PyBQs, PyShells
!!    public nneighbours, neighbours, neighb_vectors
!
!    integer(kind=8)            :: ncores
!!    integer(kind=8), pointer   :: nneighbours(:)
!!    integer(kind=8), pointer   :: neighbours(:,:)
!    integer(kind=8), parameter :: MAXCONN = 14
!!    real(kind=8)   , pointer   :: neighb_vectors(:,:)
!    type(dictType) , pointer   :: PyFrag, PyCluster, PyOptions, PyCell, PyField, PyBQs, PyShells


!    type, bind(c) :: frag_pass
        ! This type must match the frag_pass C struct on the python side
!        integer(c_int) :: nat
!        type(c_ptr)    :: ncores, nshells 
!    end type frag_pass
!    type :: fragment
     ! This is our Fortran fragment data structure
!        integer(c_int), dimension(:)  ,pointer  :: ncores
!        integer(c_int), dimension(:)  ,pointer  :: nshells
     !   real(c_double), dimension(:,:),pointer  :: pos
     !   real(c_double), dimension(:)  ,pointer  :: charge
!    end type fragment


!    contains


