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 |
|
} |