--- branches/mmeineke/OOPSE/libmdtools/Thermo.cpp 2003/03/21 17:42:12 377 +++ trunk/OOPSE/libmdtools/Thermo.cpp 2003/04/03 19:58:24 445 @@ -10,7 +10,14 @@ using namespace std; #include "Thermo.hpp" #include "SRI.hpp" #include "Integrator.hpp" +#include "simError.h" +#ifdef IS_MPI +#define __C +#include "mpiSimulation.hpp" +#endif // is_mpi + + #define BASE_SEED 123456789 Thermo::Thermo( SimInfo* the_entry_plug ) { @@ -77,30 +84,43 @@ double Thermo::getPotential(){ double Thermo::getPotential(){ + double potential_local; double potential; - double potential_global; int el, nSRI; - SRI** sris; + Molecule* molecules; - sris = entry_plug->sr_interactions; + molecules = entry_plug->molecules; nSRI = entry_plug->n_SRI; + potential_local = 0.0; potential = 0.0; - potential_global = 0.0; - potential += entry_plug->lrPot; + potential_local += entry_plug->lrPot; - for( el=0; elget_potential(); + for( el=0; eln_mol; el++ ){ + potential_local += molecules[el].getPotential(); } +#ifdef IS_MPI + /* + std::cerr << "node " << worldRank << ": before LONG RANGE pot = " << entry_plug->lrPot + << "; pot_local = " << potential_local + << "; pot = " << potential << "\n"; + */ +#endif + // Get total potential for entire system from MPI. #ifdef IS_MPI - MPI::COMM_WORLD.Allreduce(&potential,&potential_global,1,MPI_DOUBLE,MPI_SUM); - potential = potential_global; - + MPI::COMM_WORLD.Allreduce(&potential_local,&potential,1,MPI_DOUBLE,MPI_SUM); +#else + potential = potential_local; #endif // is_mpi +#ifdef IS_MPI + /* + std::cerr << "node " << worldRank << ": after pot = " << potential << "\n"; + */ +#endif + return potential; } @@ -116,20 +136,28 @@ double Thermo::getTemperature(){ const double kb = 1.9872179E-3; // boltzman's constant in kcal/(mol K) double temperature; + int ndf_local, ndf; - int ndf = 3 * entry_plug->n_atoms + 3 * entry_plug->n_oriented - - entry_plug->n_constraints - 3; + ndf_local = 3 * entry_plug->n_atoms + 3 * entry_plug->n_oriented + - entry_plug->n_constraints; +#ifdef IS_MPI + MPI::COMM_WORLD.Allreduce(&ndf_local,&ndf,1,MPI_INT,MPI_SUM); +#else + ndf = ndf_local; +#endif + + ndf = ndf - 3; + temperature = ( 2.0 * this->getKinetic() ) / ( ndf * kb ); return temperature; } double Thermo::getPressure(){ + // returns pressure 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 conv_Pa_atm = 9.901E-6; // convert Pa -> atm -// const double conv_internal_Pa = 1.661E-7; //convert amu/(fs^2 A) -> Pa -// const double conv_A_m = 1.0E-10; //convert A -> m - return 0.0; } @@ -140,13 +168,12 @@ void Thermo::velocitize() { double jx, jy, jz; int i, vr, vd; // velocity randomizer loop counters double vdrift[3]; - double mtot = 0.0; double vbar; const double kb = 8.31451e-7; // kb in amu, angstroms, fs, etc. double av2; double kebar; - int ndf; // number of degrees of freedom - int ndfRaw; // the raw number of degrees of freedom + 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; @@ -160,9 +187,22 @@ void Thermo::velocitize() { n_oriented = entry_plug->n_oriented; n_constraints = entry_plug->n_constraints; + // Raw degrees of freedom that we have to set + ndfRaw_local = 3 * entry_plug->n_atoms + 3 * entry_plug->n_oriented; - ndfRaw = 3 * n_atoms + 3 * n_oriented; - ndf = ndfRaw - n_constraints - 3; + // Degrees of freedom that can contain kinetic energy + ndf_local = 3 * entry_plug->n_atoms + 3 * entry_plug->n_oriented + - entry_plug->n_constraints; + +#ifdef IS_MPI + MPI::COMM_WORLD.Allreduce(&ndf_local,&ndf,1,MPI_INT,MPI_SUM); + MPI::COMM_WORLD.Allreduce(&ndfRaw_local,&ndfRaw,1,MPI_INT,MPI_SUM); +#else + ndfRaw = ndfRaw_local; + ndf = ndf_local; +#endif + ndf = ndf - 3; + kebar = kb * temperature * (double)ndf / ( 2.0 * (double)ndfRaw ); for(vr = 0; vr < n_atoms; vr++){ @@ -171,7 +211,7 @@ void Thermo::velocitize() { av2 = 2.0 * kebar / atoms[vr]->getMass(); vbar = sqrt( av2 ); - + // vbar = sqrt( 8.31451e-7 * temperature / atoms[vr]->getMass() ); // picks random velocities from a gaussian distribution @@ -185,35 +225,20 @@ void Thermo::velocitize() { atoms[vr]->set_vy( vy ); atoms[vr]->set_vz( vz ); } + + // 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. - - mtot = 0.0; - vdrift[0] = 0.0; - vdrift[1] = 0.0; - vdrift[2] = 0.0; - for(vd = 0; vd < n_atoms; vd++){ - - vdrift[0] += atoms[vd]->get_vx() * atoms[vd]->getMass(); - vdrift[1] += atoms[vd]->get_vy() * atoms[vd]->getMass(); - vdrift[2] += atoms[vd]->get_vz() * atoms[vd]->getMass(); - - mtot += atoms[vd]->getMass(); - } - - for (vd = 0; vd < 3; vd++) { - vdrift[vd] = vdrift[vd] / mtot; - } - 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]; @@ -246,3 +271,47 @@ void Thermo::velocitize() { } } } + +void Thermo::getCOMVel(double vdrift[3]){ + + double mtot, mtot_local; + double vdrift_local[3]; + int vd, n_atoms; + Atom** atoms; + + // 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 = 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++){ + + 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(); + + mtot_local += atoms[vd]->getMass(); + } + +#ifdef IS_MPI + MPI::COMM_WORLD.Allreduce(&mtot_local,&mtot,1,MPI_DOUBLE,MPI_SUM); + MPI::COMM_WORLD.Allreduce(vdrift_local,vdrift,3,MPI_DOUBLE,MPI_SUM); +#else + mtot = mtot_local; + for(vd = 0; vd < 3; vd++) { + vdrift[vd] = vdrift_local[vd]; + } +#endif + + for (vd = 0; vd < 3; vd++) { + vdrift[vd] = vdrift[vd] / mtot; + } + +} +