applications/staticProps/BOPofR.cpp

/*
 * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
 *
 * The University of Notre Dame grants you ("Licensee") a
 * non-exclusive, royalty free, license to use, modify and
 * redistribute this software in source and binary code form, provided
 * that the following conditions are met:
 *
 * 1. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 *
 * 2. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the
 *    distribution.
 *
 * This software is provided "AS IS," without a warranty of any
 * kind. All express or implied conditions, representations and
 * warranties, including any implied warranty of merchantability,
 * fitness for a particular purpose or non-infringement, are hereby
 * excluded.  The University of Notre Dame and its licensors shall not
 * be liable for any damages suffered by licensee as a result of
 * using, modifying or distributing the software or its
 * derivatives. In no event will the University of Notre Dame or its
 * licensors be liable for any lost revenue, profit or data, or for
 * direct, indirect, special, consequential, incidental or punitive
 * damages, however caused and regardless of the theory of liability,
 * arising out of the use of or inability to use software, even if the
 * University of Notre Dame has been advised of the possibility of
 * such damages.
 *
 * SUPPORT OPEN SCIENCE!  If you use OpenMD or its source code in your
 * research, please cite the appropriate papers when you publish your
 * work.  Good starting points are:
 *                                                                      
 * [1]  Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).             
 * [2]  Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).          
 * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).          
 * [4] Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010).
 * [4] , Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011). *
 *  Created by J. Daniel Gezelter on 09/26/06.
 *  @author  J. Daniel Gezelter
 *  @version $Id$
 *
 */
 
#include "applications/staticProps/BOPofR.hpp"
#include "utils/simError.h"
#include "io/DumpReader.hpp"
#include "primitives/Molecule.hpp"
#include "utils/NumericConstant.hpp"
#include "math/Wigner3jm.hpp"
#include "brains/Thermo.hpp"

using namespace MATPACK;
namespace OpenMD {

  BOPofR::BOPofR(SimInfo* info, const std::string& filename, const std::string& sele, double rCut,
                 int nbins, RealType len) : StaticAnalyser(info, filename), selectionScript_(sele), evaluator_(info), seleMan_(info){
    
    setOutputName(getPrefix(filename) + ".bo");

    evaluator_.loadScriptString(sele);
    if (!evaluator_.isDynamic()) {
      seleMan_.setSelectionSet(evaluator_.evaluate());
    }

    // Set up cutoff radius and order of the Legendre Polynomial:

    rCut_ = rCut;
    nBins_ = nbins;
    len_ = len;
        
    deltaR_ = len_/nBins_;
    RCount_.resize(nBins_);
    WofR_.resize(nBins_);
    QofR_.resize(nBins_);

        for (int i = 0; i < nBins_; i++){
                RCount_[i] = 0;
                WofR_[i] = 0;
                QofR_[i] = 0;
        }
        
    // Make arrays for Wigner3jm
    RealType* THRCOF = new RealType[2*lMax_+1];
    // Variables for Wigner routine
    RealType lPass, m1Pass, m2m, m2M;
    int error, mSize;
    mSize = 2*lMax_+1;

    for (int l = 0; l <= lMax_; l++) {
      lPass = (RealType)l;
      for (int m1 = -l; m1 <= l; m1++) {
        m1Pass = (RealType)m1;

        std::pair<int,int> lm = std::make_pair(l, m1);
        
        // Zero work array
        for (int ii = 0; ii < 2*l + 1; ii++){
          THRCOF[ii] = 0.0;
        }

        // Get Wigner coefficients
        Wigner3jm(lPass, lPass, lPass, 
                  m1Pass, m2m, m2M, 
                  THRCOF, mSize, error);
        
        m2Min[lm] = (int)floor(m2m);
        m2Max[lm] = (int)floor(m2M);
        
        for (int mmm = 0; mmm <= (int)(m2M - m2m); mmm++) {
          w3j[lm].push_back(THRCOF[mmm]);
        }
      }
    }

    delete [] THRCOF;
    THRCOF = NULL;      
        
  }
  
