!
! 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 .
!
!
!-------------------------------------------------------------------------
SUBROUTINE latgen_lib(ibrav,celldm,a1,a2,a3,omega, ierr, errormsg)
  !-----------------------------------------------------------------------
  !     sets up the crystallographic vectors a1, a2, and a3.
  !
  !     ibrav is the structure index:
  !       1  cubic P (sc)                8  orthorhombic P
  !       2  cubic F (fcc)               9  1-face (C) centered orthorhombic
  !       3  cubic I (bcc)              10  all face centered orthorhombic
  !       4  hexagonal and trigonal P   11  body centered orthorhombic
  !       5  trigonal R, 3-fold axis c  12  monoclinic P (unique axis: c)
  !       6  tetragonal P (st)          13  one face (base) centered monoclinic
  !       7  tetragonal I (bct)         14  triclinic P
  !     Also accepted:
  !       0  "free" structure          -12  monoclinic P (unique axis: b)
  !      -3  cubic bcc with a more symmetric choice of axis
  !      -5  trigonal R, threefold axis along (111)
  !      -9  alternate description for base centered orthorhombic
  !     -13  one face (base) centered monoclinic (unique axis: b)
  !      91  1-face (A) centered orthorombic
  !
  !     celldm are parameters which fix the shape of the unit cell
  !     omega is the unit-cell volume
  !
  !     NOTA BENE: all axis sets are right-handed
  !     Boxes for US PPs do not work properly with left-handed axis
  !
  USE kinds, ONLY: DP
  IMPLICIT NONE
  INTEGER, INTENT(in) :: ibrav
  real(DP), INTENT(inout) :: celldm(6)
  real(DP), INTENT(inout) :: a1(3), a2(3), a3(3)
  real(DP), INTENT(out) :: omega
  !
  character(len=*),INTENT(out) :: errormsg
  integer,INTENT(out) :: ierr
  !
  real(DP), PARAMETER:: sr2 = 1.414213562373d0, &
                        sr3 = 1.732050807569d0
  INTEGER :: i,j,k,l,iperm,ir
  real(DP) :: term, cbya, s, term1, term2, singam, sen
  !
  ! pre-set to zero, in case we quit because of error
  Omega = 0._dp
  ierr = 0
  errormsg = ''
  !
  !  user-supplied lattice vectors
  !
  IF (ibrav == 0) THEN
     IF (sqrt( a1(1)**2 + a1(2)**2 + a1(3)**2 ) == 0 )  &
         THEN; errormsg='wrong at for ibrav=0'; ierr= 1; RETURN; ENDIF
     IF (sqrt( a2(1)**2 + a2(2)**2 + a2(3)**2 ) == 0 )  &
         THEN; errormsg='wrong at for ibrav=0'; ierr= 2; RETURN; ENDIF
     IF (sqrt( a3(1)**2 + a3(2)**2 + a3(3)**2 ) == 0 )  &
         THEN; errormsg='wrong at for ibrav=0'; ierr= 3; RETURN; ENDIF

     IF ( celldm(1) /= 0.D0 ) THEN
     !
     ! ... input at are in units of alat => convert them to a.u.
     !
         a1(:) = a1(:) * celldm(1)
         a2(:) = a2(:) * celldm(1)
         a3(:) = a3(:) * celldm(1)
     ELSE
     !
     ! ... input at are in atomic units: define celldm(1) from a1
     !
         celldm(1) = sqrt( a1(1)**2 + a1(2)**2 + a1(3)**2 )
     ENDIF
     !
  ELSE
     a1(:) = 0.d0
     a2(:) = 0.d0
     a3(:) = 0.d0
  ENDIF
  !
  IF (celldm (1) <= 0.d0) THEN; errormsg='wrong celldm(1)'; ierr= abs(ibrav); RETURN; ENDIF
  !
  !  index of bravais lattice supplied
  !
  IF (ibrav == 1) THEN
     !
     !     simple cubic lattice
     !
     a1(1)=celldm(1)
     a2(2)=celldm(1)
     a3(3)=celldm(1)
     !
  ELSEIF (ibrav == 2) THEN
     !
     !     fcc lattice
     !
     term=celldm(1)/2.d0
     a1(1)=-term
     a1(3)=term
     a2(2)=term
     a2(3)=term
     a3(1)=-term
     a3(2)=term
     !
  ELSEIF (abs(ibrav) == 3) THEN
     !
     !     bcc lattice
     !
     term=celldm(1)/2.d0
     DO ir=1,3
        a1(ir)=term
        a2(ir)=term
        a3(ir)=term
     ENDDO
     IF ( ibrav < 0 ) THEN
        a1(1)=-a1(1)
        a2(2)=-a2(2)
        a3(3)=-a3(3)
     ELSE
        a2(1)=-a2(1)
        a3(1)=-a3(1)
        a3(2)=-a3(2)
     ENDIF
     !
  ELSEIF (ibrav == 4) THEN
     !
     !     hexagonal lattice
     !
     IF (celldm (3) <= 0.d0) THEN; errormsg='wrong celldm(3)'; ierr= ibrav; RETURN; ENDIF
     !
     cbya=celldm(3)
     a1(1)=celldm(1)
     a2(1)=-celldm(1)/2.d0
     a2(2)=celldm(1)*sr3/2.d0
     a3(3)=celldm(1)*cbya
     !
  ELSEIF (abs(ibrav) == 5) THEN
     !
     !     trigonal lattice
     !
     IF (celldm (4) <= -0.5_dp .or. celldm (4) >= 1.0_dp) &
             THEN; errormsg='wrong celldm(4)'; ierr=5; RETURN; ENDIF
     !
     term1=sqrt(1.0_dp + 2.0_dp*celldm(4))
     term2=sqrt(1.0_dp - celldm(4))
     !
     IF ( ibrav == 5) THEN
        !     threefold axis along c (001)
        a2(2)=sr2*celldm(1)*term2/sr3
        a2(3)=celldm(1)*term1/sr3
        a1(1)=celldm(1)*term2/sr2
        a1(2)=-a1(1)/sr3
        a1(3)= a2(3)
        a3(1)=-a1(1)
        a3(2)= a1(2)
        a3(3)= a2(3)
     ELSEIF ( ibrav == -5) THEN
        !     threefold axis along (111)
        ! Notice that in the cubic limit (alpha=90, celldm(4)=0, term1=term2=1)
        ! does not yield the x,y,z axis, but an equivalent rotated triplet:
        !    a/3 (-1,2,2), a/3 (2,-1,2), a/3 (2,2,-1)
        ! If you prefer the x,y,z axis as cubic limit, you should modify the
        ! definitions of a1(1) and a1(2) as follows:'
        !    a1(1) = celldm(1)*(term1+2.0_dp*term2)/3.0_dp
        !    a1(2) = celldm(1)*(term1-term2)/3.0_dp
        ! (info by G. Pizzi and A. Cepellotti)
        !
        a1(1) = celldm(1)*(term1-2.0_dp*term2)/3.0_dp
        a1(2) = celldm(1)*(term1+term2)/3.0_dp
        a1(3) = a1(2)
        a2(1) = a1(3)
        a2(2) = a1(1)
        a2(3) = a1(2)
        a3(1) = a1(2)
        a3(2) = a1(3)
        a3(3) = a1(1)
     ENDIF
  ELSEIF (ibrav == 6) THEN
     !
     !     tetragonal lattice
     !
     IF (celldm (3) <= 0.d0) THEN; errormsg='wrong celldm(3)'; ierr=6; RETURN; ENDIF
     !
     cbya=celldm(3)
     a1(1)=celldm(1)
     a2(2)=celldm(1)
     a3(3)=celldm(1)*cbya
     !
  ELSEIF (ibrav == 7) THEN
     !
     !     body centered tetragonal lattice
     !
     IF (celldm (3) <= 0.d0) THEN; errormsg='wrong celldm(3)'; ierr=7; RETURN; ENDIF
     !
     cbya=celldm(3)
     a2(1)=celldm(1)/2.d0
     a2(2)=a2(1)
     a2(3)=cbya*celldm(1)/2.d0
     a1(1)= a2(1)
     a1(2)=-a2(1)
     a1(3)= a2(3)
     a3(1)=-a2(1)
     a3(2)=-a2(1)
     a3(3)= a2(3)
     !
  ELSEIF (ibrav == 8) THEN
     !
     !     Simple orthorhombic lattice
     !
     IF (celldm (2) <= 0.d0) THEN; errormsg='wrong celldm(2)'; ierr=8; RETURN; ENDIF
     IF (celldm (3) <= 0.d0) THEN; errormsg='wrong celldm(3)'; ierr=8; RETURN; ENDIF
     !
     a1(1)=celldm(1)
     a2(2)=celldm(1)*celldm(2)
     a3(3)=celldm(1)*celldm(3)
     !
  ELSEIF ( abs(ibrav) == 9) THEN
     !
     !     One face (base) centered orthorhombic lattice  (C type)
     !
     IF (celldm (2) <= 0.d0) THEN; errormsg='wrong celldm(2)'; ierr=9; RETURN; ENDIF
     IF (celldm (3) <= 0.d0) THEN; errormsg='wrong celldm(3)'; ierr=9; RETURN; ENDIF
     !
     IF ( ibrav == 9 ) THEN
        !   old PWscf description
        a1(1) = 0.5d0 * celldm(1)
        a1(2) = a1(1) * celldm(2)
        a2(1) = - a1(1)
        a2(2) = a1(2)
     ELSE
        !   alternate description
        a1(1) = 0.5d0 * celldm(1)
        a1(2) =-a1(1) * celldm(2)
        a2(1) = a1(1)
        a2(2) =-a1(2)
     ENDIF
     a3(3) = celldm(1) * celldm(3)
     !
  ELSEIF ( ibrav == 91 ) THEN
     !
     !     One face (base) centered orthorhombic lattice  (A type)
     !
     IF (celldm (2) <= 0.d0) THEN; errormsg='wrong celldm(2)'; ierr=91; RETURN; ENDIF
     IF (celldm (3) <= 0.d0) THEN; errormsg='wrong celldm(3)'; ierr=91; RETURN; ENDIF
     !
     a1(1) = celldm(1)
     a2(2) = celldm(1) * celldm(2) * 0.5_DP
     a2(3) = - celldm(1) * celldm(3) * 0.5_DP
     a3(2) = a2(2)
     a3(3) = - a2(3)
     !
  ELSEIF (ibrav == 10) THEN
     !
     !     All face centered orthorhombic lattice
     !
     IF (celldm (2) <= 0.d0) THEN; errormsg='wrong celldm(2)'; ierr=10; RETURN; ENDIF
     IF (celldm (3) <= 0.d0) THEN; errormsg='wrong celldm(3)'; ierr=10; RETURN; ENDIF
     !
     a2(1) = 0.5d0 * celldm(1)
     a2(2) = a2(1) * celldm(2)
     a1(1) = a2(1)
     a1(3) = a2(1) * celldm(3)
     a3(2) = a2(1) * celldm(2)
     a3(3) = a1(3)
     !
  ELSEIF (ibrav == 11) THEN
     !
     !     Body centered orthorhombic lattice
     !
     IF (celldm (2) <= 0.d0) THEN; errormsg='wrong celldm(2)'; ierr=11; RETURN; ENDIF
     IF (celldm (3) <= 0.d0) THEN; errormsg='wrong celldm(3)'; ierr=11; RETURN; ENDIF
     !
     a1(1) = 0.5d0 * celldm(1)
     a1(2) = a1(1) * celldm(2)
     a1(3) = a1(1) * celldm(3)
     a2(1) = - a1(1)
     a2(2) = a1(2)
     a2(3) = a1(3)
     a3(1) = - a1(1)
     a3(2) = - a1(2)
     a3(3) = a1(3)
     !
  ELSEIF (ibrav == 12) THEN
     !
     !     Simple monoclinic lattice, unique (i.e. orthogonal to a) axis: c
     !
     IF (celldm (2) <= 0.d0)   THEN; errormsg='wrong celldm(2)'; ierr=12; RETURN; ENDIF
     IF (celldm (3) <= 0.d0)   THEN; errormsg='wrong celldm(3)'; ierr=12; RETURN; ENDIF
     IF (abs(celldm(4))>=1.d0) THEN; errormsg='wrong celldm(4)'; ierr=12; RETURN; ENDIF
     !
     sen=sqrt(1.d0-celldm(4)**2)
     a1(1)=celldm(1)
     a2(1)=celldm(1)*celldm(2)*celldm(4)
     a2(2)=celldm(1)*celldm(2)*sen
     a3(3)=celldm(1)*celldm(3)
     !
  ELSEIF (ibrav ==-12) THEN
     !
     !     Simple monoclinic lattice, unique axis: b (more common)
     !
     IF (celldm (2) <= 0.d0)   THEN; errormsg='wrong celldm(2)'; ierr=12; RETURN; ENDIF
     IF (celldm (3) <= 0.d0)   THEN; errormsg='wrong celldm(3)'; ierr=12; RETURN; ENDIF
     IF (abs(celldm(5))>=1.d0) THEN; errormsg='wrong celldm(5)'; ierr=12; RETURN; ENDIF
     !
     sen=sqrt(1.d0-celldm(5)**2)
     a1(1)=celldm(1)
     a2(2)=celldm(1)*celldm(2)
     a3(1)=celldm(1)*celldm(3)*celldm(5)
     a3(3)=celldm(1)*celldm(3)*sen
     !
  ELSEIF (ibrav == 13) THEN
     !
     !     One face centered monoclinic lattice unique axis c
     !
     IF (celldm (2) <= 0.d0)   THEN; errormsg='wrong celldm(2)'; ierr=13; RETURN; ENDIF
     IF (celldm (3) <= 0.d0)   THEN; errormsg='wrong celldm(3)'; ierr=13; RETURN; ENDIF
     IF (abs(celldm(4))>=1.d0) THEN; errormsg='wrong celldm(4)'; ierr=13; RETURN; ENDIF
     !
     sen = sqrt( 1.d0 - celldm(4) ** 2 )
     a1(1) = 0.5d0 * celldm(1)
     a1(3) =-a1(1) * celldm(3)
     a2(1) = celldm(1) * celldm(2) * celldm(4)
     a2(2) = celldm(1) * celldm(2) * sen
     a3(1) = a1(1)
     a3(3) =-a1(3)
  ELSEIF (ibrav == -13) THEN
     errormsg='BEWARE: axis for ibrav=-13 changed, see documentation!'
     !
     !     One face centered monoclinic lattice unique axis b
     !
     IF (celldm (2) <= 0.d0)   THEN; errormsg='wrong celldm(2)'; ierr=13; RETURN; ENDIF
     IF (celldm (3) <= 0.d0)   THEN; errormsg='wrong celldm(3)'; ierr=13; RETURN; ENDIF
     IF (abs(celldm(5))>=1.d0) THEN; errormsg='wrong celldm(5)'; ierr=13; RETURN; ENDIF
     !
     sen = sqrt( 1.d0 - celldm(5) ** 2 )
     a1(1) = 0.5d0 * celldm(1)
     a1(2) = a1(1) * celldm(2)
     a2(1) =-a1(1)
     a2(2) = a1(2)
     a3(1) = celldm(1) * celldm(3) * celldm(5)
     a3(3) = celldm(1) * celldm(3) * sen
     !
  ELSEIF (ibrav == 14) THEN
     !
     !     Triclinic lattice
     !
     IF (celldm (2) <= 0.d0)   THEN; errormsg='wrong celldm(2)'; ierr=14; RETURN; ENDIF
     IF (celldm (3) <= 0.d0)   THEN; errormsg='wrong celldm(3)'; ierr=14; RETURN; ENDIF
     IF (abs(celldm(4))>=1.d0) THEN; errormsg='wrong celldm(4)'; ierr=14; RETURN; ENDIF
     IF (abs(celldm(5))>=1.d0) THEN; errormsg='wrong celldm(5)'; ierr=14; RETURN; ENDIF
     IF (abs(celldm(6))>=1.d0) THEN; errormsg='wrong celldm(6)'; ierr=14; RETURN; ENDIF
     !
     singam=sqrt(1.d0-celldm(6)**2)
     term= (1.d0+2.d0*celldm(4)*celldm(5)*celldm(6)             &
          -celldm(4)**2-celldm(5)**2-celldm(6)**2)
     IF (term < 0.d0) THEN; errormsg='celldm do not make sense, check your data'; ierr=14; RETURN; ENDIF
     term= sqrt(term/(1.d0-celldm(6)**2))
     a1(1)=celldm(1)
     a2(1)=celldm(1)*celldm(2)*celldm(6)
     a2(2)=celldm(1)*celldm(2)*singam
     a3(1)=celldm(1)*celldm(3)*celldm(5)
     a3(2)=celldm(1)*celldm(3)*(celldm(4)-celldm(5)*celldm(6))/singam
     a3(3)=celldm(1)*celldm(3)*term
     !
  ELSE IF (ibrav /= 0) THEN
     !
     errormsg='nonexistent bravais lattice'
     ierr=ABS(ibrav)
     RETURN
     !
  ENDIF
  !
  !  calculate unit-cell volume omega
  !
  CALL volume (1.0_dp, a1, a2, a3, omega)
  !
  RETURN
  !
