--- trunk/OOPSE/libmdtools/NPTf.cpp 2003/07/10 17:10:56 588 +++ trunk/OOPSE/libmdtools/NPTf.cpp 2003/07/15 19:56:08 617 @@ -1,3 +1,4 @@ +#include #include "Atom.hpp" #include "SRI.hpp" #include "AbstractClasses.hpp" @@ -26,7 +27,7 @@ NPTf::NPTf ( SimInfo *theInfo, ForceFields* the_ff): chi = 0.0; for(i = 0; i < 3; i++) - for (j = 0; j < 3; j_++) + for (j = 0; j < 3; j++) eta[i][j] = 0.0; have_tau_thermostat = 0; @@ -37,16 +38,20 @@ 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]; + double bigScale, smallScale, offDiagMax; tt2 = tauThermostat * tauThermostat; tb2 = tauBarostat * tauBarostat; @@ -62,10 +67,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 +84,34 @@ 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 ); + for (j = 0; j < 3; j++ ) + pos[j] += dt * (vel[j] + sc[j]); - 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]); + atoms[i]->setPos( pos ); if( atoms[i]->isDirectional() ){ @@ -118,64 +119,62 @@ 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) - - + + bigScale = 1.0; + smallScale = 1.0; + offDiagMax = 0.0; + 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]; @@ -187,24 +186,59 @@ void NPTf::moveA() { // Taylor expansion for the exponential truncated at second order: scaleMat[i][j] += dt*eta[i][j] + 0.5*dt*dt*eta2ij; + if (i != j) + if (fabs(scaleMat[i][j]) > offDiagMax) + offDiagMax = fabs(scaleMat[i][j]); + } + + if (scaleMat[i][i] > bigScale) bigScale = scaleMat[i][i]; + if (scaleMat[i][i] < smallScale) smallScale = scaleMat[i][i]; } - - info->getBoxM(hm); - info->matMul3(hm, scaleMat, hmnew); - info->setBoxM(hmnew); + if ((bigScale > 1.1) || (smallScale < 0.9)) { + sprintf( painCave.errMsg, + "NPTf error: Attempting a Box scaling of more than 10 percent.\n" + " Check your tauBarostat, as it is probably too small!\n\n" + " scaleMat = [%lf\t%lf\t%lf]\n" + " [%lf\t%lf\t%lf]\n" + " [%lf\t%lf\t%lf]\n", + scaleMat[0][0],scaleMat[0][1],scaleMat[0][2], + scaleMat[1][0],scaleMat[1][1],scaleMat[1][2], + scaleMat[2][0],scaleMat[2][1],scaleMat[2][2]); + painCave.isFatal = 1; + simError(); + } else if (offDiagMax > 0.1) { + sprintf( painCave.errMsg, + "NPTf error: Attempting an off-diagonal Box scaling of more than 10 percent.\n" + " Check your tauBarostat, as it is probably too small!\n\n" + " scaleMat = [%lf\t%lf\t%lf]\n" + " [%lf\t%lf\t%lf]\n" + " [%lf\t%lf\t%lf]\n", + scaleMat[0][0],scaleMat[0][1],scaleMat[0][2], + scaleMat[1][0],scaleMat[1][1],scaleMat[1][2], + scaleMat[2][0],scaleMat[2][1],scaleMat[2][2]); + painCave.isFatal = 1; + simError(); + } else { + info->getBoxM(hm); + info->matMul3(hm, scaleMat, hmnew); + info->setBoxM(hmnew); + } + } 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 +272,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 ); + + } } }