--- trunk/OOPSE/libmdtools/NVT.cpp 2003/09/19 20:00:27 778 +++ trunk/OOPSE/libmdtools/NVT.cpp 2003/10/29 00:19:10 837 @@ -6,7 +6,7 @@ #include "Thermo.hpp" #include "ReadWrite.hpp" #include "Integrator.hpp" -#include "simError.h" +#include "simError.h" // Basic thermostating via Hoover, Phys.Rev.A, 1985, Vol. 31 (5) 1695-1697 @@ -14,12 +14,36 @@ template NVT::NVT ( SimInfo *theInfo, F template NVT::NVT ( SimInfo *theInfo, ForceFields* the_ff): T( theInfo, the_ff ) { + GenericData* data; + DoubleData * chiValue; + DoubleData * integralOfChidtValue; + + chiValue = NULL; + integralOfChidtValue = NULL; + chi = 0.0; have_tau_thermostat = 0; have_target_temp = 0; have_chi_tolerance = 0; integralOfChidt = 0.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(); + } + oldVel = new double[3*nAtoms]; oldJi = new double[3*nAtoms]; } @@ -30,7 +54,7 @@ template void NVT::moveA() { } template void NVT::moveA() { - + int i, j; DirectionalAtom* dAtom; double Tb[3], ji[3]; @@ -42,7 +66,7 @@ template void NVT::moveA() { // We need the temperature at time = t for the chi update below: instTemp = tStats->getTemperature(); - + for( i=0; igetVel( vel ); @@ -60,34 +84,34 @@ template void NVT::moveA() { atoms[i]->setVel( vel ); 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++) + for (j=0; j < 3; j++) ji[j] += dt2 * (Tb[j] * eConvert - ji[j]*chi); - + this->rotationPropagation( dAtom, ji ); - + dAtom->setJ( ji ); - } + } } - + if (nConstrained){ constrainA(); } - // Finally, evolve chi a half step (just like a velocity) using + // Finally, evolve chi a half step (just like a velocity) using // temperature at time t, not time t+dt/2 chi += dt2 * ( instTemp / targetTemp - 1.0) / (tauThermostat*tauThermostat); @@ -130,40 +154,40 @@ template void NVT::moveB( void ){ // do the iteration: for (k=0; k < 4; k++) { - + instTemp = tStats->getTemperature(); // evolve chi another half step using the temperature at t + dt/2 prevChi = chi; - chi = oldChi + dt2 * ( instTemp / targetTemp - 1.0) / + chi = oldChi + dt2 * ( instTemp / targetTemp - 1.0) / (tauThermostat*tauThermostat); - + for( i=0; igetFrc( frc ); atoms[i]->getVel(vel); - + mass = atoms[i]->getMass(); - + // 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 - oldVel[3*i + j]*chi); - + 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 ); } } @@ -174,12 +198,12 @@ template void NVT::moveB( void ){ if (fabs(prevChi - chi) <= chiTolerance) break; } - + integralOfChidt += dt2*chi; } template void NVT::resetIntegrator( void ){ - + chi = 0.0; integralOfChidt = 0.0; } @@ -189,10 +213,10 @@ template int NVT::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, "NVT error: You can't use the NVT integrator without a targetTemp!\n" @@ -201,9 +225,9 @@ template int NVT::readyCheck() { simError(); return -1; } - + // We must set tauThermostat. - + if (!have_tau_thermostat) { sprintf( painCave.errMsg, "NVT error: If you use the constant temperature\n" @@ -211,7 +235,7 @@ template int NVT::readyCheck() { painCave.isFatal = 1; simError(); return -1; - } + } if (!have_chi_tolerance) { sprintf( painCave.errMsg, @@ -220,9 +244,9 @@ template int NVT::readyCheck() { have_chi_tolerance = 1; painCave.isFatal = 0; simError(); - } + } - return 1; + return 1; } @@ -234,19 +258,30 @@ template double NVT::getConservedQuanti double thermostat_kinetic; double thermostat_potential; - fkBT = (double)(info->getNDF() ) * kB * targetTemp; - + fkBT = (double)(info->getNDF() ) * kB * targetTemp; + Energy = tStats->getTotalE(); - thermostat_kinetic = fkBT* tauThermostat * tauThermostat * chi * chi / + thermostat_kinetic = fkBT* tauThermostat * tauThermostat * chi * chi / (2.0 * eConvert); thermostat_potential = fkBT * integralOfChidt / eConvert; conservedQuantity = Energy + thermostat_kinetic + thermostat_potential; - - cerr << info->getTime() << "\t" << Energy << "\t" << thermostat_kinetic << + + cerr << info->getTime() << "\t" << Energy << "\t" << thermostat_kinetic << "\t" << thermostat_potential << "\t" << conservedQuantity << endl; - return conservedQuantity; + return conservedQuantity; } + +template string NVT::getAdditionalParameters(void){ + string parameters; + const int BUFFERSIZE = 2000; // size of the read buffer + char buffer[BUFFERSIZE]; + + sprintf(buffer,"\t%g\t%g;", chi, integralOfChidt); + parameters += buffer; + + return parameters; +}