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 |
+ |
#ifndef IS_MPI |
55 |
+ |
if (simParams->haveSeed()) { |
56 |
+ |
seedValue = simParams->getSeed(); |
57 |
+ |
randNumGen_ = new MTRand(seedValue); |
58 |
+ |
}else { |
59 |
+ |
randNumGen_ = new MTRand(); |
60 |
+ |
} |
61 |
+ |
#else |
62 |
+ |
if (simParams->haveSeed()) { |
63 |
+ |
seedValue = simParams->getSeed(); |
64 |
+ |
randNumGen_ = new ParallelRandNumGen(seedValue); |
65 |
+ |
}else { |
66 |
+ |
randNumGen_ = new ParallelRandNumGen(); |
67 |
+ |
} |
68 |
+ |
#endif |
69 |
+ |
} |
70 |
+ |
|
71 |
+ |
Velocitizer::~Velocitizer() { |
72 |
+ |
delete randNumGen_; |
73 |
+ |
} |
74 |
+ |
|
75 |
|
void Velocitizer::velocitize(double temperature) { |
76 |
|
Vector3d aVel; |
77 |
|
Vector3d aJ; |
91 |
|
Molecule * mol; |
92 |
|
StuntDouble * integrableObject; |
93 |
|
|
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 |
94 |
|
|
95 |
+ |
|
96 |
|
kebar = kb * temperature * info_->getNdfRaw() / (2.0 * info_->getNdf()); |
97 |
|
|
98 |
|
for( mol = info_->beginMolecule(i); mol != NULL; |
110 |
|
// centered on vbar |
111 |
|
|
112 |
|
for( int k = 0; k < 3; k++ ) { |
113 |
< |
aVel[k] = vbar * randNumGen.randNorm(0.0, 1.0); |
113 |
> |
aVel[k] = vbar * randNumGen_->randNorm(0.0, 1.0); |
114 |
|
} |
115 |
|
|
116 |
|
integrableObject->setVel(aVel); |
125 |
|
|
126 |
|
aJ[l] = 0.0; |
127 |
|
vbar = sqrt(2.0 * kebar * I(m, m)); |
128 |
< |
aJ[m] = vbar * randNumGen.randNorm(0.0, 1.0); |
128 |
> |
aJ[m] = vbar * randNumGen_->randNorm(0.0, 1.0); |
129 |
|
vbar = sqrt(2.0 * kebar * I(n, n)); |
130 |
< |
aJ[n] = vbar * randNumGen.randNorm(0.0, 1.0); |
130 |
> |
aJ[n] = vbar * randNumGen_->randNorm(0.0, 1.0); |
131 |
|
} else { |
132 |
|
for( int k = 0; k < 3; k++ ) { |
133 |
|
vbar = sqrt(2.0 * kebar * I(k, k)); |
134 |
< |
aJ[k] = vbar * randNumGen.randNorm(0.0, 1.0); |
134 |
> |
aJ[k] = vbar *randNumGen_->randNorm(0.0, 1.0); |
135 |
|
} |
136 |
|
} // else isLinear |
137 |
|
|