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)) {

     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;
  int nBonds;
  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();
        }
        
    }      

      /* 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 = -lNumber_; m <= lNumber_; m++){
                sphericalHarmonic.setM(m);
                QBar_lm(m) += sphericalHarmonic.getValueAt(theta,phi);
            }
            nBonds++;
                }  
      }
      
     /*Normalize Qbar*/
     for (int m = -lNumber_;m <= lNumber_; m++){
      QBar_lm(m) = QBar_lm(m)/nBonds;   
     }
      
    
  }
  
  /*Normalize by number of frames*/
    for (int m = -lNumber_;m <= lNumber_; m++){
      QBar_lm(m) = QBar_lm(m)/nFrames;  
     }
     
     
     /* Find second order invariant Q_l*/
     
       for (int m = -lNumber_;m <= lNumber_; m++){
         QSq_l += pow(QBar_lm(m),2);    
        }
     Q_l = sqrt((4*NumericConstant::PI/lNumber_+1)*QSq_l);
     
     /* Find Third Order Invariant W_l*/
     for (int m = -lNumber_;m<= lNumber_;m++){
        
        
     }
     
     
  writeOrderParameter();
  
}

void BondOrderParameter::initalizeHistogram() {
    std::fill(histogram_.begin(), histogram_.end(), 0);
  }

 void BondOrderParameter::processHistogram() {

    int nPairs = getNPairs();
    RealType volume = info_->getSnapshotManager()->getCurrentSnapshot()->getVolume();
    RealType pairDensity = nPairs /volume * 2.0;
    RealType pairConstant = ( 4.0 * NumericConstant::PI * pairDensity ) / 3.0;

    for(int i = 0 ; i < histogram_.size(); ++i){

      RealType rLower = i * deltaR_;
      RealType rUpper = rLower + deltaR_;
      RealType volSlice = ( rUpper * rUpper * rUpper ) - ( rLower * rLower * rLower );
      RealType nIdeal = volSlice * pairConstant;

      avgGofr_[i] += histogram_[i] / nIdeal;    
    }

  }

  void BondOrderParameter::collectHistogram(StuntDouble* sd1, StuntDouble* sd2) {

    if (sd1 == sd2) {
      return;
    }
    
    Vector3d pos1 = sd1->getPos();
    Vector3d pos2 = sd2->getPos();
    Vector3d r12 = pos2 - pos1;
    currentSnapshot_->wrapVector(r12);

    RealType distance = r12.length();

    if (distance < len_) {
      int whichBin = distance / deltaR_;
      histogram_[whichBin] += 2;
    }
  }


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:	1001
Committed:	Fri Jul 7 14:18:36 2006 UTC (19 years, 3 months ago) by chuckv
File size:	8156 byte(s)
Log Message:	Adding more changes to orientational bond order parameter.
#	User	Rev	Content
1	chuckv	994	/*
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	chuckv	1001
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	chuckv	995	#include "applications/staticProps/BondOrderParameter.hpp"
50	chuckv	994	#include "utils/simError.h"
51			#include "io/DumpReader.hpp"
52			#include "primitives/Molecule.hpp"
53			#include "utils/NumericConstant.hpp"
54	chuckv	1001	#include "math/RealSphericalHarmonic.hpp"
55	chuckv	994	namespace oopse {
56
57
58	chuckv	995	BondOrderParameter::BondOrderParameter(SimInfo* info, const std::string& filename, const std::string& sele1,
59			const std::string& sele2, double rCut, int lNumber)
60	chuckv	994	: StaticAnalyser(info, filename),
61	chuckv	995	selectionScript1_(sele1), evaluator1_(info),
62			seleMan1_(info){
63	chuckv	994
64	chuckv	1001	setOutputName(getPrefix(filename) + ".obo");
65	chuckv	994
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	chuckv	1001	/* Set up cutoff radius and type of order parameter we are calcuating*/
80			lNumber_ = lNumber;
81			rCut_ = rCut;
82			mSize_ = 2*lNumber_+1;
83	chuckv	994
84			int i;
85			int j;
86			StuntDouble* sd1;
87			StuntDouble* sd2;
88	chuckv	995	for (sd1 = seleMan1_.beginSelected(i), sd2 = seleMan1_.beginSelected(j);
89	chuckv	994	sd1 != NULL && sd2 != NULL;
90			sd1 = seleMan1_.nextSelected(i), sd2 = seleMan2_.nextSelected(j)) {
91
92			sdPairs_.push_back(std::make_pair(sd1, sd2));
93			}
94
95
96			}
97
98	chuckv	995	void BondOrderParameter::process() {
99	chuckv	994	Molecule* mol;
100			RigidBody* rb;
101			SimInfo::MoleculeIterator mi;
102			Molecule::RigidBodyIterator rbIter;
103	chuckv	1001	RealType theta;
104			RealType phi;
105			RealType r;
106			RealType dist;
107			RealType* QBar_lm;
108			int nBonds;
109			RealSphericalHarmonic sphericalHarmonic;
110	chuckv	994
111	chuckv	1001
112	chuckv	994	DumpReader reader(info_, dumpFilename_);
113			int nFrames = reader.getNFrames();
114
115	chuckv	1001	/Set the l for the spherical harmonic, it doesn't change/
116			sphericalHarmonic.setL(lNumber_);
117
118	chuckv	994	for (int i = 0; i < nFrames; i += step_) {
119			reader.readFrame(i);
120			currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
121	chuckv	1001	nBonds = 0;
122	chuckv	994
123			for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
124			//change the positions of atoms which belong to the rigidbodies
125			for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) {
126			rb->updateAtoms();
127			}
128
129			}
130
131	chuckv	1001	/* Calculate "bonds" and build Q_lm(r) where Q_lm = Y_lm(theta(r),phi(r)) */
132	chuckv	994	for (std::vector<std::pair<StuntDouble, StuntDouble> >::iterator j = sdPairs_.begin(); j != sdPairs_.end(); ++j) {
133			Vector3d vec = j->first->getPos() - j->second->getPos();
134			currentSnapshot_->wrapVector(vec);
135	chuckv	1001	/* The spherical harmonics are wrt any arbitray coordiate sysetm,
136			* we choose standard spherical coordinates */
137			r = sqrt(pow(vec.x(),2)+pow(vec.y(),2)+pow(vec.z(),2));
138
139			/* Check to see if neighbor is in bond cuttoff*/
140			if (r<rCut_){
141			theta = atan(vec.y()/vec.x());
142			phi = acos(vec.z()/r);
143			for(int m = -lNumber_; m <= lNumber_; m++){
144			sphericalHarmonic.setM(m);
145			QBar_lm(m) += sphericalHarmonic.getValueAt(theta,phi);
146			}
147			nBonds++;
148			}
149	chuckv	994	}
150
151	chuckv	1001	/Normalize Qbar/
152			for (int m = -lNumber_;m <= lNumber_; m++){
153			QBar_lm(m) = QBar_lm(m)/nBonds;
154			}
155	chuckv	994
156	chuckv	1001
157			}
158
159			/Normalize by number of frames/
160			for (int m = -lNumber_;m <= lNumber_; m++){
161			QBar_lm(m) = QBar_lm(m)/nFrames;
162			}
163
164
165
166			/* Find second order invariant Q_l*/
167
168			for (int m = -lNumber_;m <= lNumber_; m++){
169			QSq_l += pow(QBar_lm(m),2);
170			}
171			Q_l = sqrt((4NumericConstant::PI/lNumber_+1)QSq_l);
172
173			/* Find Third Order Invariant W_l*/
174			for (int m = -lNumber_;m<= lNumber_;m++){
175
176
177			}
178
179
180			writeOrderParameter();
181
182			}
183	chuckv	994
184	chuckv	1001	void BondOrderParameter::initalizeHistogram() {
185			std::fill(histogram_.begin(), histogram_.end(), 0);
186			}
187	chuckv	994
188	chuckv	1001	void BondOrderParameter::processHistogram() {
189	chuckv	994
190	chuckv	1001	int nPairs = getNPairs();
191			RealType volume = info_->getSnapshotManager()->getCurrentSnapshot()->getVolume();
192			RealType pairDensity = nPairs /volume * 2.0;
193			RealType pairConstant = ( 4.0 * NumericConstant::PI * pairDensity ) / 3.0;
194	chuckv	994
195	chuckv	1001	for(int i = 0 ; i < histogram_.size(); ++i){
196
197			RealType rLower = i * deltaR_;
198			RealType rUpper = rLower + deltaR_;
199			RealType volSlice = ( rUpper * rUpper * rUpper ) - ( rLower * rLower * rLower );
200			RealType nIdeal = volSlice * pairConstant;
201
202			avgGofr_[i] += histogram_[i] / nIdeal;
203			}
204
205			}
206
207			void BondOrderParameter::collectHistogram(StuntDouble* sd1, StuntDouble* sd2) {
208
209			if (sd1 == sd2) {
210			return;
211			}
212	chuckv	994
213	chuckv	1001	Vector3d pos1 = sd1->getPos();
214			Vector3d pos2 = sd2->getPos();
215			Vector3d r12 = pos2 - pos1;
216			currentSnapshot_->wrapVector(r12);
217
218			RealType distance = r12.length();
219
220			if (distance < len_) {
221			int whichBin = distance / deltaR_;
222			histogram_[whichBin] += 2;
223			}
224	chuckv	994	}
225
226
227	chuckv	1001
228
229
230
231	chuckv	995	void BondOrderParameter::writeOrderParameter() {
232	chuckv	994
233			std::ofstream os(getOutputFileName().c_str());
234			os << "#radial distribution function\n";
235			os<< "#selection1: (" << selectionScript1_ << ")\t";
236			os << "selection2: (" << selectionScript2_ << ")\n";
237			os << "#p2\tdirector_x\tdirector_y\tdiretor_z\tangle(degree)\n";
238
239			for (std::size_t i = 0; i < orderParams_.size(); ++i) {
240			os << orderParams_[i].p2 << "\t"
241			<< orderParams_[i].director[0] << "\t"
242			<< orderParams_[i].director[1] << "\t"
243			<< orderParams_[i].director[2] << "\t"
244			<< orderParams_[i].angle << "\n";
245
246			}
247
248			}
249
250			}
251