--- trunk/src/brains/Thermo.cpp 2004/09/24 04:16:43 2 +++ branches/development/src/brains/Thermo.cpp 2012/06/05 17:48:40 1733 @@ -1,450 +1,569 @@ +/* + * 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] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010). + * [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011). + */ + #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" +#include "types/MultipoleAdapter.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::getElectronicTemperature() { + SimInfo::MoleculeIterator miter; + std::vector::iterator iiter; + Molecule* mol; + Atom* atom; + RealType cvel; + RealType cmass; + RealType kinetic = 0.0; + RealType kinetic_global = 0.0; + + for (mol = info_->beginMolecule(miter); mol != NULL; mol = info_->nextMolecule(miter)) { + for (atom = mol->beginFluctuatingCharge(iiter); atom != NULL; + atom = mol->nextFluctuatingCharge(iiter)) { + cmass = atom->getChargeMass(); + cvel = atom->getFlucQVel(); + + kinetic += cmass * cvel * cvel; + + } + } + +#ifdef IS_MPI - total = this->getKinetic() + this->getPotential(); - return total; -} + MPI_Allreduce(&kinetic, &kinetic_global, 1, MPI_REALTYPE, MPI_SUM, + MPI_COMM_WORLD); + kinetic = kinetic_global; -double Thermo::getTemperature(){ +#endif //is_mpi - const double kb = 1.9872156E-3; // boltzman's constant in kcal/(mol K) - double temperature; + kinetic = kinetic * 0.5; + return ( 2.0 * kinetic) / (info_->getNFluctuatingCharges()* PhysicalConstants::kb ); + } - temperature = ( 2.0 * this->getKinetic() ) / ((double)info->ndf * kb ); - return temperature; -} -double Thermo::getVolume() { - return info->boxVol; -} -double Thermo::getPressure() { + RealType Thermo::getVolume() { + Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot(); + return curSnapshot->getVolume(); + } - // 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() { - 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; - return pressure; -} + Mat3x3d tensor; + RealType pressure; -double Thermo::getPressureX() { + tensor = getPressureTensor(); - // Relies on the calculation of the full molecular pressure tensor - - const double p_convert = 1.63882576e8; - double press[3][3]; - double pressureX; + pressure = PhysicalConstants::pressureConvert * (tensor(0, 0) + tensor(1, 1) + tensor(2, 2)) / 3.0; - this->getPressureTensor(press); + return pressure; + } - pressureX = p_convert * press[0][0]; + RealType Thermo::getPressure(int direction) { - return pressureX; -} + // Relies on the calculation of the full molecular pressure tensor -double Thermo::getPressureY() { + + Mat3x3d tensor; + RealType pressure; - // Relies on the calculation of the full molecular pressure tensor - - const double p_convert = 1.63882576e8; - double press[3][3]; - double pressureY; + tensor = getPressureTensor(); - this->getPressureTensor(press); + pressure = PhysicalConstants::pressureConvert * tensor(direction, direction); - pressureY = p_convert * press[1][1]; + return pressure; + } - return pressureY; -} + 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); -double Thermo::getPressureZ() { + 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)) { - // Relies on the calculation of the full molecular pressure tensor - - const double p_convert = 1.63882576e8; - double press[3][3]; - double pressureZ; + 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 - this->getPressureTensor(press); + RealType volume = this->getVolume(); + Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot(); + Mat3x3d stressTensor = curSnapshot->getStressTensor(); - pressureZ = p_convert * press[2][2]; + pressureTensor = (p_global + + PhysicalConstants::energyConvert * stressTensor)/volume; + + return pressureTensor; + } - return pressureZ; -} + 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(); -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 + 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); - const double e_convert = 4.184e-4; + // grab the simulation box dipole moment if specified + if (info_->getCalcBoxDipole()){ + Vector3d totalDipole = getBoxDipole(); + stat[Stats::BOX_DIPOLE_X] = totalDipole(0); + stat[Stats::BOX_DIPOLE_Y] = totalDipole(1); + stat[Stats::BOX_DIPOLE_Z] = totalDipole(2); + } - double molmass, volume; - double vcom[3]; - double p_local[9], p_global[9]; - int i, j, k; + Globals* simParams = info_->getSimParams(); + // grab the heat flux if desired + if (simParams->havePrintHeatFlux()) { + if (simParams->getPrintHeatFlux()){ + Vector3d heatFlux = getHeatFlux(); + stat[Stats::HEATFLUX_X] = heatFlux(0); + stat[Stats::HEATFLUX_Y] = heatFlux(1); + stat[Stats::HEATFLUX_Z] = heatFlux(2); + } + } - for (i=0; i < 9; i++) { - p_local[i] = 0.0; - p_global[i] = 0.0; - } + if (simParams->haveTaggedAtomPair() && + simParams->havePrintTaggedPairDistance()) { + if ( simParams->getPrintTaggedPairDistance()) { + + std::pair tap = simParams->getTaggedAtomPair(); + Vector3d pos1, pos2, rab; - // use velocities of integrableObjects and their masses: +#ifdef IS_MPI + std::cerr << "tap = " << tap.first << " " << tap.second << std::endl; - for (i=0; i < info->integrableObjects.size(); i++) { + int mol1 = info_->getGlobalMolMembership(tap.first); + int mol2 = info_->getGlobalMolMembership(tap.second); + std::cerr << "mols = " << mol1 << " " << mol2 << std::endl; - 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]); + int proc1 = info_->getMolToProc(mol1); + int proc2 = info_->getMolToProc(mol2); - } + std::cerr << " procs = " << proc1 << " " <getVolume(); + RealType data[3]; + if (proc1 == worldRank) { + StuntDouble* sd1 = info_->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); + } - - 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; + 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 + 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(); + } } + + /**@todo need refactorying*/ + //Conserved Quantity is set by integrator and time is set by setTime + } -} -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; - } + Vector3d Thermo::getBoxDipole() { + Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot(); + SimInfo::MoleculeIterator miter; + std::vector::iterator aiter; + Molecule* mol; + Atom* atom; + RealType charge; + RealType moment(0.0); + Vector3d ri(0.0); + Vector3d dipoleVector(0.0); + Vector3d nPos(0.0); + Vector3d pPos(0.0); + RealType nChg(0.0); + RealType pChg(0.0); + int nCount = 0; + int pCount = 0; - 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++){ + RealType chargeToC = 1.60217733e-19; + RealType angstromToM = 1.0e-10; + RealType debyeToCm = 3.33564095198e-30; - // uses equipartition theory to solve for vbar in angstrom/fs + for (mol = info_->beginMolecule(miter); mol != NULL; + mol = info_->nextMolecule(miter)) { - av2 = 2.0 * kebar / info->integrableObjects[vr]->getMass(); - vbar = sqrt( av2 ); + for (atom = mol->beginAtom(aiter); atom != NULL; + atom = mol->nextAtom(aiter)) { + + if (atom->isCharge() ) { + charge = 0.0; + GenericData* data = atom->getAtomType()->getPropertyByName("Charge"); + if (data != NULL) { - // picks random velocities from a gaussian distribution - // centered on vbar + charge = (dynamic_cast(data))->getData(); + charge *= chargeToC; - for (j=0; j<3; j++) - aVel[j] = vbar * gaussStream->getGaussian(); + ri = atom->getPos(); + currSnapshot->wrapVector(ri); + ri *= angstromToM; + + if (charge < 0.0) { + nPos += ri; + nChg -= charge; + nCount++; + } else if (charge > 0.0) { + pPos += ri; + pChg += charge; + pCount++; + } + } + } + + MultipoleAdapter ma = MultipoleAdapter(atom->getAtomType()); + if (ma.isDipole() ) { + Vector3d u_i = atom->getElectroFrame().getColumn(2); + moment = ma.getDipoleMoment(); + moment *= debyeToCm; + dipoleVector += u_i * moment; + } + } + } - info->integrableObjects[vr]->setVel( aVel ); - - if(info->integrableObjects[vr]->isDirectional()){ + +#ifdef IS_MPI + RealType pChg_global, nChg_global; + int pCount_global, nCount_global; + Vector3d pPos_global, nPos_global, dipVec_global; - info->integrableObjects[vr]->getI( I ); + MPI_Allreduce(&pChg, &pChg_global, 1, MPI_REALTYPE, MPI_SUM, + MPI_COMM_WORLD); + pChg = pChg_global; + MPI_Allreduce(&nChg, &nChg_global, 1, MPI_REALTYPE, MPI_SUM, + MPI_COMM_WORLD); + nChg = nChg_global; + MPI_Allreduce(&pCount, &pCount_global, 1, MPI_INTEGER, MPI_SUM, + MPI_COMM_WORLD); + pCount = pCount_global; + MPI_Allreduce(&nCount, &nCount_global, 1, MPI_INTEGER, MPI_SUM, + MPI_COMM_WORLD); + nCount = nCount_global; + MPI_Allreduce(pPos.getArrayPointer(), pPos_global.getArrayPointer(), 3, + MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD); + pPos = pPos_global; + MPI_Allreduce(nPos.getArrayPointer(), nPos_global.getArrayPointer(), 3, + MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD); + nPos = nPos_global; + MPI_Allreduce(dipoleVector.getArrayPointer(), + dipVec_global.getArrayPointer(), 3, + MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD); + dipoleVector = dipVec_global; +#endif //is_mpi - if (info->integrableObjects[vr]->isLinear()) { + // first load the accumulated dipole moment (if dipoles were present) + Vector3d boxDipole = dipoleVector; + // now include the dipole moment due to charges + // use the lesser of the positive and negative charge totals + RealType chg_value = nChg <= pChg ? nChg : pChg; + + // find the average positions + if (pCount > 0 && nCount > 0 ) { + pPos /= pCount; + nPos /= nCount; + } - l= info->integrableObjects[vr]->linearAxis(); - m = (l+1)%3; - n = (l+2)%3; + // dipole is from the negative to the positive (physics notation) + boxDipole += (pPos - nPos) * chg_value; - 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 - + return boxDipole; } - // Get the Center of Mass drift velocity. + // Returns the Heat Flux Vector for the system + Vector3d Thermo::getHeatFlux(){ + Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot(); + SimInfo::MoleculeIterator miter; + std::vector::iterator iiter; + Molecule* mol; + StuntDouble* integrableObject; + RigidBody::AtomIterator ai; + Atom* atom; + Vector3d vel; + Vector3d angMom; + Mat3x3d I; + int i; + int j; + int k; + RealType mass; - getCOMVel(vdrift); - - // Corrects for the center of mass drift. - // sums all the momentum and divides by total mass. + Vector3d x_a; + RealType kinetic; + RealType potential; + RealType eatom; + RealType AvgE_a_ = 0; + // Convective portion of the heat flux + Vector3d heatFluxJc = V3Zero; - for(vd = 0; vd < nobj; vd++){ - - info->integrableObjects[vd]->getVel(aVel); - - for (j=0; j < 3; j++) - aVel[j] -= vdrift[j]; + /* Calculate convective portion of the heat flux */ + for (mol = info_->beginMolecule(miter); mol != NULL; + mol = info_->nextMolecule(miter)) { + + for (integrableObject = mol->beginIntegrableObject(iiter); + integrableObject != NULL; + integrableObject = mol->nextIntegrableObject(iiter)) { - info->integrableObjects[vd]->setVel( aVel ); - } + 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(); -void Thermo::getCOMVel(double vdrift[3]){ + 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); + } + } - double mtot, mtot_local; - double aVel[3], amass; - double vdrift_local[3]; - int vd, j; - int nobj; + potential = 0.0; - nobj = info->integrableObjects.size(); + if (integrableObject->isRigidBody()) { + RigidBody* rb = dynamic_cast(integrableObject); + for (atom = rb->beginAtom(ai); atom != NULL; + atom = rb->nextAtom(ai)) { + potential += atom->getParticlePot(); + } + } else { + potential = integrableObject->getParticlePot(); + cerr << "ppot = " << potential << "\n"; + } - 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 ); + potential *= PhysicalConstants::energyConvert; // amu A^2/fs^2 + // The potential may not be a 1/2 factor + eatom = (kinetic + potential)/2.0; // amu A^2/fs^2 + heatFluxJc[0] += eatom*vel[0]; // amu A^3/fs^3 + heatFluxJc[1] += eatom*vel[1]; // amu A^3/fs^3 + heatFluxJc[2] += eatom*vel[2]; // amu A^3/fs^3 + } + } - for(j = 0; j < 3; j++) - vdrift_local[j] += aVel[j] * amass; - - mtot_local += amass; - } + std::cerr << "Heat flux heatFluxJc is: " << heatFluxJc << std::endl; + /* The J_v vector is reduced in fortan so everyone has the global + * Jv. Jc is computed over the local atoms and must be reduced + * among all processors. + */ #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]; - } + MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &heatFluxJc[0], 3, MPI::REALTYPE, + MPI::SUM); #endif - for (vd = 0; vd < 3; vd++) { - vdrift[vd] = vdrift[vd] / mtot; - } - -} + // (kcal/mol * A/fs) * conversion => (amu A^3)/fs^3 -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++){ + Vector3d heatFluxJv = currSnapshot->getConductiveHeatFlux() * + PhysicalConstants::energyConvert; - amass = info->integrableObjects[i]->getMass(); - info->integrableObjects[i]->getPos( aPos ); - - for(j = 0; j < 3; j++) - COM_local[j] += aPos[j] * amass; + std::cerr << "Heat flux Jc is: " << heatFluxJc << std::endl; + std::cerr << "Heat flux Jv is: " << heatFluxJv << std::endl; - mtot_local += amass; + // Correct for the fact the flux is 1/V (Jc + Jv) + return (heatFluxJv + heatFluxJc) / this->getVolume(); // amu / fs^3 } - -#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