applications/staticProps/BondOrderParameter.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/BondOrderParameter.hpp"
#include "utils/simError.h"
#include "io/DumpReader.hpp"
#include "primitives/Molecule.hpp"
#include "utils/NumericConstant.hpp"
#include "math/Wigner3jm.hpp"

using namespace MATPACK;
namespace OpenMD {

  BondOrderParameter::BondOrderParameter(SimInfo* info, 
                                         const std::string& filename, 
                                         const std::string& sele,
                                         double rCut, int nbins) : 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;
    Qcount_.resize(lMax_+1);
    Wcount_.resize(lMax_+1);

    // Q can take values from 0 to 1

    MinQ_ = 0.0;
    MaxQ_ = 1.1;
    deltaQ_ = (MaxQ_ - MinQ_) / nbins;

    // W_6 for icosahedral clusters is 11 / sqrt(4199) = 0.169754, so we'll
    // use values for MinW_ and MaxW_ that are slightly larger than this:

    MinW_ = -1.1;
    MaxW_ = 1.1;
    deltaW_ = (MaxW_ - MinW_) / nbins;

    // 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;
  }
  
  BondOrderParameter::~BondOrderParameter() {
    Q_histogram_.clear();
    W_histogram_.clear();
    for (int l = 0; l <= lMax_; l++) {
      for (int m = -l; m <= l; m++) {
        w3j[std::make_pair(l,m)].clear();
      }
    }
    w3j.clear();
    m2Min.clear();
    m2Max.clear();
  }
  
  void BondOrderParameter::initializeHistogram() {
    for (int bin = 0; bin < nBins_; bin++) {
      for (int l = 0; l <= lMax_; l++) {
        Q_histogram_[std::make_pair(bin,l)] = 0;
        W_histogram_[std::make_pair(bin,l)] = 0;
      }
    }
  }

  void BondOrderParameter::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;
    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;

    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();
      
      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;
          }
        }
        
        // 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 = sd->getPos() - 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);
        
        Nbonds += nBonds;
        for (int l = 0; l <= lMax_;  l++) {
          for (int m = -l; m <= l; m++) {
            QBar[std::make_pair(l,m)] += (RealType)nBonds*q[std::make_pair(l,m)];
          }
        }
      }
    }
      
    // Normalize Qbar2
    for (int l = 0; l <= lMax_; l++) {
      for (int m = -l; m <= l; m++){
        QBar[std::make_pair(l,m)] /= Nbonds;
      }
    }
    
    // Find second order invariant Q_l
    
    for (int l = 0; l <= lMax_; l++) {
      Q2[l] = 0.0;
      for (int m = -l; m <= l; m++){
        Q2[l] += norm(QBar[std::make_pair(l,m)]);
      }
      Q[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] * QBar[lm] * 
            QBar[std::make_pair(l,m2)] * QBar[std::make_pair(l,m3)];
        }
      }
      
      W_hat[l] = W[l] / pow(Q2[l], RealType(1.5));
    }
    
    writeOrderParameter(Q, W_hat);    
  }

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

    for (int l = 0; l <= lMax_; l++) {
      if (q[l] >= MinQ_ && q[l] < MaxQ_) {
        int qbin = (q[l] - MinQ_) / deltaQ_;
        Q_histogram_[std::make_pair(qbin,l)] += 1;
        Qcount_[l]++;      
      } else {
        sprintf( painCave.errMsg,
                 "q_l value outside reasonable range\n");
        painCave.severity = OPENMD_ERROR;
        painCave.isFatal = 1;
        simError();  
      }
    }

    for (int l = 0; l <= lMax_; l++) {
      if (real(what[l]) >= MinW_ && real(what[l]) < MaxW_) {
        int wbin = (real(what[l]) - MinW_) / deltaW_;
        W_histogram_[std::make_pair(wbin,l)] += 1;
        Wcount_[l]++;      
      } else {
        sprintf( painCave.errMsg,
                 "Re[w_hat] value (%lf) outside reasonable range\n", real(what[l]));
        painCave.severity = OPENMD_ERROR;
        painCave.isFatal = 1;
        simError();  
      }
    }

  }  


  void BondOrderParameter::writeOrderParameter(std::vector<RealType> Q, 
                                               std::vector<ComplexType> What) {
    
    std::ofstream osq((getOutputFileName() + "q").c_str());

    if (osq.is_open()) {
      
      osq << "# Bond Order Parameters\n";
      osq << "# selection: (" << selectionScript_ << ")\n";
      osq << "# \n";
      for (int l = 0; l <= lMax_; l++) {
        osq << "# <Q_" << l << ">: " << Q[l] << "\n";
      }
      // Normalize by number of frames and write it out:
      for (int i = 0; i < nBins_; ++i) {
        RealType Qval = MinQ_ + (i + 0.5) * deltaQ_;               
        osq << Qval;
        for (int l = 0; l <= lMax_; l++) {

          osq << "\t" << (RealType)Q_histogram_[std::make_pair(i,l)]/(RealType)Qcount_[l]/deltaQ_;
        }
        osq << "\n";
      }

      osq.close();

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

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

    if (osw.is_open()) {
      osw << "# Bond Order Parameters\n";
      osw << "# selection: (" << selectionScript_ << ")\n";
      osw << "# \n";
      for (int l = 0; l <= lMax_; l++) {
        osw << "# <W_" << l << ">: " << real(What[l]) << "\t" << imag(What[l]) << "\n";
      }
      // Normalize by number of frames and write it out:
      for (int i = 0; i < nBins_; ++i) {
        RealType Wval = MinW_ + (i + 0.5) * deltaW_;               
        osw << Wval;
        for (int l = 0; l <= lMax_; l++) {

          osw << "\t" << (RealType)W_histogram_[std::make_pair(i,l)]/(RealType)Wcount_[l]/deltaW_;
        }
        osw << "\n";
      }

      osw.close();
    } else {
      sprintf(painCave.errMsg, "BondOrderParameter: 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:	13891 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	1038	/*
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	1038	* 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	1038	* 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	gezelter	1054	*
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	gezelter	1054	* Created by J. Daniel Gezelter on 09/26/06.
42			* @author J. Daniel Gezelter
43	gezelter	1442	* @version $Id$
44	gezelter	1054	*
45	chuckv	1038	*/
46	gezelter	1039
47	chuckv	1038	#include "applications/staticProps/BondOrderParameter.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	gezelter	1041
54	gezelter	1782	using namespace MATPACK;
55	gezelter	1390	namespace OpenMD {
56	chuckv	1038
57	gezelter	1039	BondOrderParameter::BondOrderParameter(SimInfo* info,
58			const std::string& filename,
59			const std::string& sele,
60	gezelter	1052	double rCut, int nbins) : StaticAnalyser(info, filename), selectionScript_(sele), evaluator_(info), seleMan_(info){
61	gezelter	1039
62	gezelter	1041	setOutputName(getPrefix(filename) + ".bo");
63	chuckv	1038
64	gezelter	1039	evaluator_.loadScriptString(sele);
65			if (!evaluator_.isDynamic()) {
66			seleMan_.setSelectionSet(evaluator_.evaluate());
67	chuckv	1038	}
68
69	gezelter	1039	// Set up cutoff radius and order of the Legendre Polynomial:
70
71	chuckv	1038	rCut_ = rCut;
72	gezelter	1051	nBins_ = nbins;
73			Qcount_.resize(lMax_+1);
74			Wcount_.resize(lMax_+1);
75	gezelter	1048
76			// Q can take values from 0 to 1
77
78			MinQ_ = 0.0;
79	gezelter	1051	MaxQ_ = 1.1;
80	gezelter	1048	deltaQ_ = (MaxQ_ - MinQ_) / nbins;
81
82			// W_6 for icosahedral clusters is 11 / sqrt(4199) = 0.169754, so we'll
83			// use values for MinW_ and MaxW_ that are slightly larger than this:
84
85	gezelter	1100	MinW_ = -1.1;
86			MaxW_ = 1.1;
87	gezelter	1048	deltaW_ = (MaxW_ - MinW_) / nbins;
88	gezelter	1053
89			// Make arrays for Wigner3jm
90	gezelter	1782	RealType* THRCOF = new RealType[2*lMax_+1];
91	gezelter	1053	// Variables for Wigner routine
92	gezelter	1782	RealType lPass, m1Pass, m2m, m2M;
93	gezelter	1053	int error, mSize;
94			mSize = 2*lMax_+1;
95
96			for (int l = 0; l <= lMax_; l++) {
97	gezelter	1782	lPass = (RealType)l;
98	gezelter	1053	for (int m1 = -l; m1 <= l; m1++) {
99	gezelter	1782	m1Pass = (RealType)m1;
100	gezelter	1053
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	gezelter	1100
108	gezelter	1053	// Get Wigner coefficients
109	gezelter	1782	Wigner3jm(lPass, lPass, lPass,
110			m1Pass, m2m, m2M,
111			THRCOF, mSize, error);
112	gezelter	1610
113	gezelter	1053	m2Min[lm] = (int)floor(m2m);
114			m2Max[lm] = (int)floor(m2M);
115
116	gezelter	1100	for (int mmm = 0; mmm <= (int)(m2M - m2m); mmm++) {
117	gezelter	1053	w3j[lm].push_back(THRCOF[mmm]);
118			}
119			}
120			}
121	gezelter	1054	delete [] THRCOF;
122			THRCOF = NULL;
123	gezelter	1039	}
124	gezelter	1053
125	gezelter	1041	BondOrderParameter::~BondOrderParameter() {
126	gezelter	1048	Q_histogram_.clear();
127			W_histogram_.clear();
128	gezelter	1054	for (int l = 0; l <= lMax_; l++) {
129			for (int m = -l; m <= l; m++) {
130			w3j[std::make_pair(l,m)].clear();
131			}
132			}
133			w3j.clear();
134			m2Min.clear();
135			m2Max.clear();
136	gezelter	1041	}
137	gezelter	1054
138	jmichalk	1785	void BondOrderParameter::initializeHistogram() {
139	gezelter	1051	for (int bin = 0; bin < nBins_; bin++) {
140			for (int l = 0; l <= lMax_; l++) {
141			Q_histogram_[std::make_pair(bin,l)] = 0;
142			W_histogram_[std::make_pair(bin,l)] = 0;
143			}
144			}
145	gezelter	1048	}
146
147	gezelter	1039	void BondOrderParameter::process() {
148			Molecule* mol;
149			Atom* atom;
150			RigidBody* rb;
151	gezelter	1043	int myIndex;
152	gezelter	1039	SimInfo::MoleculeIterator mi;
153			Molecule::RigidBodyIterator rbIter;
154			Molecule::AtomIterator ai;
155			StuntDouble* sd;
156	gezelter	1043	Vector3d vec;
157	gezelter	1042	RealType costheta;
158	gezelter	1039	RealType phi;
159			RealType r;
160	gezelter	1051	std::map<std::pair<int,int>,ComplexType> q;
161			std::vector<RealType> q_l;
162	gezelter	1052	std::vector<RealType> q2;
163			std::vector<ComplexType> w;
164			std::vector<ComplexType> w_hat;
165	gezelter	1051	std::map<std::pair<int,int>,ComplexType> QBar;
166			std::vector<RealType> Q2;
167			std::vector<RealType> Q;
168			std::vector<ComplexType> W;
169			std::vector<ComplexType> W_hat;
170	gezelter	1048	int nBonds, Nbonds;
171	gezelter	1042	SphericalHarmonic sphericalHarmonic;
172	gezelter	1782	int i;
173	gezelter	1043
174	gezelter	1039	DumpReader reader(info_, dumpFilename_);
175			int nFrames = reader.getNFrames();
176	gezelter	1041	frameCounter_ = 0;
177	chuckv	1038
178	gezelter	1051	q_l.resize(lMax_+1);
179	gezelter	1052	q2.resize(lMax_+1);
180			w.resize(lMax_+1);
181			w_hat.resize(lMax_+1);
182
183	gezelter	1051	Q2.resize(lMax_+1);
184			Q.resize(lMax_+1);
185			W.resize(lMax_+1);
186			W_hat.resize(lMax_+1);
187	gezelter	1094	Nbonds = 0;
188	gezelter	1051
189	gezelter	1039	for (int istep = 0; istep < nFrames; istep += step_) {
190			reader.readFrame(istep);
191	gezelter	1041	frameCounter_++;
192	gezelter	1039	currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
193
194			if (evaluator_.isDynamic()) {
195			seleMan_.setSelectionSet(evaluator_.evaluate());
196			}
197	chuckv	1038
198	gezelter	1039	// update the positions of atoms which belong to the rigidbodies
199	chuckv	1038
200	gezelter	1039	for (mol = info_->beginMolecule(mi); mol != NULL;
201			mol = info_->nextMolecule(mi)) {
202			for (rb = mol->beginRigidBody(rbIter); rb != NULL;
203			rb = mol->nextRigidBody(rbIter)) {
204			rb->updateAtoms();
205			}
206	gezelter	1048	}
207
208	gezelter	1039	// outer loop is over the selected StuntDoubles:
209	chuckv	1038
210	gezelter	1039	for (sd = seleMan_.beginSelected(i); sd != NULL;
211			sd = seleMan_.nextSelected(i)) {
212	chuckv	1038
213	gezelter	1043	myIndex = sd->getGlobalIndex();
214	gezelter	1048	nBonds = 0;
215	gezelter	1051
216			for (int l = 0; l <= lMax_; l++) {
217			for (int m = -l; m <= l; m++) {
218			q[std::make_pair(l,m)] = 0.0;
219			}
220	gezelter	1048	}
221	gezelter	1039
222			// inner loop is over all other atoms in the system:
223
224			for (mol = info_->beginMolecule(mi); mol != NULL;
225			mol = info_->nextMolecule(mi)) {
226			for (atom = mol->beginAtom(ai); atom != NULL;
227			atom = mol->nextAtom(ai)) {
228	chuckv	1038
229	gezelter	1043	if (atom->getGlobalIndex() != myIndex) {
230	chuckv	1038
231	gezelter	1043	vec = sd->getPos() - atom->getPos();
232	gezelter	1078
233			if (usePeriodicBoundaryConditions_)
234			currentSnapshot_->wrapVector(vec);
235	gezelter	1042
236	gezelter	1043	// Calculate "bonds" and build Q_lm(r) where
237			// Q_lm = Y_lm(theta(r),phi(r))
238			// The spherical harmonics are wrt any arbitrary coordinate
239			// system, we choose standard spherical coordinates
240
241			r = vec.length();
242
243			// Check to see if neighbor is in bond cutoff
244
245			if (r < rCut_) {
246			costheta = vec.z() / r;
247			phi = atan2(vec.y(), vec.x());
248	gezelter	1051
249			for (int l = 0; l <= lMax_; l++) {
250			sphericalHarmonic.setL(l);
251			for(int m = -l; m <= l; m++){
252			sphericalHarmonic.setM(m);
253			q[std::make_pair(l,m)] += sphericalHarmonic.getValueAt(costheta, phi);
254	gezelter	1610
255	gezelter	1051	}
256	gezelter	1043	}
257			nBonds++;
258			}
259			}
260	gezelter	1039	}
261			}
262	gezelter	1051
263
264	gezelter	1052	for (int l = 0; l <= lMax_; l++) {
265			q2[l] = 0.0;
266			for (int m = -l; m <= l; m++){
267	gezelter	1782	q[std::make_pair(l,m)] /= (RealType)nBonds;
268	gezelter	1610
269	gezelter	1052	q2[l] += norm(q[std::make_pair(l,m)]);
270			}
271	gezelter	1094	q_l[l] = sqrt(q2[l] * 4.0 * NumericConstant::PI / (RealType)(2*l + 1));
272	gezelter	1052	}
273	gezelter	1094
274	gezelter	1052	// Find Third Order Invariant W_l
275
276			for (int l = 0; l <= lMax_; l++) {
277			w[l] = 0.0;
278			for (int m1 = -l; m1 <= l; m1++) {
279	gezelter	1053	std::pair<int,int> lm = std::make_pair(l, m1);
280	gezelter	1100	for (int mmm = 0; mmm <= (m2Max[lm] - m2Min[lm]); mmm++) {
281	gezelter	1053	int m2 = m2Min[lm] + mmm;
282			int m3 = -m1-m2;
283			w[l] += w3j[lm][mmm] * q[lm] *
284			q[std::make_pair(l,m2)] * q[std::make_pair(l,m3)];
285	gezelter	1052	}
286			}
287
288	gezelter	1782	w_hat[l] = w[l] / pow(q2[l], RealType(1.5));
289	gezelter	1052	}
290
291			collectHistogram(q_l, w_hat);
292
293	gezelter	1048	Nbonds += nBonds;
294	gezelter	1051	for (int l = 0; l <= lMax_; l++) {
295			for (int m = -l; m <= l; m++) {
296	gezelter	1094	QBar[std::make_pair(l,m)] += (RealType)nBonds*q[std::make_pair(l,m)];
297	gezelter	1051	}
298	gezelter	1048	}
299			}
300	gezelter	1047	}
301	gezelter	1051
302	gezelter	1047	// Normalize Qbar2
303	gezelter	1051	for (int l = 0; l <= lMax_; l++) {
304			for (int m = -l; m <= l; m++){
305			QBar[std::make_pair(l,m)] /= Nbonds;
306			}
307	gezelter	1039	}
308
309	gezelter	1047	// Find second order invariant Q_l
310	gezelter	1039
311	gezelter	1051	for (int l = 0; l <= lMax_; l++) {
312			Q2[l] = 0.0;
313			for (int m = -l; m <= l; m++){
314			Q2[l] += norm(QBar[std::make_pair(l,m)]);
315			}
316			Q[l] = sqrt(Q2[l] * 4.0 * NumericConstant::PI / (RealType)(2*l + 1));
317	gezelter	1039	}
318	gezelter	1047
319			// Find Third Order Invariant W_l
320
321	gezelter	1051	for (int l = 0; l <= lMax_; l++) {
322			W[l] = 0.0;
323			for (int m1 = -l; m1 <= l; m1++) {
324	gezelter	1053	std::pair<int,int> lm = std::make_pair(l, m1);
325	gezelter	1100	for (int mmm = 0; mmm <= (m2Max[lm] - m2Min[lm]); mmm++) {
326	gezelter	1053	int m2 = m2Min[lm] + mmm;
327			int m3 = -m1-m2;
328			W[l] += w3j[lm][mmm] * QBar[lm] *
329			QBar[std::make_pair(l,m2)] * QBar[std::make_pair(l,m3)];
330	gezelter	1051	}
331	gezelter	1047	}
332	gezelter	1041
333	gezelter	1782	W_hat[l] = W[l] / pow(Q2[l], RealType(1.5));
334	gezelter	1039	}
335	gezelter	1047
336	gezelter	1051	writeOrderParameter(Q, W_hat);
337	gezelter	1047	}
338	chuckv	1038
339	gezelter	1052	void BondOrderParameter::collectHistogram(std::vector<RealType> q,
340			std::vector<ComplexType> what) {
341	chuckv	1038
342	gezelter	1051	for (int l = 0; l <= lMax_; l++) {
343			if (q[l] >= MinQ_ && q[l] < MaxQ_) {
344			int qbin = (q[l] - MinQ_) / deltaQ_;
345			Q_histogram_[std::make_pair(qbin,l)] += 1;
346			Qcount_[l]++;
347			} else {
348			sprintf( painCave.errMsg,
349			"q_l value outside reasonable range\n");
350	gezelter	1390	painCave.severity = OPENMD_ERROR;
351	gezelter	1051	painCave.isFatal = 1;
352			simError();
353			}
354	gezelter	1048	}
355
356	gezelter	1052	for (int l = 0; l <= lMax_; l++) {
357			if (real(what[l]) >= MinW_ && real(what[l]) < MaxW_) {
358			int wbin = (real(what[l]) - MinW_) / deltaW_;
359			W_histogram_[std::make_pair(wbin,l)] += 1;
360			Wcount_[l]++;
361			} else {
362			sprintf( painCave.errMsg,
363	gezelter	1094	"Re[w_hat] value (%lf) outside reasonable range\n", real(what[l]));
364	gezelter	1390	painCave.severity = OPENMD_ERROR;
365	gezelter	1052	painCave.isFatal = 1;
366			simError();
367			}
368			}
369
370	gezelter	1048	}
371
372
373	gezelter	1052	void BondOrderParameter::writeOrderParameter(std::vector<RealType> Q,
374			std::vector<ComplexType> What) {
375	gezelter	1051
376	gezelter	1052	std::ofstream osq((getOutputFileName() + "q").c_str());
377
378			if (osq.is_open()) {
379	gezelter	1041
380	gezelter	1052	osq << "# Bond Order Parameters\n";
381			osq << "# selection: (" << selectionScript_ << ")\n";
382			osq << "# \n";
383	gezelter	1051	for (int l = 0; l <= lMax_; l++) {
384	gezelter	1052	osq << "# <Q_" << l << ">: " << Q[l] << "\n";
385	gezelter	1051	}
386	gezelter	1048	// Normalize by number of frames and write it out:
387	gezelter	1051	for (int i = 0; i < nBins_; ++i) {
388			RealType Qval = MinQ_ + (i + 0.5) * deltaQ_;
389	gezelter	1052	osq << Qval;
390	gezelter	1051	for (int l = 0; l <= lMax_; l++) {
391	gezelter	1094
392			osq << "\t" << (RealType)Q_histogram_[std::make_pair(i,l)]/(RealType)Qcount_[l]/deltaQ_;
393	gezelter	1051	}
394	gezelter	1052	osq << "\n";
395	gezelter	1048	}
396
397	gezelter	1052	osq.close();
398	gezelter	1047
399	gezelter	1041	} else {
400			sprintf(painCave.errMsg, "BondOrderParameter: unable to open %s\n",
401	gezelter	1052	(getOutputFileName() + "q").c_str());
402	gezelter	1041	painCave.isFatal = 1;
403			simError();
404			}
405	gezelter	1052
406			std::ofstream osw((getOutputFileName() + "w").c_str());
407
408			if (osw.is_open()) {
409			osw << "# Bond Order Parameters\n";
410			osw << "# selection: (" << selectionScript_ << ")\n";
411			osw << "# \n";
412			for (int l = 0; l <= lMax_; l++) {
413	gezelter	1090	osw << "# <W_" << l << ">: " << real(What[l]) << "\t" << imag(What[l]) << "\n";
414	gezelter	1052	}
415			// Normalize by number of frames and write it out:
416			for (int i = 0; i < nBins_; ++i) {
417			RealType Wval = MinW_ + (i + 0.5) * deltaW_;
418			osw << Wval;
419			for (int l = 0; l <= lMax_; l++) {
420	gezelter	1094
421			osw << "\t" << (RealType)W_histogram_[std::make_pair(i,l)]/(RealType)Wcount_[l]/deltaW_;
422	gezelter	1052	}
423			osw << "\n";
424			}
425
426			osw.close();
427			} else {
428			sprintf(painCave.errMsg, "BondOrderParameter: unable to open %s\n",
429			(getOutputFileName() + "w").c_str());
430			painCave.isFatal = 1;
431			simError();
432			}
433
434	gezelter	1041	}
435	gezelter	1039	}