| 48 |
|
#include "primitives/Molecule.hpp" |
| 49 |
|
#include "utils/NumericConstant.hpp" |
| 50 |
|
#include "utils/PhysicalConstants.hpp" |
| 51 |
< |
#include "math/Polynomial.hpp" |
| 51 |
> |
#include "math/Eigenvalue.hpp" |
| 52 |
|
|
| 53 |
|
namespace OpenMD { |
| 54 |
|
|
| 58 |
|
: SequentialAnalyzer(info, filename), selectionScript_(sele), |
| 59 |
|
evaluator_(info), seleMan_(info), solidZ_(solidZ), |
| 60 |
|
threshDens_(threshDens), nRBins_(nrbins), nZBins_(nzbins) { |
| 61 |
< |
|
| 61 |
> |
|
| 62 |
|
setOutputName(getPrefix(filename) + ".ca2"); |
| 63 |
|
|
| 64 |
|
evaluator_.loadScriptString(sele); |
| 77 |
|
Mat3x3d hmat = info_->getSnapshotManager()->getCurrentSnapshot()->getHmat(); |
| 78 |
|
RealType len = std::min(hmat(0, 0), hmat(1, 1)); |
| 79 |
|
RealType zLen = hmat(2,2); |
| 80 |
+ |
|
| 81 |
|
RealType dr = len / (RealType) nRBins_; |
| 82 |
|
RealType dz = zLen / (RealType) nZBins_; |
| 83 |
|
|
| 84 |
|
std::vector<std::vector<RealType> > histo; |
| 85 |
|
histo.resize(nRBins_); |
| 85 |
– |
for (int i = 0 ; i < nRBins_; ++i) { |
| 86 |
– |
histo[i].resize(nZBins_); |
| 87 |
– |
} |
| 86 |
|
for (unsigned int i = 0; i < histo.size(); ++i){ |
| 87 |
+ |
histo[i].resize(nZBins_); |
| 88 |
|
std::fill(histo[i].begin(), histo[i].end(), 0.0); |
| 89 |
|
} |
| 90 |
|
|
| 114 |
|
for (sd = seleMan_.beginSelected(i); sd != NULL; |
| 115 |
|
sd = seleMan_.nextSelected(i)) { |
| 116 |
|
pos = sd->getPos() - com; |
| 117 |
< |
|
| 117 |
> |
|
| 118 |
> |
// r goes from zero upwards |
| 119 |
|
r = sqrt(pow(pos.x(), 2) + pow(pos.y(), 2)); |
| 120 |
< |
z = pos.z() - solidZ_; |
| 121 |
< |
|
| 120 |
> |
// z is possibly symmetric around 0 |
| 121 |
> |
z = pos.z(); |
| 122 |
> |
|
| 123 |
|
int whichRBin = int(r / dr); |
| 124 |
< |
int whichZBin = int(z/ dz); |
| 124 |
> |
int whichZBin = int( (zLen/2.0 + z) / dz); |
| 125 |
|
|
| 126 |
< |
if ((r <= len) && (z <= zLen)) |
| 126 |
> |
if ((whichRBin < nRBins_) && (whichZBin >= 0) && (whichZBin < nZBins_)) |
| 127 |
|
histo[whichRBin][whichZBin] += sd->getMass(); |
| 128 |
|
|
| 129 |
|
} |
| 134 |
|
RealType rU = rL + dr; |
| 135 |
|
RealType volSlice = NumericConstant::PI * dz * (( rU*rU ) - ( rL*rL )); |
| 136 |
|
|
| 137 |
< |
for (unsigned int j = 0; j < histo[i].size(); ++j){ |
| 137 |
> |
for (unsigned int j = 0; j < histo[i].size(); ++j) { |
| 138 |
|
histo[i][j] *= PhysicalConstants::densityConvert / volSlice; |
| 139 |
|
} |
| 140 |
|
} |
| 141 |
|
|
| 142 |
< |
for (unsigned int i = 0; i < histo.size(); ++i) { |
| 143 |
< |
RealType ther = dr * (i + 0.5); |
| 144 |
< |
for(unsigned int j = 0; j < histo[i].size(); ++j) { |
| 145 |
< |
if (histo[i][j] <= threshDens_) { |
| 146 |
< |
RealType thez = dz * (j + 0.5); |
| 147 |
< |
cerr << ther << "\t" << thez << "\n"; |
| 148 |
< |
break; |
| 142 |
> |
std::vector<Vector<RealType, 2> > points; |
| 143 |
> |
points.clear(); |
| 144 |
> |
|
| 145 |
> |
for (unsigned int j = 0; j < nZBins_; ++j) { |
| 146 |
> |
|
| 147 |
> |
// The z coordinates were measured relative to the selection |
| 148 |
> |
// center of mass. However, we're interested in the elevation |
| 149 |
> |
// above the solid surface. Also, the binning was done around |
| 150 |
> |
// zero with enough bins to cover the zLength of the box: |
| 151 |
> |
|
| 152 |
> |
RealType thez = com.z() - solidZ_ - zLen/2.0 + dz * (j + 0.5); |
| 153 |
> |
bool aboveThresh = false; |
| 154 |
> |
bool foundThresh = false; |
| 155 |
> |
int rloc = 0; |
| 156 |
> |
|
| 157 |
> |
for (unsigned int i = 0; i < nRBins_; ++i) { |
| 158 |
> |
RealType ther = dr * (i + 0.5); |
| 159 |
> |
if (histo[i][j] >= threshDens_) aboveThresh = true; |
| 160 |
> |
|
| 161 |
> |
if (aboveThresh && (histo[i][j] <= threshDens_)) { |
| 162 |
> |
rloc = i; |
| 163 |
> |
foundThresh = true; |
| 164 |
> |
aboveThresh = false; |
| 165 |
|
} |
| 166 |
+ |
|
| 167 |
|
} |
| 168 |
+ |
if (foundThresh) { |
| 169 |
+ |
Vector<RealType,2> point; |
| 170 |
+ |
point[0] = dr*(rloc+0.5); |
| 171 |
+ |
point[1] = thez; |
| 172 |
+ |
points.push_back( point ); |
| 173 |
+ |
} |
| 174 |
|
} |
| 175 |
|
|
| 176 |
< |
// values_.push_back( acos(maxct)*(180.0/M_PI) ); |
| 176 |
> |
int numPoints = points.size(); |
| 177 |
> |
|
| 178 |
> |
// Compute the average of the data points. |
| 179 |
> |
Vector<RealType, 2> average = points[0]; |
| 180 |
> |
int i0; |
| 181 |
> |
for (i0 = 1; i0 < numPoints; ++i0) { |
| 182 |
> |
average += points[i0]; |
| 183 |
> |
} |
| 184 |
> |
RealType invNumPoints = ((RealType)1)/(RealType)numPoints; |
| 185 |
> |
average *= invNumPoints; |
| 186 |
|
|
| 187 |
+ |
DynamicRectMatrix<RealType> mat(4, 4); |
| 188 |
+ |
int row, col; |
| 189 |
+ |
for (row = 0; row < 4; ++row) { |
| 190 |
+ |
for (col = 0; col < 4; ++col){ |
| 191 |
+ |
mat(row,col) = 0.0; |
| 192 |
+ |
} |
| 193 |
+ |
} |
| 194 |
+ |
for (int i = 0; i < numPoints; ++i) { |
| 195 |
+ |
RealType x = points[i][0]; |
| 196 |
+ |
RealType y = points[i][1]; |
| 197 |
+ |
RealType x2 = x*x; |
| 198 |
+ |
RealType y2 = y*y; |
| 199 |
+ |
RealType xy = x*y; |
| 200 |
+ |
RealType r2 = x2+y2; |
| 201 |
+ |
RealType xr2 = x*r2; |
| 202 |
+ |
RealType yr2 = y*r2; |
| 203 |
+ |
RealType r4 = r2*r2; |
| 204 |
+ |
|
| 205 |
+ |
mat(0,1) += x; |
| 206 |
+ |
mat(0,2) += y; |
| 207 |
+ |
mat(0,3) += r2; |
| 208 |
+ |
mat(1,1) += x2; |
| 209 |
+ |
mat(1,2) += xy; |
| 210 |
+ |
mat(1,3) += xr2; |
| 211 |
+ |
mat(2,2) += y2; |
| 212 |
+ |
mat(2,3) += yr2; |
| 213 |
+ |
mat(3,3) += r4; |
| 214 |
+ |
} |
| 215 |
+ |
mat(0,0) = (RealType)numPoints; |
| 216 |
+ |
|
| 217 |
+ |
for (row = 0; row < 4; ++row) { |
| 218 |
+ |
for (col = 0; col < row; ++col) { |
| 219 |
+ |
mat(row,col) = mat(col,row); |
| 220 |
+ |
} |
| 221 |
+ |
} |
| 222 |
+ |
|
| 223 |
+ |
for (row = 0; row < 4; ++row) { |
| 224 |
+ |
for (col = 0; col < 4; ++col) { |
| 225 |
+ |
mat(row,col) *= invNumPoints; |
| 226 |
+ |
} |
| 227 |
+ |
} |
| 228 |
+ |
|
| 229 |
+ |
JAMA::Eigenvalue<RealType> eigensystem(mat); |
| 230 |
+ |
DynamicRectMatrix<RealType> evects(4, 4); |
| 231 |
+ |
DynamicVector<RealType> evals(4); |
| 232 |
+ |
|
| 233 |
+ |
eigensystem.getRealEigenvalues(evals); |
| 234 |
+ |
eigensystem.getV(evects); |
| 235 |
+ |
|
| 236 |
+ |
DynamicVector<RealType> evector = evects.getColumn(0); |
| 237 |
+ |
RealType inv = ((RealType)1)/evector[3]; // beware zero divide |
| 238 |
+ |
RealType coeff[3]; |
| 239 |
+ |
for (row = 0; row < 3; ++row) { |
| 240 |
+ |
coeff[row] = inv*evector[row]; |
| 241 |
+ |
} |
| 242 |
+ |
|
| 243 |
+ |
Vector<RealType, 2> center; |
| 244 |
+ |
|
| 245 |
+ |
center[0] = -((RealType)0.5)*coeff[1]; |
| 246 |
+ |
center[1] = -((RealType)0.5)*coeff[2]; |
| 247 |
+ |
RealType radius = sqrt(fabs(center[0]*center[0] + center[1]*center[1] |
| 248 |
+ |
- coeff[0])); |
| 249 |
+ |
RealType ev0 = fabs(evals[0]); |
| 250 |
+ |
|
| 251 |
+ |
int i1; |
| 252 |
+ |
for (i1 = 0; i1 < 100; ++i1) { |
| 253 |
+ |
// Update the iterates. |
| 254 |
+ |
Vector<RealType, 2> current = center; |
| 255 |
+ |
|
| 256 |
+ |
// Compute average L, dL/da, dL/db. |
| 257 |
+ |
RealType lenAverage = (RealType)0; |
| 258 |
+ |
Vector<RealType, 2> derLenAverage = Vector<RealType, 2>(0.0); |
| 259 |
+ |
for (i0 = 0; i0 < numPoints; ++i0) { |
| 260 |
+ |
Vector<RealType, 2> diff = points[i0] - center; |
| 261 |
+ |
RealType length = diff.length(); |
| 262 |
+ |
if (length > 1e-6) { |
| 263 |
+ |
lenAverage += length; |
| 264 |
+ |
RealType invLength = ((RealType)1)/length; |
| 265 |
+ |
derLenAverage -= invLength*diff; |
| 266 |
+ |
} |
| 267 |
+ |
} |
| 268 |
+ |
lenAverage *= invNumPoints; |
| 269 |
+ |
derLenAverage *= invNumPoints; |
| 270 |
+ |
|
| 271 |
+ |
center = average + lenAverage*derLenAverage; |
| 272 |
+ |
radius = lenAverage; |
| 273 |
+ |
|
| 274 |
+ |
Vector<RealType, 2> diff = center - current; |
| 275 |
+ |
if (fabs(diff[0]) <= 1e-6 && fabs(diff[1]) <= 1e-6) { |
| 276 |
+ |
break; |
| 277 |
+ |
} |
| 278 |
+ |
} |
| 279 |
+ |
|
| 280 |
+ |
RealType zCen = center[1]; |
| 281 |
+ |
RealType rDrop = radius; |
| 282 |
+ |
RealType ca; |
| 283 |
+ |
|
| 284 |
+ |
if (fabs(zCen) > rDrop) { |
| 285 |
+ |
ca = 180.0; |
| 286 |
+ |
} else { |
| 287 |
+ |
|
| 288 |
+ |
if (zCen >= 0.0) { |
| 289 |
+ |
ca = 90.0 + asin(zCen/rDrop)*(180.0/M_PI); |
| 290 |
+ |
} else { |
| 291 |
+ |
ca = 90 - asin(zCen/rDrop)*(180.0/M_PI); |
| 292 |
+ |
} |
| 293 |
+ |
} |
| 294 |
+ |
|
| 295 |
+ |
values_.push_back( ca ); |
| 296 |
+ |
|
| 297 |
|
} |
| 298 |
|
} |
| 299 |
|
|
