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((4*NumericConstant::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 |
|