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. Acknowledgement of the program authors must be made in any
 *    publication of scientific results based in part on use of the
 *    program.  An acceptable form of acknowledgement is citation of
 *    the article in which the program was described (Matthew
 *    A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
 *    J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
 *    Parallel Simulation Engine for Molecular Dynamics,"
 *    J. Comput. Chem. 26, pp. 252-271 (2005))
 *
 * 2. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 *
 * 3. 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.
 */


/* Creates orientational bond order parameters as outlined by
 *     Bond-orientaional order in liquids and glasses, Steinhart,Nelson,Ronchetti
 *     Phys Rev B, 28,784,1983
 * 
 */

#include "applications/staticProps/BondOrderParameter.hpp"
#include "utils/simError.h"
#include "io/DumpReader.hpp"
#include "primitives/Molecule.hpp"
#include "utils/NumericConstant.hpp"
#include "math/RealSphericalHarmonic.hpp"
namespace oopse {


  BondOrderParameter::BondOrderParameter(SimInfo* info, const std::string& filename, const std::string& sele1,
                                         const std::string& sele2, double rCut, int lNumber)
    : StaticAnalyser(info, filename),
      selectionScript1_(sele1), evaluator1_(info),
      seleMan1_(info){

    setOutputName(getPrefix(filename) + ".obo");

    evaluator1_.loadScriptString(sele1);
    evaluator2_.loadScriptString(sele2);

    if (!evaluator1_.isDynamic()) {
      seleMan1_.setSelectionSet(evaluator1_.evaluate());
    }else {
      sprintf( painCave.errMsg,
               "--sele1 must be static selection\n");
      painCave.severity = OOPSE_ERROR;
      painCave.isFatal = 1;
      simError();
    }

    /* Set up cutoff radius and type of order parameter we are calcuating*/
    lNumber_ = lNumber;
    rCut_ = rCut;
    mSize_ = 2*lNumber_+1;

    int i;
    int j;
    StuntDouble* sd1;
    StuntDouble* sd2;
    for (sd1 = seleMan1_.beginSelected(i), sd2 = seleMan1_.beginSelected(j);
         sd1 != NULL && sd2 != NULL;
         sd1 = seleMan1_.nextSelected(i), sd2 = seleMan2_.nextSelected(j)) {
      for (sd2 = seleMan1_.beginSelected(j),sd2
             sdPairs_.push_back(std::make_pair(sd1, sd2));
           }


    }

    void BondOrderParameter::process
      () {
      Molecule* mol;
      RigidBody* rb;
      SimInfo::MoleculeIterator mi;
      Molecule::RigidBodyIterator rbIter;
      RealType theta;
      RealType phi;
      RealType r;
      RealType dist;
      RealType* QBar_lm;
      RealType QSq_l;
      int nBonds;
      int m, m_index;
      RealSphericalHarmonic sphericalHarmonic;


      DumpReader reader(info_, dumpFilename_);
      int nFrames = reader.getNFrames();

      /*Set the l for the spherical harmonic, it doesn't change*/
      sphericalHarmonic.setL(lNumber_);

      for (int i = 0; i < nFrames; i += step_) {
        reader.readFrame(i);
        currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
        nBonds = 0;

        for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
          //change the positions of atoms which belong to the rigidbodies
          for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) {
            rb->updateAtoms();
          }

        }


        /* Setup QBar */
        QBar_lm = new double[mSize_];

        /* Calculate "bonds" and build Q_lm(r) where Q_lm = Y_lm(theta(r),phi(r)) */
        for (std::vector<std::pair<StuntDouble*, StuntDouble*> >::iterator j = sdPairs_.begin(); j != sdPairs_.end(); ++j) {
          Vector3d vec = j->first->getPos() - j->second->getPos();
          currentSnapshot_->wrapVector(vec);
          /* The spherical harmonics are wrt any arbitray coordiate sysetm,
           * we choose standard spherical coordinates */
          r = sqrt(pow(vec.x(),2)+pow(vec.y(),2)+pow(vec.z(),2));

          /* Check to see if neighbor is in bond cuttoff*/
          if (r<rCut_){
            theta = atan(vec.y()/vec.x());
            phi = acos(vec.z()/r);
            for(int m_index = 0; m_index < mSize_; m_index++){
              sphericalHarmonic.setM(m_index-lNumber_);
              QBar_lm(m_index) += sphericalHarmonic.getValueAt(theta,phi);
            }
            nBonds++;
          }
        }

        /*Normalize Qbar by number of Bonds*/
        for ( int m_index = 0;m_index < mSize_; m_index++){
          QBar_lm(m_index) = QBar_lm(m_index)/nBonds;
        }


      }

