| 76 |
|
// Q can take values from 0 to 1 |
| 77 |
|
|
| 78 |
|
MinQ_ = 0.0; |
| 79 |
< |
MaxQ_ = 1.0; |
| 79 |
> |
MaxQ_ = 3.0; |
| 80 |
|
deltaQ_ = (MaxQ_ - MinQ_) / nbins; |
| 81 |
|
Q_histogram_.resize(nbins); |
| 82 |
|
|
| 114 |
|
RealType phi; |
| 115 |
|
RealType r; |
| 116 |
|
RealType dist; |
| 117 |
< |
std::map<int, ComplexType> QBar_lm; |
| 117 |
> |
std::map<int,ComplexType> QBar_lm; |
| 118 |
|
RealType QSq_l; |
| 119 |
|
RealType Q_l; |
| 120 |
|
ComplexType W_l; |
| 208 |
|
// Normalize Qbar2 |
| 209 |
|
for (int m = -lNumber_;m <= lNumber_; m++){ |
| 210 |
|
QBar_lm[m] /= nBonds; |
| 211 |
– |
std::cout << "m = " << m << " QBLM = " << QBar_lm[m] << "\n"; |
| 211 |
|
} |
| 212 |
|
|
| 213 |
|
// Find second order invariant Q_l |
| 216 |
|
for (int m = -lNumber_; m <= lNumber_; m++){ |
| 217 |
|
QSq_l += norm(QBar_lm[m]); |
| 218 |
|
} |
| 219 |
< |
std::cout << "qsq_l = " << QSq_l << "\n"; |
| 219 |
> |
|
| 220 |
|
Q_l = sqrt(QSq_l * 4.0 * NumericConstant::PI / |
| 221 |
|
(2.0*(RealType)lNumber_ + 1.0)); |
| 222 |
|
|
| 226 |
|
for (int m1 = -lNumber_; m1 <= lNumber_; m1++) { |
| 227 |
|
// Zero work array |
| 228 |
|
for (int ii = 0; ii < mSize_; ii++){ |
| 229 |
< |
THRCOF[i] = 0.0; |
| 229 |
> |
THRCOF[ii] = 0.0; |
| 230 |
|
} |
| 231 |
|
// Get Wigner coefficients |
| 232 |
|
m1Pass = (double)m1; |
| 233 |
< |
Wigner3jm(&l_, &l_, &l_, &m1Pass, &m2Min, &m2Max, THRCOF, &mSize_, &error); |
| 234 |
< |
for (int m_index = 1; m_index < (int)(m2Max - m2Min-1.0); m_index++) { |
| 235 |
< |
m2 = floor(m2Min) + m_index - 1; |
| 233 |
> |
|
| 234 |
> |
Wigner3jm(&l_, &l_, &l_, |
| 235 |
> |
&m1Pass, &m2Min, &m2Max, |
| 236 |
> |
THRCOF, &mSize_, &error); |
| 237 |
> |
|
| 238 |
> |
for (int mmm = 0; mmm < (int)(m2Max - m2Min); mmm++) { |
| 239 |
> |
m2 = (int)floor(m2Min) + mmm; |
| 240 |
|
m3 = -m1-m2; |
| 241 |
< |
W_l += THRCOF[m_index]*QBar_lm[m1]*QBar_lm[m2]*QBar_lm[m3]; |
| 241 |
> |
W_l += THRCOF[mmm] * QBar_lm[m1] * QBar_lm[m2] * QBar_lm[m3]; |
| 242 |
|
} |
| 243 |
|
} |
| 244 |
|
|
| 305 |
|
// Normalize by number of frames and write it out: |
| 306 |
|
for (int i = 0; i < Q_histogram_.size(); ++i) { |
| 307 |
|
RealType Qval = MinQ_ + (i + 0.5) * deltaQ_; |
| 308 |
< |
osq << Qval << "\t" << Q_histogram_[i] / frameCounter_ << "\n"; |
| 308 |
> |
osq << Qval << "\t" << (RealType)Q_histogram_[i] / (RealType)Qcount_ << "\n"; |
| 309 |
|
} |
| 310 |
|
|
| 311 |
|
osq.close(); |
| 331 |
|
// Normalize by number of frames and write it out: |
| 332 |
|
for (int i = 0; i < W_histogram_.size(); ++i) { |
| 333 |
|
RealType Wval = MinW_ + (i + 0.5) * deltaW_; |
| 334 |
< |
osw << Wval << "\t" << W_histogram_[i] / frameCounter_ << "\n"; |
| 334 |
> |
osw << Wval << "\t" << (RealType)W_histogram_[i] / (RealType)Wcount_ << "\n"; |
| 335 |
|
} |
| 336 |
|
|
| 337 |
|
osw.close(); |