--- trunk/OOPSE/libmdtools/Integrator.cpp 2003/06/20 20:29:36 561 +++ trunk/OOPSE/libmdtools/Integrator.cpp 2003/06/23 21:24:03 563 @@ -137,6 +137,9 @@ void Integrator::checkConstraints( void ){ constrainedA[i] = temp_con[i].get_a(); constrainedB[i] = temp_con[i].get_b(); constrainedDsqr[i] = temp_con[i].get_dsqr(); + + cerr << "constraint " << constrainedA[i] << " <-> " << constrainedB[i] + << " => " << constrainedDsqr[i] << "\n"; } @@ -402,22 +405,6 @@ void Integrator::preMove( void ){ if( nConstrained ){ -// if( oldAtoms != nAtoms ){ - -// // save oldAtoms to check for lode balanceing later on. - -// oldAtoms = nAtoms; - -// delete[] moving; -// delete[] moved; -// delete[] oldPos; - -// moving = new int[nAtoms]; -// moved = new int[nAtoms]; - -// oldPos = new double[nAtoms*3]; -// } - for(i=0; i<(nAtoms*3); i++) oldPos[i] = pos[i]; } } @@ -428,7 +415,7 @@ void Integrator::constrainA(){ int done; double pxab, pyab, pzab; double rxab, ryab, rzab; - int a, b; + int a, b, ax, ay, az, bx, by, bz; double rma, rmb; double dx, dy, dz; double rpab; @@ -455,49 +442,62 @@ void Integrator::constrainA(){ a = constrainedA[i]; b = constrainedB[i]; - + + ax = (a*3) + 0; + ay = (a*3) + 1; + az = (a*3) + 2; + + bx = (b*3) + 0; + by = (b*3) + 1; + bz = (b*3) + 2; + + if( moved[a] || moved[b] ){ - pxab = pos[3*a+0] - pos[3*b+0]; - pyab = pos[3*a+1] - pos[3*b+1]; - pzab = pos[3*a+2] - pos[3*b+2]; + pxab = pos[ax] - pos[bx]; + pyab = pos[ay] - pos[by]; + pzab = pos[az] - pos[bz]; //periodic boundary condition pxab = pxab - info->box_x * copysign(1, pxab) - * int( fabs(pxab) / info->box_x + 0.5); + * (int)( fabs(pxab / info->box_x) + 0.5); pyab = pyab - info->box_y * copysign(1, pyab) - * int( fabs(pyab) / info->box_y + 0.5); + * (int)( fabs(pyab / info->box_y) + 0.5); pzab = pzab - info->box_z * copysign(1, pzab) - * int( fabs(pzab) / info->box_z + 0.5); + * (int)( fabs(pzab / info->box_z) + 0.5); pabsq = pxab * pxab + pyab * pyab + pzab * pzab; rabsq = constrainedDsqr[i]; diffsq = pabsq - rabsq; // the original rattle code from alan tidesley - if (fabs(diffsq) > tol*rabsq*2) { - rxab = oldPos[3*a+0] - oldPos[3*b+0]; - ryab = oldPos[3*a+1] - oldPos[3*b+1]; - rzab = oldPos[3*a+2] - oldPos[3*b+2]; + if (fabs(diffsq) > (tol*rabsq*2)) { + rxab = oldPos[ax] - oldPos[bx]; + ryab = oldPos[ay] - oldPos[by]; + rzab = oldPos[az] - oldPos[bz]; rxab = rxab - info->box_x * copysign(1, rxab) - * int( fabs(rxab) / info->box_x + 0.5); + * (int)( fabs(rxab / info->box_x) + 0.5); ryab = ryab - info->box_y * copysign(1, ryab) - * int( fabs(ryab) / info->box_y + 0.5); + * (int)( fabs(ryab / info->box_y) + 0.5); rzab = rzab - info->box_z * copysign(1, rzab) - * int( fabs(rzab) / info->box_z + 0.5); + * (int)( fabs(rzab / info->box_z) + 0.5); rpab = rxab * pxab + ryab * pyab + rzab * pzab; rpabsq = rpab * rpab; if (rpabsq < (rabsq * -diffsq)){ + + cerr << "rpabsq = " << rpabsq << ", rabsq = " << rabsq + << ", -diffsq = " << -diffsq << "\n"; + #ifdef IS_MPI a = atoms[a]->getGlobalIndex(); b = atoms[b]->getGlobalIndex(); #endif //is_mpi sprintf( painCave.errMsg, - "Constraint failure in constrainA at atom %d and %d\n.", + "Constraint failure in constrainA at atom %d and %d.\n", a, b ); painCave.isFatal = 1; simError(); @@ -511,25 +511,25 @@ void Integrator::constrainA(){ dy = ryab * gab; dz = rzab * gab; - pos[3*a+0] += rma * dx; - pos[3*a+1] += rma * dy; - pos[3*a+2] += rma * dz; + pos[ax] += rma * dx; + pos[ay] += rma * dy; + pos[az] += rma * dz; - pos[3*b+0] -= rmb * dx; - pos[3*b+1] -= rmb * dy; - pos[3*b+2] -= rmb * dz; + pos[bx] -= rmb * dx; + pos[by] -= rmb * dy; + pos[bz] -= rmb * dz; dx = dx / dt; dy = dy / dt; dz = dz / dt; - vel[3*a+0] += rma * dx; - vel[3*a+1] += rma * dy; - vel[3*a+2] += rma * dz; + vel[ax] += rma * dx; + vel[ay] += rma * dy; + vel[az] += rma * dz; - vel[3*b+0] -= rmb * dx; - vel[3*b+1] -= rmb * dy; - vel[3*b+2] -= rmb * dz; + vel[bx] -= rmb * dx; + vel[by] -= rmb * dy; + vel[bz] -= rmb * dz; moving[a] = 1; moving[b] = 1; @@ -545,6 +545,7 @@ void Integrator::constrainA(){ } iteration++; + cerr << "iterainA = " << iteration << "\n"; } if( !done ){ @@ -564,7 +565,7 @@ void Integrator::constrainB( void ){ int done; double vxab, vyab, vzab; double rxab, ryab, rzab; - int a, b; + int a, b, ax, ay, az, bx, by, bz; double rma, rmb; double dx, dy, dz; double rabsq, pabsq, rvab; @@ -586,22 +587,30 @@ void Integrator::constrainB( void ){ a = constrainedA[i]; b = constrainedB[i]; + ax = 3*a +0; + ay = 3*a +1; + az = 3*a +2; + + bx = 3*b +0; + by = 3*b +1; + bz = 3*b +2; + if( moved[a] || moved[b] ){ - vxab = vel[3*a+0] - vel[3*b+0]; - vyab = vel[3*a+1] - vel[3*b+1]; - vzab = vel[3*a+2] - vel[3*b+2]; + vxab = vel[ax] - vel[bx]; + vyab = vel[ay] - vel[by]; + vzab = vel[az] - vel[bz]; - rxab = pos[3*a+0] - pos[3*b+0]; - ryab = pos[3*a+1] - pos[3*b+1]; - rzab = pos[3*a+2] - pos[3*b+2]; + rxab = pos[ax] - pos[bx]; + ryab = pos[ay] - pos[by]; + rzab = pos[az] - pos[bz]; rxab = rxab - info->box_x * copysign(1, rxab) - * int( fabs(rxab) / info->box_x + 0.5); + * (int)( fabs(rxab / info->box_x) + 0.5); ryab = ryab - info->box_y * copysign(1, ryab) - * int( fabs(ryab) / info->box_y + 0.5); + * (int)( fabs(ryab / info->box_y) + 0.5); rzab = rzab - info->box_z * copysign(1, rzab) - * int( fabs(rzab) / info->box_z + 0.5); + * (int)( fabs(rzab / info->box_z) + 0.5); rma = 1.0 / atoms[a]->getMass(); rmb = 1.0 / atoms[b]->getMass(); @@ -616,13 +625,13 @@ void Integrator::constrainB( void ){ dy = ryab * gab; dz = rzab * gab; - vel[3*a+0] += rma * dx; - vel[3*a+1] += rma * dy; - vel[3*a+2] += rma * dz; + vel[ax] += rma * dx; + vel[ay] += rma * dy; + vel[az] += rma * dz; - vel[3*b+0] -= rmb * dx; - vel[3*b+1] -= rmb * dy; - vel[3*b+2] -= rmb * dz; + vel[bx] -= rmb * dx; + vel[by] -= rmb * dy; + vel[bz] -= rmb * dz; moving[a] = 1; moving[b] = 1;