--- trunk/OOPSE/libmdtools/Thermo.cpp 2003/04/08 12:44:18 475 +++ trunk/OOPSE/libmdtools/Thermo.cpp 2003/07/15 14:45:09 608 @@ -33,10 +33,10 @@ double Thermo::getKinetic(){ double Thermo::getKinetic(){ const double e_convert = 4.184E-4; // convert kcal/mol -> (amu A^2)/fs^2 - double vx2, vy2, vz2; - double kinetic, v_sqr; - int kl; - double jx2, jy2, jz2; // the square of the angular momentums + double kinetic; + double amass; + double aVel[3], aJ[3], I[3][3]; + int j, kl; DirectionalAtom *dAtom; @@ -51,24 +51,23 @@ double Thermo::getKinetic(){ kinetic = 0.0; kinetic_global = 0.0; for( kl=0; kl < n_atoms; kl++ ){ + + atoms[kl]->getVel(aVel); + amass = atoms[kl]->getMass(); + + for (j=0; j < 3; j++) + kinetic += amass * aVel[j] * aVel[j]; - vx2 = atoms[kl]->get_vx() * atoms[kl]->get_vx(); - vy2 = atoms[kl]->get_vy() * atoms[kl]->get_vy(); - vz2 = atoms[kl]->get_vz() * atoms[kl]->get_vz(); - - v_sqr = vx2 + vy2 + vz2; - kinetic += atoms[kl]->getMass() * v_sqr; - if( atoms[kl]->isDirectional() ){ dAtom = (DirectionalAtom *)atoms[kl]; + + dAtom->getJ( aJ ); + dAtom->getI( I ); - jx2 = dAtom->getJx() * dAtom->getJx(); - jy2 = dAtom->getJy() * dAtom->getJy(); - jz2 = dAtom->getJz() * dAtom->getJz(); + for (j=0; j<3; j++) + kinetic += aJ[j]*aJ[j] / I[j][j]; - kinetic += (jx2 / dAtom->getIxx()) + (jy2 / dAtom->getIyy()) - + (jz2 / dAtom->getIzz()); } } #ifdef IS_MPI @@ -134,53 +133,106 @@ double Thermo::getPressure(){ return temperature; } -double Thermo::getPressure(){ - // returns pressure in units amu*fs^-2*Ang^-1 +double Thermo::getEnthalpy() { + + const double e_convert = 4.184E-4; // convert kcal/mol -> (amu A^2)/fs^2 + double u, p, v; + double press[3][3]; + + u = this->getTotalE(); + + this->getPressureTensor(press); + p = (press[0][0] + press[1][1] + press[2][2]) / 3.0; + + v = this->getVolume(); + + return (u + (p*v)/e_convert); +} + +double Thermo::getVolume() { + + return entry_plug->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; + + this->getPressureTensor(press); + + pressure = p_convert * (press[0][0] + press[1][1] + press[2][2]) / 3.0; + + return pressure; +} + + +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 - const double convert = 4.184e-4; - double molmass; + const double e_convert = 4.184e-4; + + double molmass, volume; double vcom[3]; - double p_local, p_sum, p_mol, virial; - double theBox[3]; - double* tau; - int i, nMols; + double p_local[9], p_global[9]; + int i, j, k, l, 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: - p_local = 0.0; + 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 += (vcom[0]*vcom[0] + vcom[1]*vcom[1] + vcom[2]*vcom[2]) * molmass; + + 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_sum,1,MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(p_local,p_global,9,MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); #else - p_sum = p_local; + for (i=0; i<9; i++) { + p_global[i] = p_local[i]; + } #endif // is_mpi - virial = tau[0] + tau[4] + tau[8]; - entry_plug->getBox(theBox); + volume = entry_plug->boxVol; - p_mol = (p_sum - virial*convert) / (3.0 * theBox[0] * theBox[1]* theBox[2]); - - return p_mol; + 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; + } + } } void Thermo::velocitize() { double x,y; - double vx, vy, vz; - double jx, jy, jz; - int i, vr, vd; // velocity randomizer loop counters + double aVel[3], aJ[3], I[3][3]; + int i, j, vr, vd; // velocity randomizer loop counters double vdrift[3]; double vbar; const double kb = 8.31451e-7; // kb in amu, angstroms, fs, etc. @@ -214,13 +266,11 @@ void Thermo::velocitize() { // picks random velocities from a gaussian distribution // centered on vbar - vx = vbar * gaussStream->getGaussian(); - vy = vbar * gaussStream->getGaussian(); - vz = vbar * gaussStream->getGaussian(); + for (j=0; j<3; j++) + aVel[j] = vbar * gaussStream->getGaussian(); + + atoms[vr]->setVel( aVel ); - atoms[vr]->set_vx( vx ); - atoms[vr]->set_vy( vy ); - atoms[vr]->set_vz( vz ); } // Get the Center of Mass drift velocity. @@ -232,17 +282,12 @@ void Thermo::velocitize() { for(vd = 0; vd < n_atoms; vd++){ - vx = atoms[vd]->get_vx(); - vy = atoms[vd]->get_vy(); - vz = atoms[vd]->get_vz(); - - vx -= vdrift[0]; - vy -= vdrift[1]; - vz -= vdrift[2]; + atoms[vd]->getVel(aVel); - atoms[vd]->set_vx(vx); - atoms[vd]->set_vy(vy); - atoms[vd]->set_vz(vz); + for (j=0; j < 3; j++) + aVel[j] -= vdrift[j]; + + atoms[vd]->setVel( aVel ); } if( n_oriented ){ @@ -251,19 +296,17 @@ void Thermo::velocitize() { if( atoms[i]->isDirectional() ){ dAtom = (DirectionalAtom *)atoms[i]; + dAtom->getI( I ); + + for (j = 0 ; j < 3; j++) { - vbar = sqrt( 2.0 * kebar * dAtom->getIxx() ); - jx = vbar * gaussStream->getGaussian(); + vbar = sqrt( 2.0 * kebar * I[j][j] ); + aJ[j] = vbar * gaussStream->getGaussian(); - vbar = sqrt( 2.0 * kebar * dAtom->getIyy() ); - jy = vbar * gaussStream->getGaussian(); - - vbar = sqrt( 2.0 * kebar * dAtom->getIzz() ); - jz = vbar * gaussStream->getGaussian(); - - dAtom->setJx( jx ); - dAtom->setJy( jy ); - dAtom->setJz( jz ); + } + + dAtom->setJ( aJ ); + } } } @@ -272,8 +315,9 @@ void Thermo::getCOMVel(double vdrift[3]){ void Thermo::getCOMVel(double vdrift[3]){ double mtot, mtot_local; + double aVel[3], amass; double vdrift_local[3]; - int vd, n_atoms; + int vd, n_atoms, j; Atom** atoms; // We are very careless here with the distinction between n_atoms and n_local @@ -289,11 +333,13 @@ void Thermo::getCOMVel(double vdrift[3]){ for(vd = 0; vd < n_atoms; vd++){ - vdrift_local[0] += atoms[vd]->get_vx() * atoms[vd]->getMass(); - vdrift_local[1] += atoms[vd]->get_vy() * atoms[vd]->getMass(); - vdrift_local[2] += atoms[vd]->get_vz() * atoms[vd]->getMass(); + amass = atoms[vd]->getMass(); + atoms[vd]->getVel( aVel ); + + for(j = 0; j < 3; j++) + vdrift_local[j] += aVel[j] * amass; - mtot_local += atoms[vd]->getMass(); + mtot_local += amass; } #ifdef IS_MPI