END SUBROUTINE latgen_lib
!
!-------------------------------------------------------------------------
SUBROUTINE at2celldm (ibrav,alat,a1,a2,a3,celldm)
  !-----------------------------------------------------------------------
  !
  !     Returns celldm parameters computed from lattice vectors a1,a2,a3 
  !     a1, a2, a3 are in "alat" units
  !     If Bravais lattice index ibrav=0, only celldm(1) is set to alat
  !     See latgen for definition of celldm and lattice vectors.
  !     a1, a2, a3, ibrav, alat are not modified
  !
  USE kinds, ONLY: DP
  IMPLICIT NONE
  INTEGER, INTENT(in) :: ibrav
  REAL(DP), INTENT(in) :: alat, a1(3), a2(3), a3(3)
  REAL(DP), INTENT(out) :: celldm(6)
  !
  celldm = 0.d0
  !
  SELECT CASE  ( ibrav )
  CASE (0)
     celldm(1) = 1.0_dp
  CASE (1)
     celldm(1) = sqrt( dot_product (a1,a1) )
  CASE (2)
     celldm(1) = sqrt( dot_product (a1,a1) * 2.0_dp )
  CASE (3,-3)
     celldm(1) = sqrt( dot_product (a1,a1) / 3.0_dp ) * 2.0_dp
  CASE (4)
     celldm(1) = sqrt( dot_product (a1,a1) )
     celldm(3) = sqrt( dot_product (a3,a3) ) / celldm(1)
  CASE (5, -5 )
     celldm(1) = sqrt( dot_product (a1,a1) )
     celldm(4) = dot_product(a1(:),a2(:)) / celldm(1) &
                 / sqrt( dot_product( a2,a2) )
  CASE (6)
     celldm(1)= sqrt( dot_product (a1,a1) )
     celldm(3)= sqrt( dot_product (a3,a3) ) / celldm(1)
  CASE (7)
     celldm(1) = abs(a1(1))*2.0_dp
     celldm(3) = abs(a1(3)/a1(1))
  CASE (8)
     celldm(1) = sqrt( dot_product (a1,a1) )
     celldm(2) = sqrt( dot_product (a2,a2) ) / celldm(1)
     celldm(3) = sqrt( dot_product (a3,a3) ) / celldm(1)
  CASE (9, -9 )
     celldm(1) = abs(a1(1))*2.0_dp
     celldm(2) = abs(a2(2))*2.0_dp/celldm(1)
     celldm(3) = abs(a3(3))/celldm(1)
  CASE (91 )
     celldm(1) = sqrt( dot_product (a1,a1) )
     celldm(2) = abs (a2(2))*2.0_dp/celldm(1)
     celldm(3) = abs (a3(3))*2.0_dp/celldm(1)
  CASE (10)
     celldm(1) = abs(a1(1))*2.0_dp
     celldm(2) = abs(a2(2))*2.0_dp/celldm(1)
     celldm(3) = abs(a3(3))*2.0_dp/celldm(1)
  CASE (11)
     celldm(1) = abs(a1(1))*2.0_dp
     celldm(2) = abs(a1(2))*2.0_dp/celldm(1)
     celldm(3) = abs(a1(3))*2.0_dp/celldm(1)
  CASE (12,-12)
     celldm(1) = sqrt( dot_product (a1,a1) )
     celldm(2) = sqrt( dot_product(a2(:),a2(:)) ) / celldm(1)
     celldm(3) = sqrt( dot_product(a3(:),a3(:)) ) / celldm(1)
     IF ( ibrav == 12 ) THEN
        celldm(4) = dot_product(a1(:),a2(:)) / celldm(1) / &
             sqrt(dot_product(a2(:),a2(:)))
     ELSE
        celldm(5) = dot_product(a1(:),a3(:)) / celldm(1) / &
             sqrt(dot_product(a3(:),a3(:)))
     ENDIF
  CASE (13)
     celldm(1) = abs(a1(1))*2.0_dp
     celldm(2) = sqrt( dot_product(a2(:),a2(:))) / celldm(1)
     celldm(3) = abs (a1(3)/a1(1))
     celldm(4) = a2(1)/a1(1)/celldm(2)/2.0_dp
     !     SQRT(DOT_PRODUCT(a1(:),a1(:)) * DOT_PRODUCT(a2(:),a2(:)))
     !celldm(4) = DOT_PRODUCT(a1(:),a2(:)) / &
     !     SQRT(DOT_PRODUCT(a1(:),a1(:)) * DOT_PRODUCT(a2(:),a2(:)))
  CASE (-13)
     celldm(1) = abs(a1(1))*2.0_dp
     celldm(2) = abs (a2(2)/a2(1))
     celldm(3) = sqrt( dot_product(a3(:),a3(:))) / celldm(1)
     celldm(5) = a3(1)/a1(1)/celldm(3)/2.0_dp
     !celldm(5) = DOT_PRODUCT(a1(:),a3(:)) / &
     !     SQRT(DOT_PRODUCT(a1(:),a1(:)) * DOT_PRODUCT(a3(:),a3(:)))
  CASE (14)
     celldm(1) = sqrt(dot_product(a1(:),a1(:)))
     celldm(2) = sqrt( dot_product(a2(:),a2(:))) / celldm(1)
     celldm(3) = sqrt( dot_product(a3(:),a3(:))) / celldm(1)
     celldm(4) = dot_product(a3(:),a2(:))/sqrt(dot_product(a2(:),a2(:)) * &
          dot_product(a3(:),a3(:)))
     celldm(5) = dot_product(a3(:),a1(:)) / celldm(1) / &
          sqrt( dot_product(a3(:),a3(:)))
     celldm(6) = dot_product(a1(:),a2(:)) / celldm(1) / &
          sqrt(dot_product(a2(:),a2(:)))
  CASE DEFAULT
     CALL infomsg('at2celldm', 'wrong ibrav?')
  END SELECT
  !
  celldm(1) = celldm(1) * alat
  !
