csettle.c
1 
/*


2 
* $Id: csettle.c,v 1.19.2.1 2009/05/28 09:08:28 hess Exp $

3 
*

4 
* This source code is part of

5 
*

6 
* G R O M A C S

7 
*

8 
* GROningen MAchine for Chemical Simulations

9 
*

10 
* VERSION 3.2.0

11 
* Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.

12 
* Copyright (c) 19912000, University of Groningen, The Netherlands.

13 
* Copyright (c) 20012004, The GROMACS development team,

14 
* check out http://www.gromacs.org for more information.

15 

16 
* This program is free software; you can redistribute it and/or

17 
* modify it under the terms of the GNU General Public License

18 
* as published by the Free Software Foundation; either version 2

19 
* of the License, or (at your option) any later version.

20 
*

21 
* If you want to redistribute modifications, please consider that

22 
* scientific software is very special. Version control is crucial 

23 
* bugs must be traceable. We will be happy to consider code for

24 
* inclusion in the official distribution, but derived work must not

25 
* be called official GROMACS. Details are found in the README & COPYING

26 
* files  if they are missing, get the official version at www.gromacs.org.

27 
*

28 
* To help us fund GROMACS development, we humbly ask that you cite

29 
* the papers on the package  you can find them in the top README file.

30 
*

31 
* For more info, check our website at http://www.gromacs.org

32 
*

33 
* And Hey:

34 
* GROwing Monsters And Cloning Shrimps

35 
*/

36 
#ifdef HAVE_CONFIG_H

37 
#include <config.h> 
38 
#endif

39  
40 
#include <math.h> 
41 
#include <stdio.h> 
42 
#include "vec.h" 
43 
#include "constr.h" 
44 
#include "gmx_fatal.h" 
45  
46 
#ifdef DEBUG

47 
static void check_cons(FILE *fp,char *title,real x[],int OW1,int HW2,int HW3) 
48 
{ 
49 
rvec dOH1,dOH2,dHH; 
50 
int m;

51 

52 
for(m=0; (m<DIM); m++) { 
53 
dOH1[m]=x[OW1+m]x[HW2+m]; 
54 
dOH2[m]=x[OW1+m]x[HW3+m]; 
55 
dHH[m]=x[HW2+m]x[HW3+m]; 
56 
} 
57 
fprintf(fp,"%10s, OW1=%3d, HW2=%3d, HW3=%3d, dOH1: %8.3f, dOH2: %8.3f, dHH: %8.3f\n",

58 
title,OW1/DIM,HW2/DIM,HW3/DIM,norm(dOH1),norm(dOH2),norm(dHH)); 
59 
} 
60 
#endif

61  
62  
63 
void settle_proj(FILE *fp,int nsettle, t_iatom iatoms[],rvec x[], 
64 
real dOH,real dHH,real invmO,real invmH, 
65 
rvec *der,rvec *derp, 
66 
bool bCalcVir,tensor rmdder)

