--- trunk/OOPSE/libmdtools/NPTi.cpp 2003/07/22 19:54:52 645 +++ trunk/OOPSE/libmdtools/NPTi.cpp 2003/09/16 20:02:11 767 @@ -9,7 +9,11 @@ #include "Integrator.hpp" #include "simError.h" +#ifdef IS_MPI +#include "mpiSimulation.hpp" +#endif + // Basic isotropic thermostating and barostating via the Melchionna // modification of the Hoover algorithm: // @@ -25,15 +29,36 @@ template NPTi::NPTi ( SimInfo *theInfo, { chi = 0.0; eta = 0.0; + integralOfChidt = 0.0; have_tau_thermostat = 0; have_tau_barostat = 0; have_target_temp = 0; have_target_pressure = 0; + have_chi_tolerance = 0; + have_eta_tolerance = 0; + have_pos_iter_tolerance = 0; + + oldPos = new double[3*nAtoms]; + oldVel = new double[3*nAtoms]; + oldJi = new double[3*nAtoms]; +#ifdef IS_MPI + Nparticles = mpiSim->getTotAtoms(); +#else + Nparticles = theInfo->n_atoms; +#endif + } +template NPTi::~NPTi() { + delete[] oldPos; + delete[] oldVel; + delete[] oldJi; +} + template void NPTi::moveA() { - - int i, j; + + //new version of NPTi + int i, j, k; DirectionalAtom* dAtom; double Tb[3], ji[3]; double A[3][3], I[3][3]; @@ -43,6 +68,7 @@ template void NPTi::moveA() { double rj[3]; double instaTemp, instaPress, instaVol; double tt2, tb2, scaleFactor; + double COM[3]; tt2 = tauThermostat * tauThermostat; tb2 = tauBarostat * tauBarostat; @@ -50,38 +76,29 @@ template void NPTi::moveA() { 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)); - + tStats->getCOM(COM); + + //evolve velocity half step 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]; + // velocity half step (use chi from previous step here): + vel[j] += dt2 * ((frc[j] / mass ) * eConvert - vel[j]*(chi + eta)); + } 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 ); @@ -122,10 +139,46 @@ template void NPTi::moveA() { dAtom->setJ( ji ); dAtom->setA( A ); - } + } + } + // evolve chi and eta half step + + chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2; + eta += dt2 * ( instaVol * (instaPress - targetPressure) / (p_convert*NkBT*tb2)); + + //calculate the integral of chidt + integralOfChidt += dt2*chi; + + //save the old positions + for(i = 0; i < nAtoms; i++){ + atoms[i]->getPos(pos); + for(j = 0; j < 3; j++) + oldPos[i*3 + j] = pos[j]; } + + //the first estimation of r(t+dt) is equal to r(t) + + for(k = 0; k < 4; k ++){ + for(i =0 ; i < nAtoms; i++){ + + atoms[i]->getVel(vel); + atoms[i]->getPos(pos); + + for(j = 0; j < 3; j++) + rj[j] = (oldPos[i*3 + j] + pos[j])/2 - COM[j]; + + for(j = 0; j < 3; j++) + pos[j] = oldPos[i*3 + j] + dt*(vel[j] + eta*rj[j]); + + atoms[i]->setPos( pos ); + + } + + } + + // Scale the box after all the positions have been moved: scaleFactor = exp(dt*eta); @@ -139,68 +192,122 @@ template void NPTi::moveA() { painCave.isFatal = 1; simError(); } else { - info->scaleBox(exp(dt*eta)); - } + info->scaleBox(scaleFactor); + } } template void NPTi::moveB( void ){ - - int i, j; + + //new version of NPTi + int i, j, k; DirectionalAtom* dAtom; double Tb[3], ji[3]; double vel[3], frc[3]; double mass; - double instaTemp, instaPress, instaVol; + double instTemp, instPress, instVol; double tt2, tb2; + double oldChi, prevChi; + double oldEta, preEta; 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)); - + // Set things up for the iteration: + + oldChi = chi; + oldEta = eta; + for( i=0; igetVel( vel ); - atoms[i]->getFrc( frc ); - mass = atoms[i]->getMass(); + for (j=0; j < 3; j++) + oldVel[3*i + j] = vel[j]; - // 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); + for (j=0; j < 3; j++) + oldJi[3*i + j] = ji[j]; - dAtom->setJ( ji ); } } + + // do the iteration: + + instVol = tStats->getVolume(); + + for (k=0; k < 4; k++) { + + instTemp = tStats->getTemperature(); + instPress = tStats->getPressure(); + + // evolve chi another half step using the temperature at t + dt/2 + + prevChi = chi; + chi = oldChi + dt2 * ( instTemp / targetTemp - 1.0) / + (tauThermostat*tauThermostat); + + preEta = eta; + eta = oldEta + dt2 * ( instVol * (instPress - targetPressure) / + (p_convert*NkBT*tb2)); + + + for( i=0; igetFrc( frc ); + atoms[i]->getVel(vel); + + mass = atoms[i]->getMass(); + + // velocity half step + for (j=0; j < 3; j++) + vel[j] = oldVel[3*i+j] + dt2 * ((frc[j] / mass ) * eConvert - oldVel[3*i + 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 ); + + for (j=0; j < 3; j++) + ji[j] = oldJi[3*i + j] + dt2 * (Tb[j] * eConvert - oldJi[3*i+j]*chi); + + dAtom->setJ( ji ); + } + } + + if (fabs(prevChi - chi) <= chiTolerance && fabs(preEta -eta) <= etaTolerance) + break; + } + + //calculate integral of chida + integralOfChidt += dt2*chi; + + } +template void NPTi::resetIntegrator() { + chi = 0.0; + eta = 0.0; +} + template int NPTi::readyCheck() { + + //check parent's readyCheck() first + if (T::readyCheck() == -1) + return -1; // First check to see if we have a target temperature. // Not having one is fatal. @@ -247,9 +354,68 @@ template int NPTi::readyCheck() { return -1; } + if (!have_chi_tolerance) { + sprintf( painCave.errMsg, + "NPTi warning: setting chi tolerance to 1e-6\n"); + chiTolerance = 1e-6; + have_chi_tolerance = 1; + painCave.isFatal = 0; + simError(); + } + + if (!have_eta_tolerance) { + sprintf( painCave.errMsg, + "NPTi warning: setting eta tolerance to 1e-6\n"); + etaTolerance = 1e-6; + have_eta_tolerance = 1; + painCave.isFatal = 0; + simError(); + } // We need NkBT a lot, so just set it here: - NkBT = (double)info->ndf * kB * targetTemp; + NkBT = (double)Nparticles * kB * targetTemp; + fkBT = (double)info->ndf * kB * targetTemp; return 1; } + +template double NPTi::getConservedQuantity(void){ + + double conservedQuantity; + double tb2; + double eta2; + double E_NPT; + double U; + double TS; + double PV; + double extra; + + U = tStats->getTotalE(); + + TS = fkBT * + (integralOfChidt + tauThermostat * tauThermostat * chi * chi / 2.0) / eConvert; + + PV = (targetPressure * tStats->getVolume() / p_convert) / eConvert; + + tb2 = tauBarostat * tauBarostat; + eta2 = eta * eta; + + + extra = ((double)info->ndfTrans * kB * targetTemp * tb2 * eta2 / 2.0) / eConvert; + + cout.width(8); + cout.precision(8); + + + cout << info->getTime() << "\t" + << chi << "\t" + << eta << "\t" + << U << "\t" + << TS << "\t" + << PV << "\t" + << extra << "\t" + << U+TS+PV+extra << endl; + + conservedQuantity = U+TS+PV+extra; + return conservedQuantity; +}