43 |
|
#include "math/SquareMatrix3.hpp" |
44 |
|
#include "primitives/Molecule.hpp" |
45 |
|
#include "primitives/StuntDouble.hpp" |
46 |
< |
#include "math/MersenneTwister.hpp" |
46 |
> |
|
47 |
|
namespace oopse { |
48 |
|
|
49 |
+ |
Velocitizer::Velocitizer(SimInfo* info) { |
50 |
+ |
|
51 |
+ |
int seedValue; |
52 |
+ |
Globals * simParams = info->getSimParams(); |
53 |
+ |
|
54 |
+ |
if (simParams->haveSeed()) { |
55 |
+ |
seedValue = simParams->getSeed(); |
56 |
+ |
randNumGen_ = new OOPSERandNumGen(seedValue); |
57 |
+ |
}else { |
58 |
+ |
randNumGen_ = new OOPSERandNumGen(); |
59 |
+ |
} |
60 |
+ |
|
61 |
+ |
} |
62 |
+ |
|
63 |
+ |
Velocitizer::~Velocitizer() { |
64 |
+ |
delete randNumGen_; |
65 |
+ |
} |
66 |
+ |
|
67 |
|
void Velocitizer::velocitize(double temperature) { |
68 |
|
Vector3d aVel; |
69 |
|
Vector3d aJ; |
83 |
|
Molecule * mol; |
84 |
|
StuntDouble * integrableObject; |
85 |
|
|
68 |
– |
|
69 |
– |
#ifndef IS_MPI |
70 |
– |
MTRand randNumGen(info_->getSeed()); |
71 |
– |
#else |
72 |
– |
int nProcessors; |
73 |
– |
MPI_Comm_size(MPI_COMM_WORLD, &nProcessors); |
74 |
– |
MTRand randNumGen(info_->getSeed(), nProcessors, worldRank); |
75 |
– |
#endif |
86 |
|
|
87 |
+ |
|
88 |
|
kebar = kb * temperature * info_->getNdfRaw() / (2.0 * info_->getNdf()); |
89 |
|
|
90 |
|
for( mol = info_->beginMolecule(i); mol != NULL; |
102 |
|
// centered on vbar |
103 |
|
|
104 |
|
for( int k = 0; k < 3; k++ ) { |
105 |
< |
aVel[k] = vbar * randNumGen.randNorm(0.0, 1.0); |
105 |
> |
aVel[k] = vbar * randNumGen_->randNorm(0.0, 1.0); |
106 |
|
} |
107 |
|
|
108 |
|
integrableObject->setVel(aVel); |
117 |
|
|
118 |
|
aJ[l] = 0.0; |
119 |
|
vbar = sqrt(2.0 * kebar * I(m, m)); |
120 |
< |
aJ[m] = vbar * randNumGen.randNorm(0.0, 1.0); |
120 |
> |
aJ[m] = vbar * randNumGen_->randNorm(0.0, 1.0); |
121 |
|
vbar = sqrt(2.0 * kebar * I(n, n)); |
122 |
< |
aJ[n] = vbar * randNumGen.randNorm(0.0, 1.0); |
122 |
> |
aJ[n] = vbar * randNumGen_->randNorm(0.0, 1.0); |
123 |
|
} else { |
124 |
|
for( int k = 0; k < 3; k++ ) { |
125 |
|
vbar = sqrt(2.0 * kebar * I(k, k)); |
126 |
< |
aJ[k] = vbar * randNumGen.randNorm(0.0, 1.0); |
126 |
> |
aJ[k] = vbar *randNumGen_->randNorm(0.0, 1.0); |
127 |
|
} |
128 |
|
} // else isLinear |
129 |
|
|