1      subroutine shake_chain(n,indx,nb,
2     >                       tol,maxit,
3     >                       dsq,mass,
4     >                       r2,r1,gablambda)
5      implicit none
6      integer n,indx(*),nb
7      real*8  tol
8      integer maxit
9      real*8  dsq(*),mass(*)
10      real*8  r2(3,*), r1(3,*)
11      real*8 gablambda
12
13
14        REAL*8       RXI(n), RYI(n), RZI(n)
15        REAL*8       PXI(n), PYI(n), PZI(n)
16        LOGICAL     MOVING(n)
17        LOGICAL     MOVED(n)
18
19        LOGICAL     DONE
20        INTEGER     IT, A, B, I,J
21        REAL*8      PXAB, PYAB, PZAB, PABSQ
22        REAL*8      RXAB, RYAB, RZAB, RABSQ, DIFFSQ, RPAB
23        REAL*8      GAB, DX, DY, DZ, TOL2
24        REAL*8      RPTOL, RMA, RMB
25        PARAMETER ( RPTOL = 1.0E-6 )
26
27c        real*8 ua(3,3),ub(3,3),volume
28c        real*8 c1,c2,c3
29
30c*     **** external functions ****
31c      real*8   lattice_unita
32c      external lattice_unita
33
34
35*     ***** Determine the unit lattice vectors and distances ******
36c      do j=1,3
37c      do i=1,3
38c         ua(i,j) = lattice_unita(i,j)
39c      end do
40c      end do
41c      ub(1,1) = ua(2,2)*ua(3,3) - ua(3,2)*ua(2,3)
42c      ub(2,1) = ua(3,2)*ua(1,3) - ua(1,2)*ua(3,3)
43c      ub(3,1) = ua(1,2)*ua(2,3) - ua(2,2)*ua(1,3)
44c      ub(1,2) = ua(2,3)*ua(3,1) - ua(3,3)*ua(2,1)
45c      ub(2,2) = ua(3,3)*ua(1,1) - ua(1,3)*ua(3,1)
46c      ub(3,2) = ua(1,3)*ua(2,1) - ua(2,3)*ua(1,1)
47c      ub(1,3) = ua(2,1)*ua(3,2) - ua(3,1)*ua(2,2)
48c      ub(2,3) = ua(3,1)*ua(1,2) - ua(1,1)*ua(3,2)
49c      ub(3,3) = ua(1,1)*ua(2,2) - ua(2,1)*ua(1,2)
50c      volume = ua(1,1)*ub(1,1)
51c     >       + ua(2,1)*ub(2,1)
52c     >       + ua(3,1)*ub(3,1)
53c      volume = 1.0d0/volume
54c      call dscal(9,volume,ub,1)
55
56
57
58        TOL2   = 2.0 * tol
59        !TOL2   = 1.0d-15
60
61
62        gablambda = 0.0d0
63        do A = 1, n
64              RXI(A) = r1(1,indx(A))
65              RYI(A) = r1(2,indx(A))
66              RZI(A) = r1(3,indx(A))
67              PXI(A) = r2(1,indx(A))
68              PYI(A) = r2(2,indx(A))
69              PZI(A) = r2(3,indx(A))
70
71              MOVING(A) = .FALSE.
72              MOVED(A)  = .TRUE.
73        end do
74
75           IT = 0
76           DONE = .FALSE.
77
78C       ** BEGIN ITERATIVE LOOP **
79
801000       IF ( ( .NOT. DONE ) .AND. ( IT .LE. MAXIT ) ) THEN
81
82              DONE = .TRUE.
83
84              DO 300 A = 1, NB
85
86                 B = A + 1
87                 IF ( B .GT. N ) B = 1
88
89                 IF ( MOVED(A) .OR. MOVED(B) ) THEN
90
91                    PXAB = PXI(A) - PXI(B)
92                    PYAB = PYI(A) - PYI(B)
93                    PZAB = PZI(A) - PZI(B)
94                    call lattice_min_difference(PXAB,PYAB,PZAB)
95c                    c1 = PXAB*ub(1,1) + PYAB*ub(2,1) + PZAB*ub(3,1)
96c                    c2 = PXAB*ub(1,2) + PYAB*ub(2,2) + PZAB*ub(3,2)
97c                    c3 = PXAB*ub(1,3) + PYAB*ub(2,3) + PZAB*ub(3,3)
98c                    c1 = c1 - ANINT(c1)
99c                    c2 = c2 - ANINT(c2)
100c                    c3 = c3 - ANINT(c3)
101c                    PXAB = ua(1,1)*c1 + ua(1,2)*c2 + ua(1,3)*c3
102c                    PYAB = ua(2,1)*c1 + ua(2,2)*c2 + ua(2,3)*c3
103c                    PZAB = ua(3,1)*c1 + ua(3,2)*c2 + ua(3,3)*c3
104
105
106                    PABSQ  = PXAB ** 2 + PYAB ** 2 + PZAB ** 2
107                    RABSQ  = DSQ(A)
108                    DIFFSQ = RABSQ - PABSQ
109
110
111                    IF ( ABS(DIFFSQ) .GT. ( RABSQ * TOL2 ) ) THEN
112
113                       RXAB = RXI(A) - RXI(B)
114                       RYAB = RYI(A) - RYI(B)
115                       RZAB = RZI(A) - RZI(B)
116                       call lattice_min_difference(RXAB,RYAB,RZAB)
117c                       c1 = RXAB*ub(1,1) + RYAB*ub(2,1) + RZAB*ub(3,1)
118c                       c2 = RXAB*ub(1,2) + RYAB*ub(2,2) + RZAB*ub(3,2)
119c                       c3 = RXAB*ub(1,3) + RYAB*ub(2,3) + RZAB*ub(3,3)
120c                       c1 = c1 - ANINT(c1)
121c                       c2 = c2 - ANINT(c2)
122c                       c3 = c3 - ANINT(c3)
123c                       RXAB = ua(1,1)*c1 + ua(1,2)*c2 + ua(1,3)*c3
124c                       RYAB = ua(2,1)*c1 + ua(2,2)*c2 + ua(2,3)*c3
125c                       RZAB = ua(3,1)*c1 + ua(3,2)*c2 + ua(3,3)*c3
126
127                       RPAB = RXAB * PXAB + RYAB * PYAB + RZAB * PZAB
128
129                       IF ( RPAB .LT. ( RABSQ * RPTOL ) ) THEN
130
131                          STOP 'CONSTRAINT FAILURE '
132
133                       ENDIF
134
135                       RMA = 1.0 / mass(A)
136                       RMB = 1.0 / mass(B)
137                       GAB = DIFFSQ / ( 2.0 * ( RMA + RMB ) * RPAB )
138                       DX  = RXAB * GAB
139                       DY  = RYAB * GAB
140                       DZ  = RZAB * GAB
141                       gablambda = gablambda + GAB
142
143                       PXI(A) = PXI(A) + RMA * DX
144                       PYI(A) = PYI(A) + RMA * DY
145                       PZI(A) = PZI(A) + RMA * DZ
146                       PXI(B) = PXI(B) - RMB * DX
147                       PYI(B) = PYI(B) - RMB * DY
148                       PZI(B) = PZI(B) - RMB * DZ
149
150                       MOVING(A) = .TRUE.
151                       MOVING(B) = .TRUE.
152                       DONE = .FALSE.
153
154                    ENDIF
155
156                 ENDIF
157
158300           CONTINUE
159
160              DO 400 A = 1, n
161                 MOVED(A) = MOVING(A)
162                 MOVING(A) = .FALSE.
163400           CONTINUE
164
165              IT = IT + 1
166              GOTO 1000
167
168           ENDIF
169
170C       ** END ITERATIVE LOOP **
171
172           IF ( .NOT. DONE ) THEN
173
174              WRITE(*,'('' TOO MANY CONSTRAINT ITERATIONS '')')
175              STOP
176
177           ENDIF
178
179           do a = 1, n
180              r2(1,indx(a)) = PXI(a)
181              r2(2,indx(a)) = PYI(a)
182              r2(3,indx(a)) = PZI(a)
183              r1(1,indx(a)) = RXI(a)
184              r1(2,indx(a)) = RYI(a)
185              r1(3,indx(a)) = RZI(a)
186           end do
187
188        RETURN
189        END
190c $Id$
191