--- trunk/OOPSE/libmdtools/Thermo.cpp 2003/03/27 20:12:15 423 +++ trunk/OOPSE/libmdtools/Thermo.cpp 2003/04/10 16:22:00 486 @@ -4,12 +4,12 @@ using namespace std; #ifdef IS_MPI #include -#include #endif //is_mpi #include "Thermo.hpp" #include "SRI.hpp" #include "Integrator.hpp" +#include "simError.h" #ifdef IS_MPI #define __C @@ -72,7 +72,8 @@ double Thermo::getKinetic(){ } } #ifdef IS_MPI - MPI::COMM_WORLD.Allreduce(&kinetic,&kinetic_global,1,MPI_DOUBLE,MPI_SUM); + MPI_Allreduce(&kinetic,&kinetic_global,1,MPI_DOUBLE, + MPI_SUM, MPI_COMM_WORLD); kinetic = kinetic_global; #endif //is_mpi @@ -86,25 +87,33 @@ double Thermo::getPotential(){ double potential_local; double potential; int el, nSRI; - SRI** sris; + Molecule* molecules; - sris = entry_plug->sr_interactions; + molecules = entry_plug->molecules; nSRI = entry_plug->n_SRI; potential_local = 0.0; + potential = 0.0; potential_local += entry_plug->lrPot; for( el=0; eln_mol; el++ ){ - potential_local += entry_plug->molecules[el]->get_potential(); + potential_local += molecules[el].getPotential(); } // Get total potential for entire system from MPI. #ifdef IS_MPI - MPI::COMM_WORLD.Allreduce(&potential_local,&potential,1,MPI_DOUBLE,MPI_SUM); + MPI_Allreduce(&potential_local,&potential,1,MPI_DOUBLE, + MPI_SUM, MPI_COMM_WORLD); #else potential = potential_local; #endif // is_mpi +#ifdef IS_MPI + /* + std::cerr << "node " << worldRank << ": after pot = " << potential << "\n"; + */ +#endif + return potential; } @@ -120,30 +129,106 @@ double Thermo::getTemperature(){ const double kb = 1.9872179E-3; // boltzman's constant in kcal/(mol K) double temperature; - int ndf_local, ndf; - ndf_local = 3 * entry_plug->n_atoms + 3 * entry_plug->n_oriented - - entry_plug->n_constraints; + temperature = ( 2.0 * this->getKinetic() ) / ((double)entry_plug->ndf * kb ); + return temperature; +} -#ifdef IS_MPI - MPI::COMM_WORLD.Allreduce(&ndf_local,&ndf,1,MPI_INT,MPI_SUM); -#else - ndf = ndf_local; -#endif +double Thermo::getEnthalpy() { - ndf = ndf - 3; + const double e_convert = 4.184E-4; // convert kcal/mol -> (amu A^2)/fs^2 + double u, p, v; + double press[9]; + + u = this->getTotalE(); + + this->getPressureTensor(press); + p = (press[0] + press[4] + press[8]) / 3.0; + + v = this->getVolume(); + + return (u + (p*v)/e_convert); +} + +double Thermo::getVolume() { + double theBox[3]; + + entry_plug->getBox(theBox); + return (theBox[0] * theBox[1] * theBox[2]); +} + +double Thermo::getPressure() { + // returns the pressure in units of atm + // Relies on the calculation of the full molecular pressure tensor - temperature = ( 2.0 * this->getKinetic() ) / ( ndf * kb ); - return temperature; + const double p_convert = 1.63882576e8; + double press[9]; + double pressure; + + this->getPressureTensor(press); + + pressure = p_convert * (press[0] + press[4] + press[8]) / 3.0; + + return pressure; } -double Thermo::getPressure(){ -// const double conv_Pa_atm = 9.901E-6; // convert Pa -> atm -// const double conv_internal_Pa = 1.661E-7; //convert amu/(fs^2 A) -> Pa -// const double conv_A_m = 1.0E-10; //convert A -> m +void Thermo::getPressureTensor(double press[9]){ + // 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 - return 0.0; + const double e_convert = 4.184e-4; + + double molmass, volume; + double vcom[3]; + double p_local[9], p_global[9]; + double theBox[3]; + //double* tau; + int i, nMols; + Molecule* molecules; + + nMols = entry_plug->n_mol; + molecules = entry_plug->molecules; + //tau = entry_plug->tau; + + // use velocities of molecular centers of mass and molecular masses: + for (i=0; i < 9; i++) { + p_local[i] = 0.0; + p_global[i] = 0.0; + } + + for (i=0; i < nMols; i++) { + molmass = molecules[i].getCOMvel(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); +#else + for (i=0; i<9; i++) { + p_global[i] = p_local[i]; + } +#endif // is_mpi + + entry_plug->getBox(theBox); + + volume = theBox[0] * theBox[1] * theBox[2]; + + for(i=0; i<9; i++) { + press[i] = (p_global[i] - entry_plug->tau[i]*e_convert) / volume; + } } void Thermo::velocitize() { @@ -157,8 +242,6 @@ void Thermo::velocitize() { const double kb = 8.31451e-7; // kb in amu, angstroms, fs, etc. double av2; double kebar; - int ndf, ndf_local; // number of degrees of freedom - int ndfRaw, ndfRaw_local; // the raw number of degrees of freedom int n_atoms; Atom** atoms; DirectionalAtom* dAtom; @@ -172,31 +255,16 @@ void Thermo::velocitize() { n_oriented = entry_plug->n_oriented; n_constraints = entry_plug->n_constraints; - // Raw degrees of freedom that we have to set - ndfRaw_local = 3 * entry_plug->n_atoms + 3 * entry_plug->n_oriented; - - // Degrees of freedom that can contain kinetic energy - ndf_local = 3 * entry_plug->n_atoms + 3 * entry_plug->n_oriented - - entry_plug->n_constraints; + kebar = kb * temperature * (double)entry_plug->ndf / + ( 2.0 * (double)entry_plug->ndfRaw ); -#ifdef IS_MPI - MPI::COMM_WORLD.Allreduce(&ndf_local,&ndf,1,MPI_INT,MPI_SUM); - MPI::COMM_WORLD.Allreduce(&ndfRaw_local,&ndfRaw,1,MPI_INT,MPI_SUM); -#else - ndfRaw = ndfRaw_local; - ndf = ndf_local; -#endif - ndf = ndf - 3; - - kebar = kb * temperature * (double)ndf / ( 2.0 * (double)ndfRaw ); - for(vr = 0; vr < n_atoms; vr++){ // uses equipartition theory to solve for vbar in angstrom/fs av2 = 2.0 * kebar / atoms[vr]->getMass(); vbar = sqrt( av2 ); - + // vbar = sqrt( 8.31451e-7 * temperature / atoms[vr]->getMass() ); // picks random velocities from a gaussian distribution @@ -245,7 +313,7 @@ void Thermo::velocitize() { vbar = sqrt( 2.0 * kebar * dAtom->getIyy() ); jy = vbar * gaussStream->getGaussian(); - + vbar = sqrt( 2.0 * kebar * dAtom->getIzz() ); jz = vbar * gaussStream->getGaussian(); @@ -285,8 +353,8 @@ void Thermo::getCOMVel(double vdrift[3]){ } #ifdef IS_MPI - MPI::COMM_WORLD.Allreduce(&mtot_local,&mtot,1,MPI_DOUBLE,MPI_SUM); - MPI::COMM_WORLD.Allreduce(vdrift_local,vdrift,3,MPI_DOUBLE,MPI_SUM); + 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++) {