! This file is copied and modified from QUANTUM ESPRESSO
! Kun Cao, Henry Lambert, Feliciano Giustino
 
!
! Copyright (C) 2011 Quantum ESPRESSO group
! This file is distributed under the terms of the
! GNU General Public License. See the file `License'
! in the root directory of the present distribution,
! or http://www.gnu.org/copyleft/gpl.txt .
!
!
!----------------------------------------------------------------------
MODULE dfile_autoname
!----------------------------------------------------------------------
  USE kinds,     ONLY : DP
  !
  PUBLIC :: dfile_name, dfile_generate_name, dfile_get_qlist
  !
  PRIVATE
  CHARACTER(len=12),PARAMETER :: dfile_directory_basename='.dfile_dir'
  !
CONTAINS
!----------------------------------------------------------------------
FUNCTION dfile_directory_file(basename, prefix)
  !----------------------------------------------------------------------
  IMPLICIT NONE
  CHARACTER(len=*),INTENT(in) :: basename
  CHARACTER(len=*),INTENT(in) :: prefix
  CHARACTER(len=512) :: dfile_directory_file
  dfile_directory_file = TRIM(prefix)//"."// &
                           TRIM(basename)//dfile_directory_basename
  RETURN
  !----------------------------------------------------------------------
END FUNCTION dfile_directory_file
!----------------------------------------------------------------------
!
!----------------------------------------------------------------------
FUNCTION open_dfile_directory(basename, prefix)
  !----------------------------------------------------------------------
  USE io_files, ONLY : find_free_unit
  IMPLICIT NONE
  CHARACTER(len=*),INTENT(in) :: basename
  CHARACTER(len=*),INTENT(in) :: prefix   ! directory where to operate
  INTEGER :: open_dfile_directory
  INTEGER :: ios
  CHARACTER(len=256) :: filename
!  LOGICAL :: exst
  !
  filename = dfile_directory_file(basename, prefix)
  !print*, "opening dir:", TRIM(filename)
  open_dfile_directory = find_free_unit()
  !
!  INQUIRE( FILE = TRIM(filename), EXIST = exst )
  !IF(.not.exst) print*, "does not exist: >",TRIM(filename),"<"
#ifdef __XLF
  OPEN(UNIT  = open_dfile_directory, &
       ACCESS= 'sequential',           &
       POSITION='append',              &
       FILE  = TRIM(filename),         &
       FORM  ='formatted', status='unknown', iostat=ios)
#else
  OPEN(UNIT  = open_dfile_directory, &
       ACCESS= 'append',               &
       FILE  = TRIM(filename),         &
       FORM  ='formatted', status='unknown', iostat=ios)
#endif
  !
  IF(ios/=0) CALL errore('open_dfile_directory','Cannot open: '//TRIM(filename),ABS(ios))
  !
  RETURN
  !----------------------------------------------------------------------
END FUNCTION open_dfile_directory
!----------------------------------------------------------------------
!
!----------------------------------------------------------------------
FUNCTION scan_dfile_directory(iunit, xq, at, found, equiv) 
  !----------------------------------------------------------------------
  IMPLICIT NONE
  CHARACTER(len=256) :: scan_dfile_directory
  !
  REAL(DP),INTENT(in) :: xq(3)
  REAL(DP),INTENT(in) :: at(3,3)
  INTEGER,INTENT(in)  :: iunit
  LOGICAL,INTENT(out) :: found
  LOGICAL,INTENT(in),OPTIONAL :: equiv ! if .false. only look for exactly q
                                       ! if .true. any q+G is ok (default)
  !
  INTEGER  :: ios, iq_
  REAL(DP) :: xp(3), aq(3), ap(3)
  CHARACTER(len=256) :: xp_name
  REAL(DP),PARAMETER :: gam(3) = (/ 0._dp, 0._dp, 0._dp /), accept = 1.e-5_dp
  !
  LOGICAL :: equiv_
  LOGICAL,EXTERNAL :: eqvect
  !
  found=.false.
  scan_dfile_directory = ''
  !
  equiv_ = .true.
  IF(present(equiv)) equiv_ = equiv
  ! xq in crystal coordinates:
  aq = xq
  CALL cryst_to_cart (1,aq,at,-1)
  !
  REWIND(iunit)
  ios=0
  !
  SCAN_FILE : &
  DO WHILE(ios==0)
      READ(iunit,*,iostat=ios) xp, ap, iq_, xp_name
    ! ap = xp
    ! CALL cryst_to_cart (1,ap,at,-1)
    !
    IF (equiv_) THEN
      IF (eqvect(aq,ap,gam,accept) .and. ios==0) THEN
        found=.true.
        scan_dfile_directory = TRIM(ADJUSTL(xp_name))
        EXIT SCAN_FILE
      ENDIF
    ELSE
      IF ( ALL(ABS(ap-aq)<accept) ) THEN
        found=.true.
        scan_dfile_directory = TRIM(ADJUSTL(xp_name))
        EXIT SCAN_FILE
      ENDIF
    ENDIF
  ENDDO SCAN_FILE
  !
  RETURN
  !----------------------------------------------------------------------
END FUNCTION scan_dfile_directory
!----------------------------------------------------------------------
!
!----------------------------------------------------------------------
FUNCTION dfile_name(xq, at, name, prefix, generate, index_q, equiv)
  !----------------------------------------------------------------------
  ! automatically generate a name for fildrho file
  USE io_global,    ONLY : ionode
  IMPLICIT NONE
  ! function:
  CHARACTER(len=256) :: dfile_name
  ! input variables:
  INTEGER, INTENT(in)         :: index_q      ! index of the q-point  
  REAL(DP),INTENT(in)         :: xq(3)    ! the q point in cartesian axes
  REAL(DP),INTENT(in)         :: at(3,3)  ! the lattice vectors, to transform the q to crystal coords
  CHARACTER(len=*),INTENT(in) :: prefix   ! directory where to operate
  CHARACTER(len=*),INTENT(in) :: name     ! input fildrho
  LOGICAL,INTENT(in)          :: generate ! make a new name if not found
  LOGICAL,INTENT(in),OPTIONAL :: equiv    ! accept an equivalent point q+G instead of q
                                          ! (only used when generate=.false.)
  !
  INTEGER :: iunit = -1, ios
  LOGICAL :: found
  CHARACTER(len=256) :: basename
  REAL(DP) :: aq(3)
  !
  ! Only ionode returns a meaningful name, as only ionode should do the i/o
  IF(.not.ionode) THEN
    dfile_name = ' '
    RETURN
  ENDIF
  !
  IF(name(1:5) /= 'auto:') THEN
    dfile_name = name
    RETURN
  ENDIF
  !
  basename = TRIM(name(6:))
  !
  
  iunit = open_dfile_directory(basename, prefix)
  rewind(iunit)
 
  dfile_name = scan_dfile_directory(iunit, xq, at, found, equiv)

  CLOSE(iunit)
  
  !
  ! Return here if point was found
  !IF(found) print*, "xq found as ", TRIM(dfile_name)
  IF(found) RETURN
  !
  IF(.not.generate) THEN
    WRITE(*,'(7x,"Error: ",3f12.6)') xq
    WRITE(*,'(7x,"Error: ",a,2x,a)') TRIM(name), TRIM(prefix)
    CALL errore('dfile_name','Requested q vector not found @ '//TRIM(basename), 1)
  ENDIF
  !
  ! Make up a new name
  dfile_name = TRIM(dfile_generate_name(xq, at, basename))
  !
  ! Append the new name to the list
  iunit = open_dfile_directory(basename, prefix)
  aq = xq
  CALL cryst_to_cart (1,aq,at,-1)
  !
  WRITE(iunit,*,iostat=ios) xq, aq, index_q, TRIM(dfile_name)

  IF(ios/=0) CALL errore('dfile_name','Cannot write dfile_directory',1)
  CLOSE(iunit)
  !
  RETURN
  !----------------------------------------------------------------------
END FUNCTION dfile_name
!----------------------------------------------------------------------
!
!----------------------------------------------------------------------
SUBROUTINE dfile_get_qlist(xqs, nqs, name, prefix, index_q)
  !----------------------------------------------------------------------
  ! automatically generate a name for fildrho file
  USE io_global,    ONLY : ionode
  IMPLICIT NONE
  ! input variables:
  INTEGER,INTENT(in)          :: nqs        ! max number of points
  INTEGER                     :: index_q(nqs)    ! index of irreducible q
  REAL(DP),INTENT(out)        :: xqs(3,nqs) ! the q point in cartesian axes
  CHARACTER(len=*),INTENT(in) :: prefix     ! directory where to operate
  CHARACTER(len=*),INTENT(in) :: name       ! input fildrho
  !
  INTEGER :: iunit = -1, ios, iq
  REAL(DP) :: xq_dummy(3)
  CHARACTER(len=256) :: basename
  !
  ! Only ionode scans for the filename, and does NOT broadcast. The broadcast
  ! must be done outside, because here we do not know if all processors are
  ! calling this subroutine!
  !
  IF (.not. ionode) THEN
   xqs = 0._dp
   RETURN
  ENDIF
  !
  IF(name(1:5) == 'auto:') THEN
    basename = TRIM(name(6:))
  ELSE
    basename = TRIM(name)
  ENDIF
  !
  iunit = open_dfile_directory(basename, prefix)
  REWIND(iunit)
  !
  GET_Q_LOOP : &
  DO iq = 1,nqs
     READ(iunit,*,iostat=ios) xqs(:,iq), xq_dummy(:),index_q(iq)
     IF(ios/=0) THEN
       CALL errore('dfile_get_qlist', 'Error while reading q point', iq)
     ENDIF
  ENDDO &
  GET_Q_LOOP
  !
  CLOSE(iunit)
  !
  RETURN
  !----------------------------------------------------------------------
END SUBROUTINE dfile_get_qlist
!----------------------------------------------------------------------
!
!----------------------------------------------------------------------
FUNCTION dfile_generate_name(xq, at, name)
  !----------------------------------------------------------------------
  ! automatically generate a name for fildrho file
  IMPLICIT NONE
  ! function:
  CHARACTER(len=256) :: dfile_generate_name
  ! input variables:
  REAL(DP),INTENT(in)         :: xq(3)   ! the q point in cartesian axes
  REAL(DP),INTENT(in)         :: at(3,3) ! lattice vectors
  CHARACTER(len=*),INTENT(in) :: name    ! input fildrho
  ! local variables:
  REAL(DP) :: aq(3) ! xq in crystal axes
  !
  ! take xq to crystalline coordinates
  aq = xq
  CALL cryst_to_cart (1,aq,at,-1)
  !
  WRITE(dfile_generate_name, '(a,".",a,"_",a,"_",a)') TRIM(name), &
        TRIM(real2frac(aq(1))), TRIM(real2frac(aq(2))), TRIM(real2frac(aq(3)))
  !
  dfile_generate_name = TRIM(dfile_generate_name)
  !
  RETURN
  !
  !----------------------------------------------------------------------
END FUNCTION dfile_generate_name
!----------------------------------------------------------------------
!
! Convert a real number to a string containing fraction, using 'o' as a division sign, es.
! 0    --> "0"
! -4   --> "-4"
! 0.25 --> "1o4"
! -1.66666666667 -> "-5/3"
! 
!----------------------------------------------------------------------
FUNCTION real2frac(r) RESULT (f)
  !----------------------------------------------------------------------
  IMPLICIT NONE
  REAL(DP),INTENT(in) :: r
  CHARACTER(len=64) :: f
  !
  INTEGER :: d, n
  INTEGER,PARAMETER :: max_denominator = 48000
  REAL(DP),PARAMETER :: accept = 1.d-6
  CHARACTER(len=64) :: nc,dc
  !
  IF(max_denominator*accept*20>1._dp) &
    CALL errore('real2frac', 'incompatible parameters', 2)
  ! Threat zero and integers separately:
  IF (ABS(r)<accept) THEN
    f = '0'
    RETURN
  ELSE IF (ABS(r-NINT(r))<accept) THEN
    WRITE(nc, '(i16)') NINT(r)
    f = TRIM(ADJUSTL(nc))
    RETURN
  ENDIF
  !
  ! If the number is not an integer, recompose the fraction
  DO d = 1, max_denominator+1
    IF( ABS(r*d-NINT(r*d)) < accept ) EXIT
  ENDDO
  !
!  IF (d > max_denominator) CALL errore('real2frac', 'not a fraction', 1)
  ! 
  IF (d > max_denominator) THEN
     WRITE(*, '("WARNING from real2frac:",e25.15," is not a fraction, falling back to hex." )') r
     WRITE(f,'(Z64)') r
     f='0x'//TRIM(ADJUSTL(f))
     RETURN
  ENDIF
  !
  n = NINT(r*d)
  !
  WRITE(nc, '(i16)') n
  WRITE(dc, '(i16)') d
  !
  f = TRIM(ADJUSTL(nc))//'o'//TRIM(ADJUSTL(dc))
  !
  RETURN
  !
  !----------------------------------------------------------------------
END FUNCTION real2frac
!----------------------------------------------------------------------
!
!----------------------------------------------------------------------
END MODULE dfile_autoname
!----------------------------------------------------------------------