67 
{ 
68 
/* Settle for projection out constraint components

69 
* of derivatives of the coordinates.

70 
* Berk Hess 2008110

71 
*/

72  
73 
/* Initialized data */

74 
static bool bFirst=TRUE; 
75 
static real imO,imH,invdOH,invdHH;

76 
static matrix invmat;

77  
78 
real imOn,imHn; 
79 
matrix mat; 
80 
int i,m,m2,ow1,hw2,hw3;

81 
rvec roh2,roh3,rhh,dc,fc; 
82  
83 
if (bFirst) {

84 
if (fp)

85 
fprintf(fp,"Going to use settle for derivatives (%d waters)\n",nsettle);

86  
87 
imO = invmO; 
88 
imH = invmH; 
89 
/* We normalize the inverse masses with imO for the matrix inversion.

90 
* so we can keep using masses of almost zero for frozen particles,

91 
* without running out of the float range in m_inv.

92 
*/

93 
imOn = 1;

94 
imHn = imH/imO; 
95  
96 
/* Construct the constraint coupling matrix */

97 
mat[0][0] = imOn + imHn; 
98 
mat[0][1] = imOn*(1  0.5*dHH*dHH/(dOH*dOH)); 
99 
mat[0][2] = imHn*0.5*dHH/dOH; 
100 
mat[1][1] = mat[0][0]; 
101 
mat[1][2] = mat[0][2]; 
102 
mat[2][2] = imHn + imHn; 
103 
mat[1][0] = mat[0][1]; 
104 
mat[2][0] = mat[0][2]; 
105 
mat[2][1] = mat[1][2]; 
106  
107 
m_inv(mat,invmat); 
108  
109 
msmul(invmat,1/imO,invmat);

110  
111 
invdOH = 1/dOH;

112 
invdHH = 1/dHH;

113  
114 
bFirst = FALSE; 
115 
} 
116  
117 
#ifdef PRAGMAS

118 
#pragma ivdep

119 
#endif

120 
for (i=0; i<nsettle; i++) { 
121 
ow1 = iatoms[i*2+1]; 
122 
hw2 = ow1 + 1;

123 
hw3 = ow1 + 2;

124 

125 
for(m=0; m<DIM; m++) 
126 
roh2[m] = (x[ow1][m]  x[hw2][m])*invdOH; 
127 
for(m=0; m<DIM; m++) 
128 
roh3[m] = (x[ow1][m]  x[hw3][m])*invdOH; 
129 
for(m=0; m<DIM; m++) 
130 
rhh [m] = (x[hw2][m]  x[hw3][m])*invdHH; 
131 
/* 18 flops */

132  
133 
/* Determine the projections of der on the bonds */

134 
clear_rvec(dc); 
135 
for(m=0; m<DIM; m++) 
136 
dc[0] += (der[ow1][m]  der[hw2][m])*roh2[m];

137 
for(m=0; m<DIM; m++) 
138 
dc[1] += (der[ow1][m]  der[hw3][m])*roh3[m];

139 
for(m=0; m<DIM; m++) 
140 
dc[2] += (der[hw2][m]  der[hw3][m])*rhh [m];

141 
/* 27 flops */

142  
143 
/* Determine the correction for the three bonds */

144 
mvmul(invmat,dc,fc); 
145 
/* 15 flops */

146  
147 
/* Subtract the corrections from derp */

148 
for(m=0; m<DIM; m++) { 
149 
derp[ow1][m] = imO*( fc[0]*roh2[m] + fc[1]*roh3[m]); 
150 
derp[hw2][m] = imH*(fc[0]*roh2[m] + fc[2]*rhh [m]); 
151 
derp[hw3][m] = imH*(fc[1]*roh3[m]  fc[2]*rhh [m]); 
152 
} 
153 
/* 45 flops */

154  
155 
if (bCalcVir) {

156 
/* Determining r x m der is easy,

157 
* since fc contains the mass weighted corrections for der.

158 
*/

159 
for(m=0; m<DIM; m++) { 
160 
for(m2=0; m2<DIM; m2++) { 
161 
rmdder[m][m2] += 
162 
dOH*roh2[m]*roh2[m2]*fc[0] +

163 
dOH*roh3[m]*roh3[m2]*fc[1] +

164 
dHH*rhh [m]*rhh [m2]*fc[2];

165 
} 
166 
} 
167 
} 
168 
} 
169 
} 
170  
171  
172 
/* Our local shake routine to be used when settle breaks down due to a zero determinant */