  BOPofR::~BOPofR() {
          /*
        std::cerr << "Freeing stuff" << std::endl;
    for (int l = 0; l <= lMax_; l++) {
      for (int m = -l; m <= l; m++) {
        w3j[std::make_pair(l,m)].clear();
      }
    }
        std::cerr << "w3j made free...." << std::endl;
   for (int bin = 0; bin < nBins_; bin++) {
                QofR_[bin].clear();
                WofR_[bin].clear();
                RCount_[bin].clear();
    }
        std::cout << "R arrays made free...." << std::endl;
   w3j.clear();
    m2Min.clear();
    m2Max.clear();
    RCount_.clear();
    WofR_.clear();
    QofR_.clear();
 */
  }

  
  void BOPofR::initializeHistogram() {
        for (int i = 0; i < nBins_; i++){
                RCount_[i] = 0;
                WofR_[i] = 0;
                QofR_[i] = 0;
        }
  }


  void BOPofR::process() {
    Molecule* mol;
    Atom* atom;
    RigidBody* rb;
    int myIndex;
    SimInfo::MoleculeIterator mi;
    Molecule::RigidBodyIterator rbIter;
    Molecule::AtomIterator ai;
    StuntDouble* sd;
    Vector3d vec;
    RealType costheta;
    RealType phi;
    RealType r;
    Vector3d rCOM;
    RealType distCOM;
    Vector3d pos;
    Vector3d CenterOfMass;
    std::map<std::pair<int,int>,ComplexType> q;
    std::vector<RealType> q_l;
    std::vector<RealType> q2;
    std::vector<ComplexType> w;
    std::vector<ComplexType> w_hat;
    std::map<std::pair<int,int>,ComplexType> QBar;
    std::vector<RealType> Q2;
    std::vector<RealType> Q;
    std::vector<ComplexType> W;
    std::vector<ComplexType> W_hat;
    int nBonds, Nbonds;
    SphericalHarmonic sphericalHarmonic;
    int i;
    
    DumpReader reader(info_, dumpFilename_);    
    int nFrames = reader.getNFrames();
    frameCounter_ = 0;

    Thermo thermo(info_);

    q_l.resize(lMax_+1);
    q2.resize(lMax_+1);
    w.resize(lMax_+1);
    w_hat.resize(lMax_+1);

    Q2.resize(lMax_+1);
    Q.resize(lMax_+1);
    W.resize(lMax_+1);
    W_hat.resize(lMax_+1);
    Nbonds = 0;

    for (int istep = 0; istep < nFrames; istep += step_) {
      reader.readFrame(istep);
      frameCounter_++;
      currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
      CenterOfMass = thermo.getCom();
      if (evaluator_.isDynamic()) {
        seleMan_.setSelectionSet(evaluator_.evaluate());
      }

      // update the positions of atoms which belong to the rigidbodies

      for (mol = info_->beginMolecule(mi); mol != NULL; 
           mol = info_->nextMolecule(mi)) {
        for (rb = mol->beginRigidBody(rbIter); rb != NULL; 
             rb = mol->nextRigidBody(rbIter)) {
          rb->updateAtoms();
        }        
      }           
            
      // outer loop is over the selected StuntDoubles:

      for (sd = seleMan_.beginSelected(i); sd != NULL; 
           sd = seleMan_.nextSelected(i)) {

        myIndex = sd->getGlobalIndex();
        
        nBonds = 0;
        
        for (int l = 0; l <= lMax_; l++) {
          for (int m = -l; m <= l; m++) {
            q[std::make_pair(l,m)] = 0.0;
          }
        }
        pos = sd->getPos();
        rCOM = CenterOfMass - pos;
        if (usePeriodicBoundaryConditions_) 
          currentSnapshot_->wrapVector(rCOM);
        distCOM = rCOM.length();

        // inner loop is over all other atoms in the system:
        
        for (mol = info_->beginMolecule(mi); mol != NULL; 
             mol = info_->nextMolecule(mi)) {
          for (atom = mol->beginAtom(ai); atom != NULL; 
               atom = mol->nextAtom(ai)) {

            if (atom->getGlobalIndex() != myIndex) {
              vec = pos - atom->getPos();       

              if (usePeriodicBoundaryConditions_) 
                currentSnapshot_->wrapVector(vec);
              
              // Calculate "bonds" and build Q_lm(r) where 
              //      Q_lm = Y_lm(theta(r),phi(r))                
              // The spherical harmonics are wrt any arbitrary coordinate
              // system, we choose standard spherical coordinates 
              
              r = vec.length();
              
              // Check to see if neighbor is in bond cutoff 
              
              if (r < rCut_) { 
                costheta = vec.z() / r; 
                phi = atan2(vec.y(), vec.x());

                for (int l = 0; l <= lMax_; l++) {
                  sphericalHarmonic.setL(l);
                  for(int m = -l; m <= l; m++){
                    sphericalHarmonic.setM(m);
                    q[std::make_pair(l,m)] += sphericalHarmonic.getValueAt(costheta, phi);
                  }
                }
                nBonds++;
              }  
            }
          }
        }
        
        
        for (int l = 0; l <= lMax_; l++) {
          q2[l] = 0.0;
          for (int m = -l; m <= l; m++){
            q[std::make_pair(l,m)] /= (RealType)nBonds;            
            q2[l] += norm(q[std::make_pair(l,m)]);
          }
          q_l[l] = sqrt(q2[l] * 4.0 * NumericConstant::PI / (RealType)(2*l + 1));
        }
        
        // Find Third Order Invariant W_l
    
        for (int l = 0; l <= lMax_; l++) {
          w[l] = 0.0;
          for (int m1 = -l; m1 <= l; m1++) {
            std::pair<int,int> lm = std::make_pair(l, m1);
            for (int mmm = 0; mmm <= (m2Max[lm] - m2Min[lm]); mmm++) {
              int m2 = m2Min[lm] + mmm;
              int m3 = -m1-m2;
              w[l] += w3j[lm][mmm] * q[lm] * 
                q[std::make_pair(l,m2)] *  q[std::make_pair(l,m3)];
            }
          }
          
          w_hat[l] = w[l] / pow(q2[l], RealType(1.5));
        }

        collectHistogram(q_l, w_hat, distCOM);
                
//              printf( "%s  %18.10g %18.10g %18.10g %18.10g \n", sd->getType().c_str(),pos[0],pos[1],pos[2],real(w_hat[6]));

      }
    }
        
    writeOrderParameter();    
  }