END SUBROUTINE at2celldm
!
FUNCTION at2ibrav (a1, a2, a3) RESULT (ibrav)
  !
  !     Returns ibrav from lattice vectors if recognized, 0 otherwise
  !
  USE kinds, ONLY: dp
  IMPLICIT NONE
  REAL(dp), INTENT (in) :: a1(3), a2(3), a3(3)
  REAL(dp) :: v1, v2, v3, cosab, cosac, cosbc
  !
  INTEGER :: ibrav
  !
  v1 = sqrt( dot_product( a1,a1 ) )
  v2 = sqrt( dot_product( a2,a2 ) )
  v3 = sqrt( dot_product( a3,a3 ) )
  cosbc = dot_product(a2,a3)/v2/v3
  cosac = dot_product(a1,a3)/v1/v3
  cosab = dot_product(a1,a2)/v1/v2
  !
  ! Assume triclinic if nothing suitable found
  !
  ibrav = 14
  IF ( eqq(v1,v2) .and. eqq(v1,v3) ) THEN
     ! Case: a=b=c
     IF (eqq(cosab,cosac) .and. eqq(cosab,cosbc)) THEN
        ! Case: alpha = beta = gamma
        IF ( eqq(cosab,0.0_dp) ) THEN
           ! Cubic P - ibrav=1
           ibrav = 1
        ELSEIF ( eqq(cosab,0.5_dp) ) THEN
           ! Cubic F - ibrav=2
           ibrav = 2
        ELSEIF ( eqq(cosab,-1.0_dp/3.0_dp) ) THEN
           ! Cubic I - ibrav=-3
           ibrav = -3
        ELSE
           IF ( eqq(abs(a1(3)),abs(a2(3))) .and. eqq(abs(a2(3)),abs(a3(3))) ) THEN
              ! Trigonal 001 axis
              ibrav =5
           ELSE
              ! Trigonal, 111 axis
              ibrav =-5
           ENDIF
           !
        ENDIF
        !
     ELSEIF ( eqq(cosab,cosac) .and. neqq(cosab,cosbc) ) THEN
        IF ( eqq(abs(a1(1)),abs(a1(2))) .and. &
             eqq(abs(a2(1)),abs(a2(2))) ) THEN
           ! Tetragonal I
           ibrav = 7
        ELSE
           ! Cubic I - ibrav=3
           ibrav = 3
        ENDIF
     ELSEIF ( eqq(cosab,-cosac) .and. eqq(cosab,cosbc) .and. &
          eqq(cosab,1.0_dp/3.0_dp) ) THEN
        ! Cubic I - ibrav=3
        ibrav = 3
     ELSEIF ( eqq(abs(a1(1)),abs(a2(1))) .and. &
               eqq(abs(a2(2)),abs(a2(2))) ) THEN
        ! Orthorhombic body-centered
        ibrav = 11
     ENDIF
     !
  ELSEIF ( eqq(v1,v2) .and. neqq(v1,v3) ) THEN
     ! Case: a=b/=c
     IF ( eqq(cosab,0.0_dp) .and. eqq(cosac,0.0_dp) .and. eqq(cosbc,0.0_dp) ) THEN
        ! Case: alpha = beta = gamma = 90
        ! Simple tetragonal
        ibrav = 6
     ELSEIF ( eqq(cosab,-0.5_dp) .and. eqq(cosac,0.0_dp) .and. eqq(cosbc,0.0_dp) ) THEN
        ! Case: alpha = 120, beta = gamma = 90 => simple hexagonal
        ! Simple hexagonal
        ibrav = 4
     ELSEIF ( eqq(cosac,0.0_dp) .and. eqq(cosbc,0.0_dp) ) THEN
        ! Orthorhombic bco
        IF ( eqq(a1(1),a2(1)) .and. eqq(a1(2),-a2(2))) THEN
           ibrav = -9
        ELSEIF ( eqq(a1(1),-a2(1)) .and. eqq(a1(2),a2(2))) THEN
           ibrav = 9
        ENDIF
     ELSEIF ( eqq(cosac,-cosbc) ) THEN
        ! bco (unique axis b)
        ibrav =-13
     ENDIF

  ELSEIF ( eqq(v1,v3) .and. neqq(v1,v2) ) THEN
     ! Case: a=c/=b
     ! Monoclinic bco (unique axis c)
     ibrav = 13
     !
  ELSEIF ( eqq(v2,v3) .and. neqq(v1,v2) ) THEN
     ! Case: a/=b=c
     ! Orthorhombic 1-face bco
     ibrav = 91
     !
  ELSEIF ( neqq(v1,v2) .and. neqq(v1,v3) .and. neqq(v2,v3) ) THEN
     ! Case: a/=b/=c
     IF ( eqq(cosab,0.0_dp) .and. eqq(cosac,0.0_dp) .and. eqq(cosbc,0.0_dp) ) THEN
        ! Case: alpha = beta = gamma = 90
        ! Orthorhombic P
        ibrav = 8
     ELSEIF ( neqq(cosab,0.0_dp) .and. eqq(cosac,0.0_dp) .and. eqq(cosbc,0.0_dp) ) THEN
        ! Case: alpha /= 90,  beta = gamma = 90
        ! Monoclinic P, unique axis c
        ibrav = 12
     ELSEIF ( eqq(cosab,0.0_dp) .and. neqq(cosac,0.0_dp) .and. eqq(cosbc,0.0_dp) ) THEN
        ! Case: beta /= 90, alpha = gamma = 90
        ! Monoclinic P, unique axis b
        ibrav =-12
     ELSEIF ( neqq(cosab,0.0_dp) .and. neqq(cosac,0.0_dp) .and. neqq(cosbc,0.0_dp) ) THEN
        ! Case: alpha /= 90, beta /= 90, gamma /= 90
        IF ( eqq(abs(a1(1)),abs(a2(1))) .and. eqq(abs(a1(3)),abs(a3(3))) .and. &
             eqq(abs(a2(2)),abs(a3(2))) ) THEN
        ! Orthorhombic F
           ibrav = 10
        ELSE
           ! Triclinic
           ibrav = 14
        ENDIF
     ENDIF
     !
  ENDIF
  !
