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

subroutine frac_to_cart(ndim, nsites, ierr)
  !
  ! Converts from fractional coordinates to cartesian coordinates
  !
  ! Limitation: assumes no symmetry (i.e. the fractional coordinates describe 
  ! the full primitive unit cell, not the asymmetric unit)
  !
  implicit none

  integer, intent(in) :: ndim, nsites
  integer, intent(out) :: ierr
  integer :: isite

  double precision :: a, b, c, alpha, beta, gamma, dum
  double precision, dimension(3, 3) :: vect
  double precision, dimension(3, nsites) :: frac_coords, cart_coords

  integer, parameter :: periodic_frag = 1, cluster_frag = 2

  ierr = 0

  ! Get unit cell parameters from fragment
  ! a, b, c in atomic units
  ! alpha, beta, gamma in radians
  call get_unit_cell(periodic_frag, a, b, c, alpha, beta, gamma, ierr)
  if (ierr .ne. 0) return

  ! Call construct routines to generate the cell vectors
  vect = 0.0d0
  select case (ndim)
  case (1)
     vect(1, 1) = a
  case (2)
     call lattice_cryst_2(a, b, gamma, vect(1:2, 1:2), dum)
  case (3)
     call lattice_cryst_3(a, b, c, alpha, beta, gamma, vect, dum)
  case default
     write(6,*) "frac_to_cart: error, ndim=", ndim
     ierr = 1
     return
  end select

  ! Get fractional coordinates from fragment
  ! This is received as a single block of coordinates containing
  ! all atoms, point charges (bqs) and shells. There is no need to 
  ! differentiate between these as the same transformation is applied to all 
  ! sites.


  ! YL FIXME: get_frac_coords() is yet to implement
  call get_frac_coords(periodic_frag, frac_coords, ierr)  
  if (ierr .ne. 0) return

  ! Convert fractional coordinates to cartesian
  do isite = 1, nsites
     ! Only the periodic dimensions are fractional
     select case (ndim)
     case (1)
        cart_coords(1, isite) = frac_coords(1, isite) * a
        cart_coords(2, isite) = frac_coords(2, isite)
        cart_coords(3, isite) = frac_coords(3, isite)
     case (2)
        cart_coords(1, isite) = vect(1,1) * frac_coords(1, isite) + vect(2,1) * frac_coords(2, isite) 
        cart_coords(2, isite) = vect(1,2) * frac_coords(1, isite) + vect(2,2) * frac_coords(2, isite) 
        cart_coords(3, isite) = frac_coords(3, isite)
     case (3)
        cart_coords(1, isite) = vect(1,1) * frac_coords(1, isite) &
             + vect(2,1) * frac_coords(2, isite) + vect(3,1) * frac_coords(3, isite)
        cart_coords(2, isite) = vect(1,2) * frac_coords(1, isite) &
             + vect(2,2) * frac_coords(2, isite) + vect(3,2) * frac_coords(3, isite)
        cart_coords(3, isite) = vect(1,3) * frac_coords(1, isite) &
             + vect(2,3) * frac_coords(2, isite) + vect(3,3) * frac_coords(3, isite)
     case default
        write(6,*) "frac_to_cart: error, ndim=", ndim
        ierr = 1
        return
     end select
  enddo

  ! Look for equivalent sites closer to the origin
  call find_equiv_coords(ndim, nsites, vect, cart_coords, ierr)
  if (ierr .ne. 0) return

  ! Give cartesian coordinates to fragment
  call set_cart_unit_cell(periodic_frag, vect, cart_coords, ierr)
  if (ierr .ne. 0) return

end subroutine frac_to_cart