  void BOPofR::collectHistogram(std::vector<RealType> q, 
                                std::vector<ComplexType> what, RealType distCOM) {

    if ( distCOM < len_){
      // Figure out where this distance goes...
      int whichBin = distCOM / deltaR_;
      RCount_[whichBin]++;

      if(real(what[6]) < -0.15){                                
                WofR_[whichBin]++;
      }
          if(q[6] > 0.5){
        QofR_[whichBin]++;
          }
    }  

  }

  void BOPofR::writeOrderParameter() {
    
    std::ofstream osq((getOutputFileName() + "qr").c_str());

    if (osq.is_open()) {
      
      // Normalize by number of frames and write it out:
          
      for (int i = 0; i < nBins_; ++i) {
        RealType Rval = (i + 0.5) * deltaR_;               
        osq << Rval;
                if (RCount_[i] == 0){
                        osq << "\t" << 0;
                        osq << "\n";
                }else{
                osq << "\t" << (RealType)QofR_[i]/(RealType)RCount_[i];        
                osq << "\n";
                }
      }

      osq.close();

    } else {
      sprintf(painCave.errMsg, "BOPofR: unable to open %s\n", 
              (getOutputFileName() + "q").c_str());
      painCave.isFatal = 1;
      simError();  
    }

    std::ofstream osw((getOutputFileName() + "wr").c_str());

    if (osw.is_open()) {
      // Normalize by number of frames and write it out:
      for (int i = 0; i < nBins_; ++i) {
        RealType Rval = deltaR_ * (i + 0.5);               
        osw << Rval;
                if (RCount_[i] == 0){
                        osw << "\t" << 0;
                        osw << "\n";
                }else{
                osw << "\t" << (RealType)WofR_[i]/(RealType)RCount_[i];
                osw << "\n";
                }
      }

      osw.close();
    } else {
      sprintf(painCave.errMsg, "BOPofR: unable to open %s\n", 
              (getOutputFileName() + "w").c_str());
      painCave.isFatal = 1;
      simError();  
    
    }
        
  }
}
Revision:	1785
Committed:	Wed Aug 22 18:43:27 2012 UTC (12 years, 8 months ago) by jmichalk
File size:	11624 byte(s)
Log Message:	Trunk: The changes in this commit are confined to applications/staticProps and for the most part deal with a misspelling of initialize. The one other change took place in StaticProps.cpp and deals with the default treatment of sele2. It had previously been set to 'select all' which seems to go against what would be desired by not specifying it with regard to proper operations of many of the analysis programs ( g of r's especially)
#	User	Rev	Content
1	chuckv	1128	/*
2			* Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3			*
4			* The University of Notre Dame grants you ("Licensee") a
5			* non-exclusive, royalty free, license to use, modify and
6			* redistribute this software in source and binary code form, provided
7			* that the following conditions are met:
8			*
9	gezelter	1390	* 1. Redistributions of source code must retain the above copyright
10	chuckv	1128	* notice, this list of conditions and the following disclaimer.
11			*
12	gezelter	1390	* 2. Redistributions in binary form must reproduce the above copyright
13	chuckv	1128	* notice, this list of conditions and the following disclaimer in the
14			* documentation and/or other materials provided with the
15			* distribution.
16			*
17			* This software is provided "AS IS," without a warranty of any
18			* kind. All express or implied conditions, representations and
19			* warranties, including any implied warranty of merchantability,
20			* fitness for a particular purpose or non-infringement, are hereby
21			* excluded. The University of Notre Dame and its licensors shall not
22			* be liable for any damages suffered by licensee as a result of
23			* using, modifying or distributing the software or its
24			* derivatives. In no event will the University of Notre Dame or its
25			* licensors be liable for any lost revenue, profit or data, or for
26			* direct, indirect, special, consequential, incidental or punitive
27			* damages, however caused and regardless of the theory of liability,
28			* arising out of the use of or inability to use software, even if the
29			* University of Notre Dame has been advised of the possibility of
30			* such damages.
31			*
32	gezelter	1390	* SUPPORT OPEN SCIENCE! If you use OpenMD or its source code in your
33			* research, please cite the appropriate papers when you publish your
34			* work. Good starting points are:
35			*
36			* [1] Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).
37			* [2] Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).
38			* [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).
39	gezelter	1782	* [4] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010).
40			* [4] , Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011). *
41	chuckv	1128	* Created by J. Daniel Gezelter on 09/26/06.
42			* @author J. Daniel Gezelter
43	gezelter	1442	* @version $Id$
44	chuckv	1128	*
45			*/
46
47			#include "applications/staticProps/BOPofR.hpp"
48			#include "utils/simError.h"
49			#include "io/DumpReader.hpp"
50			#include "primitives/Molecule.hpp"
51			#include "utils/NumericConstant.hpp"
52	gezelter	1782	#include "math/Wigner3jm.hpp"
53			#include "brains/Thermo.hpp"
54	chuckv	1128
55	gezelter	1782	using namespace MATPACK;
56	gezelter	1390	namespace OpenMD {
57	chuckv	1128
58			BOPofR::BOPofR(SimInfo* info, const std::string& filename, const std::string& sele, double rCut,
59			int nbins, RealType len) : StaticAnalyser(info, filename), selectionScript_(sele), evaluator_(info), seleMan_(info){
60
61			setOutputName(getPrefix(filename) + ".bo");
62
63			evaluator_.loadScriptString(sele);
64			if (!evaluator_.isDynamic()) {
65			seleMan_.setSelectionSet(evaluator_.evaluate());
66			}
67
68			// Set up cutoff radius and order of the Legendre Polynomial:
69
70			rCut_ = rCut;
71			nBins_ = nbins;
72			len_ = len;
73
74			deltaR_ = len_/nBins_;
75			RCount_.resize(nBins_);
76			WofR_.resize(nBins_);
77			QofR_.resize(nBins_);
78	chuckv	1137
79			for (int i = 0; i < nBins_; i++){
80			RCount_[i] = 0;
81			WofR_[i] = 0;
82			QofR_[i] = 0;
83			}
84	chuckv	1128
85			// Make arrays for Wigner3jm
86	gezelter	1782	RealType* THRCOF = new RealType[2*lMax_+1];
87	chuckv	1128	// Variables for Wigner routine
88	gezelter	1782	RealType lPass, m1Pass, m2m, m2M;
89	chuckv	1128	int error, mSize;
90			mSize = 2*lMax_+1;
91
92			for (int l = 0; l <= lMax_; l++) {
93	gezelter	1782	lPass = (RealType)l;
94	chuckv	1128	for (int m1 = -l; m1 <= l; m1++) {
95	gezelter	1782	m1Pass = (RealType)m1;
96	chuckv	1128
97			std::pair<int,int> lm = std::make_pair(l, m1);
98
99			// Zero work array
100			for (int ii = 0; ii < 2*l + 1; ii++){
101			THRCOF[ii] = 0.0;
102			}
103
104			// Get Wigner coefficients
105	gezelter	1782	Wigner3jm(lPass, lPass, lPass,
106			m1Pass, m2m, m2M,
107			THRCOF, mSize, error);
108	chuckv	1128
109			m2Min[lm] = (int)floor(m2m);
110			m2Max[lm] = (int)floor(m2M);
111
112			for (int mmm = 0; mmm <= (int)(m2M - m2m); mmm++) {
113			w3j[lm].push_back(THRCOF[mmm]);
114			}
115			}
116			}
117
118			delete [] THRCOF;
119	chuckv	1137	THRCOF = NULL;
120	chuckv	1128
121			}
122
123			BOPofR::~BOPofR() {
124			/*
125			std::cerr << "Freeing stuff" << std::endl;
126			for (int l = 0; l <= lMax_; l++) {
127			for (int m = -l; m <= l; m++) {
128			w3j[std::make_pair(l,m)].clear();
129			}
130			}
131			std::cerr << "w3j made free...." << std::endl;
132			for (int bin = 0; bin < nBins_; bin++) {
133			QofR_[bin].clear();
134			WofR_[bin].clear();
135			RCount_[bin].clear();
136			}
137			std::cout << "R arrays made free...." << std::endl;
138			w3j.clear();
139			m2Min.clear();
140			m2Max.clear();
141			RCount_.clear();
142			WofR_.clear();
143			QofR_.clear();
144			*/
145			}
146
147
148	jmichalk	1785	void BOPofR::initializeHistogram() {
149	chuckv	1137	for (int i = 0; i < nBins_; i++){
150			RCount_[i] = 0;
151			WofR_[i] = 0;
152			QofR_[i] = 0;
153			}
154	chuckv	1128	}
155
156
157			void BOPofR::process() {
158			Molecule* mol;
159			Atom* atom;
160			RigidBody* rb;
161			int myIndex;
162			SimInfo::MoleculeIterator mi;
163			Molecule::RigidBodyIterator rbIter;
164			Molecule::AtomIterator ai;
165			StuntDouble* sd;
166			Vector3d vec;
167			RealType costheta;
168			RealType phi;
169			RealType r;
170			Vector3d rCOM;
171			RealType distCOM;
172			Vector3d pos;
173			Vector3d CenterOfMass;
174			std::map<std::pair<int,int>,ComplexType> q;
175			std::vector<RealType> q_l;
176			std::vector<RealType> q2;
177			std::vector<ComplexType> w;
178			std::vector<ComplexType> w_hat;
179			std::map<std::pair<int,int>,ComplexType> QBar;
180			std::vector<RealType> Q2;
181			std::vector<RealType> Q;
182			std::vector<ComplexType> W;
183			std::vector<ComplexType> W_hat;
184			int nBonds, Nbonds;
185			SphericalHarmonic sphericalHarmonic;
186	gezelter	1782	int i;
187
188	chuckv	1128	DumpReader reader(info_, dumpFilename_);
189			int nFrames = reader.getNFrames();
190			frameCounter_ = 0;
191
192	gezelter	1782	Thermo thermo(info_);
193
194	chuckv	1128	q_l.resize(lMax_+1);
195			q2.resize(lMax_+1);
196			w.resize(lMax_+1);
197			w_hat.resize(lMax_+1);
198
199			Q2.resize(lMax_+1);
200			Q.resize(lMax_+1);
201			W.resize(lMax_+1);
202			W_hat.resize(lMax_+1);
203			Nbonds = 0;
204
205			for (int istep = 0; istep < nFrames; istep += step_) {
206			reader.readFrame(istep);
207			frameCounter_++;
208			currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
209	gezelter	1782	CenterOfMass = thermo.getCom();
210	chuckv	1128	if (evaluator_.isDynamic()) {
211			seleMan_.setSelectionSet(evaluator_.evaluate());
212			}
213
214			// update the positions of atoms which belong to the rigidbodies
215
216			for (mol = info_->beginMolecule(mi); mol != NULL;
217			mol = info_->nextMolecule(mi)) {
218			for (rb = mol->beginRigidBody(rbIter); rb != NULL;
219			rb = mol->nextRigidBody(rbIter)) {
220			rb->updateAtoms();
221			}
222			}
223
224			// outer loop is over the selected StuntDoubles:
225
226			for (sd = seleMan_.beginSelected(i); sd != NULL;
227			sd = seleMan_.nextSelected(i)) {
228
229			myIndex = sd->getGlobalIndex();
230
231			nBonds = 0;
232
233			for (int l = 0; l <= lMax_; l++) {
234			for (int m = -l; m <= l; m++) {
235			q[std::make_pair(l,m)] = 0.0;
236			}
237			}
238			pos = sd->getPos();
239			rCOM = CenterOfMass - pos;
240			if (usePeriodicBoundaryConditions_)
241			currentSnapshot_->wrapVector(rCOM);
242			distCOM = rCOM.length();
243
244			// inner loop is over all other atoms in the system:
245
246			for (mol = info_->beginMolecule(mi); mol != NULL;
247			mol = info_->nextMolecule(mi)) {
248			for (atom = mol->beginAtom(ai); atom != NULL;
249			atom = mol->nextAtom(ai)) {
250
251			if (atom->getGlobalIndex() != myIndex) {
252			vec = pos - atom->getPos();
253
254			if (usePeriodicBoundaryConditions_)
255			currentSnapshot_->wrapVector(vec);
256
257			// Calculate "bonds" and build Q_lm(r) where
258			// Q_lm = Y_lm(theta(r),phi(r))
259			// The spherical harmonics are wrt any arbitrary coordinate
260			// system, we choose standard spherical coordinates
261
262			r = vec.length();
263
264			// Check to see if neighbor is in bond cutoff
265
266			if (r < rCut_) {
267			costheta = vec.z() / r;
268			phi = atan2(vec.y(), vec.x());
269
270			for (int l = 0; l <= lMax_; l++) {
271			sphericalHarmonic.setL(l);
272			for(int m = -l; m <= l; m++){
273			sphericalHarmonic.setM(m);
274			q[std::make_pair(l,m)] += sphericalHarmonic.getValueAt(costheta, phi);
275			}
276			}
277			nBonds++;
278			}
279			}
280			}
281			}
282
283
284			for (int l = 0; l <= lMax_; l++) {
285			q2[l] = 0.0;
286			for (int m = -l; m <= l; m++){
287			q[std::make_pair(l,m)] /= (RealType)nBonds;
288			q2[l] += norm(q[std::make_pair(l,m)]);
289			}
290			q_l[l] = sqrt(q2[l] * 4.0 * NumericConstant::PI / (RealType)(2*l + 1));
291			}
292
293			// Find Third Order Invariant W_l
294
295			for (int l = 0; l <= lMax_; l++) {
296			w[l] = 0.0;
297			for (int m1 = -l; m1 <= l; m1++) {
298			std::pair<int,int> lm = std::make_pair(l, m1);
299			for (int mmm = 0; mmm <= (m2Max[lm] - m2Min[lm]); mmm++) {
300			int m2 = m2Min[lm] + mmm;
301			int m3 = -m1-m2;
302			w[l] += w3j[lm][mmm] * q[lm] *
303			q[std::make_pair(l,m2)] * q[std::make_pair(l,m3)];
304			}
305			}
306
307	gezelter	1782	w_hat[l] = w[l] / pow(q2[l], RealType(1.5));
308	chuckv	1128	}
309
310			collectHistogram(q_l, w_hat, distCOM);
311	chuckv	1137
312	chuckv	1194	// printf( "%s %18.10g %18.10g %18.10g %18.10g \n", sd->getType().c_str(),pos[0],pos[1],pos[2],real(w_hat[6]));
313	chuckv	1137
314	chuckv	1128	}
315			}
316
317			writeOrderParameter();
318			}
319
320			void BOPofR::collectHistogram(std::vector<RealType> q,
321			std::vector<ComplexType> what, RealType distCOM) {
322
323			if ( distCOM < len_){
324			// Figure out where this distance goes...
325			int whichBin = distCOM / deltaR_;
326	chuckv	1137	RCount_[whichBin]++;
327
328			if(real(what[6]) < -0.15){
329			WofR_[whichBin]++;
330	chuckv	1128	}
331	chuckv	1137	if(q[6] > 0.5){
332			QofR_[whichBin]++;
333			}
334	chuckv	1128	}
335
336			}
337
338			void BOPofR::writeOrderParameter() {
339
340			std::ofstream osq((getOutputFileName() + "qr").c_str());
341
342			if (osq.is_open()) {
343
344			// Normalize by number of frames and write it out:
345
346			for (int i = 0; i < nBins_; ++i) {
347			RealType Rval = (i + 0.5) * deltaR_;
348			osq << Rval;
349	chuckv	1137	if (RCount_[i] == 0){
350			osq << "\t" << 0;
351			osq << "\n";
352			}else{
353			osq << "\t" << (RealType)QofR_[i]/(RealType)RCount_[i];
354			osq << "\n";
355			}
356	chuckv	1128	}
357
358			osq.close();
359
360			} else {
361			sprintf(painCave.errMsg, "BOPofR: unable to open %s\n",
362			(getOutputFileName() + "q").c_str());
363			painCave.isFatal = 1;
364			simError();
365			}
366
367			std::ofstream osw((getOutputFileName() + "wr").c_str());
368
369			if (osw.is_open()) {
370			// Normalize by number of frames and write it out:
371			for (int i = 0; i < nBins_; ++i) {
372			RealType Rval = deltaR_ * (i + 0.5);
373			osw << Rval;
374	chuckv	1137	if (RCount_[i] == 0){
375			osw << "\t" << 0;
376			osw << "\n";
377			}else{
378			osw << "\t" << (RealType)WofR_[i]/(RealType)RCount_[i];
379			osw << "\n";
380			}
381	chuckv	1128	}
382
383			osw.close();
384			} else {
385			sprintf(painCave.errMsg, "BOPofR: unable to open %s\n",
386			(getOutputFileName() + "w").c_str());
387			painCave.isFatal = 1;
388			simError();
389
390			}
391
392			}
393			}