MODULE casino_pp

  !
  ! All variables read from CASINO file format
  !
  ! trailing underscore means that a variable with the same name
  ! is used in module 'upf' containing variables to be written
  !
  USE upf_kinds, ONLY : dp

  CHARACTER(len=20) :: dft_
  CHARACTER(len=2)  :: psd_
  REAL(dp) :: zp_
  INTEGER nlc, nnl, lmax_, lloc, nchi, rel_
  LOGICAL :: numeric, bhstype, nlcc_
  CHARACTER(len=2), ALLOCATABLE :: els_(:)
  REAL(dp) :: zmesh
  REAL(dp) :: xmin      = -7.0_dp
  REAL(dp) :: dx        = 20.0_dp/1500.0_dp
  REAL(dp) :: tn_prefac = 0.75E-6_dp
  LOGICAL  :: tn_grid   = .true.


  REAL(dp), ALLOCATABLE::  r_(:)
  INTEGER :: mesh_

  REAL(dp), ALLOCATABLE::  vnl(:,:)
  INTEGER, ALLOCATABLE:: lchi_(:), nns_(:)
  REAL(dp), ALLOCATABLE:: chi_(:,:),  oc_(:)

CONTAINS
  !
  !     ----------------------------------------------------------
  SUBROUTINE read_casino(iunps,nofiles, waveunit)
    !     ----------------------------------------------------------
    !
    !     Reads in a CASINO tabulated pp file and it's associated
    !     awfn files. Some basic processing such as removing the
    !     r factors from the potentials is also performed.


    USE upf_kinds, ONLY : dp
    IMPLICIT NONE
    TYPE :: wavfun_list
       INTEGER :: occ,eup,edwn, nquant, lquant
       CHARACTER(len=2) :: label
       REAL(dp), ALLOCATABLE :: wavefunc(:)
       TYPE (wavfun_list), POINTER :: p

    END TYPE wavfun_list

    TYPE :: channel_list
       INTEGER :: lquant
       REAL(dp), ALLOCATABLE :: channel(:)
       TYPE (channel_list), POINTER :: p

    END TYPE channel_list


    TYPE (channel_list), POINTER :: phead
    TYPE (channel_list), POINTER :: pptr
    TYPE (channel_list), POINTER :: ptail

    TYPE (wavfun_list), POINTER :: mhead
    TYPE (wavfun_list), POINTER :: mptr
    TYPE (wavfun_list), POINTER :: mtail

    INTEGER :: iunps, nofiles, ios
    !
    LOGICAL :: groundstate, found
    CHARACTER(len=2) :: label, rellab

    INTEGER :: l, i, ir, nb, gsorbs, j,k,m,tmp, lquant, orbs, nquant
    INTEGER, ALLOCATABLE :: gs(:,:)
    INTEGER, INTENT(in) :: waveunit(nofiles)

    NULLIFY (  mhead, mptr, mtail )
    dft_ = 'HF'   !Hardcoded at the moment should eventually be HF anyway

    nlc = 0              !These two values are always 0 for numeric pps
    nnl = 0              !so lets just hard code them

    nlcc_ = .false.       !Again these two are alwas false for CASINO pps
    bhstype = .false.



    READ(iunps,'(a2,35x,a2)') rellab, psd_
    READ(iunps,*)
    IF ( rellab == 'DF' ) THEN
       rel_=1
    ELSE
       rel_=0
    ENDIF

    READ(iunps,*) zmesh,zp_  !Here we are reading zmesh (atomic #) and
    DO i=1,3                 !zp_ (pseudo charge)
       READ(iunps,*)
    ENDDO
    READ(iunps,*) lloc               !reading in lloc
    IF ( zp_<=0d0 ) &
         CALL upf_error( 'read_casino','Wrong zp ',1 )
    IF ( lloc>3.or.lloc<0 ) &
         CALL upf_error( 'read_casino','Wrong lloc ',1 )


    !
    !    compute the radial mesh
    !

    DO i=1,3
       READ(iunps,*)
    ENDDO
    READ(iunps,*) mesh_   !Reading in total no. of mesh points


    ALLOCATE(  r_(mesh_))

    READ(iunps,*)
    DO i=1,mesh_
       READ(iunps,*) r_(i)
    ENDDO


    ! Read in the different channels of V_nl
    ALLOCATE(phead)
    ptail => phead
    pptr  => phead

    ALLOCATE( pptr%channel(mesh_) )
    READ(iunps, '(15x,I1,7x)') l
    pptr%lquant=l
    READ(iunps, *)  (pptr%channel(ir),ir=1,mesh_)


    DO
       READ(iunps, '(15x,I1,7x)', IOSTAT=ios) l

       IF (ios /= 0 ) THEN
          exit
       ENDIF

       ALLOCATE(pptr%p)
       pptr=> pptr%p
       ptail=> pptr
       ALLOCATE( pptr%channel(mesh_) )
       pptr%lquant=l
       READ(iunps, *)  (pptr%channel(ir),ir=1,mesh_)

    ENDDO

    !Compute the number of channels read in.
    lmax_ =-1
    pptr => phead
    DO
       IF ( .not. associated(pptr) )exit
       lmax_=lmax_+1

       pptr =>pptr%p
    ENDDO

    ALLOCATE(vnl(mesh_,0:lmax_))
    i=0
    pptr => phead
    DO
       IF ( .not. associated(pptr) )exit
       !         lchi_(i) = pptr%lquant

       DO ir=1,mesh_
          vnl(ir,i) = pptr%channel(ir)
       ENDDO
       DEALLOCATE( pptr%channel )
       pptr =>pptr%p
       i=i+1
    ENDDO

    !Clean up the linked list (deallocate it)
    DO
       IF ( .not. associated(phead) )exit
       pptr => phead
       phead => phead%p
       DEALLOCATE( pptr )
    ENDDO

    DO l = 0, lmax_
       DO ir = 1, mesh_
          vnl(ir,l) = vnl(ir,l)/r_(ir) !Removing the factor of r CASINO has
       ENDDO
       ! correcting for possible divide by zero
       IF ( r_(1) == 0 ) THEN
          vnl(1,l) = 0
       ENDIF
    ENDDO

    ALLOCATE(mhead)
    mtail => mhead

    mptr => mhead

    NULLIFY(mtail%p)
    groundstate=.true.
    DO j=1,nofiles

       DO i=1,4

          READ(waveunit(j),*)
       ENDDO

       READ(waveunit(j),*) orbs

       IF ( groundstate ) THEN

          ALLOCATE( gs(orbs,3) )

          gs = 0
          gsorbs = orbs
       ENDIF

       DO i=1,2
          READ(waveunit(j),*)
       ENDDO

       READ(waveunit(j),*) mtail%eup, mtail%edwn
       READ(waveunit(j),*)

       DO i=1,mtail%eup+mtail%edwn
          READ(waveunit(j),*) tmp, nquant, lquant

          IF ( groundstate ) THEN
             found = .true.

             DO m=1,orbs

                IF ( (nquant==gs(m,1) .and. lquant==gs(m,2)) ) THEN
                   gs(m,3) = gs(m,3) + 1
                   exit
                ENDIF

                found = .false.

             ENDDO

             IF (.not. found ) THEN

                DO m=1,orbs

                   IF ( gs(m,1) == 0 ) THEN
                      gs(m,1) = nquant
                      gs(m,2) = lquant
                      gs(m,3) = 1

                      exit
                   ENDIF

                ENDDO

             ENDIF

          ENDIF

       ENDDO

       READ(waveunit(j),*)
       READ(waveunit(j),*)

       DO i=1,mesh_
          READ(waveunit(j),*)
       ENDDO

       DO k=1,orbs
          READ(waveunit(j),'(13x,a2)', err=300) label
          READ(waveunit(j),*) tmp, nquant, lquant

          IF ( .not. groundstate ) THEN
             found = .false.

             DO m = 1,gsorbs

                IF ( nquant == gs(m,1) .and. lquant == gs(m,2) ) THEN
                   found = .true.
                   exit
                ENDIF
             ENDDO
             mptr => mhead
             DO
                IF ( .not. associated(mptr) )exit
                IF ( nquant == mptr%nquant .and. lquant == mptr%lquant ) found = .true.
                mptr =>mptr%p
             ENDDO
             IF ( found ) THEN
                DO i=1,mesh_
                   READ(waveunit(j),*)
                ENDDO

                CYCLE
             ENDIF
          ENDIF
            IF ( allocated(mtail%wavefunc) ) THEN
                ALLOCATE(mtail%p)
                mtail=>mtail%p
                NULLIFY(mtail%p)
                ALLOCATE( mtail%wavefunc(mesh_) )
             ELSE
                ALLOCATE( mtail%wavefunc(mesh_) )
             ENDIF
             mtail%label = label
             mtail%nquant = nquant
             mtail%lquant = lquant


             READ(waveunit(j), *, err=300) (mtail%wavefunc(ir),ir=1,mesh_)
          ENDDO
          groundstate = .false.
       ENDDO

       nchi =0
       mptr => mhead
       DO
          IF ( .not. associated(mptr) )exit
          nchi=nchi+1

          mptr =>mptr%p
       ENDDO

       ALLOCATE(lchi_(nchi), els_(nchi), nns_(nchi))
       ALLOCATE(oc_(nchi))
       ALLOCATE(chi_(mesh_,nchi))
       oc_ = 0

       !Sort out the occupation numbers
       DO i=1,gsorbs
          oc_(i)=gs(i,3)
       ENDDO
       DEALLOCATE( gs )

       i=1
       mptr => mhead
       DO
          IF ( .not. associated(mptr) )exit
          nns_(i) = mptr%nquant
          lchi_(i) = mptr%lquant
          els_(i) = mptr%label

          DO ir=1,mesh_

             chi_(ir:,i) = mptr%wavefunc(ir)
          ENDDO
          DEALLOCATE( mptr%wavefunc )
          mptr =>mptr%p
          i=i+1
       ENDDO

       !Clean up the linked list (deallocate it)
       DO
          IF ( .not. associated(mhead) )exit
          mptr => mhead
          mhead => mhead%p
          DEALLOCATE( mptr )
       ENDDO


       !     ----------------------------------------------------------
       WRITE (0,'(a)') 'Pseudopotential successfully read'
       !     ----------------------------------------------------------
       RETURN

300    CALL upf_error('read_casino','pseudo file is empty or wrong',1)

     END SUBROUTINE read_casino

     !     ----------------------------------------------------------
     SUBROUTINE convert_casino(upf_out)
       !     ----------------------------------------------------------
       USE upf_kinds,    ONLY : dp
       USE pseudo_types, ONLY : pseudo_upf

       IMPLICIT NONE

       TYPE(pseudo_upf), INTENT(inout)       :: upf_out

       REAL(dp), ALLOCATABLE :: aux(:)
       REAL(dp) :: vll
       INTEGER :: kkbeta, l, iv, ir, i, nb

       !
       upf_out%nv       = "2.0.1"
       upf_out%tvanp    = .false.
       upf_out%tpawp    = .false.
       upf_out%tcoulombp= .false.
       upf_out%has_so   = .false.
       upf_out%has_wfc  = .false.
       upf_out%has_gipaw= .false.
       upf_out%paw_as_gipaw = .false.
       !
       WRITE(upf_out%generated, '("From a Trail & Needs tabulated &
            &PP for CASINO")')
       WRITE(upf_out%author,'("unknown")')
       WRITE(upf_out%date,'("unknown")')
       upf_out%comment = 'Info: automatically converted from CASINO &
            &Tabulated format'

       IF (rel_== 0) THEN
          upf_out%rel = 'no'
       ELSEIF (rel_==1 ) THEN
          upf_out%rel = 'scalar'
       ELSE
          upf_out%rel = 'full'
       ENDIF

       IF (xmin == 0 ) THEN
          xmin= log(zmesh * r_(2) )
       ENDIF

       ! Allocate and assign the radial grid

       upf_out%mesh  = mesh_
       upf_out%zmesh = zmesh
       upf_out%dx    = dx
       upf_out%xmin  = xmin

       ALLOCATE(upf_out%rab(upf_out%mesh))
       ALLOCATE(  upf_out%r(upf_out%mesh))

       upf_out%r = r_
       DEALLOCATE( r_ )

       upf_out%rmax = maxval(upf_out%r)


       !
       ! subtract out the local part from the different
       ! potential channels
       !

       DO l = 0, lmax_
          IF ( l/=lloc ) vnl(:,l) = vnl(:,l) - vnl(:,lloc)
       ENDDO

       ALLOCATE (upf_out%vloc(upf_out%mesh))
       upf_out%vloc(:) = vnl(:,lloc)


       ! Compute the derivatives of the grid. The Trail and Needs
       ! grids use r(i) = (tn_prefac / zmesh)*( exp(i*dx) - 1 ) so
       ! must be treated differently to standard QE grids.

       IF ( tn_grid ) THEN
          DO ir = 1, upf_out%mesh
             upf_out%rab(ir) = dx * ( upf_out%r(ir) + tn_prefac / zmesh )
          ENDDO
       ELSE
          DO ir = 1, upf_out%mesh
             upf_out%rab(ir) = dx  * upf_out%r(ir)
          ENDDO
       ENDIF


       !
       !    compute the atomic charges
       !
       ALLOCATE (upf_out%rho_at(upf_out%mesh))
       upf_out%rho_at(:) = 0.d0

       DO nb = 1, nchi
          IF( oc_(nb)/=0.d0) THEN
             upf_out%rho_at(:) = upf_out%rho_at(:) +&
               &  oc_(nb)*chi_(:,nb)**2
          ENDIF
       ENDDO

       ! This section deals with the pseudo wavefunctions.
       ! These values are just given directly to the pseudo_upf structure
       upf_out%nwfc  = nchi

       ALLOCATE( upf_out%oc(upf_out%nwfc), upf_out%epseu(upf_out%nwfc) )
       ALLOCATE( upf_out%lchi(upf_out%nwfc), upf_out%nchi(upf_out%nwfc) )
       ALLOCATE( upf_out%els(upf_out%nwfc) )
       ALLOCATE( upf_out%rcut_chi(upf_out%nwfc) )
       ALLOCATE( upf_out%rcutus_chi (upf_out%nwfc) )

       DO i=1, upf_out%nwfc
          upf_out%nchi(i)  = nns_(i)
          upf_out%lchi(i)  = lchi_(i)
          upf_out%rcut_chi(i)  = 0.0d0
          upf_out%rcutus_chi(i)= 0.0d0
          upf_out%oc (i)   = oc_(i)
          upf_out%els(i) = els_(i)
          upf_out%epseu(i) = 0.0d0
       ENDDO
       DEALLOCATE (lchi_, oc_, nns_)

       upf_out%psd = psd_
       upf_out%typ = 'NC'
       upf_out%nlcc = nlcc_
       upf_out%zp = zp_
       upf_out%etotps = 0.0d0
       upf_out%ecutrho=0.0d0
       upf_out%ecutwfc=0.0d0
       upf_out%lloc=lloc

       IF ( lmax_ == lloc) THEN
          upf_out%lmax = lmax_-1
       ELSE
          upf_out%lmax = lmax_
       ENDIF
       upf_out%nbeta = lmax_

       ALLOCATE ( upf_out%els_beta(upf_out%nbeta) )
       ALLOCATE ( upf_out%rcut(upf_out%nbeta) )
       ALLOCATE ( upf_out%rcutus(upf_out%nbeta) )

       upf_out%rcut=0.0d0
       upf_out%rcutus=0.0d0
       upf_out%dft =dft_


       IF (upf_out%nbeta > 0) THEN

          ALLOCATE(upf_out%kbeta(upf_out%nbeta), upf_out%lll(upf_out%nbeta))
          upf_out%kkbeta=upf_out%mesh
          DO ir = 1,upf_out%mesh
             IF ( upf_out%r(ir) > upf_out%rmax ) THEN
                upf_out%kkbeta=ir
                exit
             ENDIF
          ENDDO

          ! make sure kkbeta is odd as required for simpson
          IF(mod(upf_out%kkbeta,2) == 0) upf_out%kkbeta=upf_out%kkbeta-1
          upf_out%kbeta(:) = upf_out%kkbeta
          ALLOCATE(aux(upf_out%kkbeta))
          ALLOCATE(upf_out%beta(upf_out%mesh,upf_out%nbeta))
          ALLOCATE(upf_out%dion(upf_out%nbeta,upf_out%nbeta))

          upf_out%dion(:,:) =0.d0

          iv=0
          DO i=1,upf_out%nwfc
             l=upf_out%lchi(i)
             IF (l/=upf_out%lloc) THEN
                iv=iv+1
                upf_out%els_beta(iv)=upf_out%els(i)
                upf_out%lll(iv)=l
                DO ir=1,upf_out%kkbeta

                   upf_out%beta(ir,iv)=chi_(ir,i)*vnl(ir,l)
                   aux(ir) = chi_(ir,i)**2*vnl(ir,l)

                ENDDO
                CALL simpson(upf_out%kkbeta,aux,upf_out%rab,vll)
                upf_out%dion(iv,iv) = 1.0d0/vll
             ENDIF

             IF(iv >= upf_out%nbeta) exit  ! skip additional pseudo wfns
          ENDDO


          DEALLOCATE (vnl, aux)

          !
          !   redetermine ikk2
          !
          DO iv=1,upf_out%nbeta
             upf_out%kbeta(iv)=upf_out%kkbeta
             DO ir = upf_out%kkbeta,1,-1
                IF ( abs(upf_out%beta(ir,iv)) > 1.d-12 ) THEN
                   upf_out%kbeta(iv)=ir
                   exit
                ENDIF
             ENDDO
          ENDDO
       ENDIF

       ALLOCATE (upf_out%chi(upf_out%mesh,upf_out%nwfc))
       upf_out%chi = chi_
       DEALLOCATE (chi_)

       RETURN
     END SUBROUTINE convert_casino


     SUBROUTINE write_casino_tab(upf_in, fileout)

       USE pseudo_types, ONLY : pseudo_upf

       IMPLICIT NONE

       CHARACTER(LEN=*), INTENT(in)       :: fileout
       TYPE(pseudo_upf), INTENT(in)       :: upf_in
       INTEGER :: i, lp1, unout_

       INTEGER, EXTERNAL :: atomic_number

       OPEN ( NEWUNIT=unout_, FILE=TRIM(fileout), ACTION = 'WRITE') 

       WRITE(unout_,*) "Converted Pseudopotential in REAL space for ", upf_in%psd
       WRITE(unout_,*) "Atomic number and pseudo-charge"
       WRITE(unout_,"(I3,F8.2)") atomic_number( upf_in%psd ),upf_in%zp
       WRITE(unout_,*) "Energy units (rydberg/hartree/ev):"
       WRITE(unout_,*) "rydberg"
       WRITE(unout_,*) "Angular momentum of local component (0=s,1=p,2=d..)"
       WRITE(unout_,"(I2)") upf_in%lloc
       WRITE(unout_,*) "NLRULE override (1) VMC/DMC (2) config gen (0 ==> &
            &input/default VALUE)"
       WRITE(unout_,*) "0 0"
       WRITE(unout_,*) "Number of grid points"
       WRITE(unout_,*) upf_in%mesh
       WRITE(unout_,*) "R(i) in atomic units"
       WRITE(unout_, "(T4,E22.15)") upf_in%r(:)

       lp1 = size ( vnl, 2 )
       DO i=1,lp1
          WRITE(unout_, "(A,I1,A)") 'r*potential (L=',i-1,') in Ry'
          WRITE(unout_, "(T4,E22.15)") vnl(:,i)
       ENDDO
       CLOSE (unout_) 
       DEALLOCATE(vnl)

     END SUBROUTINE write_casino_tab

     SUBROUTINE conv_upf2casino(upf_in)

       USE pseudo_types, ONLY : pseudo_upf


       IMPLICIT NONE

       TYPE(pseudo_upf), INTENT(in)       :: upf_in
       INTEGER :: i, l, channels

       REAL(dp), PARAMETER :: offset=1E-20_dp
       !This is an offset added to the wavefunctions to
       !eliminate any divide by zeros that may be caused by
       !zeroed wavefunction terms.

       IF (upf_in%typ /= 'NC') THEN
          WRITE(0,*) ''
          WRITE(0,*) 'WRONG PSEUDOPOTENTIAL!'
          WRITE(0,*) 'Only norm-conserving pps can be used in CASINO!'
          STOP
       ENDIF
       
       WRITE(0,*) "Number of grid points: ", upf_in%mesh
       WRITE(0,*) "Number of KB projectors: ", upf_in%nbeta
       WRITE(0,*) "Channel(s) of KB projectors: ", upf_in%lll
       WRITE(0,*) "Number of channels to be re-constructed: ", upf_in%nbeta+1

       channels=upf_in%nbeta+1
       ALLOCATE ( vnl(upf_in%mesh,channels) )

       !Set up the local component of each channel
       DO i=1,channels
          vnl(:,i)=upf_in%r(:)*upf_in%vloc(:)
       ENDDO


       DO i=1,upf_in%nbeta
          l=upf_in%lll(i)+1

          !Check if any wfc components have been zeroed
          !and apply the offset IF they have

          IF ( minval(abs(upf_in%chi(:,l))) /= 0 ) THEN
             vnl(:,l)= (upf_in%beta(:,l)/(upf_in%chi(:,l)) &
                  *upf_in%r(:)) + vnl(:,l)
          ELSE
             WRITE(0,"(A,ES10.3,A)") 'Applying ',offset , ' offset to &
                  &wavefunction to avoid divide by zero'
             vnl(:,l)= (upf_in%beta(:,l)/(upf_in%chi(:,l)+offset) &
                  *upf_in%r(:)) + vnl(:,l)
          ENDIF

       ENDDO

     END SUBROUTINE conv_upf2casino

END MODULE casino_pp