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:	1038
Committed:	Mon Sep 18 21:31:23 2006 UTC (18 years, 7 months ago) by chuckv
Original Path:	*trunk/src/applications/staticProps/BondOrderParameter.cpp*
File size:	7493 byte(s)
Log Message:	Latest Version
#	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			* 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