subroutine centre_cart_coords(ndim, nsites, ierr)
  !
  ! Wrapper around routine to find equivalent cartesian sites closer
  ! to the origin
  !
  implicit none

  integer, intent(in) :: ndim, nsites
  integer, intent(out) :: ierr
  integer :: isite

  double precision, dimension(3, 3) :: vect
  double precision, dimension(3, nsites) :: cart_coords

  integer, parameter :: periodic_frag = 1, cluster_frag = 2

  ierr = 0

  ! Get unit cell vectors from fragment
  vect = 0.0d0
  call get_cell_vectors(periodic_frag, vect, ierr)
  if (ierr .ne. 0) return  

  ! Get cartesian coordinates from fragment
  ! This is received as a single block of coordinates containing
  ! all atoms, point charges (bqs) and shells. There is no need to 
  ! differentiate between these as the same transformation is applied to all 
  ! sites.
  ! Note we read and write cart coords directly as the shells have 
  ! not yet been assigned to relative positions
  call get_cart_coords(periodic_frag, 0, cart_coords, ierr)  
  if (ierr .ne. 0) return

  ! Look for equivalent sites closer to the origin
  call find_equiv_coords(ndim, nsites, vect, cart_coords, ierr)
  if (ierr .ne. 0) return

  ! Give cartesian coordinates to fragment
  call set_cart_coords(periodic_frag, 0, cart_coords, ierr)
  if (ierr .ne. 0) return

end subroutine centre_cart_coords


subroutine find_equiv_coords(ndim, nsites, vect, cart_coords, ierr)
  !
  ! Following xtal, look for equivalent sites closer to the origin
  ! Note that we do not worry here about keeping shells associated with cores
  ! In c_create, this is handled afterwards by a call to FRAG_assign_shells_to_atoms
  !
  implicit none

  integer, intent(in) :: ndim, nsites
  double precision, dimension(3, 3), intent(in) :: vect
  double precision, dimension(3, nsites), intent(inout) :: cart_coords 
  integer, intent(out) :: ierr

  integer :: isite
  integer :: ap, bp, cp, loopa, loopb, loopc

  double precision, parameter :: search_tol = 1.0d-4
  double precision :: transx, transy, transz, bestx, besty, bestz, best_r2, r2

  ierr = 0

  ap = 0
  bp = 0
  cp = 0
  ! Only search over periodic dimensions
  select case (ndim)
  case (1)
     ap = 1
  case (2)
     ap = 1
     bp = 1
  case (3)
     ap = 1
     bp = 1
     cp = 1
  case default
     write(6,*) "find_equiv_coords: error, ndim=", ndim
     ierr = 1
     return
  end select
  do isite = 1, nsites
     ! Refinement - calculate the untranslated position first so that there is 
     ! no unecessary switching to equivalent coordinates.
     ! Note the results may now not exactly reproduce the old xtal code
     ! (though they will always be equivalent).
     bestx = cart_coords(1, isite)
     besty = cart_coords(2, isite)
     bestz = cart_coords(3, isite)
     r2 = bestx*bestx + besty*besty + bestz*bestz
     best_r2 = r2 - search_tol
     ! Now check possible translated positions
     do loopa = -ap, ap
        do loopb = -bp, bp
           do loopc = -cp, cp
              transx = cart_coords(1, isite) + loopa*vect(1, 1) + loopb*vect(2, 1) + loopc*vect(3, 1)
              transy = cart_coords(2, isite) + loopa*vect(1, 2) + loopb*vect(2, 2) + loopc*vect(3, 2)
              transz = cart_coords(3, isite) + loopa*vect(1, 3) + loopb*vect(2, 3) + loopc*vect(3, 3)
              r2 = transx*transx + transy*transy + transz*transz
              if (r2 .lt. best_r2) then
                 bestx = transx
                 besty = transy
                 bestz = transz
                 ! Tolerance term here is again following xtal. 
                 best_r2 = r2 - search_tol
              end if
           end do ! loopc
        end do ! loopb
     end do ! loopa
     cart_coords(1, isite) = bestx
     cart_coords(2, isite) = besty
     cart_coords(3, isite) = bestz
  end do ! isite

end subroutine find_equiv_coords
