--- trunk/OOPSE/libmdtools/NPT.cpp 2003/10/28 16:03:37 829 +++ trunk/OOPSE/libmdtools/NPT.cpp 2003/10/29 00:19:10 837 @@ -7,7 +7,7 @@ #include "Thermo.hpp" #include "ReadWrite.hpp" #include "Integrator.hpp" -#include "simError.h" +#include "simError.h" #ifdef IS_MPI #include "mpiSimulation.hpp" @@ -18,15 +18,22 @@ // modification of the Hoover algorithm: // // Melchionna, S., Ciccotti, G., and Holian, B. L., 1993, -// Molec. Phys., 78, 533. +// Molec. Phys., 78, 533. // // and -// +// // Hoover, W. G., 1986, Phys. Rev. A, 34, 2499. template NPT::NPT ( SimInfo *theInfo, ForceFields* the_ff): T( theInfo, the_ff ) { + GenericData* data; + DoubleData * chiValue; + DoubleData * integralOfChidtValue; + + chiValue = NULL; + integralOfChidtValue = NULL; + chi = 0.0; integralOfChidt = 0.0; have_tau_thermostat = 0; @@ -37,6 +44,23 @@ template NPT::NPT ( SimInfo *theInfo, F have_eta_tolerance = 0; have_pos_iter_tolerance = 0; + // retrieve chi and integralOfChidt from simInfo + data = info->getProperty(CHIVALUE_ID); + if(data){ + chiValue = dynamic_cast(data); + } + + data = info->getProperty(INTEGRALOFCHIDT_ID); + if(data){ + integralOfChidtValue = dynamic_cast(data); + } + + // chi and integralOfChidt should appear by pair + if(chiValue && integralOfChidtValue){ + chi = chiValue->getData(); + integralOfChidt = integralOfChidtValue->getData(); + } + oldPos = new double[3*nAtoms]; oldVel = new double[3*nAtoms]; oldJi = new double[3*nAtoms]; @@ -55,7 +79,7 @@ template void NPT::moveA() { } template void NPT::moveA() { - + //new version of NPT int i, j, k; DirectionalAtom* dAtom; @@ -69,9 +93,9 @@ template void NPT::moveA() { tStats->getPressureTensor( press ); instaPress = p_convert * (press[0][0] + press[1][1] + press[2][2]) / 3.0; instaVol = tStats->getVolume(); - + tStats->getCOM(COM); - + //evolve velocity half step for( i=0; i void NPT::moveA() { getVelScaleA( sc, vel ); for (j=0; j < 3; j++) { - + // velocity half step (use chi from previous step here): vel[j] += dt2 * ((frc[j] / mass ) * eConvert - sc[j]); - + } 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++) + for (j=0; j < 3; j++) ji[j] += dt2 * (Tb[j] * eConvert - ji[j]*chi); - + this->rotationPropagation( dAtom, ji ); - + dAtom->setJ( ji ); - } + } } // evolve chi and eta half step - + evolveChiA(); evolveEtaA(); @@ -127,9 +151,9 @@ template void NPT::moveA() { for(j = 0; j < 3; j++) oldPos[i*3 + j] = pos[j]; } - - //the first estimation of r(t+dt) is equal to r(t) - + + //the first estimation of r(t+dt) is equal to r(t) + for(k = 0; k < 5; k ++){ for(i =0 ; i < nAtoms; i++){ @@ -138,33 +162,33 @@ template void NPT::moveA() { atoms[i]->getPos(pos); this->getPosScale( pos, COM, i, sc ); - + for(j = 0; j < 3; j++) pos[j] = oldPos[i*3 + j] + dt*(vel[j] + sc[j]); atoms[i]->setPos( pos ); } - + if (nConstrained){ constrainA(); } } - + // Scale the box after all the positions have been moved: - + this->scaleSimBox(); } template void NPT::moveB( void ){ - + //new version of NPT int i, j, k; DirectionalAtom* dAtom; double Tb[3], ji[3], sc[3]; double vel[3], frc[3]; double mass; - + // Set things up for the iteration: for( i=0; i void NPT::moveB( void ){ // do the iteration: instaVol = tStats->getVolume(); - + for (k=0; k < 4; k++) { - + instaTemp = tStats->getTemperature(); instaPress = tStats->getPressure(); @@ -199,42 +223,42 @@ template void NPT::moveB( void ){ this->evolveChiB(); this->evolveEtaB(); - + for( i=0; igetFrc( frc ); atoms[i]->getVel(vel); - + mass = atoms[i]->getMass(); - + getVelScaleB( sc, i ); // velocity half step - for (j=0; j < 3; j++) + for (j=0; j < 3; j++) vel[j] = oldVel[3*i+j] + dt2 * ((frc[j] / mass ) * eConvert - sc[j]); - + atoms[i]->setVel( vel ); - + if( atoms[i]->isDirectional() ){ dAtom = (DirectionalAtom *)atoms[i]; - - // get and convert the torque to body frame - + + // get and convert the torque to body frame + dAtom->getTrq( Tb ); - dAtom->lab2Body( Tb ); - - for (j=0; j < 3; j++) + 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 (nConstrained){ constrainB(); - } - + } + if ( this->chiConverged() && this->etaConverged() ) break; } @@ -255,14 +279,14 @@ template void NPT::evolveChiB() { } template void NPT::evolveChiB() { - + prevChi = chi; chi = oldChi + dt2 * ( instaTemp / targetTemp - 1.0) / tt2; } template bool NPT::chiConverged() { - - return ( fabs( prevChi - chi ) <= chiTolerance ); + + return ( fabs( prevChi - chi ) <= chiTolerance ); } template int NPT::readyCheck() { @@ -270,10 +294,10 @@ template int NPT::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. - + + // First check to see if we have a target temperature. + // Not having one is fatal. + if (!have_target_temp) { sprintf( painCave.errMsg, "NPT error: You can't use the NPT integrator\n" @@ -293,9 +317,9 @@ template int NPT::readyCheck() { simError(); return -1; } - + // We must set tauThermostat. - + if (!have_tau_thermostat) { sprintf( painCave.errMsg, "NPT error: If you use the NPT\n" @@ -303,10 +327,10 @@ template int NPT::readyCheck() { painCave.isFatal = 1; simError(); return -1; - } + } // We must set tauBarostat. - + if (!have_tau_barostat) { sprintf( painCave.errMsg, "NPT error: If you use the NPT\n" @@ -314,7 +338,7 @@ template int NPT::readyCheck() { painCave.isFatal = 1; simError(); return -1; - } + } if (!have_chi_tolerance) { sprintf( painCave.errMsg, @@ -323,7 +347,7 @@ template int NPT::readyCheck() { have_chi_tolerance = 1; painCave.isFatal = 0; simError(); - } + } if (!have_eta_tolerance) { sprintf( painCave.errMsg, @@ -332,18 +356,18 @@ template int NPT::readyCheck() { have_eta_tolerance = 1; painCave.isFatal = 0; simError(); - } - + } + // We need NkBT a lot, so just set it here: This is the RAW number // of particles, so no subtraction or addition of constraints or // orientational degrees of freedom: - + NkBT = (double)Nparticles * kB * targetTemp; - + // fkBT is used because the thermostat operates on more degrees of freedom // than the barostat (when there are particles with orientational degrees // of freedom). ndf = 3 * (n_atoms + n_oriented -1) - n_constraint - nZcons - + fkBT = (double)info->ndf * kB * targetTemp; tt2 = tauThermostat * tauThermostat;