--- trunk/OOPSE/libmdtools/Integrator.cpp 2003/06/20 20:29:36 561 +++ trunk/OOPSE/libmdtools/Integrator.cpp 2003/07/02 21:26:55 572 @@ -137,6 +137,7 @@ 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(); + } @@ -256,7 +257,7 @@ void Integrator::integrate( void ){ } - dumpOut->writeFinal(); + dumpOut->writeFinal(currTime); delete dumpOut; delete statOut; @@ -293,19 +294,19 @@ void Integrator::moveA( void ){ double ji[3]; double angle; + + for( i=0; igetMass() ) * eConvert; // position whole step - for( j=atomIndex; j<(atomIndex+3); j++ ) - pos[j] += dt * vel[j]; - - + for( j=atomIndex; j<(atomIndex+3); j++ ) pos[j] += dt * vel[j]; + if( atoms[i]->isDirectional() ){ dAtom = (DirectionalAtom *)atoms[i]; @@ -402,22 +403,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]; } } @@ -426,9 +411,9 @@ void Integrator::constrainA(){ int i,j,k; int done; - double pxab, pyab, pzab; - double rxab, ryab, rzab; - int a, b; + double pab[3]; + double rab[3]; + int a, b, ax, ay, az, bx, by, bz; double rma, rmb; double dx, dy, dz; double rpab; @@ -444,8 +429,7 @@ void Integrator::constrainA(){ moving[i] = 0; moved[i] = 1; } - - + iteration = 0; done = 0; while( !done && (iteration < maxIteration )){ @@ -455,49 +439,51 @@ 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]; + pab[0] = pos[ax] - pos[bx]; + pab[1] = pos[ay] - pos[by]; + pab[2] = pos[az] - pos[bz]; - //periodic boundary condition - pxab = pxab - info->box_x * copysign(1, pxab) - * int( fabs(pxab) / info->box_x + 0.5); - pyab = pyab - info->box_y * copysign(1, pyab) - * int( fabs(pyab) / info->box_y + 0.5); - pzab = pzab - info->box_z * copysign(1, pzab) - * int( fabs(pzab) / info->box_z + 0.5); - - pabsq = pxab * pxab + pyab * pyab + pzab * pzab; + //periodic boundary condition + + info->wrapVector( pab ); + + pabsq = pab[0] * pab[0] + pab[1] * pab[1] + pab[2] * pab[2]; + rabsq = constrainedDsqr[i]; - diffsq = pabsq - rabsq; + diffsq = rabsq - pabsq; // 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]; - - rxab = rxab - info->box_x * copysign(1, rxab) - * int( fabs(rxab) / info->box_x + 0.5); - ryab = ryab - info->box_y * copysign(1, ryab) - * int( fabs(ryab) / info->box_y + 0.5); - rzab = rzab - info->box_z * copysign(1, rzab) - * int( fabs(rzab) / info->box_z + 0.5); + if (fabs(diffsq) > (tol*rabsq*2)) { + rab[0] = oldPos[ax] - oldPos[bx]; + rab[1] = oldPos[ay] - oldPos[by]; + rab[2] = oldPos[az] - oldPos[bz]; - rpab = rxab * pxab + ryab * pyab + rzab * pzab; + info->wrapVector( rab ); + + rpab = rab[0] * pab[0] + rab[1] * pab[1] + rab[2] * pab[2]; + rpabsq = rpab * rpab; if (rpabsq < (rabsq * -diffsq)){ + #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(); @@ -505,31 +491,32 @@ void Integrator::constrainA(){ rma = 1.0 / atoms[a]->getMass(); rmb = 1.0 / atoms[b]->getMass(); - + gab = diffsq / ( 2.0 * ( rma + rmb ) * rpab ); - dx = rxab * gab; - dy = ryab * gab; - dz = rzab * gab; + + dx = rab[0] * gab; + dy = rab[1] * gab; + dz = rab[2] * 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; @@ -563,8 +550,8 @@ void Integrator::constrainB( void ){ int i,j,k; int done; double vxab, vyab, vzab; - double rxab, ryab, rzab; - int a, b; + double rab[3]; + int a, b, ax, ay, az, bx, by, bz; double rma, rmb; double dx, dy, dz; double rabsq, pabsq, rvab; @@ -581,48 +568,53 @@ void Integrator::constrainB( void ){ iteration = 0; while( !done && (iteration < maxIteration ) ){ + done = 1; + for(i=0; ibox_x * copysign(1, rxab) - * int( fabs(rxab) / info->box_x + 0.5); - ryab = ryab - info->box_y * copysign(1, ryab) - * int( fabs(ryab) / info->box_y + 0.5); - rzab = rzab - info->box_z * copysign(1, rzab) - * int( fabs(rzab) / info->box_z + 0.5); - + info->wrapVector( rab ); + rma = 1.0 / atoms[a]->getMass(); rmb = 1.0 / atoms[b]->getMass(); - rvab = rxab * vxab + ryab * vyab + rzab * vzab; + rvab = rab[0] * vxab + rab[1] * vyab + rab[2] * vzab; gab = -rvab / ( ( rma + rmb ) * constrainedDsqr[i] ); if (fabs(gab) > tol) { - dx = rxab * gab; - dy = ryab * gab; - dz = rzab * gab; + dx = rab[0] * gab; + dy = rab[1] * gab; + dz = rab[2] * 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;