| 51 |
|
#include "io/DumpReader.hpp" |
| 52 |
|
#include "primitives/Molecule.hpp" |
| 53 |
|
#include "utils/NumericConstant.hpp" |
| 54 |
< |
#include "math/RealSphericalHarmonic.hpp" |
| 54 |
> |
#include "math/SphericalHarmonic.hpp" |
| 55 |
|
|
| 56 |
|
namespace oopse { |
| 57 |
|
|
| 108 |
|
Molecule::RigidBodyIterator rbIter; |
| 109 |
|
Molecule::AtomIterator ai; |
| 110 |
|
StuntDouble* sd; |
| 111 |
< |
RealType theta; |
| 111 |
> |
RealType costheta; |
| 112 |
|
RealType phi; |
| 113 |
|
RealType r; |
| 114 |
|
RealType dist; |
| 115 |
< |
std::map<int, RealType> QBar_lm; |
| 115 |
> |
std::map<int, ComplexType> QBar_lm; |
| 116 |
|
RealType QSq_l; |
| 117 |
|
RealType Q_l; |
| 118 |
|
int nBonds; |
| 119 |
< |
RealSphericalHarmonic sphericalHarmonic; |
| 119 |
> |
SphericalHarmonic sphericalHarmonic; |
| 120 |
|
int i, j; |
| 121 |
|
|
| 122 |
|
// Set the l for the spherical harmonic, it doesn't change |
| 180 |
|
// Check to see if neighbor is in bond cutoff |
| 181 |
|
|
| 182 |
|
if (r < rCut_) { |
| 183 |
< |
theta = atan2(vec.y(), vec.x()); |
| 184 |
< |
phi = acos(vec.z()/r); |
| 183 |
> |
costheta = vec.z() / r; |
| 184 |
> |
phi = atan2(vec.y(), vec.x()); |
| 185 |
> |
|
| 186 |
|
for(int m = -lNumber_; m <= lNumber_; m++){ |
| 187 |
|
sphericalHarmonic.setM(m); |
| 188 |
< |
QBar_lm[m] += sphericalHarmonic.getValueAt(theta,phi); |
| 188 |
> |
QBar_lm[m] += sphericalHarmonic.getValueAt(costheta,phi); |
| 189 |
|
} |
| 190 |
|
nBonds++; |
| 191 |
|
} |
| 192 |
|
} |
| 193 |
|
} |
| 194 |
|
|
| 195 |
< |
// Normalize Qbar |
| 195 |
> |
// Normalize Qbar2 |
| 196 |
|
for (int m = -lNumber_;m <= lNumber_; m++){ |
| 197 |
|
QBar_lm[m] /= nBonds; |
| 198 |
|
} |
| 201 |
|
|
| 202 |
|
QSq_l = 0.0; |
| 203 |
|
for (int m = -lNumber_; m <= lNumber_; m++){ |
| 204 |
< |
QSq_l += pow(QBar_lm[m], 2); |
| 204 |
> |
QSq_l += norm(QBar_lm[m]); |
| 205 |
|
} |
| 206 |
|
Q_l = sqrt(QSq_l*(4.0 * NumericConstant::PI / (2.0*(RealType)lNumber_ + 1))); |
| 207 |
|
|
| 214 |
|
double m1Pass, m2Min, m2Max; |
| 215 |
|
int error, m1, m2, m3; |
| 216 |
|
|
| 217 |
< |
RealType W_l; |
| 218 |
< |
RealType W_l_hat; |
| 217 |
> |
ComplexType W_l; |
| 218 |
> |
ComplexType W_l_hat; |
| 219 |
|
W_l = 0.0; |
| 220 |
|
for (int m1 = -lNumber_; m1 <= lNumber_; m1++) { |
| 221 |
|
// Zero work array |
| 236 |
|
|
| 237 |
|
// accumulate histogram data for Q_l and W_l_hat: |
| 238 |
|
|
| 239 |
< |
collectHistogram(Q_l, W_l_hat); |
| 239 |
> |
collectHistogram(Q_l, real(W_l_hat)); |
| 240 |
|
|
| 241 |
|
} |
| 242 |
|
} |