CONTAINS
  !
  LOGICAL FUNCTION eqq (x,y)
    REAL(dp), INTENT (in) :: x,y
    eqq = abs(x-y) < 1.0d-5
  END FUNCTION eqq
  LOGICAL FUNCTION neqq (x,y)
    REAL(dp), INTENT (in) :: x,y
    neqq= abs(x-y) > 1.0d-5
  END FUNCTION neqq
  !
END FUNCTION at2ibrav
!
SUBROUTINE abc2celldm ( ibrav, a,b,c,cosab,cosac,cosbc, celldm )
  !
  !  returns internal parameters celldm from crystallographics ones
  !
  USE kinds,     ONLY: dp
  USE constants, ONLY: bohr_radius_angs
  IMPLICIT NONE
  !
  INTEGER,  INTENT (in) :: ibrav
  REAL(DP), INTENT (in) :: a,b,c, cosab, cosac, cosbc
  REAL(DP), INTENT (out) :: celldm(6)
  !
  IF (a <= 0.0_dp) CALL errore('abc2celldm','incorrect lattice parameter (a)',1)
  IF (b <  0.0_dp) CALL errore('abc2celldm','incorrect lattice parameter (b)',1)
  IF (c <  0.0_dp) CALL errore('abc2celldm','incorrect lattice parameter (c)',1)
  IF ( abs (cosab) > 1.0_dp) CALL errore('abc2celldm', &
                   'incorrect lattice parameter (cosab)',1)
  IF ( abs (cosac) > 1.0_dp) CALL errore('abc2celldm', &
                   'incorrect lattice parameter (cosac)',1)
  IF ( abs (cosbc) > 1.0_dp) CALL errore('abc2celldm', &
       'incorrect lattice parameter (cosbc)',1)
  !
  celldm(1) = a / bohr_radius_angs
  celldm(2) = b / a
  celldm(3) = c / a
  !
  IF ( ibrav == 14 .or. ibrav == 0 ) THEN
     !
     ! ... triclinic lattice
     !
     celldm(4) = cosbc
     celldm(5) = cosac
     celldm(6) = cosab
     !
  ELSEIF ( ibrav ==-12 .or. ibrav ==-13 ) THEN
     !
     ! ... monoclinic P or base centered lattice, unique axis b
     !
     celldm(4) = 0.0_dp
     celldm(5) = cosac
     celldm(6) = 0.0_dp
     !
  ELSEIF ( ibrav ==-5 .or. ibrav ==5 .or. ibrav ==12 .or. ibrav ==13 ) THEN
     !
     ! ... trigonal and monoclinic lattices, unique axis c
     !
     celldm(4) = cosab
     celldm(5) = 0.0_dp
     celldm(6) = 0.0_dp
     !
  ELSE
     !
     celldm(4) = 0.0_dp
     celldm(5) = 0.0_dp
     celldm(6) = 0.0_dp
     !
  ENDIF
  !
