1c
2c $Id$
3c
4
5*     ********************************************
6*     *                                          *
7*     *           shake_bonddiff                 *
8*     *                                          *
9*     ********************************************
10*
11*   This routine handles the bond difference constraint.
12*
13*   sigma = |Ra-Rb|^2 - |Rc-Rd|^2 - gamma
14*
15*   Entry - na - number of atoms (na==3) or (na==4)
16*           indx(*) - indexes of the na atoms a,b,c,d
17*           tol - tolerance of lagrange multiplier iteration
18*           maxit - number iterations to solve for Lagrange multiplier
19*           gamma - constraint value
20*           mass(*) - masses of the atoms a,b,c,d
21*           r2(3,*) - partially updated positions
22*           r1(3,*) - previous positions
23*           x - lagrange multiplier constraint
24*
25      subroutine shake_bonddiff(na,indx,
26     >                       tol,maxit,
27     >                       gamma,mass,
28     >                       r2,r1,x)
29      implicit none
30      integer na
31      integer indx(4)
32      real*8  tol
33      integer maxit
34      real*8  gamma,mass(4)
35      real*8  r2(3,*), r1(3,*)
36      real*8  x
37
38      logical done
39      integer it
40      real*8 Rab(3),Rcd(3),tRab(3),tRcd(3),mabinv,mcdinv
41      real*8 Rac(3),Rbc(3),tRac(3),tRbc(3),macinv,mbcinv
42      real*8 mainv,mbinv,mcinv
43      real*8 A,B,C,xold,dsc,y1,y2
44
45      if (na.eq.3) then
46         Rac(1) = r1(1,indx(1)) - r1(1,indx(3))
47         Rac(2) = r1(2,indx(1)) - r1(2,indx(3))
48         Rac(3) = r1(3,indx(1)) - r1(3,indx(3))
49         call lattice_min_difference(Rac(1),Rac(2),Rac(3))
50
51         Rbc(1) = r1(1,indx(2)) - r1(1,indx(3))
52         Rbc(2) = r1(2,indx(2)) - r1(2,indx(3))
53         Rbc(3) = r1(3,indx(2)) - r1(3,indx(3))
54         call lattice_min_difference(Rbc(1),Rbc(2),Rbc(3))
55
56         tRac(1) = r2(1,indx(1)) - r2(1,indx(3))
57         tRac(2) = r2(2,indx(1)) - r2(2,indx(3))
58         tRac(3) = r2(3,indx(1)) - r2(3,indx(3))
59         call lattice_min_difference(tRac(1),tRac(2),tRac(3))
60
61         tRbc(1) = r2(1,indx(2)) - r2(1,indx(3))
62         tRbc(2) = r2(2,indx(2)) - r2(2,indx(3))
63         tRbc(3) = r2(3,indx(2)) - r2(3,indx(3))
64         call lattice_min_difference(tRbc(1),tRbc(2),tRbc(3))
65
66         macinv = 1.0d0/mass(1) + 1.0d0/mass(3)
67         mbcinv = 1.0d0/mass(2) + 1.0d0/mass(3)
68         mainv  = 1.0d0/mass(1)
69         mbinv  = 1.0d0/mass(2)
70         mcinv  = 1.0d0/mass(3)
71
72         A = (macinv**2 - mcinv**2)*(Rac(1)**2+Rac(2)**2+Rac(3)**2)
73     >     + (mcinv**2 - mbcinv**2)*(Rbc(1)**2+Rbc(2)**2+Rbc(3)**2)
74     >     + 2.0d0*mcinv*(mbinv+mainv)
75     >       *(Rac(1)*Rbc(1)+Rac(2)*Rbc(2)+Rac(3)*Rbc(3))
76
77         B = 2.0d0*macinv*(tRac(1)*Rac(1)+tRac(2)*Rac(2)+tRac(3)*Rac(3))
78     >     - 2.0d0*mcinv *(tRac(1)*Rbc(1)+tRac(2)*Rbc(2)+tRac(3)*Rbc(3))
79     >     + 2.0d0*mbcinv*(tRbc(1)*Rbc(1)+tRbc(2)*Rbc(2)+tRbc(3)*Rbc(3))
80     >     - 2.0d0*mcinv *(tRbc(1)*Rac(1)+tRbc(2)*Rac(2)+tRbc(3)*Rac(3))
81
82         C = (tRac(1)**2 + tRac(2)**2 + tRac(3)**2)
83     >     - (tRbc(1)**2 + tRbc(2)**2 + tRbc(3)**2) - gamma
84
85
86         if (dabs(A).gt.1.0d-12) then
87            x = -C
88            it = 0
89            done = .false.
90            do while (.not. done)
91               xold = x
92               x    = -C + (1.0d0-B)*x - A*x*x
93               it = it + 1
94               done = (dabs(x-xold).lt.tol).or.(it.gt.maxit)
95            end do
96
97            dsc = B*B-4.0d0*A*C
98            if (dsc.ge.0.0d0) then
99               y1 = (-B + dsqrt(dsc))/(2.0d0*A)
100               y2 = (-B - dsqrt(dsc))/(2.0d0*A)
101               if (dabs(x-y1).lt.1.0d-1) then
102                  x = y1
103               else if (dabs(x-y2).lt.1.0d-1) then
104                  x = y2
105               end if
106            end if
107         else
108            x = -C/B
109            it = 0
110         end if
111
112         r2(1,indx(1)) = r2(1,indx(1)) + x*Rac(1)*mainv
113         r2(2,indx(1)) = r2(2,indx(1)) + x*Rac(2)*mainv
114         r2(3,indx(1)) = r2(3,indx(1)) + x*Rac(3)*mainv
115
116         r2(1,indx(2)) = r2(1,indx(2)) - x*Rbc(1)*mbinv
117         r2(2,indx(2)) = r2(2,indx(2)) - x*Rbc(2)*mbinv
118         r2(3,indx(2)) = r2(3,indx(2)) - x*Rbc(3)*mbinv
119
120         r2(1,indx(3)) = r2(1,indx(3)) + x*(Rbc(1)-Rac(1))*mcinv
121         r2(2,indx(3)) = r2(2,indx(3)) + x*(Rbc(2)-Rac(2))*mcinv
122         r2(3,indx(3)) = r2(3,indx(3)) + x*(Rbc(3)-Rac(3))*mcinv
123
124c         tRac(1) = r2(1,indx(1)) - r2(1,indx(3))
125c         tRac(2) = r2(2,indx(1)) - r2(2,indx(3))
126c         tRac(3) = r2(3,indx(1)) - r2(3,indx(3))
127c
128c         tRbc(1) = r2(1,indx(2)) - r2(1,indx(3))
129c         tRbc(2) = r2(2,indx(2)) - r2(2,indx(3))
130c         tRbc(3) = r2(3,indx(2)) - r2(3,indx(3))
131c
132c         C = (tRac(1)**2 + tRac(2)**2 + tRac(3)**2)
133c     >     - (tRbc(1)**2 + tRbc(2)**2 + tRbc(3)**2) - gamma
134c
135c         write(*,*) "x,y1,y2,constraint=",x,y1,y2,C
136
137      else
138         Rab(1) = r1(1,indx(1)) - r1(1,indx(2))
139         Rab(2) = r1(2,indx(1)) - r1(2,indx(2))
140         Rab(3) = r1(3,indx(1)) - r1(3,indx(2))
141         call lattice_min_difference(Rab(1),Rab(2),Rab(3))
142
143         Rcd(1) = r1(1,indx(3)) - r1(1,indx(4))
144         Rcd(2) = r1(2,indx(3)) - r1(2,indx(4))
145         Rcd(3) = r1(3,indx(3)) - r1(3,indx(4))
146         call lattice_min_difference(Rcd(1),Rcd(2),Rcd(3))
147
148         tRab(1) = r2(1,indx(1)) - r2(1,indx(2))
149         tRab(2) = r2(2,indx(1)) - r2(2,indx(2))
150         tRab(3) = r2(3,indx(1)) - r2(3,indx(2))
151         call lattice_min_difference(tRab(1),tRab(2),tRab(3))
152
153         tRcd(1) = r2(1,indx(3)) - r2(1,indx(4))
154         tRcd(2) = r2(2,indx(3)) - r2(2,indx(4))
155         tRcd(3) = r2(3,indx(3)) - r2(3,indx(4))
156         call lattice_min_difference(tRcd(1),tRcd(2),tRcd(3))
157
158         mabinv = 1.0d0/mass(1) + 1.0d0/mass(2)
159         mcdinv = 1.0d0/mass(3) + 1.0d0/mass(4)
160
161         A = (mabinv**2)*(Rab(1)**2+Rab(2)**2+Rab(3)**2)
162     >     - (mcdinv**2)*(Rcd(1)**2+Rcd(2)**2+Rcd(3)**2)
163
164         B = 2.0d0*mabinv*(tRab(1)*Rab(1)+tRab(2)*Rab(2)+tRab(3)*Rab(3))
165     >     + 2.0d0*mcdinv*(tRcd(1)*Rcd(1)+tRcd(2)*Rcd(2)+tRcd(3)*Rcd(3))
166
167         C = (tRab(1)**2 + tRab(2)**2 + tRab(3)**2)
168     >     - (tRcd(1)**2 + tRcd(2)**2 + tRcd(3)**2) - gamma
169
170         if (dabs(A).gt.1.0d-12) then
171            x = -C
172            it = 0
173            done = .false.
174            do while (.not. done)
175               xold = x
176               x    = -C + (1.0d0-B)*x - A*x*x
177               it = it + 1
178               done = (dabs(x-xold).lt.tol).or.(it.gt.maxit)
179            end do
180
181            dsc = B*B-4.0d0*A*C
182            if (dsc.ge.0.0d0) then
183               y1 = (-B + dsqrt(dsc))/(2.0d0*A)
184               y2 = (-B - dsqrt(dsc))/(2.0d0*A)
185               if (dabs(x-y1).lt.1.0d-1) then
186                  x = y1
187               else if (dabs(x-y2).lt.1.0d-1) then
188                  x = y2
189               end if
190            end if
191         else
192            x = -C/B
193            it = 0
194         end if
195
196         r2(1,indx(1)) = r2(1,indx(1)) + x*Rab(1)/mass(1)
197         r2(2,indx(1)) = r2(2,indx(1)) + x*Rab(2)/mass(1)
198         r2(3,indx(1)) = r2(3,indx(1)) + x*Rab(3)/mass(1)
199
200         r2(1,indx(2)) = r2(1,indx(2)) - x*Rab(1)/mass(2)
201         r2(2,indx(2)) = r2(2,indx(2)) - x*Rab(2)/mass(2)
202         r2(3,indx(2)) = r2(3,indx(2)) - x*Rab(3)/mass(2)
203
204         r2(1,indx(3)) = r2(1,indx(3)) - x*Rcd(1)/mass(3)
205         r2(2,indx(3)) = r2(2,indx(3)) - x*Rcd(2)/mass(3)
206         r2(3,indx(3)) = r2(3,indx(3)) - x*Rcd(3)/mass(3)
207
208         r2(1,indx(4)) = r2(1,indx(4)) + x*Rcd(1)/mass(4)
209         r2(2,indx(4)) = r2(2,indx(4)) + x*Rcd(2)/mass(4)
210         r2(3,indx(4)) = r2(3,indx(4)) + x*Rcd(3)/mass(4)
211
212      end if
213
214
215      return
216      end
217