--- trunk/OOPSE/libmdtools/NPTf.cpp 2003/07/10 22:15:53 590 +++ trunk/OOPSE/libmdtools/NPTf.cpp 2003/07/14 22:38:13 600 @@ -37,15 +37,18 @@ void NPTf::moveA() { void NPTf::moveA() { - int i,j,k; - int atomIndex, aMatIndex; + int i, j, k; DirectionalAtom* dAtom; - double Tb[3]; - double ji[3]; - double ri[3], vi[3], sc[3]; - double instaTemp, instaVol; - double tt2, tb2, eta2ij; - double angle; + double Tb[3], ji[3]; + double A[3][3], I[3][3]; + double angle, mass; + double vel[3], pos[3], frc[3]; + + double rj[3]; + double instaTemp, instaPress, instaVol; + double tt2, tb2; + double sc[3]; + double eta2ij; double press[3][3], vScale[3][3], hm[3][3], hmnew[3][3], scaleMat[3][3]; tt2 = tauThermostat * tauThermostat; @@ -62,10 +65,10 @@ void NPTf::moveA() { for (i = 0; i < 3; i++ ) { for (j = 0; j < 3; j++ ) { if (i == j) { - + eta[i][j] += dt2 * instaVol * (press[i][j] - targetPressure/p_convert) / (NkBT*tb2); - + vScale[i][j] = eta[i][j] + chi; } else { @@ -79,38 +82,32 @@ void NPTf::moveA() { } for( i=0; igetVel( vel ); + atoms[i]->getPos( pos ); + atoms[i]->getFrc( frc ); + + mass = atoms[i]->getMass(); // velocity half step + + info->matVecMul3( vScale, vel, sc ); - vi[0] = vel[atomIndex]; - vi[1] = vel[atomIndex+1]; - vi[2] = vel[atomIndex+2]; - - info->matVecMul3( vScale, vi, sc ); - - vi[0] += dt2 * ((frc[atomIndex] /atoms[i]->getMass())*eConvert - sc[0]); - vi[1] += dt2 * ((frc[atomIndex+1]/atoms[i]->getMass())*eConvert - sc[1]); - vi[2] += dt2 * ((frc[atomIndex+2]/atoms[i]->getMass())*eConvert - sc[2]); + for (j = 0; j < 3; j++) { + vel[j] += dt2 * ((frc[j] / mass) * eConvert - sc[j]); + rj[j] = pos[j]; + } - vel[atomIndex] = vi[0]; - vel[atomIndex+1] = vi[1]; - vel[atomIndex+2] = vi[2]; + atoms[i]->setVel( vel ); // position whole step - ri[0] = pos[atomIndex]; - ri[1] = pos[atomIndex+1]; - ri[2] = pos[atomIndex+2]; + info->wrapVector(rj); - info->wrapVector(ri); + info->matVecMul3( eta, rj, sc ); - info->matVecMul3( eta, ri, sc ); - - pos[atomIndex] += dt * (vel[atomIndex] + sc[0]); - pos[atomIndex+1] += dt * (vel[atomIndex+1] + sc[1]); - pos[atomIndex+2] += dt * (vel[atomIndex+2] + sc[2]); + for (j = 0; j < 3; j++ ) + pos[j] += dt * (vel[j] + sc[j]); if( atoms[i]->isDirectional() ){ @@ -118,64 +115,59 @@ void NPTf::moveA() { // get and convert the torque to body frame - Tb[0] = dAtom->getTx(); - Tb[1] = dAtom->getTy(); - Tb[2] = dAtom->getTz(); - + dAtom->getTrq( Tb ); dAtom->lab2Body( Tb ); // get the angular momentum, and propagate a half step - ji[0] = dAtom->getJx(); - ji[1] = dAtom->getJy(); - ji[2] = dAtom->getJz(); + dAtom->getJ( ji ); + + for (j=0; j < 3; j++) + ji[j] += dt2 * (Tb[j] * eConvert - ji[j]*chi); - ji[0] += dt2 * (Tb[0] * eConvert - ji[0]*chi); - ji[1] += dt2 * (Tb[1] * eConvert - ji[1]*chi); - ji[2] += dt2 * (Tb[2] * eConvert - ji[2]*chi); - // use the angular velocities to propagate the rotation matrix a // full time step - + + dAtom->getA(A); + dAtom->getI(I); + // rotate about the x-axis - angle = dt2 * ji[0] / dAtom->getIxx(); - this->rotate( 1, 2, angle, ji, &Amat[aMatIndex] ); - + angle = dt2 * ji[0] / I[0][0]; + 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] ); + angle = dt2 * ji[1] / I[1][1]; + 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] ); + angle = dt * ji[2] / I[2][2]; + 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] ); + angle = dt2 * ji[1] / I[1][1]; + 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] ); + angle = dt2 * ji[0] / I[0][0]; + this->rotate( 1, 2, angle, ji, A ); - dAtom->setJx( ji[0] ); - dAtom->setJy( ji[1] ); - dAtom->setJz( ji[2] ); - } - + dAtom->setJ( ji ); + dAtom->setA( A ); + } } - + // Scale the box after all the positions have been moved: - + // Use a taylor expansion for eta products: Hmat = Hmat . exp(dt * etaMat) // Hmat = Hmat . ( Ident + dt * etaMat + dt^2 * etaMat*etaMat / 2) - - + + for(i=0; i<3; i++){ for(j=0; j<3; j++){ - + // Calculate the matrix Product of the eta array (we only need // the ij element right now): - + eta2ij = 0.0; for(k=0; k<3; k++){ eta2ij += eta[i][k] * eta[k][j]; @@ -186,10 +178,10 @@ void NPTf::moveA() { if (i == j) scaleMat[i][j] = 1.0; // Taylor expansion for the exponential truncated at second order: scaleMat[i][j] += dt*eta[i][j] + 0.5*dt*dt*eta2ij; - + } } - + info->getBoxM(hm); info->matMul3(hm, scaleMat, hmnew); info->setBoxM(hmnew); @@ -197,14 +189,16 @@ void NPTf::moveB( void ){ } void NPTf::moveB( void ){ - int i,j, k; - int atomIndex; + + int i, j; DirectionalAtom* dAtom; - double Tb[3]; - double ji[3]; - double vi[3], sc[3]; - double instaTemp, instaVol; + double Tb[3], ji[3]; + double vel[3], frc[3]; + double mass; + + double instaTemp, instaPress, instaVol; double tt2, tb2; + double sc[3]; double press[3][3], vScale[3][3]; tt2 = tauThermostat * tauThermostat; @@ -238,51 +232,41 @@ void NPTf::moveB( void ){ } for( i=0; igetVel( vel ); + atoms[i]->getFrc( frc ); + + mass = atoms[i]->getMass(); + // velocity half step + + info->matVecMul3( vScale, vel, sc ); - vi[0] = vel[atomIndex]; - vi[1] = vel[atomIndex+1]; - vi[2] = vel[atomIndex+2]; - - info->matVecMul3( vScale, vi, sc ); - - vi[0] += dt2 * ((frc[atomIndex] /atoms[i]->getMass())*eConvert - sc[0]); - vi[1] += dt2 * ((frc[atomIndex+1]/atoms[i]->getMass())*eConvert - sc[1]); - vi[2] += dt2 * ((frc[atomIndex+2]/atoms[i]->getMass())*eConvert - sc[2]); + for (j = 0; j < 3; j++) { + vel[j] += dt2 * ((frc[j] / mass) * eConvert - sc[j]); + } - vel[atomIndex] = vi[0]; - vel[atomIndex+1] = vi[1]; - vel[atomIndex+2] = vi[2]; + atoms[i]->setVel( vel ); if( atoms[i]->isDirectional() ){ - + dAtom = (DirectionalAtom *)atoms[i]; - + // get and convert the torque to body frame - Tb[0] = dAtom->getTx(); - Tb[1] = dAtom->getTy(); - Tb[2] = dAtom->getTz(); - + dAtom->getTrq( Tb ); dAtom->lab2Body( Tb ); - // get the angular momentum, and complete the angular momentum - // half step + // get the angular momentum, and propagate a half step - ji[0] = dAtom->getJx(); - ji[1] = dAtom->getJy(); - ji[2] = dAtom->getJz(); + dAtom->getJ( ji ); - ji[0] += dt2 * (Tb[0] * eConvert - ji[0]*chi); - ji[1] += dt2 * (Tb[1] * eConvert - ji[1]*chi); - ji[2] += dt2 * (Tb[2] * eConvert - ji[2]*chi); + for (j=0; j < 3; j++) + ji[j] += dt2 * (Tb[j] * eConvert - ji[j]*chi); - dAtom->setJx( ji[0] ); - dAtom->setJy( ji[1] ); - dAtom->setJz( ji[2] ); - } + dAtom->setJ( ji ); + + } } }