--- trunk/OOPSE/libmdtools/Thermo.cpp 2003/04/03 20:21:54 447 +++ trunk/OOPSE/libmdtools/Thermo.cpp 2004/08/23 15:11:36 1452 @@ -1,4 +1,4 @@ -#include +#include #include using namespace std; @@ -10,73 +10,76 @@ using namespace std; #include "SRI.hpp" #include "Integrator.hpp" #include "simError.h" +#include "MatVec3.h" +#include "ConstraintManager.hpp" +#include "Mat3x3d.hpp" #ifdef IS_MPI #define __C #include "mpiSimulation.hpp" #endif // is_mpi +inline double roundMe( double x ){ + return ( x >= 0 ) ? floor( x + 0.5 ) : ceil( x - 0.5 ); +} -#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 ); + + cpIter = info->consMan->createPairIterator(); } Thermo::~Thermo(){ delete gaussStream; + delete cpIter; } 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 i, j, k, kl; - DirectionalAtom *dAtom; - - int n_atoms; double kinetic_global; - Atom** atoms; - + vector integrableObjects = info->integrableObjects; - n_atoms = entry_plug->n_atoms; - atoms = entry_plug->atoms; - kinetic = 0.0; kinetic_global = 0.0; - for( kl=0; kl < n_atoms; kl++ ){ - 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(); + for (kl=0; klgetVel(aVel); + amass = integrableObjects[kl]->getMass(); - v_sqr = vx2 + vy2 + vz2; - kinetic += atoms[kl]->getMass() * v_sqr; + for(j=0; j<3; j++) + kinetic += amass*aVel[j]*aVel[j]; - if( atoms[kl]->isDirectional() ){ - - dAtom = (DirectionalAtom *)atoms[kl]; - - jx2 = dAtom->getJx() * dAtom->getJx(); - jy2 = dAtom->getJy() * dAtom->getJy(); - jz2 = dAtom->getJz() * dAtom->getJz(); - - kinetic += (jx2 / dAtom->getIxx()) + (jy2 / dAtom->getIyy()) - + (jz2 / dAtom->getIzz()); - } + 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; @@ -89,14 +92,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(); } @@ -108,12 +111,6 @@ double Thermo::getPotential(){ potential = potential_local; #endif // is_mpi -#ifdef IS_MPI - /* - std::cerr << "node " << worldRank << ": after pot = " << potential << "\n"; - */ -#endif - return potential; } @@ -127,96 +124,199 @@ 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; - int ndf_local, ndf; + + temperature = ( 2.0 * this->getKinetic() ) / ((double)info->ndf * kb ); + return temperature; +} + +double Thermo::getVolume() { + + return info->boxVol; +} + +double Thermo::getPressure() { + + // Relies on the calculation of the full molecular pressure tensor - ndf_local = 3 * entry_plug->n_atoms + 3 * entry_plug->n_oriented - - entry_plug->n_constraints; + const double p_convert = 1.63882576e8; + double press[3][3]; + double pressure; -#ifdef IS_MPI - MPI_Allreduce(&ndf_local,&ndf,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD); -#else - ndf = ndf_local; -#endif + this->getPressureTensor(press); - ndf = ndf - 3; + 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 - temperature = ( 2.0 * this->getKinetic() ) / ( ndf * kb ); - return temperature; + const double p_convert = 1.63882576e8; + double press[3][3]; + double pressureX; + + this->getPressureTensor(press); + + pressureX = p_convert * press[0][0]; + + return pressureX; } -double Thermo::getPressure(){ - // returns pressure in units amu*fs^-2*Ang^-1 +double Thermo::getPressureY() { + + // 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::getPressureZ() { + + // Relies on the calculation of the full molecular pressure tensor + + const double p_convert = 1.63882576e8; + double press[3][3]; + double pressureZ; + + this->getPressureTensor(press); + + pressureZ = p_convert * press[2][2]; + + return pressureZ; +} + + +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 - return 0.0; + const double e_convert = 4.184e-4; + + double molmass, volume; + double vcom[3]; + double p_local[9], p_global[9]; + int i, j, k; + + for (i=0; i < 9; i++) { + p_local[i] = 0.0; + p_global[i] = 0.0; + } + + // use velocities of integrableObjects and their masses: + + for (i=0; i < info->integrableObjects.size(); i++) { + + molmass = info->integrableObjects[i]->getMass(); + + info->integrableObjects[i]->getVel(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 + + volume = this->getVolume(); + + + + 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; + } + } } 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, 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 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; double temperature; - int n_oriented; - int n_constraints; + int nobj; - 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; + nobj = info->integrableObjects.size(); - // 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; + temperature = info->target_temp; -#ifdef IS_MPI - MPI_Allreduce(&ndf_local,&ndf,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD); - MPI_Allreduce(&ndfRaw_local,&ndfRaw,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD); -#else - ndfRaw = ndfRaw_local; - ndf = ndf_local; -#endif - ndf = ndf - 3; - - kebar = kb * temperature * (double)ndf / ( 2.0 * (double)ndfRaw ); + 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 - vx = vbar * gaussStream->getGaussian(); - vy = vbar * gaussStream->getGaussian(); - vz = vbar * gaussStream->getGaussian(); + for (j=0; j<3; j++) + aVel[j] = vbar * gaussStream->getGaussian(); + + info->integrableObjects[vr]->setVel( aVel ); + + if(info->integrableObjects[vr]->isDirectional()){ - atoms[vr]->set_vx( vx ); - atoms[vr]->set_vy( vy ); - atoms[vr]->set_vz( vz ); + 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. @@ -226,70 +326,42 @@ 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++){ - vx = atoms[vd]->get_vx(); - vy = atoms[vd]->get_vy(); - vz = atoms[vd]->get_vz(); - - vx -= vdrift[0]; - vy -= vdrift[1]; - vz -= vdrift[2]; + info->integrableObjects[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]; + + info->integrableObjects[vd]->setVel( aVel ); } - if( n_oriented ){ - - for( i=0; iisDirectional() ){ - - dAtom = (DirectionalAtom *)atoms[i]; - vbar = sqrt( 2.0 * kebar * dAtom->getIxx() ); - jx = 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 ); - } - } - } } void Thermo::getCOMVel(double vdrift[3]){ double mtot, mtot_local; + double aVel[3], amass; double vdrift_local[3]; - int vd, n_atoms; - 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 = entry_plug->n_atoms; - atoms = entry_plug->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++){ - 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 = info->integrableObjects[vd]->getMass(); + info->integrableObjects[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 @@ -308,3 +380,220 @@ 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, j; + int nobj; + + mtot_local = 0.0; + COM_local[0] = 0.0; + COM_local[1] = 0.0; + COM_local[2] = 0.0; + + nobj = info->integrableObjects.size(); + for(i = 0; i < nobj; i++){ + + amass = info->integrableObjects[i]->getMass(); + info->integrableObjects[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; + } +} + +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 ); + } +} + +void Thermo::removeAngularMomentum(){ + Vector3d vcom; + Vector3d qcom; + Vector3d pos; + Vector3d vel; + double mass; + double xx; + double yy; + double zz; + double xy; + double xz; + double yz; + Vector3d localAngMom; + Vector3d angMom; + Vector3d omega; + vector integrableObjects; + double localInertiaVec[9]; + double inertiaVec[9]; + vector qMinusQCom; + vector vMinusVCom; + Mat3x3d inertiaMat; + Mat3x3d inverseInertiaMat; + + integrableObjects = info->integrableObjects; + qMinusQCom.resize(integrableObjects.size()); + vMinusVCom.resize(integrableObjects.size()); + + getCOM(qcom.vec); + getCOMVel(vcom.vec); + + //initialize components for inertia tensor + xx = 0.0; + yy = 0.0; + zz = 0.0; + xy = 0.0; + xz = 0.0; + yz = 0.0; + + //build components of Inertia tensor + // + // [ Ixx -Ixy -Ixz ] + // J = | -Iyx Iyy -Iyz | + // [ -Izx -Iyz Izz ] + //See Fowles and Cassidy Chapter 9 or Goldstein Chapter 5 + for(size_t i = 0; i < integrableObjects.size(); i++){ + integrableObjects[i]->getPos(pos.vec); + integrableObjects[i]->getVel(vel.vec); + mass = integrableObjects[i]->getMass(); + + qMinusQCom[i] = pos - qcom; + info->wrapVector(qMinusQCom[i].vec); + + vMinusVCom[i] = vel - vcom; + + //compute moment of inertia coefficents + xx += qMinusQCom[i].x * qMinusQCom[i].x * mass; + yy += qMinusQCom[i].y * qMinusQCom[i].y * mass; + zz += qMinusQCom[i].z * qMinusQCom[i].z * mass; + + // compute products of inertia + xy += qMinusQCom[i].x * qMinusQCom[i].y * mass; + xz += qMinusQCom[i].x * qMinusQCom[i].z * mass; + yz += qMinusQCom[i].y * qMinusQCom[i].z * mass; + + localAngMom += crossProduct(qMinusQCom[i] , vMinusVCom[i] ) * mass; + + } + + localInertiaVec[0] =yy+zz; + localInertiaVec[1] = -xy; + localInertiaVec[2] = -xz; + localInertiaVec[3] = -xy; + localInertiaVec[4] = xx+zz; + localInertiaVec[5] = -yz; + localInertiaVec[6] = -xz; + localInertiaVec[7] = -yz; + localInertiaVec[8] = xx+yy; + + //Sum and distribute inertia and angmom arrays +#ifdef MPI + + MPI_Allreduce(localInertiaVec, inertiaVec, 9, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); + + MPI_Allreduce(localAngMom.vec, angMom.vec, 3, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); + + inertiaMat.element[0][0] = inertiaVec[0]; + inertiaMat.element[0][1] = inertiaVec[1]; + inertiaMat.element[0][2] = inertiaVec[2]; + + inertiaMat.element[1][0] = inertiaVec[3]; + inertiaMat.element[1][1] = inertiaVec[4]; + inertiaMat.element[1][2] = inertiaVec[5]; + + inertiaMat.element[2][0] = inertiaVec[6]; + inertiaMat.element[2][1] = inertiaVec[7]; + inertiaMat.element[2][2] = inertiaVec[8]; + +#else + + inertiaMat.element[0][0] = localInertiaVec[0]; + inertiaMat.element[0][1] = localInertiaVec[1]; + inertiaMat.element[0][2] = localInertiaVec[2]; + + inertiaMat.element[1][0] = localInertiaVec[3]; + inertiaMat.element[1][1] = localInertiaVec[4]; + inertiaMat.element[1][2] = localInertiaVec[5]; + + inertiaMat.element[2][0] = localInertiaVec[6]; + inertiaMat.element[2][1] = localInertiaVec[7]; + inertiaMat.element[2][2] = localInertiaVec[8]; + + angMom = localAngMom; +#endif + + //invert the moment of inertia tensor by LU-decomposition / backsolving: + + inverseInertiaMat = inertiaMat.inverse(); + + //calculate the angular velocities: omega = I^-1 . L + + omega = inverseInertiaMat * angMom; + + //subtract out center of mass velocity and angular momentum from + //particle velocities + + for(size_t i = 0; i < integrableObjects.size(); i++){ + vel = vMinusVCom[i] - crossProduct(omega, qMinusQCom[i]); + integrableObjects[i]->setVel(vel.vec); + } +} + +double Thermo::getConsEnergy(){ + ConstraintPair* consPair; + double totConsEnergy; + double bondLen2; + double dist; + double lamda; + + totConsEnergy = 0; + + for(cpIter->first(); !cpIter->isEnd(); cpIter->next()){ + consPair = cpIter->currentItem(); + bondLen2 = consPair->getBondLength2(); + lamda = consPair->getLamda(); + //dist = consPair->getDistance(); + + //totConsEnergy += lamda * (dist*dist - bondLen2); + } + + return totConsEnergy; +} + +