--- trunk/OOPSE/libmdtools/NPTf.cpp 2003/07/09 01:41:11 577 +++ trunk/OOPSE/libmdtools/NPTf.cpp 2003/07/10 17:10:56 588 @@ -9,7 +9,7 @@ #include "simError.h" -// Basic isotropic thermostating and barostating via the Melchionna +// Basic non-isotropic thermostating and barostating via the Melchionna // modification of the Hoover algorithm: // // Melchionna, S., Ciccotti, G., and Holian, B. L., 1993, @@ -22,9 +22,13 @@ NPTf::NPTf ( SimInfo *theInfo, ForceFields* the_ff): NPTf::NPTf ( SimInfo *theInfo, ForceFields* the_ff): Integrator( theInfo, the_ff ) { - int i; + int i, j; chi = 0.0; - for(i = 0; i < 9; i++) eta[i] = 0.0; + + for(i = 0; i < 3; i++) + for (j = 0; j < 3; j_++) + eta[i][j] = 0.0; + have_tau_thermostat = 0; have_tau_barostat = 0; have_target_temp = 0; @@ -38,74 +42,75 @@ void NPTf::moveA() { DirectionalAtom* dAtom; double Tb[3]; double ji[3]; - double rj[3]; - double instaTemp, instaPress, instaVol; - double tt2, tb2; + double ri[3], vi[3], sc[3]; + double instaTemp, instaVol; + double tt2, tb2, eta2ij; double angle; - double press[9]; - const double p_convert = 1.63882576e8; + double press[3][3], vScale[3][3], hm[3][3], hmnew[3][3], scaleMat[3][3]; tt2 = tauThermostat * tauThermostat; tb2 = tauBarostat * tauBarostat; instaTemp = tStats->getTemperature(); tStats->getPressureTensor(press); - - for (i=0; i < 9; i++) press[i] *= p_convert; - instaVol = tStats->getVolume(); // first evolve chi a half step chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2; - - eta[0] += dt2 * instaVol * (press[0] - targetPressure) / (NkBT*tb2); - eta[1] += dt2 * instaVol * press[1] / (NkBT*tb2); - eta[2] += dt2 * instaVol * press[2] / (NkBT*tb2); - eta[3] += dt2 * instaVol * press[3] / (NkBT*tb2); - eta[4] += dt2 * instaVol * (press[4] - targetPressure) / (NkBT*tb2); - eta[5] += dt2 * instaVol * press[5] / (NkBT*tb2); - eta[6] += dt2 * instaVol * press[6] / (NkBT*tb2); - eta[7] += dt2 * instaVol * press[7] / (NkBT*tb2); - eta[8] += dt2 * instaVol * (press[8] - targetPressure) / (NkBT*tb2); - + + 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 { + + eta[i][j] += dt2 * instaVol * press[i][j] / (NkBT*tb2); + + vScale[i][j] = eta[i][j]; + + } + } + } + for( i=0; imatVecMul3( vScale, vi, sc ); - vx += dt2 * ((frc[atomIndex] /atoms[i]->getMass())*eConvert - scx); - vy += dt2 * ((frc[atomIndex+1]/atoms[i]->getMass())*eConvert - scy); - vz += dt2 * ((frc[atomIndex+2]/atoms[i]->getMass())*eConvert - scz); + 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]); - vel[atomIndex] = vx; - vel[atomIndex+1] = vy; - vel[atomIndex+2] = vz; + vel[atomIndex] = vi[0] + vel[atomIndex+1] = vi[1]; + vel[atomIndex+2] = vi[2]; // position whole step - rj[0] = pos[atomIndex]; - rj[1] = pos[atomIndex+1]; - rj[2] = pos[atomIndex+2]; + ri[0] = pos[atomIndex]; + ri[1] = pos[atomIndex+1]; + ri[2] = pos[atomIndex+2]; - info->wrapVector(rj); + info->wrapVector(ri); - scx = eta[0]*rj[0] + eta[1]*rj[1] + eta[2]*rj[2]; - scy = eta[3]*rj[0] + eta[4]*rj[1] + eta[5]*rj[2]; - scz = eta[6]*rj[0] + eta[7]*rj[1] + eta[8]*rj[2]; + info->matVecMul3( eta, ri, sc ); - pos[atomIndex] += dt * (vel[atomIndex] + scx); - pos[atomIndex+1] += dt * (vel[atomIndex+1] + scy); - pos[atomIndex+2] += dt * (vel[atomIndex+2] + scz); + 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]); if( atoms[i]->isDirectional() ){ @@ -160,51 +165,97 @@ void NPTf::moveA() { } // Scale the box after all the positions have been moved: - - - // Use a taylor expansion for eta products - - info->getBoxM(hm); - + // 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]; + } + + scaleMat[i][j] = 0.0; + // identity matrix (see above): + 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->scaleBox(exp(dt*eta)); - - + } + } + + info->getBoxM(hm); + info->matMul3(hm, scaleMat, hmnew); + info->setBoxM(hmnew); + } -void NPTi::moveB( void ){ - int i,j,k; +void NPTf::moveB( void ){ + int i,j, k; int atomIndex; DirectionalAtom* dAtom; double Tb[3]; double ji[3]; - double instaTemp, instaPress, instaVol; + double vi[3], sc[3]; + double instaTemp, instaVol; double tt2, tb2; + double press[3][3], vScale[3][3]; tt2 = tauThermostat * tauThermostat; tb2 = tauBarostat * tauBarostat; instaTemp = tStats->getTemperature(); - instaPress = tStats->getPressure(); + tStats->getPressureTensor(press); instaVol = tStats->getVolume(); - + + // first evolve chi a half step + chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2; - eta += dt2 * ( instaVol * (instaPress - targetPressure) / (NkBT*tb2)); + 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 { + + eta[i][j] += dt2 * instaVol * press[i][j] / (NkBT*tb2); + + vScale[i][j] = eta[i][j]; + + } + } + } + for( i=0; igetMass())*eConvert - - vel[j]*(chi+eta)); + 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]); + + vel[atomIndex] = vi[0] + vel[atomIndex+1] = vi[1]; + vel[atomIndex+2] = vi[2]; + if( atoms[i]->isDirectional() ){ dAtom = (DirectionalAtom *)atoms[i]; @@ -235,14 +286,14 @@ void NPTi::moveB( void ){ } } -int NPTi::readyCheck() { +int NPTf::readyCheck() { // First check to see if we have a target temperature. // Not having one is fatal. if (!have_target_temp) { sprintf( painCave.errMsg, - "NPTi error: You can't use the NPTi integrator\n" + "NPTf error: You can't use the NPTf integrator\n" " without a targetTemp!\n" ); painCave.isFatal = 1; @@ -252,7 +303,7 @@ int NPTi::readyCheck() { if (!have_target_pressure) { sprintf( painCave.errMsg, - "NPTi error: You can't use the NPTi integrator\n" + "NPTf error: You can't use the NPTf integrator\n" " without a targetPressure!\n" ); painCave.isFatal = 1; @@ -264,7 +315,7 @@ int NPTi::readyCheck() { if (!have_tau_thermostat) { sprintf( painCave.errMsg, - "NPTi error: If you use the NPTi\n" + "NPTf error: If you use the NPTf\n" " integrator, you must set tauThermostat.\n"); painCave.isFatal = 1; simError(); @@ -275,7 +326,7 @@ int NPTi::readyCheck() { if (!have_tau_barostat) { sprintf( painCave.errMsg, - "NPTi error: If you use the NPTi\n" + "NPTf error: If you use the NPTf\n" " integrator, you must set tauBarostat.\n"); painCave.isFatal = 1; simError();