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