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 cart2frac(na,xc,yc,zc,rv,xyzfrac)
9C
10C  Converts Cartesian coordinates to fractional coordinates
11C
12C  On entry :
13C
14C  na            = number of atoms
15C  xc(na)        = Cartesian X coordinate
16C  yc(na)        = Cartesian Y coordinate
17C  zc(na)        = Cartesian Z coordinate
18C  rv(3,3)       = cell vectors
19C
20C  On exit :
21C
22C  xyzfrac(3,na) = fractional X coordinate
23C
24C  Julian Gale, NRI, Curtin University, May 2004
25C
26      use precision
27
28      implicit none
29C
30C  Passed variables
31C
32      integer, intent(in)     :: na
33      real(dp),  intent(inout)  :: xc(na)
34      real(dp),  intent(inout)  :: yc(na)
35      real(dp),  intent(inout)  :: zc(na)
36      real(dp),  intent(out)    :: xyzfrac(3,na)
37      real(dp),  intent(in)     :: rv(3,3)
38C
39C  Local variables
40C
41      integer                 :: i
42      real(dp)                  :: rmat(3,3)
43C
44C  Copy lattice vectors to scratch array
45C
46      do i = 1,3
47        rmat(1,i) = rv(1,i)
48        rmat(2,i) = rv(2,i)
49        rmat(3,i) = rv(3,i)
50      enddo
51C
52C  Copy Cartesian coordinates to fractional array
53C
54      do i = 1,na
55        xyzfrac(1,i) = xc(i)
56        xyzfrac(2,i) = yc(i)
57        xyzfrac(3,i) = zc(i)
58      enddo
59C
60C  Convert cartesian coordinates to fractional for 3D case
61C
62      call gaussxyz(na,rmat,xyzfrac)
63C
64C  Place fractional coordinates in the range 0 to 1
65C
66      do i = 1,na
67        xyzfrac(1,i) = dmod(xyzfrac(1,i)+10.0d0,1.0d0)
68        xyzfrac(2,i) = dmod(xyzfrac(2,i)+10.0d0,1.0d0)
69        xyzfrac(3,i) = dmod(xyzfrac(3,i)+10.0d0,1.0d0)
70      enddo
71C
72C  Convert fractional back to Cartesian to place within cell
73C
74      do i = 1,na
75        xc(i) = xyzfrac(1,i)*rv(1,1) +
76     .          xyzfrac(2,i)*rv(1,2) +
77     .          xyzfrac(3,i)*rv(1,3)
78        yc(i) = xyzfrac(1,i)*rv(2,1) +
79     .          xyzfrac(2,i)*rv(2,2) +
80     .          xyzfrac(3,i)*rv(2,3)
81        zc(i) = xyzfrac(1,i)*rv(3,1) +
82     .          xyzfrac(2,i)*rv(3,2) +
83     .          xyzfrac(3,i)*rv(3,3)
84      enddo
85C
86      return
87      end
88C
89      subroutine gaussxyz(m,a,xf)
90C
91C  Invert matrix A by Gaussian elimination and apply to
92C  multiple vectors x
93C
94      use precision
95      use sys, only : die
96      implicit none
97C
98C  Passed variables
99C
100      integer, intent(in)     :: m
101      real(dp),  intent(inout)  :: a(3,3)
102      real(dp),  intent(inout)  :: xf(3,m)
103C
104C  Local variables
105C
106      integer                 :: i
107      integer                 :: ie
108      integer                 :: in
109      integer                 :: ix
110      integer                 :: j
111      integer                 :: k
112      integer                 :: kk
113      real(dp)                  :: delt
114      real(dp)                  :: u
115      real(dp)                  :: x
116C
117      delt = 1.0d-10
118      do k = 1,2
119        u = abs(a(k,k))
120        kk = k + 1
121        in = k
122C
123C  Search for index in of maximum pivot value
124C
125        do i = kk,3
126          if (abs(a(i,k)).gt.u) then
127            u = abs(a(i,k))
128            in = i
129          endif
130        enddo
131        if (k.ne.in) then
132C
133C  Interchange rows k and index in
134C
135          do j = k,3
136            x = a(k,j)
137            a(k,j) = a(in,j)
138            a(in,j) = x
139          enddo
140          do j = 1,m
141            x = xf(k,j)
142            xf(k,j) = xf(in,j)
143            xf(in,j) = x
144          enddo
145        endif
146C
147C  Check to see if pivot value is too small
148C
149        if (u.lt.delt) then
150          call die('Cell matrix is singular')
151        endif
152C
153C  Forward elimination step
154C
155        do j = kk,3
156          do i = kk,3
157            a(i,j) = a(i,j) - a(i,k)*a(k,j)/a(k,k)
158          enddo
159        enddo
160        do j = 1,m
161          do i = kk,3
162            xf(i,j) = xf(i,j) - a(i,k)*xf(k,j)/a(k,k)
163          enddo
164        enddo
165      enddo
166      if (abs(a(3,3)).lt.delt) then
167        call die('Cell matrix is singular')
168      endif
169C
170C  Back substitution
171C
172      do k = 1,m
173        xf(3,k) = xf(3,k)/a(3,3)
174      	do ie = 1,2
175      	  i = 3 - ie
176      	  ix = i + 1
177      	  do j = ix,3
178            xf(i,k) = xf(i,k) - xf(j,k)*a(i,j)
179          enddo
180      	  xf(i,k) = xf(i,k)/a(i,i)
181        enddo
182      enddo
183C
184      return
185      end
186