--- trunk/OOPSE/libmdtools/Thermo.cpp 2003/07/15 14:45:09 608 +++ trunk/OOPSE/libmdtools/Thermo.cpp 2003/11/06 19:11:38 853 @@ -1,4 +1,4 @@ -#include +#include #include using namespace std; @@ -16,12 +16,9 @@ using namespace std; #include "mpiSimulation.hpp" #endif // is_mpi - -#define BASE_SEED 123456789 - -Thermo::Thermo( SimInfo* the_entry_plug ) { - entry_plug = the_entry_plug; - int baseSeed = BASE_SEED; +Thermo::Thermo( SimInfo* the_info ) { + info = the_info; + int baseSeed = the_info->getSeed(); gaussStream = new gaussianSPRNG( baseSeed ); } @@ -45,8 +42,8 @@ double Thermo::getKinetic(){ Atom** atoms; - n_atoms = entry_plug->n_atoms; - atoms = entry_plug->atoms; + n_atoms = info->n_atoms; + atoms = info->atoms; kinetic = 0.0; kinetic_global = 0.0; @@ -88,14 +85,14 @@ double Thermo::getPotential(){ int el, nSRI; Molecule* molecules; - molecules = entry_plug->molecules; - nSRI = entry_plug->n_SRI; + molecules = info->molecules; + nSRI = info->n_SRI; potential_local = 0.0; potential = 0.0; - potential_local += entry_plug->lrPot; + potential_local += info->lrPot; - for( el=0; eln_mol; el++ ){ + for( el=0; eln_mol; el++ ){ potential_local += molecules[el].getPotential(); } @@ -107,12 +104,6 @@ double Thermo::getPotential(){ potential = potential_local; #endif // is_mpi -#ifdef IS_MPI - /* - std::cerr << "node " << worldRank << ": after pot = " << potential << "\n"; - */ -#endif - return potential; } @@ -126,47 +117,76 @@ double Thermo::getTemperature(){ double Thermo::getTemperature(){ - const double kb = 1.9872179E-3; // boltzman's constant in kcal/(mol K) + const double kb = 1.9872156E-3; // boltzman's constant in kcal/(mol K) double temperature; - temperature = ( 2.0 * this->getKinetic() ) / ((double)entry_plug->ndf * kb ); + temperature = ( 2.0 * this->getKinetic() ) / ((double)info->ndf * kb ); return temperature; } -double Thermo::getEnthalpy() { +double Thermo::getVolume() { - const double e_convert = 4.184E-4; // convert kcal/mol -> (amu A^2)/fs^2 - double u, p, v; + return info->boxVol; +} + +double Thermo::getPressure() { + + // Relies on the calculation of the full molecular pressure tensor + + const double p_convert = 1.63882576e8; double press[3][3]; + double pressure; - u = this->getTotalE(); + this->getPressureTensor(press); + pressure = p_convert * (press[0][0] + press[1][1] + press[2][2]) / 3.0; + + return pressure; +} + +double Thermo::getPressureX() { + + // Relies on the calculation of the full molecular pressure tensor + + const double p_convert = 1.63882576e8; + double press[3][3]; + double pressureX; + this->getPressureTensor(press); - p = (press[0][0] + press[1][1] + press[2][2]) / 3.0; - v = this->getVolume(); + pressureX = p_convert * press[0][0]; - return (u + (p*v)/e_convert); + return pressureX; } -double Thermo::getVolume() { +double Thermo::getPressureY() { - return entry_plug->boxVol; + // Relies on the calculation of the full molecular pressure tensor + + const double p_convert = 1.63882576e8; + double press[3][3]; + double pressureY; + + this->getPressureTensor(press); + + pressureY = p_convert * press[1][1]; + + return pressureY; } -double Thermo::getPressure() { +double Thermo::getPressureZ() { // Relies on the calculation of the full molecular pressure tensor const double p_convert = 1.63882576e8; double press[3][3]; - double pressure; + double pressureZ; this->getPressureTensor(press); - pressure = p_convert * (press[0][0] + press[1][1] + press[2][2]) / 3.0; + pressureZ = p_convert * press[2][2]; - return pressure; + return pressureZ; } @@ -180,12 +200,12 @@ void Thermo::getPressureTensor(double press[3][3]){ double molmass, volume; double vcom[3]; double p_local[9], p_global[9]; - int i, j, k, l, nMols; + int i, j, k, nMols; Molecule* molecules; - nMols = entry_plug->n_mol; - molecules = entry_plug->molecules; - //tau = entry_plug->tau; + nMols = info->n_mol; + molecules = info->molecules; + //tau = info->tau; // use velocities of molecular centers of mass and molecular masses: for (i=0; i < 9; i++) { @@ -217,20 +237,19 @@ void Thermo::getPressureTensor(double press[3][3]){ } #endif // is_mpi - volume = entry_plug->boxVol; + volume = this->getVolume(); for(i = 0; i < 3; i++) { for (j = 0; j < 3; j++) { k = 3*i + j; - l = 3*j + i; - press[i][j] = (p_global[k] - entry_plug->tau[l]*e_convert) / volume; + press[i][j] = (p_global[k] + info->tau[k]*e_convert) / volume; + } } } void Thermo::velocitize() { - double x,y; double aVel[3], aJ[3], I[3][3]; int i, j, vr, vd; // velocity randomizer loop counters double vdrift[3]; @@ -245,14 +264,14 @@ void Thermo::velocitize() { int n_oriented; int n_constraints; - atoms = entry_plug->atoms; - n_atoms = entry_plug->n_atoms; - temperature = entry_plug->target_temp; - n_oriented = entry_plug->n_oriented; - n_constraints = entry_plug->n_constraints; + atoms = info->atoms; + n_atoms = info->n_atoms; + temperature = info->target_temp; + n_oriented = info->n_oriented; + n_constraints = info->n_constraints; - kebar = kb * temperature * (double)entry_plug->ndf / - ( 2.0 * (double)entry_plug->ndfRaw ); + kebar = kb * temperature * (double)info->ndfRaw / + ( 2.0 * (double)info->ndf ); for(vr = 0; vr < n_atoms; vr++){ @@ -260,9 +279,7 @@ void Thermo::velocitize() { 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 // centered on vbar @@ -323,8 +340,8 @@ void Thermo::getCOMVel(double vdrift[3]){ // We are very careless here with the distinction between n_atoms and n_local // We should really fix this before someone pokes an eye out. - n_atoms = entry_plug->n_atoms; - atoms = entry_plug->atoms; + n_atoms = info->n_atoms; + atoms = info->atoms; mtot_local = 0.0; vdrift_local[0] = 0.0; @@ -358,3 +375,47 @@ void Thermo::getCOMVel(double vdrift[3]){ } +void Thermo::getCOM(double COM[3]){ + + double mtot, mtot_local; + double aPos[3], amass; + double COM_local[3]; + int i, n_atoms, j; + Atom** atoms; + + // We are very careless here with the distinction between n_atoms and n_local + // We should really fix this before someone pokes an eye out. + + n_atoms = info->n_atoms; + atoms = info->atoms; + + mtot_local = 0.0; + COM_local[0] = 0.0; + COM_local[1] = 0.0; + COM_local[2] = 0.0; + + for(i = 0; i < n_atoms; i++){ + + amass = atoms[i]->getMass(); + atoms[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; + } +}