--- trunk/OOPSE/libmdtools/Thermo.cpp 2004/05/24 21:03:30 1192 +++ trunk/OOPSE/libmdtools/Thermo.cpp 2004/08/23 15:11:36 1452 @@ -11,6 +11,8 @@ using namespace std; #include "Integrator.hpp" #include "simError.h" #include "MatVec3.h" +#include "ConstraintManager.hpp" +#include "Mat3x3d.hpp" #ifdef IS_MPI #define __C @@ -26,10 +28,13 @@ Thermo::Thermo( SimInfo* the_info ) { int baseSeed = the_info->getSeed(); gaussStream = new gaussianSPRNG( baseSeed ); + + cpIter = info->consMan->createPairIterator(); } Thermo::~Thermo(){ delete gaussStream; + delete cpIter; } double Thermo::getKinetic(){ @@ -200,37 +205,22 @@ void Thermo::getPressureTensor(double press[3][3]){ const double e_convert = 4.184e-4; double molmass, volume; - double vcom[3], pcom[3], fcom[3], scaled[3]; + double vcom[3]; double p_local[9], p_global[9]; - int i, j, k, nMols; - Molecule* molecules; + int i, j, k; - 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++) { 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); - 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]); p_local[1] += molmass * (vcom[0] * vcom[1]); @@ -241,11 +231,11 @@ void Thermo::getPressureTensor(double press[3][3]){ 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 @@ -256,6 +246,8 @@ void Thermo::getPressureTensor(double press[3][3]){ volume = this->getVolume(); + + for(i = 0; i < 3; i++) { for (j = 0; j < 3; j++) { k = 3*i + j; @@ -449,5 +441,159 @@ void Thermo::removeCOMdrift() { 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; +} + +