--- trunk/src/brains/Thermo.cpp 2004/09/24 04:16:43 2 +++ branches/development/src/brains/Thermo.cpp 2010/07/09 23:08:25 1465 @@ -1,450 +1,307 @@ +/* + * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved. + * + * The University of Notre Dame grants you ("Licensee") a + * non-exclusive, royalty free, license to use, modify and + * redistribute this software in source and binary code form, provided + * that the following conditions are met: + * + * 1. Redistributions of source code must retain the above copyright + * notice, this list of conditions and the following disclaimer. + * + * 2. Redistributions in binary form must reproduce the above copyright + * notice, this list of conditions and the following disclaimer in the + * documentation and/or other materials provided with the + * distribution. + * + * This software is provided "AS IS," without a warranty of any + * kind. All express or implied conditions, representations and + * warranties, including any implied warranty of merchantability, + * fitness for a particular purpose or non-infringement, are hereby + * excluded. The University of Notre Dame and its licensors shall not + * be liable for any damages suffered by licensee as a result of + * using, modifying or distributing the software or its + * derivatives. In no event will the University of Notre Dame or its + * licensors be liable for any lost revenue, profit or data, or for + * direct, indirect, special, consequential, incidental or punitive + * damages, however caused and regardless of the theory of liability, + * arising out of the use of or inability to use software, even if the + * University of Notre Dame has been advised of the possibility of + * such damages. + * + * SUPPORT OPEN SCIENCE! If you use OpenMD or its source code in your + * research, please cite the appropriate papers when you publish your + * work. Good starting points are: + * + * [1] Meineke, et al., J. Comp. Chem. 26, 252-271 (2005). + * [2] Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006). + * [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008). + * [4] Vardeman & Gezelter, in progress (2009). + */ + #include #include -using namespace std; #ifdef IS_MPI #include #endif //is_mpi -#include "Thermo.hpp" -#include "SRI.hpp" -#include "Integrator.hpp" -#include "simError.h" -#include "MatVec3.h" +#include "brains/Thermo.hpp" +#include "primitives/Molecule.hpp" +#include "utils/simError.h" +#include "utils/PhysicalConstants.hpp" -#ifdef IS_MPI -#define __C -#include "mpiSimulation.hpp" -#endif // is_mpi +namespace OpenMD { -inline double roundMe( double x ){ - return ( x >= 0 ) ? floor( x + 0.5 ) : ceil( x - 0.5 ); -} + RealType Thermo::getKinetic() { + SimInfo::MoleculeIterator miter; + std::vector::iterator iiter; + Molecule* mol; + StuntDouble* integrableObject; + Vector3d vel; + Vector3d angMom; + Mat3x3d I; + int i; + int j; + int k; + RealType mass; + 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)) { + + mass = integrableObject->getMass(); + 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(); -Thermo::Thermo( SimInfo* the_info ) { - info = the_info; - int baseSeed = the_info->getSeed(); - - gaussStream = new gaussianSPRNG( baseSeed ); -} + if (integrableObject->isLinear()) { + i = integrableObject->linearAxis(); + j = (i + 1) % 3; + k = (i + 2) % 3; + kinetic += angMom[j] * angMom[j] / I(j, j) + angMom[k] * angMom[k] / I(k, k); + } else { + kinetic += angMom[0]*angMom[0]/I(0, 0) + angMom[1]*angMom[1]/I(1, 1) + + angMom[2]*angMom[2]/I(2, 2); + } + } + + } + } + +#ifdef IS_MPI -Thermo::~Thermo(){ - delete gaussStream; -} + MPI_Allreduce(&kinetic, &kinetic_global, 1, MPI_REALTYPE, MPI_SUM, + MPI_COMM_WORLD); + kinetic = kinetic_global; -double Thermo::getKinetic(){ +#endif //is_mpi - const double e_convert = 4.184E-4; // convert kcal/mol -> (amu A^2)/fs^2 - double kinetic; - double amass; - double aVel[3], aJ[3], I[3][3]; - int i, j, k, kl; + kinetic = kinetic * 0.5 / PhysicalConstants::energyConvert; - double kinetic_global; - vector integrableObjects = info->integrableObjects; - - kinetic = 0.0; - kinetic_global = 0.0; + return kinetic; + } - for (kl=0; klgetVel(aVel); - amass = integrableObjects[kl]->getMass(); + RealType Thermo::getPotential() { + RealType potential = 0.0; + Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot(); + RealType shortRangePot_local = curSnapshot->statData[Stats::SHORT_RANGE_POTENTIAL] ; - for(j=0; j<3; j++) - kinetic += amass*aVel[j]*aVel[j]; + // Get total potential for entire system from MPI. - if (integrableObjects[kl]->isDirectional()){ - - integrableObjects[kl]->getJ( aJ ); - integrableObjects[kl]->getI( I ); - - if (integrableObjects[kl]->isLinear()) { - i = integrableObjects[kl]->linearAxis(); - j = (i+1)%3; - k = (i+2)%3; - kinetic += aJ[j]*aJ[j]/I[j][j] + aJ[k]*aJ[k]/I[k][k]; - } else { - for (j=0; j<3; j++) - kinetic += aJ[j]*aJ[j] / I[j][j]; - } - } - } #ifdef IS_MPI - MPI_Allreduce(&kinetic,&kinetic_global,1,MPI_DOUBLE, - MPI_SUM, MPI_COMM_WORLD); - kinetic = kinetic_global; -#endif //is_mpi - - kinetic = kinetic * 0.5 / e_convert; - return kinetic; -} + MPI_Allreduce(&shortRangePot_local, &potential, 1, MPI_REALTYPE, MPI_SUM, + MPI_COMM_WORLD); + potential += curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL]; -double Thermo::getPotential(){ - - double potential_local; - double potential; - int el, nSRI; - Molecule* molecules; +#else - molecules = info->molecules; - nSRI = info->n_SRI; + potential = shortRangePot_local + curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL]; - potential_local = 0.0; - potential = 0.0; - potential_local += info->lrPot; +#endif // is_mpi - for( el=0; eln_mol; el++ ){ - potential_local += molecules[el].getPotential(); + return potential; } - // Get total potential for entire system from MPI. -#ifdef IS_MPI - MPI_Allreduce(&potential_local,&potential,1,MPI_DOUBLE, - MPI_SUM, MPI_COMM_WORLD); -#else - potential = potential_local; -#endif // is_mpi + RealType Thermo::getTotalE() { + RealType total; - return potential; -} + total = this->getKinetic() + this->getPotential(); + return total; + } -double Thermo::getTotalE(){ + RealType Thermo::getTemperature() { + + RealType temperature = ( 2.0 * this->getKinetic() ) / (info_->getNdf()* PhysicalConstants::kb ); + return temperature; + } - double total; + RealType Thermo::getVolume() { + Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot(); + return curSnapshot->getVolume(); + } - total = this->getKinetic() + this->getPotential(); - return total; -} + RealType Thermo::getPressure() { -double Thermo::getTemperature(){ + // Relies on the calculation of the full molecular pressure tensor - const double kb = 1.9872156E-3; // boltzman's constant in kcal/(mol K) - double temperature; - temperature = ( 2.0 * this->getKinetic() ) / ((double)info->ndf * kb ); - return temperature; -} + Mat3x3d tensor; + RealType pressure; -double Thermo::getVolume() { + tensor = getPressureTensor(); - return info->boxVol; -} + pressure = PhysicalConstants::pressureConvert * (tensor(0, 0) + tensor(1, 1) + tensor(2, 2)) / 3.0; -double Thermo::getPressure() { + return pressure; + } - // Relies on the calculation of the full molecular pressure tensor - - const double p_convert = 1.63882576e8; - double press[3][3]; - double pressure; + RealType Thermo::getPressure(int direction) { - this->getPressureTensor(press); + // Relies on the calculation of the full molecular pressure tensor - pressure = p_convert * (press[0][0] + press[1][1] + press[2][2]) / 3.0; + + Mat3x3d tensor; + RealType pressure; - return pressure; -} + tensor = getPressureTensor(); -double Thermo::getPressureX() { + pressure = PhysicalConstants::pressureConvert * tensor(direction, direction); - // Relies on the calculation of the full molecular pressure tensor - - const double p_convert = 1.63882576e8; - double press[3][3]; - double pressureX; + return pressure; + } - this->getPressureTensor(press); + Mat3x3d Thermo::getPressureTensor() { + // returns pressure tensor in units amu*fs^-2*Ang^-1 + // routine derived via viral theorem description in: + // Paci, E. and Marchi, M. J.Phys.Chem. 1996, 100, 4314-4322 + Mat3x3d pressureTensor; + Mat3x3d p_local(0.0); + Mat3x3d p_global(0.0); - pressureX = p_convert * press[0][0]; + SimInfo::MoleculeIterator i; + std::vector::iterator j; + Molecule* mol; + StuntDouble* integrableObject; + for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) { + for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL; + integrableObject = mol->nextIntegrableObject(j)) { - return pressureX; -} + 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_REALTYPE, MPI_SUM, MPI_COMM_WORLD); +#else + p_global = p_local; +#endif // is_mpi -double Thermo::getPressureY() { + RealType volume = this->getVolume(); + Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot(); + Mat3x3d tau = curSnapshot->statData.getTau(); - // Relies on the calculation of the full molecular pressure tensor - - const double p_convert = 1.63882576e8; - double press[3][3]; - double pressureY; + pressureTensor = (p_global + PhysicalConstants::energyConvert* tau)/volume; + + return pressureTensor; + } - this->getPressureTensor(press); - pressureY = p_convert * press[1][1]; + void Thermo::saveStat(){ + Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot(); + Stats& stat = currSnapshot->statData; + + stat[Stats::KINETIC_ENERGY] = getKinetic(); + stat[Stats::POTENTIAL_ENERGY] = getPotential(); + stat[Stats::TOTAL_ENERGY] = stat[Stats::KINETIC_ENERGY] + stat[Stats::POTENTIAL_ENERGY] ; + stat[Stats::TEMPERATURE] = getTemperature(); + stat[Stats::PRESSURE] = getPressure(); + stat[Stats::VOLUME] = getVolume(); - return pressureY; -} + Mat3x3d tensor =getPressureTensor(); + stat[Stats::PRESSURE_TENSOR_XX] = tensor(0, 0); + stat[Stats::PRESSURE_TENSOR_XY] = tensor(0, 1); + stat[Stats::PRESSURE_TENSOR_XZ] = tensor(0, 2); + stat[Stats::PRESSURE_TENSOR_YX] = tensor(1, 0); + stat[Stats::PRESSURE_TENSOR_YY] = tensor(1, 1); + stat[Stats::PRESSURE_TENSOR_YZ] = tensor(1, 2); + stat[Stats::PRESSURE_TENSOR_ZX] = tensor(2, 0); + stat[Stats::PRESSURE_TENSOR_ZY] = tensor(2, 1); + stat[Stats::PRESSURE_TENSOR_ZZ] = tensor(2, 2); -double Thermo::getPressureZ() { - // Relies on the calculation of the full molecular pressure tensor - - const double p_convert = 1.63882576e8; - double press[3][3]; - double pressureZ; + Globals* simParams = info_->getSimParams(); - this->getPressureTensor(press); + if (simParams->haveTaggedAtomPair() && + simParams->havePrintTaggedPairDistance()) { + if ( simParams->getPrintTaggedPairDistance()) { + + std::pair tap = simParams->getTaggedAtomPair(); + Vector3d pos1, pos2, rab; - pressureZ = p_convert * press[2][2]; +#ifdef IS_MPI + std::cerr << "tap = " << tap.first << " " << tap.second << std::endl; - return pressureZ; -} + int mol1 = info_->getGlobalMolMembership(tap.first); + int mol2 = info_->getGlobalMolMembership(tap.second); + std::cerr << "mols = " << mol1 << " " << mol2 << std::endl; + int proc1 = info_->getMolToProc(mol1); + int proc2 = info_->getMolToProc(mol2); -void Thermo::getPressureTensor(double press[3][3]){ - // returns pressure tensor in units amu*fs^-2*Ang^-1 - // routine derived via viral theorem description in: - // Paci, E. and Marchi, M. J.Phys.Chem. 1996, 100, 4314-4322 + std::cerr << " procs = " << proc1 << " " <getIOIndexToIntegrableObject(tap.first); + std::cerr << " on proc " << proc1 << ", sd1 has global index= " << sd1->getGlobalIndex() << std::endl; + pos1 = sd1->getPos(); + data[0] = pos1.x(); + data[1] = pos1.y(); + data[2] = pos1.z(); + MPI_Bcast(data, 3, MPI_REALTYPE, proc1, MPI_COMM_WORLD); + } else { + MPI_Bcast(data, 3, MPI_REALTYPE, proc1, MPI_COMM_WORLD); + pos1 = Vector3d(data); + } - double molmass, volume; - double vcom[3]; - double p_local[9], p_global[9]; - int i, j, k; - for (i=0; i < 9; i++) { - p_local[i] = 0.0; - p_global[i] = 0.0; - } - - // use velocities of integrableObjects and their masses: - - for (i=0; i < info->integrableObjects.size(); i++) { - - molmass = info->integrableObjects[i]->getMass(); - - info->integrableObjects[i]->getVel(vcom); - - p_local[0] += molmass * (vcom[0] * vcom[0]); - p_local[1] += molmass * (vcom[0] * vcom[1]); - p_local[2] += molmass * (vcom[0] * vcom[2]); - p_local[3] += molmass * (vcom[1] * vcom[0]); - p_local[4] += molmass * (vcom[1] * vcom[1]); - p_local[5] += molmass * (vcom[1] * vcom[2]); - p_local[6] += molmass * (vcom[2] * vcom[0]); - p_local[7] += molmass * (vcom[2] * vcom[1]); - p_local[8] += molmass * (vcom[2] * vcom[2]); - - } - - // Get total for entire system from MPI. - -#ifdef IS_MPI - MPI_Allreduce(p_local,p_global,9,MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); + if (proc2 == worldRank) { + StuntDouble* sd2 = info_->getIOIndexToIntegrableObject(tap.second); + std::cerr << " on proc " << proc2 << ", sd2 has global index= " << sd2->getGlobalIndex() << std::endl; + pos2 = sd2->getPos(); + data[0] = pos2.x(); + data[1] = pos2.y(); + data[2] = pos2.z(); + MPI_Bcast(data, 3, MPI_REALTYPE, proc2, MPI_COMM_WORLD); + } else { + MPI_Bcast(data, 3, MPI_REALTYPE, proc2, MPI_COMM_WORLD); + pos2 = Vector3d(data); + } #else - for (i=0; i<9; i++) { - p_global[i] = p_local[i]; - } -#endif // is_mpi - - volume = this->getVolume(); - - - - for(i = 0; i < 3; i++) { - for (j = 0; j < 3; j++) { - k = 3*i + j; - press[i][j] = (p_global[k] + info->tau[k]*e_convert) / volume; + StuntDouble* at1 = info_->getIOIndexToIntegrableObject(tap.first); + StuntDouble* at2 = info_->getIOIndexToIntegrableObject(tap.second); + pos1 = at1->getPos(); + pos2 = at2->getPos(); +#endif + rab = pos2 - pos1; + currSnapshot->wrapVector(rab); + stat[Stats::TAGGED_PAIR_DISTANCE] = rab.length(); + } } - } -} - -void Thermo::velocitize() { - - double aVel[3], aJ[3], I[3][3]; - int i, j, l, m, n, vr, vd; // velocity randomizer loop counters - double vdrift[3]; - double vbar; - const double kb = 8.31451e-7; // kb in amu, angstroms, fs, etc. - double av2; - double kebar; - double temperature; - int nobj; - - if (!info->have_target_temp) { - sprintf( painCave.errMsg, - "You can't resample the velocities without a targetTemp!\n" - ); - painCave.isFatal = 1; - painCave.severity = OOPSE_ERROR; - simError(); - return; - } - - nobj = info->integrableObjects.size(); - - temperature = info->target_temp; - - kebar = kb * temperature * (double)info->ndfRaw / - ( 2.0 * (double)info->ndf ); - - for(vr = 0; vr < nobj; vr++){ - - // uses equipartition theory to solve for vbar in angstrom/fs - - av2 = 2.0 * kebar / info->integrableObjects[vr]->getMass(); - vbar = sqrt( av2 ); - - // picks random velocities from a gaussian distribution - // centered on vbar - - for (j=0; j<3; j++) - aVel[j] = vbar * gaussStream->getGaussian(); - - info->integrableObjects[vr]->setVel( aVel ); - - if(info->integrableObjects[vr]->isDirectional()){ - - info->integrableObjects[vr]->getI( I ); - - if (info->integrableObjects[vr]->isLinear()) { - - l= info->integrableObjects[vr]->linearAxis(); - m = (l+1)%3; - n = (l+2)%3; - - aJ[l] = 0.0; - vbar = sqrt( 2.0 * kebar * I[m][m] ); - aJ[m] = vbar * gaussStream->getGaussian(); - vbar = sqrt( 2.0 * kebar * I[n][n] ); - aJ[n] = vbar * gaussStream->getGaussian(); - - } else { - for (j = 0 ; j < 3; j++) { - vbar = sqrt( 2.0 * kebar * I[j][j] ); - aJ[j] = vbar * gaussStream->getGaussian(); - } - } // else isLinear - - info->integrableObjects[vr]->setJ( aJ ); - }//isDirectional - - } - - // Get the Center of Mass drift velocity. - - getCOMVel(vdrift); - - // Corrects for the center of mass drift. - // sums all the momentum and divides by total mass. - - for(vd = 0; vd < nobj; vd++){ + /**@todo need refactorying*/ + //Conserved Quantity is set by integrator and time is set by setTime - info->integrableObjects[vd]->getVel(aVel); - - for (j=0; j < 3; j++) - aVel[j] -= vdrift[j]; - - info->integrableObjects[vd]->setVel( aVel ); } -} - -void Thermo::getCOMVel(double vdrift[3]){ - - double mtot, mtot_local; - double aVel[3], amass; - double vdrift_local[3]; - int vd, j; - int nobj; - - nobj = info->integrableObjects.size(); - - mtot_local = 0.0; - vdrift_local[0] = 0.0; - vdrift_local[1] = 0.0; - vdrift_local[2] = 0.0; - - for(vd = 0; vd < nobj; vd++){ - - amass = info->integrableObjects[vd]->getMass(); - info->integrableObjects[vd]->getVel( aVel ); - - for(j = 0; j < 3; j++) - vdrift_local[j] += aVel[j] * amass; - - mtot_local += amass; - } - -#ifdef IS_MPI - MPI_Allreduce(&mtot_local,&mtot,1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); - MPI_Allreduce(vdrift_local,vdrift,3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); -#else - mtot = mtot_local; - for(vd = 0; vd < 3; vd++) { - vdrift[vd] = vdrift_local[vd]; - } -#endif - - for (vd = 0; vd < 3; vd++) { - vdrift[vd] = vdrift[vd] / mtot; - } - -} - -void Thermo::getCOM(double COM[3]){ - - double mtot, mtot_local; - double aPos[3], amass; - double COM_local[3]; - int i, j; - int nobj; - - mtot_local = 0.0; - COM_local[0] = 0.0; - COM_local[1] = 0.0; - COM_local[2] = 0.0; - - nobj = info->integrableObjects.size(); - for(i = 0; i < nobj; i++){ - - amass = info->integrableObjects[i]->getMass(); - info->integrableObjects[i]->getPos( aPos ); - - for(j = 0; j < 3; j++) - COM_local[j] += aPos[j] * amass; - - mtot_local += amass; - } - -#ifdef IS_MPI - MPI_Allreduce(&mtot_local,&mtot,1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); - MPI_Allreduce(COM_local,COM,3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); -#else - mtot = mtot_local; - for(i = 0; i < 3; i++) { - COM[i] = COM_local[i]; - } -#endif - - for (i = 0; i < 3; i++) { - COM[i] = COM[i] / mtot; - } -} - -void Thermo::removeCOMdrift() { - double vdrift[3], aVel[3]; - int vd, j, nobj; - - nobj = info->integrableObjects.size(); - - // Get the Center of Mass drift velocity. - - getCOMVel(vdrift); - - // Corrects for the center of mass drift. - // sums all the momentum and divides by total mass. - - for(vd = 0; vd < nobj; vd++){ - - info->integrableObjects[vd]->getVel(aVel); - - for (j=0; j < 3; j++) - aVel[j] -= vdrift[j]; - - info->integrableObjects[vd]->setVel( aVel ); - } -} +} //end namespace OpenMD