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 uncell(rv,a,b,c,alpha,beta,gamma,degtorad) 9C 10C Convert cell vectors to parameters 11C 12C Julian Gale, April 1992 13C 14 use precision 15 implicit none 16C 17C Passed variables 18C 19 real(dp) :: a 20 real(dp) :: b 21 real(dp) :: c 22 real(dp) :: alpha 23 real(dp) :: beta 24 real(dp) :: gamma 25 real(dp) :: degtorad 26 real(dp) :: rv(3,3) 27C 28C Local variables 29C 30 integer :: i 31 integer :: j 32 real(dp) :: temp(6) 33C 34 do i = 1,3 35 temp(i) = 0.0_dp 36 do j = 1,3 37 temp(i) = temp(i) + rv(j,i)**2 38 enddo 39 temp(i) = sqrt(temp(i)) 40 enddo 41 a = abs(temp(1)) 42 b = abs(temp(2)) 43 c = abs(temp(3)) 44 do i = 1,3 45 temp(3+i) = 0.0_dp 46 enddo 47 do j = 1,3 48 temp(4) = temp(4) + rv(j,2)*rv(j,3) 49 temp(5) = temp(5) + rv(j,1)*rv(j,3) 50 temp(6) = temp(6) + rv(j,1)*rv(j,2) 51 enddo 52 temp(4) = temp(4)/(temp(2)*temp(3)) 53 temp(5) = temp(5)/(temp(1)*temp(3)) 54 temp(6) = temp(6)/(temp(1)*temp(2)) 55 alpha = acos(temp(4))/degtorad 56 beta = acos(temp(5))/degtorad 57 gamma = acos(temp(6))/degtorad 58C 59C Avoid round off errors for 90.0 and 120.0 degrees 60C 61 if (abs(alpha-90.0).lt.0.00001) alpha = 90.0_dp 62 if (abs(alpha-120.0).lt.0.00001) alpha = 120.0_dp 63 if (abs(beta-90.0).lt.0.00001) beta = 90.0_dp 64 if (abs(beta-120.0).lt.0.00001) beta = 120.0_dp 65 if (abs(gamma-90.0).lt.0.00001) gamma = 90.0_dp 66 if (abs(gamma-120.0).lt.0.00001) gamma = 120.0_dp 67C 68 return 69 end 70C 71 subroutine cellimagelist(ucell,xvec1cell,yvec1cell,zvec1cell, 72 . ixvec1cell,iyvec1cell,izvec1cell) 73C 74C Store linear array of lattice vectors for ii/jj/kk=-1,1 75C => 27 lattice vectors. This tidies up and speeds up the 76C loops over these indices. 77C 78C Julian Gale, NRI, Curtin University, March 2004 79C 80 use precision 81 implicit none 82C 83C Passed variables 84C 85 integer :: ixvec1cell(27) 86 integer :: iyvec1cell(27) 87 integer :: izvec1cell(27) 88 real(dp) :: ucell(3,3) 89 real(dp) :: xvec1cell(27) 90 real(dp) :: yvec1cell(27) 91 real(dp) :: zvec1cell(27) 92C 93C Local variables 94C 95 integer :: ii 96 integer :: iimax 97 integer :: jj 98 integer :: kk 99 real(dp) :: xcdi 100 real(dp) :: ycdi 101 real(dp) :: zcdi 102 real(dp) :: xcdj 103 real(dp) :: ycdj 104 real(dp) :: zcdj 105 real(dp) :: xcrd 106 real(dp) :: ycrd 107 real(dp) :: zcrd 108C 109!! Why this -2 ? 110!! (and the rest below: Simply to modify 111!! the indexes) 112!! The final ordering is: 113!! 1 --> (-1,-1,-1) 114!! 2 --> (-1,-1, 0) 115!! 3 --> (-1,-1, 1) 116!! 4 --> (-1, 0,-1) 117! ... 118! So it might be clearer to do this: 119c$$$ n = 0 120c$$$ a1 => ucell(:,1) ... etc 121c$$$ do kk = -1, 1 122c$$$ do jj = -1, 1 123c$$$ do ii = -1, 1 124c$$$ n = n + 1 125c$$$ vlin(1:3) = ii*a1 + jj*a2 + kk*a3 126c$$$ xvec1(n) = vlin(1) ... etc 127c$$$ enddo 128c$$$ enddo 129c$$$ enddo 130c$$$ 131!! 132!! vector cdi = -2*a1 133 xcdi = -2.0d0*ucell(1,1) 134 ycdi = -2.0d0*ucell(2,1) 135 zcdi = -2.0d0*ucell(3,1) 136 iimax = 0 137C 138C Loop over unit cells 139C 140 do ii = -1,1 141 ! cdi = cdi + a1 142 xcdi = xcdi + ucell(1,1) 143 ycdi = ycdi + ucell(2,1) 144 zcdi = zcdi + ucell(3,1) 145 ! cdj = cdi - 2*a2 146 xcdj = xcdi - 2.0d0*ucell(1,2) 147 ycdj = ycdi - 2.0d0*ucell(2,2) 148 zcdj = zcdi - 2.0d0*ucell(3,2) 149 do jj = -1,1 150 ! cdj = cdj + a2 151 xcdj = xcdj + ucell(1,2) 152 ycdj = ycdj + ucell(2,2) 153 zcdj = zcdj + ucell(3,2) 154 ! crd = cdj - 2*a3 155 xcrd = xcdj - 2.0d0*ucell(1,3) 156 ycrd = ycdj - 2.0d0*ucell(2,3) 157 zcrd = zcdj - 2.0d0*ucell(3,3) 158 do kk = -1,1 159 iimax = iimax + 1 160 ! crd = crd + a3 161 xcrd = xcrd + ucell(1,3) 162 ycrd = ycrd + ucell(2,3) 163 zcrd = zcrd + ucell(3,3) 164 ixvec1cell(iimax) = ii 165 iyvec1cell(iimax) = jj 166 izvec1cell(iimax) = kk 167 ! Store successive points in linear arrays 168 ! vec1(iimax) = crd 169 xvec1cell(iimax) = xcrd 170 yvec1cell(iimax) = ycrd 171 zvec1cell(iimax) = zcrd 172 enddo 173 enddo 174 enddo 175C 176 return 177 end 178