END SUBROUTINE abc2celldm
!
SUBROUTINE celldm2abc ( ibrav, celldm, a,b,c,cosab,cosac,cosbc )
  !
  !  returns crystallographic parameters a,b,c from celldm
  !
  USE kinds,     ONLY: dp
  USE constants, ONLY: bohr_radius_angs
  IMPLICIT NONE
  !
  INTEGER,  INTENT (in) :: ibrav
  REAL(DP), INTENT (in) :: celldm(6)
  REAL(DP), INTENT (out) :: a,b,c, cosab, cosac, cosbc
  !
  !
  a = celldm(1) * bohr_radius_angs
  b = celldm(1)*celldm(2) * bohr_radius_angs
  c = celldm(1)*celldm(3) * bohr_radius_angs
  !
  IF ( ibrav == 14 .or. ibrav == 0 ) THEN
     !
     ! ... triclinic lattice
     !
     cosbc = celldm(4)
     cosac = celldm(5)
     cosab = celldm(6)
     !
  ELSEIF ( ibrav ==-12 .or. ibrav ==-13 ) THEN
     !
     ! ... monoclinic P or base centered lattice, unique axis b
     !
     cosab = 0.0_dp
     cosac = celldm(5)
     cosbc = 0.0_dp
     !
  ELSEIF ( ibrav ==-5 .or. ibrav ==5 .or. ibrav ==12 .or. ibrav ==13 ) THEN
     !
     ! ... trigonal and monoclinic lattices, unique axis c
     !
     cosab = celldm(4)
     cosac = 0.0_dp
     cosbc = 0.0_dp
     !
  ELSE
     cosab = 0.0_dp
     cosac = 0.0_dp
     cosbc = 0.0_dp
  ENDIF
  !
