1!
2! Copyright (C) 1996-2016	The SIESTA group
3!  This file is distributed under the terms of the
4!  GNU General Public License: see COPYING in the top directory
5!  or http://www.gnu.org/copyleft/gpl.txt.
6! See Docs/Contributors.txt for a list of contributors.
7!
8      subroutine typecell(cell,ctype,lv)
9c *******************************************************************
10c Finds out if cell is cubic, fcc or bcc
11c
12c Written by P. Ordejon, March 1999.
13C (from outcell of E. Artacho)
14c ********* INPUT ***************************************************
15c double precision cell(3,3): Lattice (supercell) vectors
16c ********* OUTPUT **************************************************
17c character        ctype:   : Type of cell (sc, fcc, bcc or none)
18c double precision lv       : Lattice vector
19c *******************************************************************
20
21      use precision, only : dp
22
23      implicit none
24
25      real(dp)         :: cell(3,3), lv
26
27      character(len=4) :: ctype
28
29c Internal variables and arrays
30
31      integer          :: iv, ix
32      real(dp)         :: cellm(3), celang(3), tol, pi
33
34      parameter (tol=1.d-4)
35
36      data pi / 3.1415926d0 /
37
38c Cell-vector modules
39
40      do iv = 1, 3
41         cellm(iv) = 0.d0
42         do ix = 1, 3
43            cellm(iv) = cellm(iv) + cell(ix,iv)*cell(ix,iv)
44         enddo
45         cellm(iv) = sqrt(cellm(iv))
46      enddo
47
48c Cell-vector angles
49
50      celang(1) = 0.d0
51      do ix = 1, 3
52         celang(1) = celang(1) + cell(ix,1)*cell(ix,2)
53      enddo
54      celang(1) = acos(celang(1)/(cellm(1)*cellm(2)))*180.d0/pi
55      celang(2) = 0.d0
56      do ix = 1, 3
57         celang(2) = celang(2) + cell(ix,1)*cell(ix,3)
58      enddo
59      celang(2) = acos(celang(2)/(cellm(1)*cellm(3)))*180.d0/pi
60      celang(3) = 0.d0
61      do ix = 1, 3
62         celang(3) = celang(3) + cell(ix,2)*cell(ix,3)
63      enddo
64      celang(3) = acos(celang(3)/(cellm(2)*cellm(3)))*180.d0/pi
65
66      ctype = 'none'
67
68C Check if lattice vectors have same length, and return if not
69
70      if (abs(cellm(1)-cellm(2)) .gt. tol) return
71      if (abs(cellm(1)-cellm(3)) .gt. tol) return
72
73C Check if angles are 90 deg, in which case the cell is cubic
74      if ((abs(celang(1) - 90.0) .lt. tol) .and.
75     .    (abs(celang(2) - 90.0) .lt. tol) .and.
76     .    (abs(celang(3) - 90.0) .lt. tol)) then
77        ctype = 'sc'
78        lv = cellm(1)
79        return
80      endif
81
82C Check if angles are 60 deg, in which case the cell is fcc
83      if ((abs(celang(1) - 60.0) .lt. tol) .and.
84     .    (abs(celang(2) - 60.0) .lt. tol) .and.
85     .    (abs(celang(3) - 60.0) .lt. tol)) then
86        ctype = 'fcc'
87        lv = sqrt(2.0d0)*cellm(1)
88        return
89      endif
90
91C Check if angles are 109.47122 deg, in which case the cell is bcc
92      if ((abs(celang(1) - acos(-1./3.)*180.0/pi) .lt. tol) .and.
93     .    (abs(celang(2) - acos(-1./3.)*180.0/pi) .lt. tol) .and.
94     .    (abs(celang(3) - acos(-1./3.)*180.0/pi) .lt. tol)) then
95        ctype = 'bcc'
96        lv = (2.0d0/sqrt(3.0d0))*cellm(1)
97        return
98      endif
99
100      return
101      end
102