--- trunk/OOPSE/libmdtools/Thermo.cpp 2003/04/09 04:06:43 483 +++ trunk/OOPSE/libmdtools/Thermo.cpp 2003/04/10 16:22:00 486 @@ -134,6 +134,29 @@ double Thermo::getPressure() { return temperature; } +double Thermo::getEnthalpy() { + + 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 @@ -161,13 +184,13 @@ void Thermo::getPressureTensor(double press[9]){ double vcom[3]; double p_local[9], p_global[9]; double theBox[3]; - double* tau; + //double* tau; int i, nMols; Molecule* molecules; nMols = entry_plug->n_mol; molecules = entry_plug->molecules; - tau = entry_plug->tau; + //tau = entry_plug->tau; // use velocities of molecular centers of mass and molecular masses: for (i=0; i < 9; i++) { @@ -204,7 +227,7 @@ void Thermo::getPressureTensor(double press[9]){ volume = theBox[0] * theBox[1] * theBox[2]; for(i=0; i<9; i++) { - press[i] = (p_global[i] - tau[i]*e_convert) / volume; + press[i] = (p_global[i] - entry_plug->tau[i]*e_convert) / volume; } }