43 |
|
#include "utils/simError.h" |
44 |
|
#include "io/DumpReader.hpp" |
45 |
|
#include "primitives/Molecule.hpp" |
46 |
< |
|
46 |
> |
#include "utils/NumericConstant.hpp" |
47 |
|
namespace oopse { |
48 |
|
|
49 |
|
|
112 |
|
for (int i = 0; i < nFrames; i += step_) { |
113 |
|
reader.readFrame(i); |
114 |
|
currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot(); |
115 |
< |
Mat3x3d orderTensor(0.0); |
115 |
> |
|
116 |
|
|
117 |
|
for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) { |
118 |
|
//change the positions of atoms which belong to the rigidbodies |
122 |
|
|
123 |
|
} |
124 |
|
|
125 |
+ |
Mat3x3d orderTensor(0.0); |
126 |
|
for (std::vector<std::pair<StuntDouble*, StuntDouble*> >::iterator j = sdPairs_.begin(); j != sdPairs_.end(); ++j) { |
127 |
|
Vector3d vec = j->first->getPos() - j->second->getPos(); |
128 |
+ |
currentSnapshot_->wrapVector(vec); |
129 |
|
vec.normalize(); |
130 |
|
orderTensor +=outProduct(vec, vec); |
131 |
|
} |
132 |
|
|
133 |
|
orderTensor /= sdPairs_.size(); |
134 |
< |
orderTensor -= 1.0/3.0 * Mat3x3d::identity(); |
134 |
> |
orderTensor -= (RealType)(1.0/3.0) * Mat3x3d::identity(); |
135 |
|
|
136 |
|
Vector3d eigenvalues; |
137 |
|
Mat3x3d eigenvectors; |
138 |
|
Mat3x3d::diagonalize(orderTensor, eigenvalues, eigenvectors); |
139 |
|
|
140 |
|
int which; |
141 |
< |
double maxEval = 0.0; |
141 |
> |
RealType maxEval = 0.0; |
142 |
|
for(int k = 0; k< 3; k++){ |
143 |
|
if(fabs(eigenvalues[k]) > maxEval){ |
144 |
|
which = k; |
145 |
|
maxEval = fabs(eigenvalues[k]); |
146 |
|
} |
147 |
|
} |
148 |
< |
double p2 = 1.5 * maxEval; |
148 |
> |
RealType p2 = 1.5 * maxEval; |
149 |
|
|
150 |
|
//the eigen vector is already normalized in SquareMatrix3::diagonalize |
151 |
|
Vector3d director = eigenvectors.getColumn(which); |
153 |
|
director.negate(); |
154 |
|
} |
155 |
|
|
156 |
< |
double angle = 0.0; |
156 |
> |
RealType angle = 0.0; |
157 |
|
for (std::vector<std::pair<StuntDouble*, StuntDouble*> >::iterator j = sdPairs_.begin(); j != sdPairs_.end(); ++j) { |
158 |
|
Vector3d vec = j->first->getPos() - j->second->getPos(); |
159 |
+ |
currentSnapshot_->wrapVector(vec); |
160 |
|
vec.normalize(); |
161 |
|
|
162 |
|
angle += acos(dot(vec, director)) ; |
163 |
|
} |
164 |
< |
angle /= sdPairs_.size(); |
164 |
> |
angle = angle / (sdPairs_.size() * NumericConstant::PI) * 180.0; |
165 |
|
|
166 |
|
OrderParam param; |
167 |
|
param.p2 = p2; |
172 |
|
|
173 |
|
} |
174 |
|
|
175 |
< |
writeOrderParam(); |
175 |
> |
writeP2(); |
176 |
|
|
177 |
|
} |
178 |
|
|
179 |
< |
void P2OrderParameter::writeOrderParam() { |
179 |
> |
void P2OrderParameter::writeP2() { |
180 |
|
|
181 |
|
std::ofstream os(getOutputFileName().c_str()); |
182 |
|
os << "#radial distribution function\n"; |