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