| 50 |
|
#include "math/ParallelRandNumGen.hpp" |
| 51 |
|
#endif |
| 52 |
|
|
| 53 |
+ |
/* Remove me after testing*/ |
| 54 |
+ |
#include <cstdio> |
| 55 |
+ |
#include <iostream> |
| 56 |
+ |
/*End remove me*/ |
| 57 |
+ |
|
| 58 |
|
namespace oopse { |
| 59 |
|
|
| 60 |
|
Velocitizer::Velocitizer(SimInfo* info) : info_(info) { |
| 83 |
|
delete randNumGen_; |
| 84 |
|
} |
| 85 |
|
|
| 86 |
< |
void Velocitizer::velocitize(double temperature) { |
| 86 |
> |
void Velocitizer::velocitize(RealType temperature) { |
| 87 |
|
Vector3d aVel; |
| 88 |
|
Vector3d aJ; |
| 89 |
|
Mat3x3d I; |
| 91 |
|
int m; |
| 92 |
|
int n; |
| 93 |
|
Vector3d vdrift; |
| 94 |
< |
double vbar; |
| 94 |
> |
RealType vbar; |
| 95 |
|
/**@todo refactory kb */ |
| 96 |
< |
const double kb = 8.31451e-7; // kb in amu, angstroms, fs, etc. |
| 97 |
< |
double av2; |
| 98 |
< |
double kebar; |
| 96 |
> |
const RealType kb = 8.31451e-7; // kb in amu, angstroms, fs, etc. |
| 97 |
> |
RealType av2; |
| 98 |
> |
RealType kebar; |
| 99 |
|
|
| 100 |
+ |
Globals * simParams = info_->getSimParams(); |
| 101 |
+ |
|
| 102 |
|
SimInfo::MoleculeIterator i; |
| 103 |
|
Molecule::IntegrableObjectIterator j; |
| 104 |
|
Molecule * mol; |
| 156 |
|
|
| 157 |
|
|
| 158 |
|
removeComDrift(); |
| 159 |
+ |
// Remove angular drift if we are not using periodic boundary conditions. |
| 160 |
+ |
if(!simParams->getUsePeriodicBoundaryConditions()) removeAngularDrift(); |
| 161 |
|
|
| 162 |
|
} |
| 163 |
|
|
| 185 |
|
|
| 186 |
|
} |
| 187 |
|
|
| 188 |
+ |
|
| 189 |
+ |
void Velocitizer::removeAngularDrift() { |
| 190 |
+ |
// Get the Center of Mass drift velocity. |
| 191 |
+ |
|
| 192 |
+ |
Vector3d vdrift; |
| 193 |
+ |
Vector3d com; |
| 194 |
+ |
|
| 195 |
+ |
info_->getComAll(com,vdrift); |
| 196 |
+ |
|
| 197 |
+ |
Mat3x3d inertiaTensor; |
| 198 |
+ |
Vector3d angularMomentum; |
| 199 |
+ |
Vector3d omega; |
| 200 |
+ |
|
| 201 |
+ |
|
| 202 |
+ |
|
| 203 |
+ |
info_->getInertiaTensor(inertiaTensor,angularMomentum); |
| 204 |
+ |
// We now need the inverse of the inertia tensor. |
| 205 |
+ |
/* |
| 206 |
+ |
std::cerr << "Angular Momentum before is " |
| 207 |
+ |
<< angularMomentum << std::endl; |
| 208 |
+ |
std::cerr << "Inertia Tensor before is " |
| 209 |
+ |
<< inertiaTensor << std::endl; |
| 210 |
+ |
*/ |
| 211 |
+ |
|
| 212 |
+ |
inertiaTensor =inertiaTensor.inverse(); |
| 213 |
+ |
/* |
| 214 |
+ |
std::cerr << "Inertia Tensor after inverse is " |
| 215 |
+ |
<< inertiaTensor << std::endl; |
| 216 |
+ |
*/ |
| 217 |
+ |
omega = inertiaTensor*angularMomentum; |
| 218 |
+ |
|
| 219 |
+ |
SimInfo::MoleculeIterator i; |
| 220 |
+ |
Molecule::IntegrableObjectIterator j; |
| 221 |
+ |
Molecule * mol; |
| 222 |
+ |
StuntDouble * integrableObject; |
| 223 |
+ |
Vector3d tempComPos; |
| 224 |
+ |
|
| 225 |
+ |
// Corrects for the center of mass angular drift. |
| 226 |
+ |
// sums all the angular momentum and divides by total mass. |
| 227 |
+ |
for( mol = info_->beginMolecule(i); mol != NULL; |
| 228 |
+ |
mol = info_->nextMolecule(i) ) { |
| 229 |
+ |
for( integrableObject = mol->beginIntegrableObject(j); |
| 230 |
+ |
integrableObject != NULL; |
| 231 |
+ |
integrableObject = mol->nextIntegrableObject(j) ) { |
| 232 |
+ |
tempComPos = integrableObject->getPos()-com; |
| 233 |
+ |
integrableObject->setVel((integrableObject->getVel() - vdrift)-cross(omega,tempComPos)); |
| 234 |
+ |
} |
| 235 |
+ |
} |
| 236 |
+ |
|
| 237 |
+ |
angularMomentum = info_->getAngularMomentum(); |
| 238 |
+ |
/* |
| 239 |
+ |
std::cerr << "Angular Momentum after is " |
| 240 |
+ |
<< angularMomentum << std::endl; |
| 241 |
+ |
*/ |
| 242 |
+ |
|
| 243 |
+ |
} |
| 244 |
+ |
|
| 245 |
+ |
|
| 246 |
+ |
|
| 247 |
+ |
|
| 248 |
|
} |
