Cleanup of various conversion routines between the different celldm, abc, etc.

Incorrect routines removed, some corner cases fixed
This commit is contained in:
Paolo Giannozzi 2018-05-28 18:31:53 +02:00
parent 71a61c005d
commit dc9c0e4289
2 changed files with 70 additions and 113 deletions

View File

@ -40,8 +40,8 @@ SUBROUTINE run_dist ( exit_status )
! and reprinted along with the lattice vectors, irrespective of
! what was provided in output - useful for checking and conversion
!
CALL cell2abc ( alat, at, a,b,c,cosab,cosac,cosbc )
CALL abc2celldm ( ibrav, a,b,c,cosab,cosac,cosbc, celldm )
CALL lat2celldm ( ibrav, alat, at(1,1), at(1,2), at(1,3), celldm )
CALL celldm2abc ( ibrav, celldm, a,b,c,cosab,cosac,cosbc )
!
WRITE(stdout,'(/,5x,"Bravais lattice index ibrav =",i3)') ibrav
WRITE(stdout,'(5X,"celldm(1) = ",f12.8," (au)",t40,"a = ",f12.8," (A)")') &

View File

@ -393,47 +393,39 @@ SUBROUTINE lat2celldm (ibrav,alat,a1,a2,a3,celldm)
SELECT CASE ( ibrav )
CASE (1:3,-3)
celldm(1) = alat
celldm(2:6) = 0.d0
CASE (4)
celldm(1) = alat
celldm(2) = 0.d0
celldm(3) = SQRT( DOT_PRODUCT(a3(:),a3(:)))/alat
celldm(4:6) = 0.d0
celldm(3) = ABS(a3(3)/a1(1))
CASE (5, -5 )
celldm(1)= alat
celldm(2:3) = 0.d0
celldm(4) = DOT_PRODUCT(a1(:),a2(:))/(alat**2)
celldm(5:6) = 0.d0
celldm(1) = alat
celldm(4) = DOT_PRODUCT(a1(:),a2(:)) / SQRT(a1(1)**2+a1(2)**2+a1(3)**2) &
/ SQRT(a2(1)**2+a2(2)**2+a2(3)**2)
CASE (6)
celldm(1)= alat
celldm(3)= SQRT( DOT_PRODUCT(a3(:),a3(:)))/alat
celldm(2)= 1.d0
celldm(4:6) = 0.d0
celldm(3)= ABS(a3(3)/a1(1))
CASE (7)
celldm(1) = alat
celldm(3) = ABS(a3(3)/a3(1))
celldm(2)= 0.d0
celldm(4:6) = 0.d0
CASE (8)
celldm(1) = alat
celldm(2) = SQRT( DOT_PRODUCT (a2(:),a2(:)))/alat
celldm(3) = SQRT( DOT_PRODUCT (a3(:),a3(:)))/alat
celldm(4:6) = 0.d0
celldm(2) = ABS(a2(2)/a1(1))
celldm(3) = ABS(a3(3)/a1(1))
CASE (9, -9 )
celldm(1) = alat
celldm(2) = ABS ( a1(2)/a1(1))
celldm(3) = ABS ( a3(3)/2.d0/a1(1))
celldm(4:6) = 0.d0
CASE (91 )
celldm(1) = alat
celldm(2) = ABS ( a2(2)/a1(1))*2.d0
celldm(3) = ABS ( a3(3)/a1(1))*2.d0
CASE (10)
celldm(1) = alat
celldm(2) = ABS ( a2(2)/a1(2))
celldm(2) = ABS ( a2(2)/a2(1))
celldm(3) = ABS ( a1(3)/a1(1))
celldm(4:6) = 0.d0
CASE (11)
celldm(1) = alat
celldm(2) = ABS(a1(2)/a1(1))
celldm(3) = ABS(a1(3)/a1(1))
celldm(4:6) = 0.d0
CASE (12, -12)
celldm(1) = alat
celldm(2) = SQRT( DOT_PRODUCT(a2(:),a2(:))/DOT_PRODUCT(a1(:),a1(:)))
@ -442,13 +434,16 @@ SUBROUTINE lat2celldm (ibrav,alat,a1,a2,a3,celldm)
SQRT(DOT_PRODUCT(a1(:),a1(:))*DOT_PRODUCT(a2(:),a2(:)))
celldm(5) = DOT_PRODUCT(a1(:),a3(:))/&
SQRT(DOT_PRODUCT(a1(:),a1(:))*DOT_PRODUCT(a3(:),a3(:)))
celldm(6) = 0.d0
CASE (13)
celldm(1) = alat
celldm(2) = SQRT( DOT_PRODUCT(a2(:),a2(:)))/(2.d0*a1(1))
celldm(3) = ABS (a3(3)/a3(1))
celldm(4) = COS( ATAN2( a2(2), a2(1) ) )
celldm(5:6) = 0.d0
CASE (-13)
celldm(1) = alat
celldm(2) = ABS (a2(2)/a2(1))
celldm(3) = SQRT( DOT_PRODUCT(a3(:),a3(:)))/(2.d0*a1(1))
celldm(5) = COS( ATAN2( a3(3), a3(1) ) )
CASE (14)
celldm(1) = alat
celldm(2) = SQRT( DOT_PRODUCT(a2(:),a2(:))/DOT_PRODUCT(a1(:),a1(:)))
@ -462,86 +457,16 @@ SUBROUTINE lat2celldm (ibrav,alat,a1,a2,a3,celldm)
CASE default
celldm(1) = 1.d0
IF (alat > 0.d0 ) celldm(1) = alat
celldm (2:6) = 0.d0
END SELECT
END SUBROUTINE lat2celldm
!
INTEGER FUNCTION abc2ibrav ( a,b,c,cosab,cosac,cosbc ) RESULT (ibrav)
!
! guess of bravais lattice from crystallographic parameters
! returns valid value for ibrav, or else ibrav=-1 if failed
!
USE kinds, ONLY: dp
IMPLICIT NONE
!
REAL(DP), INTENT (IN) :: a,b,c, cosab, cosac, cosbc
!
!
IF ( (b == 0.0_dp .AND. c == 0.0_dp) .OR. (b == a .AND. c == a) ) THEN
! Case: a=b=c
IF ( cosab == 0.0_dp .AND. cosac == 0.0_dp .AND. cosbc == 0.0_dp ) THEN
! Case: simple cubic
ibrav = 1
ELSE IF ( cosab/= 0.0_dp .AND. cosac== 0.0_dp .AND. cosbc== 0.0_dp) THEN
! Case: trigonal R with threefold axis along <111>
ibrav =-5
ELSE IF ( cosab==0.5_dp .AND. cosac == cosab .AND. cosbc == cosab ) THEN
! Case: fcc
ibrav = 2
ELSE IF ( ABS ( cosab + 1.0_dp/sqrt(3.0_dp) ) < 1.0d-6 .AND. &
cosac == cosab .AND. cosbc == cosab ) THEN
! Case: bcc with symmetric basis
ibrav =-3
ELSE
! Case: unknown (possibly wrong)
ibrav =-1
END IF
ELSE IF ( (b == 0.0_dp .AND. c > 0.0_dp) .OR. (b == a .AND. c /= a) ) THEN
! Case: a=b/=c
IF ( cosab == 0.0_dp .AND. cosac == 0.0_dp .AND. cosbc == 0.0_dp ) THEN
! Case: simple tetragonal
ibrav = 6
ELSE IF ( cosab ==-0.5_dp .AND. cosac == 0.0_dp .AND. cosbc == 0.0_dp ) THEN
! Case: simple hexagonal
ibrav = 4
ELSE
! Case: unknown (body-centered tetragonal or wrong data)
ibrav =-1
END IF
ELSE IF ( b > 0.0_dp .AND. c > 0.0_dp .AND. b /= a .AND. c /= a ) THEN
! Case: a/=b/=c
IF ( cosab == 0.0_dp .AND. cosac == 0.0_dp .AND. cosbc == 0.0_dp ) THEN
! Case: simple orthorhombic
ibrav = 8
ELSE IF ( cosab /= 0.0_dp .AND. cosac == 0.0_dp .AND. cosbc == 0.0_dp ) THEN
! Case: monoclinic P, unique axis c
ibrav = 12
ELSE IF ( cosab == 0.0_dp .AND. cosac /= 0.0_dp .AND. cosbc == 0.0_dp ) THEN
! Case: monoclinic P, unique axis b
ibrav =-12
ELSE IF ( cosab /= 0.0_dp .AND. cosac /= 0.0_dp .AND. cosbc /= 0.0_dp ) THEN
! Case: triclinic
ibrav = 14
ELSE
! Case: unknown (base-, face-, body-centered orthorombic,
! (base-centered monoclinic)
ibrav = -1
END IF
END IF
!
END FUNCTION abc2ibrav
!
SUBROUTINE abc2celldm ( ibrav, a,b,c,cosab,cosac,cosbc, celldm )
!
! returns internal parameters celldm from crystallographics ones
!
USE kinds, ONLY: dp
USE kinds, ONLY: dp
USE constants, ONLY: bohr_radius_angs
IMPLICIT NONE
!
@ -563,7 +488,7 @@ SUBROUTINE abc2celldm ( ibrav, a,b,c,cosab,cosac,cosbc, celldm )
celldm(2) = b / a
celldm(3) = c / a
!
IF ( ibrav == 14 ) THEN
IF ( ibrav == 14 .OR. ibrav == 0 ) THEN
!
! ... triclinic lattice
!
@ -575,41 +500,73 @@ SUBROUTINE abc2celldm ( ibrav, a,b,c,cosab,cosac,cosbc, celldm )
!
! ... monoclinic P or base centered lattice, unique axis b
!
celldm(4) = 0.0_dp
celldm(5) = cosac
celldm(6) = 0.0_dp
!
ELSE
ELSE IF ( 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 cell2abc ( alat, at, a,b,c,cosab,cosac,cosbc )
SUBROUTINE celldm2abc ( ibrav, celldm, a,b,c,cosab,cosac,cosbc )
!
! returns crystallographic parameters a,b,c from lattice vectors
! returns crystallographic parameters a,b,c from celldm
!
USE kinds, ONLY: dp
USE kinds, ONLY: dp
USE constants, ONLY: bohr_radius_angs
IMPLICIT NONE
!
REAL(DP), INTENT (IN) :: alat, at(3,3)
INTEGER, INTENT (IN) :: ibrav
REAL(DP), INTENT (IN) :: celldm(6)
REAL(DP), INTENT (OUT) :: a,b,c, cosab, cosac, cosbc
REAL(DP) :: norm1, norm2, norm3
!
!
norm1 = SQRT ( at(1,1)**2 + at(2,1)**2 + at(3,1)**2 )
norm2 = SQRT ( at(1,2)**2 + at(2,2)**2 + at(3,2)**2 )
norm3 = SQRT ( at(1,3)**2 + at(2,3)**2 + at(3,3)**2 )
a = celldm(1) * bohr_radius_angs
b = celldm(1)*celldm(2) * bohr_radius_angs
c = celldm(1)*celldm(3) * bohr_radius_angs
!
a = alat * norm1 / bohr_radius_angs
b = alat * norm2 / bohr_radius_angs
c = alat * norm3 / bohr_radius_angs
IF ( ibrav == 14 .OR. ibrav == 0 ) THEN
!
! ... triclinic lattice
!
cosbc = celldm(4)
cosac = celldm(5)
cosab = celldm(6)
!
ELSE IF ( 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
!
ELSE IF ( 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
!
cosab = (at(1,1)*at(1,2) + at(2,1)*at(2,2) + at(3,1)*at(3,2))/norm1/norm2
cosac = (at(1,1)*at(1,3) + at(2,1)*at(2,3) + at(3,1)*at(3,3))/norm1/norm3
cosbc = (at(1,3)*at(1,2) + at(2,3)*at(2,2) + at(3,3)*at(3,2))/norm3/norm2
!
END SUBROUTINE cell2abc
END SUBROUTINE celldm2abc