| 56 |
|
const std::string& filename, |
| 57 |
|
const std::string& sele1, |
| 58 |
|
const std::string& sele2, |
| 59 |
< |
RealType rCut, RealType voxelSize, |
| 59 |
> |
RealType rCut, |
| 60 |
> |
RealType voxelSize, |
| 61 |
|
RealType gaussWidth) |
| 62 |
|
: StaticAnalyser(info, filename), selectionScript1_(sele1), |
| 63 |
|
evaluator1_(info), seleMan1_(info), selectionScript2_(sele2), |
| 115 |
|
RealType cospsi; |
| 116 |
|
RealType Qk; |
| 117 |
|
std::vector<std::pair<RealType,StuntDouble*> > myNeighbors; |
| 118 |
< |
std::vector<std::pair<Vector3d, RealType> > qvals; |
| 119 |
< |
std::vector<std::pair<Vector3d, RealType> >::iterator qiter; |
| 118 |
> |
//std::vector<std::pair<Vector3d, RealType> > qvals; |
| 119 |
> |
//std::vector<std::pair<Vector3d, RealType> >::iterator qiter; |
| 120 |
|
int isd1; |
| 121 |
|
int isd2; |
| 122 |
|
|
| 123 |
+ |
|
| 124 |
+ |
int kMax = int(5.0 * gaussWidth_ / voxelSize_); |
| 125 |
+ |
int kSqLim = kMax*kMax; |
| 126 |
+ |
cerr << "gw = " << gaussWidth_ << " vS = " << voxelSize_ << " kMax = " << kMax << " kSqLim = " << kSqLim << "\n"; |
| 127 |
+ |
|
| 128 |
|
DumpReader reader(info_, dumpFilename_); |
| 129 |
|
int nFrames = reader.getNFrames(); |
| 130 |
|
|
| 152 |
|
} |
| 153 |
|
} |
| 154 |
|
|
| 155 |
< |
qvals.clear(); |
| 155 |
> |
//qvals.clear(); |
| 156 |
|
|
| 157 |
|
// outer loop is over the selected StuntDoubles: |
| 158 |
|
for (sd = seleMan1_.beginSelected(isd1); sd != NULL; |
| 223 |
|
if (nang > 0) { |
| 224 |
|
if (usePeriodicBoundaryConditions_) |
| 225 |
|
currentSnapshot_->wrapVector(rk); |
| 226 |
< |
qvals.push_back(std::make_pair(rk, Qk)); |
| 227 |
< |
} |
| 228 |
< |
} |
| 226 |
> |
//qvals.push_back(std::make_pair(rk, Qk)); |
| 227 |
> |
|
| 228 |
> |
Vector3d pos = rk + halfBox; |
| 229 |
|
|
| 230 |
< |
for (int i = 0; i < nBins_(0); ++i) { |
| 231 |
< |
for(int j = 0; j < nBins_(1); ++j) { |
| 232 |
< |
for(int k = 0; k < nBins_(2); ++k) { |
| 233 |
< |
Vector3d pos = Vector3d(i, j, k) * voxelSize_ - halfBox; |
| 234 |
< |
for(qiter = qvals.begin(); qiter != qvals.end(); ++qiter) { |
| 235 |
< |
Vector3d d = pos - (*qiter).first; |
| 236 |
< |
currentSnapshot_->wrapVector(d); |
| 237 |
< |
RealType denom = pow(2.0 * sqrt(M_PI) * gaussWidth_, 3); |
| 238 |
< |
RealType exponent = -dot(d,d) / pow(2.0*gaussWidth_, 2); |
| 239 |
< |
RealType weight = exp(exponent) / denom; |
| 240 |
< |
count_[i][j][k] += weight; |
| 241 |
< |
hist_[i][j][k] += weight * (*qiter).second; |
| 230 |
> |
|
| 231 |
> |
Vector3i whichVoxel(int(pos[0] / voxelSize_), |
| 232 |
> |
int(pos[1] / voxelSize_), |
| 233 |
> |
int(pos[2] / voxelSize_)); |
| 234 |
> |
|
| 235 |
> |
for (int l = -kMax; l <= kMax; l++) { |
| 236 |
> |
for (int m = -kMax; m <= kMax; m++) { |
| 237 |
> |
for (int n = -kMax; n <= kMax; n++) { |
| 238 |
> |
int kk = l*l + m*m + n*n; |
| 239 |
> |
if(kk <= kSqLim) { |
| 240 |
> |
|
| 241 |
> |
int ll = (whichVoxel[0] + l) % nBins_(0); |
| 242 |
> |
ll = ll < 0 ? nBins_(0) + ll : ll; |
| 243 |
> |
int mm = (whichVoxel[1] + m) % nBins_(1); |
| 244 |
> |
mm = mm < 0 ? nBins_(1) + mm : mm; |
| 245 |
> |
int nn = (whichVoxel[2] + n) % nBins_(2); |
| 246 |
> |
nn = nn < 0 ? nBins_(2) + nn : nn; |
| 247 |
> |
|
| 248 |
> |
Vector3d bPos = Vector3d(ll,mm,nn) * voxelSize_ - halfBox; |
| 249 |
> |
Vector3d d = bPos - rk; |
| 250 |
> |
currentSnapshot_->wrapVector(d); |
| 251 |
> |
RealType denom = pow(2.0 * sqrt(M_PI) * gaussWidth_, 3); |
| 252 |
> |
RealType exponent = -dot(d,d) / pow(2.0*gaussWidth_, 2); |
| 253 |
> |
RealType weight = exp(exponent) / denom; |
| 254 |
> |
count_[ll][mm][nn] += weight; |
| 255 |
> |
hist_[ll][mm][nn] += weight * Qk; |
| 256 |
> |
} |
| 257 |
> |
} |
| 258 |
|
} |
| 259 |
|
} |
| 260 |
|
} |
| 261 |
|
} |
| 262 |
+ |
|
| 263 |
+ |
// for (int i = 0; i < nBins_(0); ++i) { |
| 264 |
+ |
// for(int j = 0; j < nBins_(1); ++j) { |
| 265 |
+ |
// for(int k = 0; k < nBins_(2); ++k) { |
| 266 |
+ |
// Vector3d pos = Vector3d(i, j, k) * voxelSize_ - halfBox; |
| 267 |
+ |
// for(qiter = qvals.begin(); qiter != qvals.end(); ++qiter) { |
| 268 |
+ |
// Vector3d d = pos - (*qiter).first; |
| 269 |
+ |
// currentSnapshot_->wrapVector(d); |
| 270 |
+ |
// RealType denom = pow(2.0 * sqrt(M_PI) * gaussWidth_, 3); |
| 271 |
+ |
// RealType exponent = -dot(d,d) / pow(2.0*gaussWidth_, 2); |
| 272 |
+ |
// RealType weight = exp(exponent) / denom; |
| 273 |
+ |
// count_[i][j][k] += weight; |
| 274 |
+ |
// hist_[i][j][k] += weight * (*qiter).second; |
| 275 |
+ |
// } |
| 276 |
+ |
// } |
| 277 |
+ |
// } |
| 278 |
+ |
// } |
| 279 |
|
} |
| 280 |
|
writeQxyz(); |
| 281 |
|
} |
