--- trunk/OOPSE/libmdtools/Thermo.cpp 2003/09/25 19:27:15 787 +++ trunk/OOPSE/libmdtools/Thermo.cpp 2004/04/22 21:33:55 1131 @@ -1,4 +1,4 @@ -#include +#include #include using namespace std; \ No newline at end of file @@ -10,6 +10,7 @@ using namespace std; #include "SRI.hpp" #include "Integrator.hpp" #include "simError.h" +#include "MatVec3.h" #ifdef IS_MPI #define __C \ No newline at end of file @@ -33,46 +34,43 @@ double Thermo::getKinetic(){ double kinetic; double amass; double aVel[3], aJ[3], I[3][3]; - int j, kl; + int i, j, k, kl; - DirectionalAtom *dAtom; - - int n_atoms; double kinetic_global; - Atom** atoms; - + vector integrableObjects = info->integrableObjects; - n_atoms = info->n_atoms; - atoms = info->atoms; - 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]; - if( atoms[kl]->isDirectional() ){ - - dAtom = (DirectionalAtom *)atoms[kl]; + for (kl=0; klgetVel(aVel); + amass = integrableObjects[kl]->getMass(); - dAtom->getJ( aJ ); - dAtom->getI( I ); - - for (j=0; j<3; j++) - kinetic += aJ[j]*aJ[j] / I[j][j]; - - } + for(j=0; j<3; j++) + kinetic += amass*aVel[j]*aVel[j]; + + 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; \ No newline at end of file @@ -104,12 +102,6 @@ double Thermo::getPotential(){ potential = potential_local; #endif // is_mpi -#ifdef IS_MPI - /* - std::cerr << "node " << worldRank << ": after pot = " << potential << "\n"; - */ -#endif - return potential; } \ No newline at end of file @@ -125,27 +117,11 @@ double Thermo::getTemperature(){ const double kb = 1.9872156E-3; // boltzman's constant in kcal/(mol K) double temperature; - + temperature = ( 2.0 * this->getKinetic() ) / ((double)info->ndf * kb ); 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[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 info->boxVol; \ No newline at end of file @@ -220,7 +196,7 @@ void Thermo::getPressureTensor(double press[3][3]){ const double e_convert = 4.184e-4; double molmass, volume; - double vcom[3]; + double vcom[3], pcom[3], fcom[3], scaled[3]; double p_local[9], p_global[9]; int i, j, k, nMols; Molecule* molecules; \ No newline at end of file @@ -235,18 +211,33 @@ void Thermo::getPressureTensor(double press[3][3]){ p_global[i] = 0.0; } - for (i=0; i < nMols; i++) { - molmass = molecules[i].getCOMvel(vcom); + for (i=0; i < info->integrableObjects.size(); i++) { - 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]); + molmass = info->integrableObjects[i]->getMass(); + + info->integrableObjects[i]->getVel(vcom); + info->integrableObjects[i]->getPos(pcom); + info->integrableObjects[i]->getFrc(fcom); + + matVecMul3(info->HmatInv, pcom, scaled); + + for(j=0; j<3; j++) + scaled[j] -= roundMe(scaled[j]); + + // calc the wrapped real coordinates from the wrapped scaled coordinates + + matVecMul3(info->Hmat, scaled, pcom); + + p_local[0] += molmass * (vcom[0] * vcom[0]) + fcom[0]*pcom[0]*eConvert; + p_local[1] += molmass * (vcom[0] * vcom[1]) + fcom[0]*pcom[1]*eConvert; + p_local[2] += molmass * (vcom[0] * vcom[2]) + fcom[0]*pcom[2]*eConvert; + p_local[3] += molmass * (vcom[1] * vcom[0]) + fcom[1]*pcom[0]*eConvert; + p_local[4] += molmass * (vcom[1] * vcom[1]) + fcom[1]*pcom[1]*eConvert; + p_local[5] += molmass * (vcom[1] * vcom[2]) + fcom[1]*pcom[2]*eConvert; + p_local[6] += molmass * (vcom[2] * vcom[0]) + fcom[2]*pcom[0]*eConvert; + p_local[7] += molmass * (vcom[2] * vcom[1]) + fcom[2]*pcom[1]*eConvert; + p_local[8] += molmass * (vcom[2] * vcom[2]) + fcom[2]*pcom[2]*eConvert; + } // Get total for entire system from MPI. \ No newline at end of file @@ -264,7 +255,7 @@ void Thermo::getPressureTensor(double press[3][3]){ 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; + press[i][j] = p_global[k] / volume; } } \ No newline at end of file @@ -273,45 +264,64 @@ void Thermo::velocitize() { void Thermo::velocitize() { double aVel[3], aJ[3], I[3][3]; - int i, j, vr, vd; // velocity randomizer loop counters + 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; - int n_atoms; - Atom** atoms; - DirectionalAtom* dAtom; double temperature; - int n_oriented; - int n_constraints; + int nobj; - atoms = info->atoms; - n_atoms = info->n_atoms; + nobj = info->integrableObjects.size(); + temperature = info->target_temp; - n_oriented = info->n_oriented; - n_constraints = info->n_constraints; kebar = kb * temperature * (double)info->ndfRaw / ( 2.0 * (double)info->ndf ); - for(vr = 0; vr < n_atoms; vr++){ + for(vr = 0; vr < nobj; vr++){ // uses equipartition theory to solve for vbar in angstrom/fs - av2 = 2.0 * kebar / atoms[vr]->getMass(); + av2 = 2.0 * kebar / info->integrableObjects[vr]->getMass(); vbar = sqrt( av2 ); - -// vbar = sqrt( 8.31451e-7 * temperature / atoms[vr]->getMass() ); - + // picks random velocities from a gaussian distribution // centered on vbar for (j=0; j<3; j++) aVel[j] = vbar * gaussStream->getGaussian(); - atoms[vr]->setVel( aVel ); + info->integrableObjects[vr]->setVel( aVel ); + + if(info->integrableObjects[vr]->isDirectional()){ + info->integrableObjects[vr]->getI( I ); + + if (info->integrableObjects[vr]->isLinear()) { + + l= info->integrableObjects[vr]->linearAxis(); + m = (l+1)%3; + n = (l+2)%3; + + 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 + } // Get the Center of Mass drift velocity. \ No newline at end of file @@ -321,36 +331,16 @@ void Thermo::velocitize() { // Corrects for the center of mass drift. // sums all the momentum and divides by total mass. - for(vd = 0; vd < n_atoms; vd++){ + for(vd = 0; vd < nobj; vd++){ - atoms[vd]->getVel(aVel); + info->integrableObjects[vd]->getVel(aVel); for (j=0; j < 3; j++) aVel[j] -= vdrift[j]; - atoms[vd]->setVel( aVel ); + info->integrableObjects[vd]->setVel( aVel ); } - if( n_oriented ){ - - for( i=0; iisDirectional() ){ - - dAtom = (DirectionalAtom *)atoms[i]; - dAtom->getI( I ); - - for (j = 0 ; j < 3; j++) { - vbar = sqrt( 2.0 * kebar * I[j][j] ); - aJ[j] = vbar * gaussStream->getGaussian(); - - } - - dAtom->setJ( aJ ); - - } - } - } } void Thermo::getCOMVel(double vdrift[3]){ \ No newline at end of file @@ -358,24 +348,20 @@ void Thermo::getCOMVel(double vdrift[3]){ double mtot, mtot_local; double aVel[3], amass; double vdrift_local[3]; - int vd, n_atoms, j; - Atom** atoms; + int vd, j; + int nobj; - // 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. + nobj = info->integrableObjects.size(); - n_atoms = info->n_atoms; - atoms = info->atoms; - mtot_local = 0.0; vdrift_local[0] = 0.0; vdrift_local[1] = 0.0; vdrift_local[2] = 0.0; - for(vd = 0; vd < n_atoms; vd++){ + for(vd = 0; vd < nobj; vd++){ - amass = atoms[vd]->getMass(); - atoms[vd]->getVel( aVel ); + amass = info->integrableObjects[vd]->getMass(); + info->integrableObjects[vd]->getVel( aVel ); for(j = 0; j < 3; j++) vdrift_local[j] += aVel[j] * amass; \ No newline at end of file @@ -404,24 +390,19 @@ 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; + int i, j; + int nobj; - // 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++){ + + nobj = info->integrableObjects.size(); + for(i = 0; i < nobj; i++){ - amass = atoms[i]->getMass(); - atoms[i]->getPos( aPos ); + amass = info->integrableObjects[i]->getMass(); + info->integrableObjects[i]->getPos( aPos ); for(j = 0; j < 3; j++) COM_local[j] += aPos[j] * amass; \ No newline at end of file @@ -443,3 +424,27 @@ void Thermo::getCOM(double COM[3]){ 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 ); + } +} \ No newline at end of file