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:	2925
Committed:	Fri Jul 7 14:18:36 2006 UTC (18 years ago) by chuckv
File size:	8156 byte(s)
Log Message:	Adding more changes to orientational bond order parameter.
#	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
92	sdPairs_.push_back(std::make_pair(sd1, sd2));
93	}
94
95
96	}
97
98	void BondOrderParameter::process() {
99	Molecule* mol;
100	RigidBody* rb;
101	SimInfo::MoleculeIterator mi;
102	Molecule::RigidBodyIterator rbIter;
103	RealType theta;
104	RealType phi;
105	RealType r;
106	RealType dist;
107	RealType* QBar_lm;
108	int nBonds;
109	RealSphericalHarmonic sphericalHarmonic;
110
111
112	DumpReader reader(info_, dumpFilename_);
113	int nFrames = reader.getNFrames();
114
115	/Set the l for the spherical harmonic, it doesn't change/
116	sphericalHarmonic.setL(lNumber_);
117
118	for (int i = 0; i < nFrames; i += step_) {
119	reader.readFrame(i);
120	currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
121	nBonds = 0;
122
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	/* Calculate "bonds" and build Q_lm(r) where Q_lm = Y_lm(theta(r),phi(r)) */
132	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	/* 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	}
150
151	/Normalize Qbar/
152	for (int m = -lNumber_;m <= lNumber_; m++){
153	QBar_lm(m) = QBar_lm(m)/nBonds;
154	}
155
156
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
184	void BondOrderParameter::initalizeHistogram() {
185	std::fill(histogram_.begin(), histogram_.end(), 0);
186	}
187
188	void BondOrderParameter::processHistogram() {
189
190	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
195	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
213	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	}
225
226
227
228
229
230
231	void BondOrderParameter::writeOrderParameter() {
232
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