END SUBROUTINE celldm2abc

SUBROUTINE remake_cell(ibrav, alat, a1,a2,a3)
  USE kinds, ONLY : DP
  USE io_global, ONLY : stdout
  IMPLICIT NONE
  INTEGER,INTENT(in) :: ibrav
  REAL(DP),INTENT(inout) :: alat, a1(3),a2(3),a3(3)
  REAL(DP) :: e1(3), e2(3), e3(3)
  REAL(DP) :: celldm_internal(6), lat_internal, omega
  ! Better not to do the following, or it may cause problems with ibrav=0 from input
!  ibrav = at2ibrav (a(:,1), a(:,2), a(:,3))
  ! Instead, let's print a warning and do nothing:
  IF(ibrav==0)THEN
    WRITE(stdout,'(a)') "WARNING! With ibrav=0, cell_dofree='ibrav' does not have any effect. "
    RETURN
  ENDIF
  !
  CALL  at2celldm (ibrav,alat,a1, a2, a3,celldm_internal)
  WRITE(stdout,'("ibrav = ",i6)') ibrav
  WRITE(stdout,'(" celldm(1) = ",f15.8)') celldm_internal(1)
  IF( celldm_internal(2) /= 0._dp) WRITE(stdout,'(" celldm(2) = ",f15.8)') celldm_internal(2)
  IF( celldm_internal(3) /= 0._dp) WRITE(stdout,'(" celldm(3) = ",f15.8)') celldm_internal(3)
  IF( celldm_internal(4) /= 0._dp) WRITE(stdout,'(" celldm(4) = ",f15.8)') celldm_internal(4)
  IF( celldm_internal(5) /= 0._dp) WRITE(stdout,'(" celldm(5) = ",f15.8)') celldm_internal(5)
  IF( celldm_internal(6) /= 0._dp) WRITE(stdout,'(" celldm(6) = ",f15.8)') celldm_internal(6)
  !
  e1=a1
  e2=a2
  e3=a3
  CALL latgen( ibrav, celldm_internal, a1,a2,a3, omega )
  !WRITE(*, '("New lattice vectors in bohr:")')
  !WRITE(*,'(3f15.8)') e(:,1)
  !WRITE(*,'(3f15.8)') e(:,2)
  !WRITE(*,'(3f15.8)') e(:,3)
  WRITE(stdout,'("Input lattice vectors:")')
  WRITE(stdout,'(3f15.8)') e1
  WRITE(stdout,'(3f15.8)') e2
  WRITE(stdout,'(3f15.8)') e3
  WRITE(stdout,'("New lattice vectors in INITIAL alat:")')
  WRITE(stdout,'(3f15.8)') a1/alat
  WRITE(stdout,'(3f15.8)') a2/alat
  WRITE(stdout,'(3f15.8)') a3/alat
  WRITE(stdout, '("New lattice vectors in NEW alat (for information only):")')
  WRITE(stdout,'(3f15.8)') a1/celldm_internal(1)
  WRITE(stdout,'(3f15.8)') a2/celldm_internal(1)
  WRITE(stdout,'(3f15.8)') a3/celldm_internal(1)
  a1=a1/alat
  a2=a2/alat
  a3=a3/alat
  WRITE(*,'("Discrepancy in bohr = ", 3f12.6)') DSQRT(SUM((a1-e1)**2)), DSQRT(SUM((a2-e2)**2)), DSQRT(SUM((a3-e3)**2))


END SUBROUTINE

!-------------------------------------------------------------------------
SUBROUTINE latgen(ibrav,celldm,a1,a2,a3,omega)
  !-----------------------------------------------------------------------
  USE kinds, ONLY: DP
  IMPLICIT NONE
  INTEGER, INTENT(in) :: ibrav
  real(DP), INTENT(inout) :: celldm(6)
  real(DP), INTENT(inout) :: a1(3), a2(3), a3(3)
  real(DP), INTENT(out) :: omega
  !
  character(len=54) :: errormsg
  integer :: ierr

  CALL latgen_lib(ibrav,celldm,a1,a2,a3,omega, ierr, errormsg)
  IF(ierr /= 0 ) THEN
            CALL errore("latgen", errormsg, abs(ierr)) 
  ELSE
    IF(errormsg/='') CALL infomsg('latgen',errormsg)
  ENDIF
  !
END SUBROUTINE 

