1*
2* $Id$
3*
4      subroutine incell2(ni,r2,r1)
5
6*     =======================================================
7*	This routine places the ions inside the cell defined
8*	by the lattice vectors centered at zero.
9*     =======================================================
10      implicit none
11      integer ni
12      real*8  r2(3,*),r1(3,*)
13
14*     **** Local variables defined ****
15      real*8  fa1,fa2,fa3
16      real*8  a(3,3),b(3,3),volume
17      integer i,j
18
19*      **** external functions ****
20       real*8   lattice_unita
21       external lattice_unita
22
23*     ***** Determine the unit lattice vectors and distances ******
24      do j=1,3
25      do i=1,3
26        a(i,j) = lattice_unita(i,j)
27      end do
28      end do
29
30      b(1,1) = a(2,2)*a(3,3) - a(3,2)*a(2,3)
31      b(2,1) = a(3,2)*a(1,3) - a(1,2)*a(3,3)
32      b(3,1) = a(1,2)*a(2,3) - a(2,2)*a(1,3)
33      b(1,2) = a(2,3)*a(3,1) - a(3,3)*a(2,1)
34      b(2,2) = a(3,3)*a(1,1) - a(1,3)*a(3,1)
35      b(3,2) = a(1,3)*a(2,1) - a(2,3)*a(1,1)
36      b(1,3) = a(2,1)*a(3,2) - a(3,1)*a(2,2)
37      b(2,3) = a(3,1)*a(1,2) - a(1,1)*a(3,2)
38      b(3,3) = a(1,1)*a(2,2) - a(2,1)*a(1,2)
39      volume = a(1,1)*b(1,1)
40     >       + a(2,1)*b(2,1)
41     >       + a(3,1)*b(3,1)
42
43      volume = 1.0d0/volume
44      call dscal(9,volume,b,1)
45
46!$OMP DO
47      do i =1,ni
48
49*        *** Break the Ion positions into the a1, a2, and a3 components ***
50         fa1 =  b(1,1) * r2(1,i)
51     >       +  b(2,1) * r2(2,i)
52     >       +  b(3,1) * r2(3,i)
53
54         fa2 =  b(1,2) * r2(1,i)
55     >       +  b(2,2) * r2(2,i)
56     >       +  b(3,2) * r2(3,i)
57
58         fa3 =  b(1,3) * r2(1,i)
59     >       +  b(2,3) * r2(2,i)
60     >       +  b(3,3) * r2(3,i)
61
62
63*	**** Change the a1, a2 and a3  components to ****
64*	**** make the ion be in the cell             ****
65
66   23   IF (fa1 .GT. (0.5d0)) THEN
67*          WRITE (*,*) 'a1>', I, R2A1, DA(1)/2.0d0
68           r2(1,i) = r2(1,i) - lattice_unita(1,1)
69           r2(2,i) = r2(2,i) - lattice_unita(2,1)
70           r2(3,i) = r2(3,i) - lattice_unita(3,1)
71
72           r1(1,i) = r1(1,i) - lattice_unita(1,1)
73           r1(2,i) = r1(2,i) - lattice_unita(2,1)
74           r1(3,i) = r1(3,i) - lattice_unita(3,1)
75
76
77           fa1 =  b(1,1) * r2(1,i)
78     >         +  b(2,1) * r2(2,i)
79     >         +  b(3,1) * r2(3,i)
80           GO TO 23
81        ENDIF
82
83   24   IF (fa1 .LE. (-0.5d0)) THEN
84*          WRITE (*,*) 'a1<', I, R2A1, DA(1)/2.0d0
85           r2(1,i) = r2(1,i) + lattice_unita(1,1)
86           r2(2,i) = r2(2,i) + lattice_unita(2,1)
87           r2(3,i) = r2(3,i) + lattice_unita(3,1)
88
89           r1(1,i) = r1(1,i) + lattice_unita(1,1)
90           r1(2,i) = r1(2,i) + lattice_unita(2,1)
91           r1(3,i) = r1(3,i) + lattice_unita(3,1)
92
93
94           fa1 =  b(1,1) * r2(1,i)
95     >         +  b(2,1) * r2(2,i)
96     >         +  b(3,1) * r2(3,i)
97            GO TO 24
98        ENDIF
99
100   25   IF (fa2 .GT. (0.5d0)) THEN
101*          WRITE (*,*) 'a2>', I, R2A2, DA(2)/2.0d0
102           r2(1,i) = r2(1,i) - lattice_unita(1,2)
103           r2(2,i) = r2(2,i) - lattice_unita(2,2)
104           r2(3,i) = r2(3,i) - lattice_unita(3,2)
105
106           r1(1,i) = r1(1,i) - lattice_unita(1,2)
107           r1(2,i) = r1(2,i) - lattice_unita(2,2)
108           r1(3,i) = r1(3,i) - lattice_unita(3,2)
109
110
111           fa2 =  b(1,2) * r2(1,i)
112     >         +  b(2,2) * r2(2,i)
113     >         +  b(3,2) * r2(3,i)
114          GO TO 25
115        ENDIF
116
117   26   IF (fa2 .LE. (-0.5d0)) THEN
118*          WRITE (*,*) 'a2<', I, R2A2, DA(2)/2.0d0
119           r2(1,i) = r2(1,i) + lattice_unita(1,2)
120           r2(2,i) = r2(2,i) + lattice_unita(2,2)
121           r2(3,i) = r2(3,i) + lattice_unita(3,2)
122
123           r1(1,i) = r1(1,i) + lattice_unita(1,2)
124           r1(2,i) = r1(2,i) + lattice_unita(2,2)
125           r1(3,i) = r1(3,i) + lattice_unita(3,2)
126
127
128           fa2 =  b(1,2) * r2(1,i)
129     >         +  b(2,2) * r2(2,i)
130     >         +  b(3,2) * r2(3,i)
131           GO TO 26
132        ENDIF
133
134
135   27   IF (fa3 .GT. (0.5d0)) THEN
136*         WRITE (*,*) 'a3>', i, R2A3, DA(3)/2.0d0
137          r2(1,i) = r2(1,i) - lattice_unita(1,3)
138          r2(2,i) = r2(2,i) - lattice_unita(2,3)
139          r2(3,i) = r2(3,i) - lattice_unita(3,3)
140
141          r1(1,i) = r1(1,i) - lattice_unita(1,3)
142          r1(2,i) = r1(2,i) - lattice_unita(2,3)
143          r1(3,i) = r1(3,i) - lattice_unita(3,3)
144
145          fa3 =  b(1,3) * r2(1,i)
146     >        +  b(2,3) * r2(2,i)
147     >        +  b(3,3) * r2(3,i)
148           GO TO 27
149        ENDIF
150
151   28   IF (fa3 .LE. (-0.5d0)) THEN
152*         WRITE (*,*) 'a3<', I, R2A3, DA(3)/2.0d0
153          r2(1,i) = r2(1,i) + lattice_unita(1,3)
154          r2(2,i) = r2(2,i) + lattice_unita(2,3)
155          r2(3,i) = r2(3,i) + lattice_unita(3,3)
156
157          r1(1,i) = r1(1,i) + lattice_unita(1,3)
158          r1(2,i) = r1(2,i) + lattice_unita(2,3)
159          r1(3,i) = r1(3,i) + lattice_unita(3,3)
160
161          fa3 =  b(1,3) * r2(1,i)
162     >        +  b(2,3) * r2(2,i)
163     >        +  b(3,3) * r2(3,i)
164           GO TO 28
165        ENDIF
166
167      end do
168!$OMP END DO
169
170      return
171      end
172