43 |
|
#include "math/SquareMatrix3.hpp" |
44 |
|
#include "primitives/Molecule.hpp" |
45 |
|
#include "primitives/StuntDouble.hpp" |
46 |
< |
#include "math/MersenneTwister.hpp" |
46 |
> |
|
47 |
> |
#ifndef IS_MPI |
48 |
> |
#include "math/SeqRandNumGen.hpp" |
49 |
> |
#else |
50 |
> |
#include "math/ParallelRandNumGen.hpp" |
51 |
> |
#endif |
52 |
> |
|
53 |
|
namespace oopse { |
54 |
|
|
55 |
+ |
Velocitizer::Velocitizer(SimInfo* info) { |
56 |
+ |
|
57 |
+ |
int seedValue; |
58 |
+ |
Globals * simParams = info->getSimParams(); |
59 |
+ |
|
60 |
+ |
#ifndef IS_MPI |
61 |
+ |
if (simParams->haveSeed()) { |
62 |
+ |
seedValue = simParams->getSeed(); |
63 |
+ |
randNumGen_ = new SeqRandNumGen(seedValue); |
64 |
+ |
}else { |
65 |
+ |
randNumGen_ = new SeqRandNumGen(); |
66 |
+ |
} |
67 |
+ |
#else |
68 |
+ |
if (simParams->haveSeed()) { |
69 |
+ |
seedValue = simParams->getSeed(); |
70 |
+ |
randNumGen_ = new ParallelRandNumGen(seedValue); |
71 |
+ |
}else { |
72 |
+ |
randNumGen_ = new ParallelRandNumGen(); |
73 |
+ |
} |
74 |
+ |
#endif |
75 |
+ |
} |
76 |
+ |
|
77 |
+ |
Velocitizer::~Velocitizer() { |
78 |
+ |
delete randNumGen_; |
79 |
+ |
} |
80 |
+ |
|
81 |
|
void Velocitizer::velocitize(double temperature) { |
82 |
|
Vector3d aVel; |
83 |
|
Vector3d aJ; |
97 |
|
Molecule * mol; |
98 |
|
StuntDouble * integrableObject; |
99 |
|
|
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 |
100 |
|
|
101 |
+ |
|
102 |
|
kebar = kb * temperature * info_->getNdfRaw() / (2.0 * info_->getNdf()); |
103 |
|
|
104 |
|
for( mol = info_->beginMolecule(i); mol != NULL; |
116 |
|
// centered on vbar |
117 |
|
|
118 |
|
for( int k = 0; k < 3; k++ ) { |
119 |
< |
aVel[k] = vbar * randNumGen.randNorm(0.0, 1.0); |
119 |
> |
aVel[k] = vbar * randNumGen_->randNorm(0.0, 1.0); |
120 |
|
} |
121 |
|
|
122 |
|
integrableObject->setVel(aVel); |
131 |
|
|
132 |
|
aJ[l] = 0.0; |
133 |
|
vbar = sqrt(2.0 * kebar * I(m, m)); |
134 |
< |
aJ[m] = vbar * randNumGen.randNorm(0.0, 1.0); |
134 |
> |
aJ[m] = vbar * randNumGen_->randNorm(0.0, 1.0); |
135 |
|
vbar = sqrt(2.0 * kebar * I(n, n)); |
136 |
< |
aJ[n] = vbar * randNumGen.randNorm(0.0, 1.0); |
136 |
> |
aJ[n] = vbar * randNumGen_->randNorm(0.0, 1.0); |
137 |
|
} else { |
138 |
|
for( int k = 0; k < 3; k++ ) { |
139 |
|
vbar = sqrt(2.0 * kebar * I(k, k)); |
140 |
< |
aJ[k] = vbar * randNumGen.randNorm(0.0, 1.0); |
140 |
> |
aJ[k] = vbar *randNumGen_->randNorm(0.0, 1.0); |
141 |
|
} |
142 |
|
} // else isLinear |
143 |
|
|