      /*Normalize by number of frames*/
      for ( int m_index = 0;m_index < mSize_; m_index++){
        QBar_lm(m_index) = QBar_lm(m_index)/nFrames;
      }


      /* Find second order invariant Q_l*/

      for (int m_index = 0 ;m_index <= sizeM_; m++){
        QSq_l += pow(QBar_lm(m),2);
      }
      Q_l_ = sqrt((4*NumericConstant::PI/lNumber_+1)*QSq_l);

      /* Find Third Order Invariant W_l*/

      /* Make arrays for Wigner3jm */
      double* THRCOF = new double[mSize_];
      /* Variables for Wigner routine */
      double l_ = (double)lNumber_;
      double m2Min;
      double m2Max;
      int error;
      int m1;
      int m2;
      int m3;

      for (int m1 = -lNumber_;m <= lNumber_;m1++){
        /* Zero work array */
        for (i=0; i<mSize_;i++){
          THRCOF[i] = 0.0;      
        }
        /* Get wigner coefficients */
        Wigner3jm(&l_,&l_,&l_,&(double)m1,&m2Min,&m2Max&,THRCOF,&mSize_,&error);
        for (m_index=1; i<(m2Max-M2Min-1.0);m_index++){
          m2 = floor(m2Min) + m_index - 1;
          m3 = -m1-m2;
          W_l_ += THRCOF(m_index)*QBar_lm(m1+lNumber_)*QBar_lm(m2+lNumber_)*QBar_lm(m3+lNumber_);
        }
      }


      writeOrderParameter();

    }


    void BondOrderParameter::writeOrderParameter() {

      std::ofstream os(getOutputFileName().c_str());
      os << "#radial distribution function\n";
      os<< "#selection1: (" << selectionScript1_ << ")\t";
      os << "selection2: (" << selectionScript2_ << ")\n";
      os << "#p2\tdirector_x\tdirector_y\tdiretor_z\tangle(degree)\n";

      for (std::size_t i = 0; i < orderParams_.size(); ++i) {
        os <<  orderParams_[i].p2 << "\t"
           <<  orderParams_[i].director[0] << "\t"
           <<  orderParams_[i].director[1] << "\t"
           <<  orderParams_[i].director[2] << "\t"
           <<  orderParams_[i].angle << "\n";

      }

    }

  }

