ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/applications/staticProps/BondOrderParameter.cpp
(Generate patch)

Comparing trunk/src/applications/staticProps/BondOrderParameter.cpp (file contents):
Revision 1052 by gezelter, Tue Sep 26 01:30:32 2006 UTC vs.
Revision 1053 by gezelter, Tue Sep 26 14:33:56 2006 UTC

# Line 85 | Line 85 | namespace oopse {
85      MinW_ = -0.25;
86      MaxW_ = 0.25;
87      deltaW_ = (MaxW_ - MinW_) / nbins;
88 –  }
88  
89 +    // Make arrays for Wigner3jm
90 +    double* THRCOF = new double[2*lMax_+1];
91 +    // Variables for Wigner routine
92 +    double lPass, m1Pass, m2m, m2M;
93 +    int error, mSize;
94 +    mSize = 2*lMax_+1;
95 +
96 +    for (int l = 0; l <= lMax_; l++) {
97 +      lPass = (double)l;
98 +      for (int m1 = -l; m1 <= l; m1++) {
99 +        m1Pass = (double)m1;
100 +
101 +        std::pair<int,int> lm = std::make_pair(l, m1);
102 +        
103 +        // Zero work array
104 +        for (int ii = 0; ii < 2*l + 1; ii++){
105 +          THRCOF[ii] = 0.0;
106 +        }
107 +            
108 +        // Get Wigner coefficients
109 +        Wigner3jm(&lPass, &lPass, &lPass,
110 +                  &m1Pass, &m2m, &m2M,
111 +                  THRCOF, &mSize, &error);
112 +        
113 +        m2Min[lm] = (int)floor(m2m);
114 +        m2Max[lm] = (int)floor(m2M);
115 +        
116 +        for (int mmm = 0; mmm < (int)(m2M - m2m); mmm++) {
117 +          w3j[lm].push_back(THRCOF[mmm]);
118 +        }
119 +      }
120 +    }
121 +  }
122 +  
123    BondOrderParameter::~BondOrderParameter() {
124      Q_histogram_.clear();
125      W_histogram_.clear();
# Line 128 | Line 161 | namespace oopse {
161      int nBonds, Nbonds;
162      SphericalHarmonic sphericalHarmonic;
163      int i, j;
131 –    // Make arrays for Wigner3jm
132 –    double* THRCOF = new double[2*lMax_+1];
133 –    // Variables for Wigner routine
134 –    double lPass, m1Pass, m2Min, m2Max;
135 –    int error, m1, m2, m3, mSize;
136 –    mSize = 2*lMax_+1;
164  
165      DumpReader reader(info_, dumpFilename_);    
166      int nFrames = reader.getNFrames();
# Line 245 | Line 272 | namespace oopse {
272      
273          for (int l = 0; l <= lMax_; l++) {
274            w[l] = 0.0;
248 –          lPass = (double)l;
275            for (int m1 = -l; m1 <= l; m1++) {
276 <            // Zero work array
277 <            for (int ii = 0; ii < 2*l + 1; ii++){
278 <              THRCOF[ii] = 0.0;
279 <            }
280 <            // Get Wigner coefficients
281 <            m1Pass = (double)m1;
256 <            
257 <            Wigner3jm(&lPass, &lPass, &lPass,
258 <                      &m1Pass, &m2Min, &m2Max,
259 <                      THRCOF, &mSize, &error);
260 <            
261 <            for (int mmm = 0; mmm < (int)(m2Max - m2Min); mmm++) {
262 <              m2 = (int)floor(m2Min) + mmm;
263 <              m3 = -m1-m2;
264 <              w[l] += THRCOF[mmm] *
265 <                q[std::make_pair(l,m1)] *
266 <                q[std::make_pair(l,m2)] *
267 <                q[std::make_pair(l,m3)];
276 >            std::pair<int,int> lm = std::make_pair(l, m1);
277 >            for (int mmm = 0; mmm < (m2Max[lm] - m2Min[lm]); mmm++) {
278 >              int m2 = m2Min[lm] + mmm;
279 >              int m3 = -m1-m2;
280 >              w[l] += w3j[lm][mmm] * q[lm] *
281 >                q[std::make_pair(l,m2)] *  q[std::make_pair(l,m3)];
282              }
283            }
284            
# Line 303 | Line 317 | namespace oopse {
317      
318      for (int l = 0; l <= lMax_; l++) {
319        W[l] = 0.0;
306 –      lPass = (double)l;
320        for (int m1 = -l; m1 <= l; m1++) {
321 <        // Zero work array
322 <        for (int ii = 0; ii < 2*l + 1; ii++){
323 <          THRCOF[ii] = 0.0;
321 >        std::pair<int,int> lm = std::make_pair(l, m1);
322 >        for (int mmm = 0; mmm < (m2Max[lm] - m2Min[lm]); mmm++) {
323 >          int m2 = m2Min[lm] + mmm;
324 >          int m3 = -m1-m2;
325 >          W[l] += w3j[lm][mmm] * QBar[lm] *
326 >            QBar[std::make_pair(l,m2)] * QBar[std::make_pair(l,m3)];
327          }
312 –        // Get Wigner coefficients
313 –        m1Pass = (double)m1;
314 –        
315 –        Wigner3jm(&lPass, &lPass, &lPass,
316 –                  &m1Pass, &m2Min, &m2Max,
317 –                  THRCOF, &mSize, &error);
318 –        
319 –        for (int mmm = 0; mmm < (int)(m2Max - m2Min); mmm++) {
320 –          m2 = (int)floor(m2Min) + mmm;
321 –          m3 = -m1-m2;
322 –          W[l] += THRCOF[mmm] *
323 –            QBar[std::make_pair(l,m1)] *
324 –            QBar[std::make_pair(l,m2)] *
325 –            QBar[std::make_pair(l,m3)];
326 –        }
328        }
329        
330        W_hat[l] = W[l] / pow(Q2[l], 1.5);

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines