!
! Copyright (C) 2001-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 .
!
!
!---------------------------------------------------------------------
PROGRAM fhi2upf
  !---------------------------------------------------------------------
  !
  !     Convert a pseudopotential file in Fritz-Haber numerical format
  !     either ".cpi" (fhi88pp) or ".fhi" (abinit)
  !     to unified pseudopotential format (v.2)
  !     Adapted from the converter written by Andrea Ferretti
  !
  USE pseudo_types, ONLY : pseudo_upf, nullify_pseudo_upf, &
                           deallocate_pseudo_upf
  USE write_upf_v2_module, ONLY :  write_upf_v2
  !
  IMPLICIT NONE
  TYPE(pseudo_upf) :: upf
  CHARACTER(len=256) filein, fileout
  INTEGER :: ios
  !
  CALL get_file ( filein )
  IF ( TRIM(filein) == ' ') &
       CALL errore ('fhi2upf', 'usage: fhi2upf "file-to-be-converted"', 1)
  OPEN ( unit=1, file=filein, status = 'old', form='formatted', iostat=ios )
  IF ( ios /= 0) CALL errore ('cpmd2upf', 'file: '//trim(filein)//' not found', 2)
  !
  CALL read_fhi(1)
  !
  CLOSE (1)

  ! convert variables read from FHI format into those needed
  ! by the upf format - add missing quantities
  !
  CALL nullify_pseudo_upf ( upf )
  !
  CALL convert_fhi (upf)
  !
  ! write to file
  !
  fileout=trim(filein)//'.UPF'
  PRINT '(''Output PP file in UPF format :  '',a)', fileout
  OPEN(unit=2,file=fileout,status='unknown',form='formatted')
  !
  CALL write_upf_v2 (2, upf )
  !
  CLOSE (unit=2)
  CALL deallocate_pseudo_upf ( upf )
  !     ----------------------------------------------------------
  WRITE (6,"('Pseudopotential successfully written')")
  WRITE (6,"('Please review the content of the PP_INFO fields')")
  WRITE (6,"('*** Please TEST BEFORE USING !!! ***')")
  !     ----------------------------------------------------------
  !
  STOP
   STOP
END PROGRAM fhi2upf

MODULE fhi
  !
  ! All variables read from FHI file format
  !

  TYPE angular_comp
     real(8), POINTER     :: pot(:)
     real(8), POINTER     :: wfc(:)
     real(8), POINTER     :: grid(:)
     real(8)              :: amesh
     INTEGER             :: nmesh
     INTEGER             :: lcomp
  END TYPE angular_comp

  !------------------------------

  real(8) :: Zval           ! valence charge
  INTEGER      :: lmax          ! max l-component used

  LOGICAL      :: nlcc_
  real(8), ALLOCATABLE :: rho_atc(:) ! core  charge

  TYPE (angular_comp), POINTER :: comp(:)  ! PP numerical info
                                           ! (wfc, grid, potentials...)
  !------------------------------

  ! variables for the abinit header

  real(8) :: Zatom, Zion, r2well, rchrg, fchrg, qchrg
  INTEGER :: pspdat = 0, pspcod = 0 , pspxc = 0, lloc = -1, mmax = 0
  CHARACTER(len=256) :: info

END MODULE fhi
!
!     ----------------------------------------------------------
SUBROUTINE read_fhi(iunps)
  !     ----------------------------------------------------------
  !
  USE fhi
  IMPLICIT NONE
  INTEGER, PARAMETER    :: Nl=7  ! max number of l-components
  INTEGER :: iunps
  real(8) :: r, drhoc, d2rhoc
  !

  INTEGER               :: l, i, idum, mesh

  ! Start reading file

  READ(iunps,'(a)') info
  READ(info,*,iostat=i) Zval, l
  IF ( i /= 0 .or. zval <= 0.0 .or. zval > 100.0 ) THEN
     WRITE (6,'("Assuming abinit format. First line:",/,A)') TRIM(info)
     READ(iunps,*) Zatom, Zion, pspdat
     READ(iunps,*) pspcod, pspxc, lmax,lloc, mmax, r2well
     IF (pspcod /= 6) THEN
        WRITE (6,'("read_fhi: unknown PP type ",i1,"...stopping")') pspcod
        STOP
     ENDIF
     READ(iunps,*) rchrg, fchrg, qchrg
     !
     READ(iunps,*)
     READ(iunps,*)
     READ(iunps,*)
     !
     READ(iunps,*) Zval, l
     IF (abs(Zion-Zval) > 1.0d-8) THEN
        WRITE (6,'("read_fhi: Zval/Zion mismatch...stopping")')
        STOP
     ENDIF
     IF (l-1 /= lmax) THEN
        WRITE (6,'("read_fhi: lmax mismatch...stopping")')
        STOP
     ENDIF
  ELSE
     info = ' '
  ENDIF
  lmax = l - 1

  IF (lmax+1 > Nl) THEN
     WRITE (6,'("read_fhi: too many l-components...stopping")')
     STOP
  ENDIF

  DO i=1,10
     READ(iunps,*)     ! skipping 11 lines
  ENDDO

  ALLOCATE( comp(0:lmax) )

  DO l=0,lmax
     comp(l)%lcomp = l
     READ(iunps,*) comp(l)%nmesh, comp(l)%amesh
     IF (mmax > 0 .and. mmax /= comp(l)%nmesh) THEN
        WRITE (6,'("read_fhi: mismatched number of grid points...stopping")')
        STOP
     ENDIF
     IF ( l > 0) THEN
        IF (comp(l)%nmesh /= comp(0)%nmesh .or.   &
            comp(l)%amesh /= comp(0)%amesh )      THEN
           WRITE(6,'("read_fhi: different radial grids not allowed...stopping")')
           STOP
        ENDIF
     ENDIF
     mesh = comp(l)%nmesh
     ALLOCATE( comp(l)%wfc(mesh),            &      ! wave-functions
               comp(l)%pot(mesh),            &      ! potentials
               comp(l)%grid(mesh)            )      ! real space radial grid
     ! read the above quantities
     DO i=1,mesh
        READ(iunps,*) idum, comp(l)%grid(i),   &
                            comp(l)%wfc(i),    &
                            comp(l)%pot(i)
     ENDDO
  ENDDO

  nlcc_ =.false.
  ALLOCATE(rho_atc(comp(0)%nmesh))
  mesh = comp(0)%nmesh
  DO i=1,mesh
     READ(iunps,*,end=10, err=20) r, rho_atc(i), drhoc, d2rhoc
     IF ( abs( r - comp(0)%grid(i) ) > 1.d-6 ) THEN
        WRITE(6,'("read_fhi: radial grid for core charge? stopping")')
        STOP
     ENDIF
  ENDDO
  nlcc_ = .true.
  !     ----------------------------------------------------------
  WRITE (6,'(a)') 'Pseudopotential with NLCC successfully read'
  !     ----------------------------------------------------------
  RETURN
20 WRITE(6,'("read_fhi: error reading core charge, assuming no core charge")')
   WRITE(6,'("this error may be due to the presence of additional", &
 &           " lines at the end of file")')
10 CONTINUE
  !     ----------------------------------------------------------
  WRITE (6,'(a)') 'Pseudopotential without NLCC successfully read'
  !     ----------------------------------------------------------
  RETURN
  !
  STOP

END SUBROUTINE read_fhi

!     ----------------------------------------------------------
SUBROUTINE convert_fhi (upf)
  !     ----------------------------------------------------------
  !
  USE fhi
  USE pseudo_types, ONLY : pseudo_upf
  USE funct, ONLY : set_dft_from_name, get_iexch, get_icorr, get_igcx, get_igcc
  USE constants, ONLY : fpi
  !
  IMPLICIT NONE
  !
  TYPE(pseudo_upf) :: upf
  !
  real(8), PARAMETER :: rmax = 10.0d0
  real(8), ALLOCATABLE :: aux(:)
  real(8) :: vll
  CHARACTER (len=2):: label
  CHARACTER (len=2), EXTERNAL:: atom_name
  INTEGER :: l, i, ir, iv
  !
  upf%nv       = "2.0.1"
  upf%generated= "Generated using Fritz-Haber code"
  upf%author   = "Author: unknown"
  upf%date     = "Generation date: as well"
  IF (trim(info) /= ' ') THEN
     upf%comment = trim(info)
  ELSE
     upf%comment = 'Info: automatically converted from FHI format'
  ENDIF
  upf%rel = 'scalar'  ! just guessing 
  IF (nint(Zatom) > 0) THEN
     upf%psd = atom_name(nint(Zatom))
     IF (nint(Zatom) > 18) upf%rel = 'no' ! just guessing
  ELSE
     PRINT '("Atom name > ",$)'
     READ (5,'(a)') upf%psd
  ENDIF
  upf%typ = 'SL'
  upf%tvanp = .FALSE.
  upf%tpawp = .FALSE.
  upf%tcoulombp=.FALSE.
  upf%nlcc = nlcc_
  !
  IF (pspxc == 7) THEN
     upf%dft = 'SLA-PW'
  ELSE IF (pspxc == 11) THEN
     upf%dft = 'PBE'
  ELSE 
     IF (pspxc > 0) THEN
        PRINT '("DFT read from abinit file: ",i1)', pspxc
     ENDIF
     PRINT '("DFT > ",$)'
     READ (5,'(a)') upf%dft
  ENDIF
  !
  upf%zp   = Zval
  upf%etotps =0.0d0
  upf%ecutrho=0.0d0
  upf%ecutwfc=0.0d0
  !
  PRINT '("Confirm or modify l max, l loc (read:",2i3,") > ",$)', lmax, lloc
  READ (5,*) lmax, upf%lloc

  IF ( lmax == upf%lloc) THEN
     upf%lmax = lmax-1
  ELSE
     upf%lmax = lmax
  ENDIF
  upf%lmax_rho = 0
  upf%nwfc  = lmax+1
  !
  ALLOCATE( upf%els(upf%nwfc) )
  ALLOCATE( upf%oc(upf%nwfc) )
  ALLOCATE( upf%epseu(upf%nwfc) )
  ALLOCATE( upf%lchi(upf%nwfc) )
  ALLOCATE( upf%nchi(upf%nwfc) )
  ALLOCATE( upf%rcut_chi (upf%nwfc) )
  ALLOCATE( upf%rcutus_chi(upf%nwfc) )

  PRINT '("PPs in FHI format do not contain information on atomic valence (pseudo-)wavefunctions")'
  PRINT '("Provide the label and the occupancy for each atomic wavefunction used in the PP generation")'
  PRINT '("If unknown: list valence wfcts and occupancies for the atomic ground state ", &
         &"in increasing l order: s,p,d,f")'
  DO i=1, upf%nwfc
