--- trunk/OOPSE/libmdtools/Integrator.cpp 2003/06/25 21:12:14 567 +++ trunk/OOPSE/libmdtools/Integrator.cpp 2003/07/11 22:34:48 594 @@ -27,6 +27,8 @@ Integrator::Integrator( SimInfo *theInfo, ForceFields* nAtoms = info->n_atoms; + std::cerr << "integ nAtoms = " << nAtoms << "\n"; + // check for constraints constrainedA = NULL; @@ -72,9 +74,14 @@ void Integrator::checkConstraints( void ){ for(int j=0; jis_constrained(); + + std::cerr << "Is the folowing bond constrained \n"; + theArray[j]->printMe(); if(constrained){ + std::cerr << "Yes\n"; + dummy_plug = theArray[j]->get_constraint(); temp_con[nConstrained].set_a( dummy_plug->get_a() ); temp_con[nConstrained].set_b( dummy_plug->get_b() ); @@ -82,7 +89,8 @@ void Integrator::checkConstraints( void ){ nConstrained++; constrained = 0; - } + } + else std::cerr << "No.\n"; } theArray = (SRI**) molecules[i].getMyBends(); @@ -226,6 +234,9 @@ void Integrator::integrate( void ){ calcStress = 1; } + std::cerr << "calcPot = " << calcPot << "; calcStress = " + << calcStress << "\n"; + integrateStep( calcPot, calcStress ); currTime += dt; @@ -257,7 +268,7 @@ void Integrator::integrate( void ){ } - dumpOut->writeFinal(); + dumpOut->writeFinal(currTime); delete dumpOut; delete statOut; @@ -293,9 +304,9 @@ void Integrator::moveA( void ){ double Tb[3]; double ji[3]; double angle; + double A[3][3]; - for( i=0; igetMass() ) * eConvert; + std::cerr<< "MoveA vel[" << i << "] = " + << vel[atomIndex] << "\t" + << vel[atomIndex+1]<< "\t" + << vel[atomIndex+2]<< "\n"; + // position whole step for( j=atomIndex; j<(atomIndex+3); j++ ) pos[j] += dt * vel[j]; + + std::cerr<< "MoveA pos[" << i << "] = " + << pos[atomIndex] << "\t" + << pos[atomIndex+1]<< "\t" + << pos[atomIndex+2]<< "\n"; + if( atoms[i]->isDirectional() ){ dAtom = (DirectionalAtom *)atoms[i]; @@ -328,25 +350,39 @@ void Integrator::moveA( void ){ // use the angular velocities to propagate the rotation matrix a // full time step + // get the atom's rotation matrix + + A[0][0] = dAtom->getAxx(); + A[0][1] = dAtom->getAxy(); + A[0][2] = dAtom->getAxz(); + + A[1][0] = dAtom->getAyx(); + A[1][1] = dAtom->getAyy(); + A[1][2] = dAtom->getAyz(); + + A[2][0] = dAtom->getAzx(); + A[2][1] = dAtom->getAzy(); + A[2][2] = dAtom->getAzz(); + // rotate about the x-axis angle = dt2 * ji[0] / dAtom->getIxx(); - this->rotate( 1, 2, angle, ji, &Amat[aMatIndex] ); + this->rotate( 1, 2, angle, ji, A ); // rotate about the y-axis angle = dt2 * ji[1] / dAtom->getIyy(); - this->rotate( 2, 0, angle, ji, &Amat[aMatIndex] ); + this->rotate( 2, 0, angle, ji, A ); // rotate about the z-axis angle = dt * ji[2] / dAtom->getIzz(); - this->rotate( 0, 1, angle, ji, &Amat[aMatIndex] ); + this->rotate( 0, 1, angle, ji, A ); // rotate about the y-axis angle = dt2 * ji[1] / dAtom->getIyy(); - this->rotate( 2, 0, angle, ji, &Amat[aMatIndex] ); + this->rotate( 2, 0, angle, ji, A ); // rotate about the x-axis angle = dt2 * ji[0] / dAtom->getIxx(); - this->rotate( 1, 2, angle, ji, &Amat[aMatIndex] ); + this->rotate( 1, 2, angle, ji, A ); dAtom->setJx( ji[0] ); dAtom->setJy( ji[1] ); @@ -371,6 +407,12 @@ void Integrator::moveB( void ){ for( j=atomIndex; j<(atomIndex+3); j++ ) vel[j] += ( dt2 * frc[j] / atoms[i]->getMass() ) * eConvert; + std::cerr<< "MoveB vel[" << i << "] = " + << vel[atomIndex] << "\t" + << vel[atomIndex+1]<< "\t" + << vel[atomIndex+2]<< "\n"; + + if( atoms[i]->isDirectional() ){ dAtom = (DirectionalAtom *)atoms[i]; @@ -411,8 +453,8 @@ void Integrator::constrainA(){ int i,j,k; int done; - double pxab, pyab, pzab; - double rxab, ryab, rzab; + double pab[3]; + double rab[3]; int a, b, ax, ay, az, bx, by, bz; double rma, rmb; double dx, dy, dz; @@ -422,8 +464,6 @@ void Integrator::constrainA(){ double gab; int iteration; - - for( i=0; ibox_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; + info->wrapVector( pab ); + pabsq = pab[0] * pab[0] + pab[1] * pab[1] + pab[2] * pab[2]; + rabsq = constrainedDsqr[i]; diffsq = rabsq - pabsq; // the original rattle code from alan tidesley if (fabs(diffsq) > (tol*rabsq*2)) { - rxab = oldPos[ax] - oldPos[bx]; - ryab = oldPos[ay] - oldPos[by]; - rzab = oldPos[az] - oldPos[bz]; + rab[0] = oldPos[ax] - oldPos[bx]; + rab[1] = oldPos[ay] - oldPos[by]; + rab[2] = oldPos[az] - oldPos[bz]; - 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); + info->wrapVector( rab ); - rpab = rxab * pxab + ryab * pyab + rzab * pzab; + rpab = rab[0] * pab[0] + rab[1] * pab[1] + rab[2] * pab[2]; rpabsq = rpab * rpab; @@ -503,9 +534,9 @@ void Integrator::constrainA(){ 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[ax] += rma * dx; pos[ay] += rma * dy; @@ -559,7 +590,7 @@ void Integrator::constrainB( void ){ int i,j,k; int done; double vxab, vyab, vzab; - double rxab, ryab, rzab; + double rab[3]; int a, b, ax, ay, az, bx, by, bz; double rma, rmb; double dx, dy, dz; @@ -598,30 +629,24 @@ void Integrator::constrainB( void ){ vyab = vel[ay] - vel[by]; vzab = vel[az] - vel[bz]; - rxab = pos[ax] - pos[bx]; - ryab = pos[ay] - pos[by]; - rzab = pos[az] - pos[bz]; + rab[0] = pos[ax] - pos[bx]; + rab[1] = pos[ay] - pos[by]; + rab[2] = pos[az] - pos[bz]; - - 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); + 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[ax] += rma * dx; vel[ay] += rma * dy; @@ -665,7 +690,7 @@ void Integrator::rotate( int axes1, int axes2, double void Integrator::rotate( int axes1, int axes2, double angle, double ji[3], - double A[9] ){ + double A[3][3] ){ int i,j,k; double sinAngle; @@ -681,7 +706,7 @@ void Integrator::rotate( int axes1, int axes2, double for(i=0; i<3; i++){ for(j=0; j<3; j++){ - tempA[j][i] = A[3*i + j]; + tempA[j][i] = A[i][j]; } } @@ -738,9 +763,9 @@ void Integrator::rotate( int axes1, int axes2, double for(i=0; i<3; i++){ for(j=0; j<3; j++){ - A[3*j + i] = 0.0; + A[j][i] = 0.0; for(k=0; k<3; k++){ - A[3*j + i] += tempA[i][k] * rot[j][k]; + A[j][i] += tempA[i][k] * rot[j][k]; } } }