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