--- trunk/OOPSE/libmdtools/NPTi.cpp 2003/09/15 16:52:02 763 +++ trunk/OOPSE/libmdtools/NPTi.cpp 2003/09/16 20:02:11 767 @@ -56,119 +56,6 @@ template void NPTi::moveA() { } template void NPTi::moveA() { - - -// int i, j; -// DirectionalAtom* dAtom; -// 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, scaleFactor; - -// tt2 = tauThermostat * tauThermostat; -// tb2 = tauBarostat * tauBarostat; - -// instaTemp = tStats->getTemperature(); -// instaPress = tStats->getPressure(); -// instaVol = tStats->getVolume(); - -// // first evolve chi a half step - -// chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2; -// eta += dt2 * ( instaVol * (instaPress - targetPressure) / -// (p_convert*NkBT*tb2)); - -// integralOfChidt += dt2* chi; - -// for( i=0; igetVel( vel ); -// atoms[i]->getPos( pos ); -// atoms[i]->getFrc( frc ); - -// mass = atoms[i]->getMass(); - -// for (j=0; j < 3; j++) { -// vel[j] += dt2 * ((frc[j] / mass ) * eConvert - vel[j]*(chi+eta)); -// rj[j] = pos[j]; -// } - -// atoms[i]->setVel( vel ); - -// info->wrapVector(rj); - -// for (j = 0; j < 3; j++) -// pos[j] += dt * (vel[j] + eta*rj[j]); - -// atoms[i]->setPos( pos ); - -// if( atoms[i]->isDirectional() ){ - -// dAtom = (DirectionalAtom *)atoms[i]; - -// // get and convert the torque to body frame - -// dAtom->getTrq( Tb ); -// dAtom->lab2Body( Tb ); - -// // get the angular momentum, and propagate a half step - -// dAtom->getJ( ji ); - -// for (j=0; j < 3; j++) -// ji[j] += dt2 * (Tb[j] * eConvert - ji[j]*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] / I[0][0]; -// this->rotate( 1, 2, angle, ji, A ); - -// // rotate about the y-axis -// angle = dt2 * ji[1] / I[1][1]; -// this->rotate( 2, 0, angle, ji, A ); - -// // rotate about the z-axis -// angle = dt * ji[2] / I[2][2]; -// this->rotate( 0, 1, angle, ji, A); - -// // rotate about the y-axis -// angle = dt2 * ji[1] / I[1][1]; -// this->rotate( 2, 0, angle, ji, A ); - -// // rotate about the x-axis -// angle = dt2 * ji[0] / I[0][0]; -// this->rotate( 1, 2, angle, ji, A ); - -// dAtom->setJ( ji ); -// dAtom->setA( A ); -// } - -// } - -// // Scale the box after all the positions have been moved: - -// scaleFactor = exp(dt*eta); - -// if ((scaleFactor > 1.1) || (scaleFactor < 0.9)) { -// sprintf( painCave.errMsg, -// "NPTi error: Attempting a Box scaling of more than 10 percent" -// " check your tauBarostat, as it is probably too small!\n" -// " eta = %lf, scaleFactor = %lf\n", eta, scaleFactor -// ); -// painCave.isFatal = 1; -// simError(); -// } else { -// info->scaleBox(exp(dt*eta)); -// } - //new version of NPTi int i, j, k; @@ -280,12 +167,8 @@ template void NPTi::moveA() { atoms[i]->getPos(pos); for(j = 0; j < 3; j++) - rj[j] = (oldPos[i*3 + j] + pos[j])/2 - COM[j]; - + rj[j] = (oldPos[i*3 + j] + pos[j])/2 - COM[j]; - //wrapVector(r(t)) = r(t)-R0 - //info->wrapVector(rj); - for(j = 0; j < 3; j++) pos[j] = oldPos[i*3 + j] + dt*(vel[j] + eta*rj[j]); @@ -312,68 +195,9 @@ template void NPTi::moveA() { info->scaleBox(scaleFactor); } - //advance volume; - volume = volume * exp(dt*eta); } template void NPTi::moveB( void ){ - -/* - int i, j; - DirectionalAtom* dAtom; - double Tb[3], ji[3]; - double vel[3], frc[3]; - double mass; - - double instaTemp, instaPress, instaVol; - double tt2, tb2; - - tt2 = tauThermostat * tauThermostat; - tb2 = tauBarostat * tauBarostat; - - instaTemp = tStats->getTemperature(); - instaPress = tStats->getPressure(); - instaVol = tStats->getVolume(); - - chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2; - eta += dt2 * ( instaVol * (instaPress - targetPressure) / - (p_convert*NkBT*tb2)); - integralOfChidt += dt2*chi; - - for( i=0; igetVel( vel ); - atoms[i]->getFrc( frc ); - - mass = atoms[i]->getMass(); - - // velocity half step - for (j=0; j < 3; j++) - vel[j] += dt2 * ((frc[j] / mass ) * eConvert - vel[j]*(chi+eta)); - - atoms[i]->setVel( vel ); - - if( atoms[i]->isDirectional() ){ - - dAtom = (DirectionalAtom *)atoms[i]; - - // get and convert the torque to body frame - - dAtom->getTrq( Tb ); - dAtom->lab2Body( Tb ); - - // get the angular momentum, and propagate a half step - - dAtom->getJ( ji ); - - for (j=0; j < 3; j++) - ji[j] += dt2 * (Tb[j] * eConvert - ji[j]*chi); - - dAtom->setJ( ji ); - } - } - -*/ //new version of NPTi int i, j, k; @@ -566,17 +390,6 @@ template double NPTi::getConservedQuant double PV; double extra; - static double pre_U; - static double pre_TS; - static double pre_PV; - static double pre_extra; - static int hackCount = 0; - - double delta_U; - double delta_TS; - double delta_PV; - double delta_extra; - U = tStats->getTotalE(); TS = fkBT * @@ -587,21 +400,9 @@ template double NPTi::getConservedQuant tb2 = tauBarostat * tauBarostat; eta2 = eta * eta; - extra = (fkBT * tb2 * eta2 / 2.0 ) / eConvert; - /* - if(hackCount == 0){ - pre_U = U; - pre_TS =TS; - pre_PV = PV; - pre_extra =extra; - hackCount ++; - } - delta_U = U - pre_U; - delta_TS = TS - pre_TS; - delta_PV = PV - pre_PV; - delta_extra = extra - pre_extra; -*/ + extra = ((double)info->ndfTrans * kB * targetTemp * tb2 * eta2 / 2.0) / eConvert; + cout.width(8); cout.precision(8); @@ -615,19 +416,6 @@ template double NPTi::getConservedQuant << extra << "\t" << U+TS+PV+extra << endl; -/* - pre_U = U; - pre_TS =TS; - pre_PV = PV; - pre_extra =extra; - - - cout << info->getTime() << "\t" - << U << "\t" - << U+TS << "\t" - << U+TS+PV << "\t" - << U+TS+PV+extra << endl; -*/ conservedQuantity = U+TS+PV+extra; return conservedQuantity; }