--- trunk/src/brains/Thermo.cpp 2005/04/15 22:04:00 507 +++ trunk/src/brains/Thermo.cpp 2006/05/17 21:51:42 963 @@ -53,7 +53,7 @@ namespace oopse { namespace oopse { - double Thermo::getKinetic() { + RealType Thermo::getKinetic() { SimInfo::MoleculeIterator miter; std::vector::iterator iiter; Molecule* mol; @@ -64,18 +64,18 @@ namespace oopse { int i; int j; int k; - double kinetic = 0.0; - double kinetic_global = 0.0; + RealType kinetic = 0.0; + RealType kinetic_global = 0.0; for (mol = info_->beginMolecule(miter); mol != NULL; mol = info_->nextMolecule(miter)) { for (integrableObject = mol->beginIntegrableObject(iiter); integrableObject != NULL; integrableObject = mol->nextIntegrableObject(iiter)) { - - double mass = integrableObject->getMass(); + + RealType mass = integrableObject->getMass(); Vector3d vel = integrableObject->getVel(); - + kinetic += mass * (vel[0]*vel[0] + vel[1]*vel[1] + vel[2]*vel[2]); - + if (integrableObject->isDirectional()) { angMom = integrableObject->getJ(); I = integrableObject->getI(); @@ -96,7 +96,7 @@ namespace oopse { #ifdef IS_MPI - MPI_Allreduce(&kinetic, &kinetic_global, 1, MPI_DOUBLE, MPI_SUM, + MPI_Allreduce(&kinetic, &kinetic_global, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD); kinetic = kinetic_global; @@ -107,53 +107,53 @@ namespace oopse { return kinetic; } - double Thermo::getPotential() { - double potential = 0.0; + RealType Thermo::getPotential() { + RealType potential = 0.0; Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot(); - double potential_local = curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL] + - curSnapshot->statData[Stats::SHORT_RANGE_POTENTIAL] ; + RealType shortRangePot_local = curSnapshot->statData[Stats::SHORT_RANGE_POTENTIAL] ; // Get total potential for entire system from MPI. #ifdef IS_MPI - MPI_Allreduce(&potential_local, &potential, 1, MPI_DOUBLE, MPI_SUM, + MPI_Allreduce(&shortRangePot_local, &potential, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD); + potential += curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL]; #else - potential = potential_local; + potential = shortRangePot_local + curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL]; #endif // is_mpi return potential; } - double Thermo::getTotalE() { - double total; + RealType Thermo::getTotalE() { + RealType total; total = this->getKinetic() + this->getPotential(); return total; } - double Thermo::getTemperature() { + RealType Thermo::getTemperature() { - double temperature = ( 2.0 * this->getKinetic() ) / (info_->getNdf()* OOPSEConstant::kb ); + RealType temperature = ( 2.0 * this->getKinetic() ) / (info_->getNdf()* OOPSEConstant::kb ); return temperature; } - double Thermo::getVolume() { + RealType Thermo::getVolume() { Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot(); return curSnapshot->getVolume(); } - double Thermo::getPressure() { + RealType Thermo::getPressure() { // Relies on the calculation of the full molecular pressure tensor Mat3x3d tensor; - double pressure; + RealType pressure; tensor = getPressureTensor(); @@ -162,6 +162,23 @@ namespace oopse { return pressure; } + RealType Thermo::getPressure(int direction) { + + // Relies on the calculation of the full molecular pressure tensor + + + Mat3x3d tensor; + RealType pressure; + + tensor = getPressureTensor(); + + pressure = OOPSEConstant::pressureConvert * tensor(direction, direction); + + return pressure; + } + + + Mat3x3d Thermo::getPressureTensor() { // returns pressure tensor in units amu*fs^-2*Ang^-1 // routine derived via viral theorem description in: @@ -178,19 +195,19 @@ namespace oopse { for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL; integrableObject = mol->nextIntegrableObject(j)) { - double mass = integrableObject->getMass(); + RealType mass = integrableObject->getMass(); Vector3d vcom = integrableObject->getVel(); p_local += mass * outProduct(vcom, vcom); } } #ifdef IS_MPI - MPI_Allreduce(p_local.getArrayPointer(), p_global.getArrayPointer(), 9, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(p_local.getArrayPointer(), p_global.getArrayPointer(), 9, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD); #else p_global = p_local; #endif // is_mpi - double volume = this->getVolume(); + RealType volume = this->getVolume(); Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot(); Mat3x3d tau = curSnapshot->statData.getTau(); @@ -210,6 +227,12 @@ namespace oopse { stat[Stats::PRESSURE] = getPressure(); stat[Stats::VOLUME] = getVolume(); + Mat3x3d tensor =getPressureTensor(); + stat[Stats::PRESSURE_TENSOR_X] = tensor(0, 0); + stat[Stats::PRESSURE_TENSOR_Y] = tensor(1, 1); + stat[Stats::PRESSURE_TENSOR_Z] = tensor(2, 2); + + /**@todo need refactorying*/ //Conserved Quantity is set by integrator and time is set by setTime