Revision:	3005
Committed:	Mon Sep 18 21:31:23 2006 UTC (17 years, 9 months ago) by chuckv
File size:	7493 byte(s)
Log Message:	Latest Version
#	Content
1	/*
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	* 1. Acknowledgement of the program authors must be made in any
10	* publication of scientific results based in part on use of the
11	* program. An acceptable form of acknowledgement is citation of
12	* the article in which the program was described (Matthew
13	* A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
14	* J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
15	* Parallel Simulation Engine for Molecular Dynamics,"
16	* J. Comput. Chem. 26, pp. 252-271 (2005))
17	*
18	* 2. Redistributions of source code must retain the above copyright
19	* notice, this list of conditions and the following disclaimer.
20	*
21	* 3. Redistributions in binary form must reproduce the above copyright
22	* notice, this list of conditions and the following disclaimer in the
23	* documentation and/or other materials provided with the
24	* distribution.
25	*
26	* This software is provided "AS IS," without a warranty of any
27	* kind. All express or implied conditions, representations and
28	* warranties, including any implied warranty of merchantability,
29	* fitness for a particular purpose or non-infringement, are hereby
30	* excluded. The University of Notre Dame and its licensors shall not
31	* be liable for any damages suffered by licensee as a result of
32	* using, modifying or distributing the software or its
33	* derivatives. In no event will the University of Notre Dame or its
34	* licensors be liable for any lost revenue, profit or data, or for
35	* direct, indirect, special, consequential, incidental or punitive
36	* damages, however caused and regardless of the theory of liability,
37	* arising out of the use of or inability to use software, even if the
38	* University of Notre Dame has been advised of the possibility of
39	* such damages.
40	*/
41
42
43	/* Creates orientational bond order parameters as outlined by
44	* Bond-orientaional order in liquids and glasses, Steinhart,Nelson,Ronchetti
45	* Phys Rev B, 28,784,1983
46	*
47	*/
48
49	#include "applications/staticProps/BondOrderParameter.hpp"
50	#include "utils/simError.h"
51	#include "io/DumpReader.hpp"
52	#include "primitives/Molecule.hpp"
53	#include "utils/NumericConstant.hpp"
54	#include "math/RealSphericalHarmonic.hpp"
55	namespace oopse {
56
57
58	BondOrderParameter::BondOrderParameter(SimInfo* info, const std::string& filename, const std::string& sele1,
59	const std::string& sele2, double rCut, int lNumber)
60	: StaticAnalyser(info, filename),
61	selectionScript1_(sele1), evaluator1_(info),
62	seleMan1_(info){
63
64	setOutputName(getPrefix(filename) + ".obo");
65
66	evaluator1_.loadScriptString(sele1);
67	evaluator2_.loadScriptString(sele2);
68
69	if (!evaluator1_.isDynamic()) {
70	seleMan1_.setSelectionSet(evaluator1_.evaluate());
71	}else {
72	sprintf( painCave.errMsg,
73	"--sele1 must be static selection\n");
74	painCave.severity = OOPSE_ERROR;
75	painCave.isFatal = 1;
76	simError();
77	}
78
79	/* Set up cutoff radius and type of order parameter we are calcuating*/
80	lNumber_ = lNumber;
81	rCut_ = rCut;
82	mSize_ = 2*lNumber_+1;
83
84	int i;
85	int j;
86	StuntDouble* sd1;
87	StuntDouble* sd2;
88	for (sd1 = seleMan1_.beginSelected(i), sd2 = seleMan1_.beginSelected(j);
89	sd1 != NULL && sd2 != NULL;
90	sd1 = seleMan1_.nextSelected(i), sd2 = seleMan2_.nextSelected(j)) {
91	for (sd2 = seleMan1_.beginSelected(j),sd2
92	sdPairs_.push_back(std::make_pair(sd1, sd2));
93	}
94
95
96	}
97
98	void BondOrderParameter::process
99	() {
100	Molecule* mol;
101	RigidBody* rb;
102	SimInfo::MoleculeIterator mi;
103	Molecule::RigidBodyIterator rbIter;
104	RealType theta;
105	RealType phi;
106	RealType r;
107	RealType dist;
108	RealType* QBar_lm;
109	RealType QSq_l;
110	int nBonds;
111	int m, m_index;
112	RealSphericalHarmonic sphericalHarmonic;
113
114
115	DumpReader reader(info_, dumpFilename_);
116	int nFrames = reader.getNFrames();
117
118	/Set the l for the spherical harmonic, it doesn't change/
119	sphericalHarmonic.setL(lNumber_);
120
121	for (int i = 0; i < nFrames; i += step_) {
122	reader.readFrame(i);
123	currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
124	nBonds = 0;
125
126	for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
127	//change the positions of atoms which belong to the rigidbodies
128	for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) {
129	rb->updateAtoms();
130	}
131
132	}
133
134
135	/* Setup QBar */
136	QBar_lm = new double[mSize_];
137
138	/* Calculate "bonds" and build Q_lm(r) where Q_lm = Y_lm(theta(r),phi(r)) */
139	for (std::vector<std::pair<StuntDouble, StuntDouble> >::iterator j = sdPairs_.begin(); j != sdPairs_.end(); ++j) {
140	Vector3d vec = j->first->getPos() - j->second->getPos();
141	currentSnapshot_->wrapVector(vec);
142	/* The spherical harmonics are wrt any arbitray coordiate sysetm,
143	* we choose standard spherical coordinates */
144	r = sqrt(pow(vec.x(),2)+pow(vec.y(),2)+pow(vec.z(),2));
145
146	/* Check to see if neighbor is in bond cuttoff*/
147	if (r<rCut_){
148	theta = atan(vec.y()/vec.x());
149	phi = acos(vec.z()/r);
150	for(int m_index = 0; m_index < mSize_; m_index++){
151	sphericalHarmonic.setM(m_index-lNumber_);
152	QBar_lm(m_index) += sphericalHarmonic.getValueAt(theta,phi);
153	}
154	nBonds++;
155	}
156	}
157
158	/Normalize Qbar by number of Bonds/
159	for ( int m_index = 0;m_index < mSize_; m_index++){
160	QBar_lm(m_index) = QBar_lm(m_index)/nBonds;
161	}
162
163
164	}
165
166	/Normalize by number of frames/
167	for ( int m_index = 0;m_index < mSize_; m_index++){
168	QBar_lm(m_index) = QBar_lm(m_index)/nFrames;
169	}
170
171
172
173	/* Find second order invariant Q_l*/
174
175	for (int m_index = 0 ;m_index <= sizeM_; m++){
176	QSq_l += pow(QBar_lm(m),2);
177	}
178	Q_l_ = sqrt((4NumericConstant::PI/lNumber_+1)QSq_l);
179
180	/* Find Third Order Invariant W_l*/
181
182	/* Make arrays for Wigner3jm */
183	double* THRCOF = new double[mSize_];
184	/* Variables for Wigner routine */
185	double l_ = (double)lNumber_;
186	double m2Min;
187	double m2Max;
188	int error;
189	int m1;
190	int m2;
191	int m3;
192
193	for (int m1 = -lNumber_;m <= lNumber_;m1++){
194	/* Zero work array */
195	for (i=0; i<mSize_;i++){
196	THRCOF[i] = 0.0;
197	}
198	/* Get wigner coefficients */
199	Wigner3jm(&l_,&l_,&l_,&(double)m1,&m2Min,&m2Max&,THRCOF,&mSize_,&error);
200	for (m_index=1; i<(m2Max-M2Min-1.0);m_index++){
201	m2 = floor(m2Min) + m_index - 1;
202	m3 = -m1-m2;
203	W_l_ += THRCOF(m_index)QBar_lm(m1+lNumber_)QBar_lm(m2+lNumber_)*QBar_lm(m3+lNumber_);
204	}
205	}
206
207
208	writeOrderParameter();
209
210	}
211
212
213	void BondOrderParameter::writeOrderParameter() {
214
215	std::ofstream os(getOutputFileName().c_str());
216	os << "#radial distribution function\n";
217	os<< "#selection1: (" << selectionScript1_ << ")\t";
218	os << "selection2: (" << selectionScript2_ << ")\n";
219	os << "#p2\tdirector_x\tdirector_y\tdiretor_z\tangle(degree)\n";
220
221	for (std::size_t i = 0; i < orderParams_.size(); ++i) {
222	os << orderParams_[i].p2 << "\t"
223	<< orderParams_[i].director[0] << "\t"
224	<< orderParams_[i].director[1] << "\t"
225	<< orderParams_[i].director[2] << "\t"
226	<< orderParams_[i].angle << "\n";
227
228	}
229
230	}
231
232	}
233