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