Project

General

Profile

csettle.c

fixed src/mdlib/csettle.c - Berk Hess, 05/28/2009 11:11 AM

 
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) 1991-2000, University of Groningen, The Netherlands.
13
 * Copyright (c) 2001-2004, 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 2008-1-10
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)<1e-8) {
230
        iconv++;
231
      } else {
232
        rp = rijx*tx+rijy*ty+rijz*tz;
233
        if(rp<1e-8) {
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 C-settle (%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*dOH-rc*rc)/wohh;
309
    rb     = sqrt(dOH*dOH-rc*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
}