173 
static int xshake(real b4[], real after[], real dOH, real dHH, real mO, real mH) 
174 
{ 
175 
real bondsq[3];

176 
real bond[9];

177 
real invmass[3];

178 
real M2[3];

179 
int iconv;

180 
int iatom[3]={0,0,1}; 
181 
int jatom[3]={1,2,2}; 
182 
real rijx,rijy,rijz,tx,ty,tz,im,jm,acor,rp,diff; 
183 
int i,ll,ii,jj,l3,ix,iy,iz,jx,jy,jz,conv;

184  
185 
invmass[0]=1.0/mO; 
186 
invmass[1]=1.0/mH; 
187 
invmass[2]=1.0/mH; 
188  
189 
bondsq[0]=dOH*dOH;

190 
bondsq[1]=bondsq[0]; 
191 
bondsq[2]=dHH*dHH;

192 

193 
M2[0]=1.0/(2.0*(invmass[0]+invmass[1])); 
194 
M2[1]=M2[0]; 
195 
M2[2]=1.0/(2.0*(invmass[1]+invmass[2])); 
196  
197 
for(ll=0;ll<3;ll++) { 
198 
l3=3*ll;

199 
ix=3*iatom[ll];

200 
jx=3*jatom[ll];

201 
for(i=0;i<3;i++) 
202 
bond[l3+i]= b4[ix+i]  b4[jx+i]; 
203 
} 
204  
205 
for(i=0,iconv=0;i<1000 && iconv<3; i++) { 
206 
for(ll=0;ll<3;ll++) { 
207 
ii = iatom[ll]; 
208 
jj = jatom[ll]; 
209 
l3 = 3*ll;

210 
ix = 3*ii;

211 
jx = 3*jj;

212 
iy = ix+1;

213 
jy = jx+1;

214 
iz = ix+2;

215 
jz = jx+2;

216  
217 
rijx = bond[l3]; 
218 
rijy = bond[l3+1];

219 
rijz = bond[l3+2];

220  
221 

222 
tx = after[ix]after[jx]; 
223 
ty = after[iy]after[jy]; 
224 
tz = after[iz]after[jz]; 
225 

226 
rp = tx*tx+ty*ty+tz*tz; 
227 
diff = bondsq[ll]  rp; 
228  
229 
if(fabs(diff)<1e8) { 
230 
iconv++; 
231 
} else {

232 
rp = rijx*tx+rijy*ty+rijz*tz; 
233 
if(rp<1e8) { 
234 
return 1; 
235 
} 
236 
acor = diff*M2[ll]/rp; 
237 
im = invmass[ii]; 
238 
jm = invmass[jj]; 
239 
tx = rijx*acor; 
240 
ty = rijy*acor; 
241 
tz = rijz*acor; 
242 
after[ix] += tx*im; 
243 
after[iy] += ty*im; 
244 
after[iz] += tz*im; 
245 
after[jx] = tx*jm; 
246 
after[jy] = ty*jm; 
247 
after[jz] = tz*jm; 
248 
} 
249 
} 
250 
} 
251 
return 0; 
252 
} 
253  
254  
255 
void csettle(FILE *fp,int nsettle, t_iatom iatoms[],real b4[], real after[], 
256 
real dOH,real dHH,real mO,real mH, 
257 
real invdt,real *v,bool bCalcVir,tensor rmdr,int *error) 
258 
{ 
259 
/* ***************************************************************** */

260 
/* ** */

261 
/* Subroutine : setlep  reset positions of TIP3P waters ** */

262 
/* Author : Shuichi Miyamoto ** */

263 
/* Date of last update : Oct. 1, 1992 ** */

264 
/* ** */

265 
/* Reference for the SETTLE algorithm ** */

266 
/* S. Miyamoto et al., J. Comp. Chem., 13, 952 (1992). ** */

267 
/* ** */

268 
/* ***************************************************************** */

269 

270 
/* Initialized data */

271 
static bool bFirst=TRUE; 
272 
/* These three weights need have double precision. Using single precision

273 
* can result in huge velocity and pressure deviations. */

274 
static double wo,wh,wohh; 
275 
static real ra,rb,rc,rc2,rone;

276 
#ifdef DEBUG_PRES

277 
static int step = 0; 
278 
#endif

279 

280 

281 
/* Local variables */

282 
real gama, beta, alpa, xcom, ycom, zcom, al2be2, tmp, tmp2; 
283 
real axlng, aylng, azlng, trns11, trns21, trns31, trns12, trns22, 
284 
trns32, trns13, trns23, trns33, cosphi, costhe, sinphi, sinthe, 
285 
cospsi, xaksxd, yaksxd, xakszd, yakszd, zakszd, zaksxd, xaksyd, 
286 
xb0, yb0, zb0, xc0, yc0, zc0, xa1; 
287 
real ya1, za1, xb1, yb1; 
288 
real zb1, xc1, yc1, zc1, yaksyd, zaksyd, sinpsi, xa3, ya3, za3, 
289 
xb3, yb3, zb3, xc3, yc3, zc3, xb0d, yb0d, xc0d, yc0d, 
290 
za1d, xb1d, yb1d, zb1d, xc1d, yc1d, zc1d, ya2d, xb2d, yb2d, yc2d, 
291 
xa3d, ya3d, za3d, xb3d, yb3d, zb3d, xc3d, yc3d, zc3d; 
292 
real t1,t2; 
293 
real dax, day, daz, dbx, dby, dbz, dcx, dcy, dcz; 
294 
real mdax, mday, mdaz, mdbx, mdby, mdbz, mdcx, mdcy, mdcz; 
295  
296 
int doshake;

297 

298 
int i, shakeret, ow1, hw2, hw3;

299  
300 
*error=1;

301 
if (bFirst) {

302 
if (fp)

303 
fprintf(fp,"Going to use Csettle (%d waters)\n",nsettle);

304 
wo = mO; 
305 
wh = mH; 
306 
wohh = mO+2.0*mH; 
307 
rc = dHH/2.0; 
308 
ra = 2.0*wh*sqrt(dOH*dOHrc*rc)/wohh; 
309 
rb = sqrt(dOH*dOHrc*rc)ra; 
310 
rc2 = dHH; 
311 
rone = 1.0; 
312  
313 
wo /= wohh; 
314 
wh /= wohh; 
315  
316 
if (fp) {

317 
fprintf(fp,"wo = %g, wh =%g, wohh = %g, rc = %g, ra = %g\n",

318 
wo,wh,wohh,rc,ra); 
319 
fprintf(fp,"rb = %g, rc2 = %g, rone = %g, dHH = %g, dOH = %g\n",

320 
rb,rc2,rone,dHH,dOH); 
321 
} 
322  
323 
bFirst = FALSE; 
324 
} 
325 
#ifdef PRAGMAS

326 
#pragma ivdep

327 
#endif

328 
for (i = 0; i < nsettle; ++i) { 
329 
doshake = 0;

330 
/*  Step1 A1'  */

331 
ow1 = iatoms[i*2+1] * 3; 
332 
hw2 = ow1 + 3;

333 
hw3 = ow1 + 6;

334 
xb0 = b4[hw2 ]  b4[ow1]; 
335 
yb0 = b4[hw2 + 1]  b4[ow1 + 1]; 
336 
zb0 = b4[hw2 + 2]  b4[ow1 + 2]; 
337 
xc0 = b4[hw3 ]  b4[ow1]; 
338 
yc0 = b4[hw3 + 1]  b4[ow1 + 1]; 
339 
zc0 = b4[hw3 + 2]  b4[ow1 + 2]; 
340 
/* 6 flops */

341 

342 
xcom = (after[ow1 ] * wo + (after[hw2 ] + after[hw3 ]) * wh); 
343 
ycom = (after[ow1 + 1] * wo + (after[hw2 + 1] + after[hw3 + 1]) * wh); 
344 
zcom = (after[ow1 + 2] * wo + (after[hw2 + 2] + after[hw3 + 2]) * wh); 
345 
/* 12 flops */

346 

347 
xa1 = after[ow1 ]  xcom; 
348 
ya1 = after[ow1 + 1]  ycom;

349 
za1 = after[ow1 + 2]  zcom;

350 
xb1 = after[hw2 ]  xcom; 
351 
yb1 = after[hw2 + 1]  ycom;

352 
zb1 = after[hw2 + 2]  zcom;

353 
xc1 = after[hw3 ]  xcom; 
354 
yc1 = after[hw3 + 1]  ycom;

355 
zc1 = after[hw3 + 2]  zcom;

356 
/* 9 flops */

357 

358 
xakszd = yb0 * zc0  zb0 * yc0; 
359 
yakszd = zb0 * xc0  xb0 * zc0; 
360 
zakszd = xb0 * yc0  yb0 * xc0; 
361 
xaksxd = ya1 * zakszd  za1 * yakszd; 
362 
yaksxd = za1 * xakszd  xa1 * zakszd; 
363 
zaksxd = xa1 * yakszd  ya1 * xakszd; 
364 
xaksyd = yakszd * zaksxd  zakszd * yaksxd; 
365 
yaksyd = zakszd * xaksxd  xakszd * zaksxd; 
366 
zaksyd = xakszd * yaksxd  yakszd * xaksxd; 
367 
/* 27 flops */

368 

369 
axlng = invsqrt(xaksxd * xaksxd + yaksxd * yaksxd + zaksxd * zaksxd); 
370 
aylng = invsqrt(xaksyd * xaksyd + yaksyd * yaksyd + zaksyd * zaksyd); 
371 
azlng = invsqrt(xakszd * xakszd + yakszd * yakszd + zakszd * zakszd); 
372 

373 
trns11 = xaksxd * axlng; 
374 
trns21 = yaksxd * axlng; 
375 
trns31 = zaksxd * axlng; 
376 
trns12 = xaksyd * aylng; 
377 
trns22 = yaksyd * aylng; 
378 
trns32 = zaksyd * aylng; 
379 
trns13 = xakszd * azlng; 
380 
trns23 = yakszd * azlng; 
381 
trns33 = zakszd * azlng; 
382 
/* 24 flops */

383 

384 
xb0d = trns11 * xb0 + trns21 * yb0 + trns31 * zb0; 
385 
yb0d = trns12 * xb0 + trns22 * yb0 + trns32 * zb0; 
386 
xc0d = trns11 * xc0 + trns21 * yc0 + trns31 * zc0; 
387 
yc0d = trns12 * xc0 + trns22 * yc0 + trns32 * zc0; 
388 
/*

389 
xa1d = trns11 * xa1 + trns21 * ya1 + trns31 * za1;

390 
ya1d = trns12 * xa1 + trns22 * ya1 + trns32 * za1;

391 
*/

392 
za1d = trns13 * xa1 + trns23 * ya1 + trns33 * za1; 
393 
xb1d = trns11 * xb1 + trns21 * yb1 + trns31 * zb1; 
394 
yb1d = trns12 * xb1 + trns22 * yb1 + trns32 * zb1; 
395 
zb1d = trns13 * xb1 + trns23 * yb1 + trns33 * zb1; 
396 
xc1d = trns11 * xc1 + trns21 * yc1 + trns31 * zc1; 
397 
yc1d = trns12 * xc1 + trns22 * yc1 + trns32 * zc1; 
398 
zc1d = trns13 * xc1 + trns23 * yc1 + trns33 * zc1; 
399 
/* 65 flops */

400 

401 
sinphi = za1d / ra; 
402 
tmp = rone  sinphi * sinphi; 
403 
if (tmp <= 0) { 
404 
*error = i; 
405 
doshake = 1;

406 
cosphi = 0;

407 
} 
408 
else

409 
cosphi = tmp*invsqrt(tmp); 
410 
sinpsi = (zb1d  zc1d) / (rc2 * cosphi); 
411 
tmp2 = rone  sinpsi * sinpsi; 
412 
if (tmp2 <= 0) { 
413 
*error = i; 
414 
doshake = 1;

415 
cospsi = 0;

416 
} 
417 
else

418 
cospsi = tmp2*invsqrt(tmp2); 
419 
/* 46 flops */

420 

421 
if(!doshake) {

422 
ya2d = ra * cosphi; 
423 
xb2d = rc * cospsi; 
424 
t1 = rb * cosphi; 
425 
t2 = rc * sinpsi * sinphi; 
426 
yb2d = t1  t2; 
427 
yc2d = t1 + t2; 
428 
/* 7 flops */

429 

430 
/*  Step3 al,be,ga  */

431 
alpa = xb2d * (xb0d  xc0d) + yb0d * yb2d + yc0d * yc2d; 
432 
beta = xb2d * (yc0d  yb0d) + xb0d * yb2d + xc0d * yc2d; 
433 
gama = xb0d * yb1d  xb1d * yb0d + xc0d * yc1d  xc1d * yc0d; 
434 
al2be2 = alpa * alpa + beta * beta; 
435 
tmp2 = (al2be2  gama * gama); 
436 
sinthe = (alpa * gama  beta * tmp2*invsqrt(tmp2)) / al2be2; 
437 
/* 47 flops */

438 

439 
/*  Step4 A3'  */

440 
tmp2 = rone  sinthe *sinthe; 
441 
costhe = tmp2*invsqrt(tmp2); 
442 
xa3d = ya2d * sinthe; 
443 
ya3d = ya2d * costhe; 
444 
za3d = za1d; 
445 
xb3d = xb2d * costhe  yb2d * sinthe; 
446 
yb3d = xb2d * sinthe + yb2d * costhe; 
447 
zb3d = zb1d; 
448 
xc3d = xb2d * costhe  yc2d * sinthe; 
449 
yc3d = xb2d * sinthe + yc2d * costhe; 
450 
zc3d = zc1d; 
451 
/* 26 flops */

452 

453 
/*  Step5 A3  */

454 
xa3 = trns11 * xa3d + trns12 * ya3d + trns13 * za3d; 
455 
ya3 = trns21 * xa3d + trns22 * ya3d + trns23 * za3d; 
456 
za3 = trns31 * xa3d + trns32 * ya3d + trns33 * za3d; 
457 
xb3 = trns11 * xb3d + trns12 * yb3d + trns13 * zb3d; 
458 
yb3 = trns21 * xb3d + trns22 * yb3d + trns23 * zb3d; 
459 
zb3 = trns31 * xb3d + trns32 * yb3d + trns33 * zb3d; 
460 
xc3 = trns11 * xc3d + trns12 * yc3d + trns13 * zc3d; 
461 
yc3 = trns21 * xc3d + trns22 * yc3d + trns23 * zc3d; 
462 
zc3 = trns31 * xc3d + trns32 * yc3d + trns33 * zc3d; 
463 
/* 45 flops */

464 
after[ow1] = xcom + xa3; 
465 
after[ow1 + 1] = ycom + ya3;

466 
after[ow1 + 2] = zcom + za3;

467 
after[hw2] = xcom + xb3; 
468 
after[hw2 + 1] = ycom + yb3;

469 
after[hw2 + 2] = zcom + zb3;

470 
after[hw3] = xcom + xc3; 
471 
after[hw3 + 1] = ycom + yc3;

472 
after[hw3 + 2] = zcom + zc3;

473 
/* 9 flops */

474  
475 
dax = xa3  xa1; 
476 
day = ya3  ya1; 
477 
daz = za3  za1; 
478 
dbx = xb3  xb1; 
479 
dby = yb3  yb1; 
480 
dbz = zb3  zb1; 
481 
dcx = xc3  xc1; 
482 
dcy = yc3  yc1; 
483 
dcz = zc3  zc1; 
484 
/* 9 flops, counted with the virial */

485  
486 
if (v) {

487 
v[ow1] += dax*invdt; 
488 
v[ow1 + 1] += day*invdt;

489 
v[ow1 + 2] += daz*invdt;

490 
v[hw2] += dbx*invdt; 
491 
v[hw2 + 1] += dby*invdt;

492 
v[hw2 + 2] += dbz*invdt;

493 
v[hw3] += dcx*invdt; 
494 
v[hw3 + 1] += dcy*invdt;

495 
v[hw3 + 2] += dcz*invdt;

496 
/* 3*6 flops */

497 
} 
498  
499 
if (bCalcVir) {

500 
mdax = mO*dax; 
501 
mday = mO*day; 
502 
mdaz = mO*daz; 
503 
mdbx = mH*dbx; 
504 
mdby = mH*dby; 
505 
mdbz = mH*dbz; 
506 
mdcx = mH*dcx; 
507 
mdcy = mH*dcy; 
508 
mdcz = mH*dcz; 
509 
rmdr[XX][XX] = b4[ow1]*mdax + b4[hw2]*mdbx + b4[hw3]*mdcx; 
510 
rmdr[XX][YY] = b4[ow1]*mday + b4[hw2]*mdby + b4[hw3]*mdcy; 
511 
rmdr[XX][ZZ] = b4[ow1]*mdaz + b4[hw2]*mdbz + b4[hw3]*mdcz; 
512 
rmdr[YY][XX] = b4[ow1+1]*mdax + b4[hw2+1]*mdbx + b4[hw3+1]*mdcx; 
513 
rmdr[YY][YY] = b4[ow1+1]*mday + b4[hw2+1]*mdby + b4[hw3+1]*mdcy; 
514 
rmdr[YY][ZZ] = b4[ow1+1]*mdaz + b4[hw2+1]*mdbz + b4[hw3+1]*mdcz; 
515 
rmdr[ZZ][XX] = b4[ow1+2]*mdax + b4[hw2+2]*mdbx + b4[hw3+2]*mdcx; 
516 
rmdr[ZZ][YY] = b4[ow1+2]*mday + b4[hw2+2]*mdby + b4[hw3+2]*mdcy; 
517 
rmdr[ZZ][ZZ] = b4[ow1+2]*mdaz + b4[hw2+2]*mdbz + b4[hw3+2]*mdcz; 
518 
/* 3*24  9 flops */

519 
} 
520 
} else {

521 
/* If we couldn't settle this water, try a simplified iterative shake instead */

522 
if(xshake(b4+ow1,after+ow1,dOH,dHH,mO,mH)!=0) 
523 
*error=i; 
524 
} 
525 
#ifdef DEBUG

526 
check_cons(fp,"settle",after,ow1,hw2,hw3);

527 
#endif

528 
} 
529 
} 