| 96 |
|
frames_ = 0; |
| 97 |
|
nProcessed_ = nFrames/step_; |
| 98 |
|
|
| 99 |
< |
//positions_ and moBool_ are 2D arrays, need the second dimension filled as well |
| 99 |
> |
// positions_ and moBool_ are 2D arrays, need the second dimension |
| 100 |
> |
// filled as well |
| 101 |
|
for(int i = 0; i < selectionCount_; i++){ |
| 102 |
|
moBool_[i].resize(nFrames); |
| 103 |
|
positions_[i].resize(nFrames); |
| 106 |
|
int iterator; |
| 107 |
|
int index = 0; |
| 108 |
|
/* Loop over all frames storing the positions in a vec< vec<pos> > |
| 109 |
< |
* At the end, positions.length() should equal seleMan1_.size() or w/e |
| 110 |
< |
* And positions[index].length() should equal nFrames (or nFrames/istep) |
| 109 |
> |
* At the end, positions.length() should equal seleMan1_.size() or |
| 110 |
> |
* w/e And positions[index].length() should equal nFrames (or |
| 111 |
> |
* nFrames/istep) |
| 112 |
|
*/ |
| 113 |
|
for(int istep = 0; istep < nFrames; istep += step_){ |
| 114 |
|
frames_++; |
| 115 |
|
reader.readFrame(istep); |
| 116 |
|
currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot(); |
| 117 |
|
|
| 118 |
< |
for(mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)){ |
| 118 |
> |
for(mol = info_->beginMolecule(mi); mol != NULL; |
| 119 |
> |
mol = info_->nextMolecule(mi)){ |
| 120 |
|
//change the positions of atoms which belong to the rigidbodies |
| 121 |
< |
for(rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)){ |
| 121 |
> |
for(rb = mol->beginRigidBody(rbIter); rb != NULL; |
| 122 |
> |
rb = mol->nextRigidBody(rbIter)){ |
| 123 |
|
rb->updateAtoms(); |
| 124 |
|
} |
| 125 |
|
} |
| 126 |
|
|
| 127 |
< |
index = 0; //count over atoms since iterators aren't the most friendly for such plebian things |
| 128 |
< |
for(sd = seleMan1_.beginSelected(iterator); sd != NULL; sd = seleMan1_.nextSelected(iterator)){ |
| 127 |
> |
index = 0; // count over atoms since iterators aren't the most |
| 128 |
> |
// friendly for such plebian things |
| 129 |
> |
for(sd = seleMan1_.beginSelected(iterator); sd != NULL; |
| 130 |
> |
sd = seleMan1_.nextSelected(iterator)){ |
| 131 |
|
Vector3d pos = sd->getPos(); |
| 132 |
|
positions_[index][istep] = pos; |
| 133 |
|
index++; |
| 134 |
|
} |
| 135 |
|
} |
| 136 |
|
|
| 131 |
– |
|
| 137 |
|
cout << "Position Array size: " << positions_.size() << "\n"; |
| 138 |
|
cout << "Frames analyzed: " << positions_[0].size() << "\n"; |
| 139 |
|
|
| 140 |
< |
|
| 136 |
< |
for(int i = 0; i < positions_.size(); i++){ |
| 140 |
> |
for(std::size_t i = 0; i < positions_.size(); i++){ |
| 141 |
|
int frameIndex = positions_[i].size(); |
| 142 |
|
for(int j = 1; j < frameIndex; j++){ |
| 143 |
|
Vector3d posF1 = positions_[i][j-1]; |
| 156 |
|
} |
| 157 |
|
|
| 158 |
|
int mobileAtomCount = 0; |
| 159 |
< |
for(int i = 0; i < moBool_.size(); i++){ |
| 159 |
> |
for(std::size_t i = 0; i < moBool_.size(); i++){ |
| 160 |
|
int frameIndex = moBool_[i].size(); |
| 161 |
|
bool mobileAtom = false; |
| 162 |
|
for(int j = 0; j < frameIndex; j++){ |
| 163 |
|
mobileAtom = mobileAtom || moBool_[i][j]; |
| 164 |
|
} |
| 165 |
< |
moBool_[i][0] = mobileAtom; //is true if any value later in the array is true, false otherwise |
| 165 |
> |
moBool_[i][0] = mobileAtom; // is true if any value later in the |
| 166 |
> |
// array is true, false otherwise |
| 167 |
|
if(mobileAtom){ |
| 168 |
|
mobileAtomCount++; |
| 169 |
|
} |
| 171 |
|
|
| 172 |
|
cout << "Mobile atom count: " << mobileAtomCount << "\n"; |
| 173 |
|
|
| 174 |
< |
//Here I shrink the size of the arrays, why look through 3888, when you only need ~800. |
| 175 |
< |
//Additionally, all of these are mobile at some point in time, the others aren't, dead weight and |
| 176 |
< |
//memory |
| 174 |
> |
// Here I shrink the size of the arrays, why look through 3888, |
| 175 |
> |
// when you only need ~800. Additionally, all of these are mobile |
| 176 |
> |
// at some point in time, the others aren't, dead weight and |
| 177 |
> |
// memory |
| 178 |
|
positions2_.resize(mobileAtomCount); |
| 179 |
|
moBool2_.resize(mobileAtomCount); |
| 180 |
|
int pos2index = 0; |
| 181 |
< |
for(int i = 0; i < positions_.size(); i++){ |
| 181 |
> |
for(std::size_t i = 0; i < positions_.size(); i++){ |
| 182 |
|
int frameCount = positions_[i].size(); |
| 183 |
|
if(moBool_[i][0]){ |
| 184 |
|
for(int j = 0; j < frameCount; j++){ |
| 206 |
|
diffStream.open(outputFilename_.c_str()); |
| 207 |
|
diffStream << "#X&Y diffusion amounts\n"; |
| 208 |
|
diffStream << "#singleMoveDistance_: " << singleMoveDistance_ << "\n"; |
| 203 |
– |
|
| 209 |
|
diffStream << "#Number of mobile atoms: " << positions2_.size() << "\n"; |
| 205 |
– |
|
| 210 |
|
diffStream << "#time, <x(t)-x(0)>, <y(t)-y(0)>, <r(t)-r(0)>\n"; |
| 211 |
< |
for(int i = 0; i < xHist_.size(); i++){ |
| 212 |
< |
diffStream << i << ", " << xHist_[i] << ", " << yHist_[i] << ", " << rHist_[i] << "\n"; |
| 211 |
> |
|
| 212 |
> |
for(std::size_t i = 0; i < xHist_.size(); i++){ |
| 213 |
> |
diffStream << i << ", " << xHist_[i] << ", " << yHist_[i] << ", " |
| 214 |
> |
<< rHist_[i] << "\n"; |
| 215 |
|
} |
| 216 |
|
diffStream.close(); |
| 217 |
|
|
| 237 |
|
//loop over particles |
| 238 |
|
// loop over frames starting at j |
| 239 |
|
// loop over frames starting at k = j (time shift of 0) |
| 240 |
< |
for(int i = 0; i < positions2_.size(); i++){ |
| 241 |
< |
int frames = positions2_[i].size() - 1; //for counting properly, otherwise moBool2_[i][j+1] will |
| 242 |
< |
// go over |
| 240 |
> |
for(std::size_t i = 0; i < positions2_.size(); i++){ |
| 241 |
> |
int frames = positions2_[i].size() - 1; // for counting |
| 242 |
> |
// properly, otherwise |
| 243 |
> |
// moBool2_[i][j+1] will |
| 244 |
> |
// go over |
| 245 |
|
for(int j = 0; j < frames; j++){ |
| 246 |
< |
//if the particle is mobile between j and j + 1, then count it for all timeShifts |
| 246 |
> |
// if the particle is mobile between j and j + 1, then count |
| 247 |
> |
// it for all timeShifts |
| 248 |
|
if(moBool2_[i][j+1]){ |
| 249 |
< |
for(int k = j; k < positions2_[0].size(); k++){ |
| 249 |
> |
for(std::size_t k = j; k < positions2_[0].size(); k++){ |
| 250 |
|
//<x(t)-x(0)> <y(t)-y(0)> <r(t)-r(0)> |
| 251 |
< |
//The positions stored are not wrapped, thus I don't need to worry about pbc |
| 251 |
> |
//The positions stored are not wrapped, thus I don't need |
| 252 |
> |
//to worry about pbc |
| 253 |
|
//Mean square displacement |
| 254 |
|
//So I do want the squared distances |
| 255 |
|
|
| 275 |
|
} |
| 276 |
|
cout << "X, Y, R calculated\n"; |
| 277 |
|
|
| 278 |
< |
for(int i = 0; i < xHist_.size(); i++){ |
| 278 |
> |
for(std::size_t i = 0; i < xHist_.size(); i++){ |
| 279 |
|
xHist_[i] = xHist_[i]/(count_[i]); |
| 280 |
|
yHist_[i] = yHist_[i]/(count_[i]); |
| 281 |
|
rHist_[i] = rHist_[i]/(count_[i]); |