10   PRINT '("Wavefunction # ",i1,": label (e.g. 4s), occupancy > ",$)', i
     READ (5,*) label, upf%oc(i)
     READ (label(1:1),*, err=10) l
     upf%els(i)  = label
     upf%nchi(i)  = l
     IF ( label(2:2) == 's' .OR. label(2:2) == 'S') then
        l=0
     ELSE IF ( label(2:2) == 'p' .OR. label(2:2) == 'P') then
        l=1
     ELSE IF ( label(2:2) == 'd' .OR. label(2:2) == 'D') then
        l=2
     ELSE IF ( label(2:2) == 'f' .OR. label(2:2) == 'F') then
        l=3
     ELSE
        l=i-1
     END IF
     upf%lchi(i)  = l
     upf%rcut_chi(i)  = 0.0d0
     upf%rcutus_chi(i)= 0.0d0
     upf%epseu(i) = 0.0d0
  ENDDO

  upf%mesh = comp(0)%nmesh
  upf%dx   = log( comp(0)%amesh )
  upf%rmax = comp(0)%grid(upf%mesh)
  upf%xmin = log( comp(0)%grid(1)*Zatom )
  upf%zmesh= Zatom
  ALLOCATE(upf%rab(upf%mesh))
  ALLOCATE(upf%r(upf%mesh))
  upf%r(:) = comp(0)%grid
  upf%rab(:)=upf%r(:)*upf%dx

  ALLOCATE (upf%rho_atc(upf%mesh))
  IF (upf%nlcc) upf%rho_atc(:) = rho_atc(1:upf%mesh) / fpi

  ALLOCATE (upf%vloc(upf%mesh))
  ! the factor 2 converts from Hartree to Rydberg
  upf%vloc(:) = 2.d0*comp(lloc)%pot
  upf%rcloc = 0.0d0

  ALLOCATE(upf%vnl(upf%mesh,0:upf%lmax,1))
  DO l=0, upf%lmax
     upf%vnl(:,l,1) = 2.d0*comp(l)%pot(:)
  END DO

  ! calculate number of nonlocal projectors
  IF ( upf%lloc >= 0 .AND. upf%lloc <= upf%lmax ) THEN
     upf%nbeta= upf%lmax
  ELSE
     upf%nbeta= upf%lmax+1
  ENDIF

  IF (upf%nbeta > 0) THEN

     ALLOCATE(upf%els_beta(upf%nbeta) )
     ALLOCATE(upf%lll(upf%nbeta))
     ALLOCATE(upf%kbeta(upf%nbeta))
     iv=0  ! counter on beta functions
     DO i=1,upf%nwfc
        l=upf%lchi(i)
        IF (l/=upf%lloc) THEN
           iv=iv+1
           DO ir = upf%mesh,1,-1
              IF ( ABS ( upf%vnl(ir,l,1) - upf%vnl(ir,upf%lloc,1) ) > 1.0E-6 ) THEN
                 upf%kbeta(iv)=ir
                 exit
              ENDIF
           ENDDO
        ENDIF
     ENDDO
     ! the number of points used in the evaluation of integrals
     ! should be even (for simpson integration)
     DO i=1,upf%nbeta
        IF ( MOD (upf%kbeta(i),2) == 0 .AND. upf%kbeta(i) < upf%mesh) &
           upf%kbeta(i)=upf%kbeta(i)+1
     END DO
     upf%kkbeta = MAXVAL(upf%kbeta(:))
     ALLOCATE(upf%beta(upf%mesh,upf%nbeta))
     ALLOCATE(upf%dion(upf%nbeta,upf%nbeta))
     upf%beta(:,:) =0.d0
     upf%dion(:,:) =0.d0
     ALLOCATE(upf%rcut  (upf%nbeta))
     ALLOCATE(upf%rcutus(upf%nbeta))
     ALLOCATE(aux(upf%kkbeta))
     iv=0  ! counter on beta functions
     DO i=1,upf%nwfc
        l=upf%lchi(i)
        IF (l/=upf%lloc) THEN
           iv=iv+1
           upf%lll(iv)=l
           upf%els_beta(iv)=upf%els(i)
           DO ir=1,upf%kbeta(iv)
              ! the factor 2 converts from Hartree to Rydberg
              upf%beta(ir,iv) = 2.d0 * comp(l)%wfc(ir) * &
                   ( comp(l)%pot(ir) - comp(upf%lloc)%pot(ir) )
              aux(ir) = comp(l)%wfc(ir) * upf%beta(ir,iv)
           ENDDO
           upf%rcut  (iv) = upf%r(upf%kbeta(iv))
           upf%rcutus(iv) = 0.0
           CALL simpson(upf%kbeta(iv),aux,upf%rab,vll)
           upf%dion(iv,iv) = 1.0d0/vll
        ENDIF
     ENDDO
     DEALLOCATE(aux)
  ENDIF

  ALLOCATE (upf%chi(upf%mesh,upf%nwfc))
  DO i=1,upf%nwfc
     upf%chi(:,i) = comp(i-1)%wfc(:)
  ENDDO

  ALLOCATE (upf%rho_at(upf%mesh))
  upf%rho_at(:) = 0.d0
  DO i=1,upf%nwfc
     upf%rho_at(:) = upf%rho_at(:) + upf%oc(i) * upf%chi(:,i) ** 2
  ENDDO
  !     ----------------------------------------------------------
  WRITE (6,'(a)') 'Pseudopotential successfully converted'
  !     ----------------------------------------------------------
  RETURN
END SUBROUTINE convert_fhi