--- trunk/OOPSE/libmdtools/Integrator.cpp 2003/06/25 21:12:14 567 +++ trunk/OOPSE/libmdtools/Integrator.cpp 2003/07/14 21:28:54 597 @@ -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(); @@ -216,8 +224,6 @@ void Integrator::integrate( void ){ pos = Atom::getPosArray(); vel = Atom::getVelArray(); frc = Atom::getFrcArray(); - trq = Atom::getTrqArray(); - Amat = Atom::getAmatArray(); while( currTime < runTime ){ @@ -226,6 +232,8 @@ void Integrator::integrate( void ){ calcStress = 1; } + std::cerr << currTime << "\n"; + integrateStep( calcPot, calcStress ); currTime += dt; @@ -257,7 +265,7 @@ void Integrator::integrate( void ){ } - dumpOut->writeFinal(); + dumpOut->writeFinal(currTime); delete dumpOut; delete statOut; @@ -271,7 +279,7 @@ void Integrator::integrateStep( int calcPot, int calcS preMove(); moveA(); - if( nConstrained ) constrainA(); + //if( nConstrained ) constrainA(); // calc forces @@ -293,9 +301,9 @@ void Integrator::moveA( void ){ double Tb[3]; double ji[3]; double angle; + double A[3][3], At[3][3]; - for( i=0; igetMass() ) * eConvert; + // position whole step for( j=atomIndex; j<(atomIndex+3); j++ ) pos[j] += dt * vel[j]; + if( atoms[i]->isDirectional() ){ dAtom = (DirectionalAtom *)atoms[i]; @@ -316,9 +326,9 @@ void Integrator::moveA( void ){ Tb[0] = dAtom->getTx(); Tb[1] = dAtom->getTy(); Tb[2] = dAtom->getTz(); - + dAtom->lab2Body( Tb ); - + // get the angular momentum, and propagate a half step ji[0] = dAtom->getJx() + ( dt2 * Tb[0] ) * eConvert; @@ -331,7 +341,7 @@ void Integrator::moveA( void ){ // rotate about the x-axis angle = dt2 * ji[0] / dAtom->getIxx(); this->rotate( 1, 2, angle, ji, &Amat[aMatIndex] ); - + // rotate about the y-axis angle = dt2 * ji[1] / dAtom->getIyy(); this->rotate( 2, 0, angle, ji, &Amat[aMatIndex] ); @@ -351,6 +361,15 @@ void Integrator::moveA( void ){ dAtom->setJx( ji[0] ); dAtom->setJy( ji[1] ); dAtom->setJz( ji[2] ); + + std::cerr << "Amat[" << i << "]\n"; + info->printMat9( &Amat[aMatIndex] ); + + std::cerr << "ji[" << i << "]\t" + << ji[0] << "\t" + << ji[1] << "\t" + << ji[2] << "\n"; + } } @@ -359,18 +378,20 @@ void Integrator::moveB( void ){ void Integrator::moveB( void ){ int i,j,k; - int atomIndex; + int atomIndex, aMatIndex; DirectionalAtom* dAtom; double Tb[3]; double ji[3]; for( i=0; igetMass() ) * eConvert; + if( atoms[i]->isDirectional() ){ dAtom = (DirectionalAtom *)atoms[i]; @@ -381,6 +402,11 @@ void Integrator::moveB( void ){ Tb[1] = dAtom->getTy(); Tb[2] = dAtom->getTz(); + std::cerr << "TrqB[" << i << "]\t" + << Tb[0] << "\t" + << Tb[1] << "\t" + << Tb[2] << "\n"; + dAtom->lab2Body( Tb ); // get the angular momentum, and complete the angular momentum @@ -393,6 +419,15 @@ void Integrator::moveB( void ){ dAtom->setJx( ji[0] ); dAtom->setJy( ji[1] ); dAtom->setJz( ji[2] ); + + + std::cerr << "Amat[" << i << "]\n"; + info->printMat9( &Amat[aMatIndex] ); + + std::cerr << "ji[" << i << "]\t" + << ji[0] << "\t" + << ji[1] << "\t" + << ji[2] << "\n"; } } @@ -411,8 +446,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 +457,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 +527,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 +583,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 +622,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; @@ -677,11 +695,12 @@ void Integrator::rotate( int axes1, int axes2, double double tempA[3][3]; double tempJ[3]; + // initialize the tempA for(i=0; i<3; i++){ for(j=0; j<3; j++){ - tempA[j][i] = A[3*i + j]; + tempA[j][i] = A[3*i+j]; } } @@ -738,9 +757,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[3*j+i] = 0.0; for(k=0; k<3; k++){ - A[3*j + i] += tempA[i][k] * rot[j][k]; + A[3*j+i] += tempA[i][k] * rot[j][k]; } } }