| 155 |
|
bool firsttime = true; |
| 156 |
|
int junkframe = 0; |
| 157 |
|
for (int i = 0; i < nblocks; ++i) { |
| 158 |
– |
//std::cerr << "block = " << i << "\n"; |
| 158 |
|
bsMan_->loadBlock(i); |
| 159 |
|
assert(bsMan_->isBlockActive(i)); |
| 160 |
|
SnapshotBlock block1 = bsMan_->getSnapshotBlock(i); |
| 162 |
|
|
| 163 |
|
// go snapshot-by-snapshot through this block: |
| 164 |
|
Snapshot* snap = bsMan_->getSnapshot(j); |
| 165 |
< |
|
| 165 |
> |
|
| 166 |
|
// update the positions and velocities of the atoms belonging |
| 167 |
|
// to rigid bodies: |
| 168 |
< |
|
| 168 |
> |
|
| 169 |
|
updateFrame(j); |
| 170 |
< |
|
| 172 |
< |
|
| 170 |
> |
|
| 171 |
|
// do the forces: |
| 172 |
< |
|
| 173 |
< |
forceMan->calcForces(true, true); |
| 174 |
< |
|
| 172 |
> |
|
| 173 |
> |
forceMan->calcForces(true, true); |
| 174 |
> |
|
| 175 |
|
int index = 0; |
| 176 |
< |
|
| 176 |
> |
|
| 177 |
|
for (mol = info_->beginMolecule(mi); mol != NULL; |
| 178 |
< |
mol = info_->nextMolecule(mi)) { |
| 178 |
> |
mol = info_->nextMolecule(mi)) { |
| 179 |
|
for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) { |
| 180 |
|
RealType mass = atom->getMass(); |
| 181 |
|
Vector3d vel = atom->getVel(); |
| 182 |
< |
RealType kinetic = mass * (vel[0]*vel[0] + vel[1]*vel[1] + vel[2]*vel[2]); |
| 182 |
> |
RealType kinetic = mass * (vel[0]*vel[0] + vel[1]*vel[1] + |
| 183 |
> |
vel[2]*vel[2]); |
| 184 |
|
RealType eatom = (kinetic + atom->getParticlePot())/2.0; |
| 185 |
|
particleEnergies.push_back(eatom); |
| 186 |
|
if(firsttime) |
| 187 |
< |
{ |
| 188 |
< |
AvgE_a_.push_back(eatom); |
| 189 |
< |
}else{ |
| 187 |
> |
{ |
| 188 |
> |
AvgE_a_.push_back(eatom); |
| 189 |
> |
} else { |
| 190 |
|
/* We assume the the number of atoms does not change.*/ |
| 191 |
|
AvgE_a_[index] += eatom; |
| 192 |
|
} |
| 200 |
|
bsMan_->unloadBlock(i); |
| 201 |
|
} |
| 202 |
|
|
| 203 |
< |
int nframes = bsMan_->getNFrames(); |
| 204 |
< |
for (int i = 0; i < AvgE_a_.size(); i++){ |
| 205 |
< |
AvgE_a_[i] /= nframes; |
| 206 |
< |
} |
| 208 |
< |
|
| 209 |
< |
|
| 210 |
< |
|
| 211 |
< |
int frame = 0; |
| 212 |
< |
|
| 213 |
< |
// Do it again to compute G^(kappa)(t) for x,y,z |
| 214 |
< |
for (int i = 0; i < nblocks; ++i) { |
| 215 |
< |
//std::cerr << "block = " << i << "\n"; |
| 216 |
< |
bsMan_->loadBlock(i); |
| 217 |
< |
|
| 218 |
< |
assert(bsMan_->isBlockActive(i)); |
| 219 |
< |
SnapshotBlock block1 = bsMan_->getSnapshotBlock(i); |
| 220 |
< |
for (int j = block1.first; j < block1.second; ++j) { |
| 221 |
< |
|
| 222 |
< |
// go snapshot-by-snapshot through this block: |
| 223 |
< |
Snapshot* snap = bsMan_->getSnapshot(j); |
| 224 |
< |
|
| 225 |
< |
updateFrame(j); |
| 226 |
< |
particleEnergies = E_a_[frame]; |
| 227 |
< |
|
| 228 |
< |
int thisAtom = 0; |
| 229 |
< |
Vector3d G_t; |
| 230 |
< |
|
| 231 |
< |
for (mol = info_->beginMolecule(mi); mol != NULL; |
| 232 |
< |
mol = info_->nextMolecule(mi)) { |
| 233 |
< |
for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) { |
| 234 |
< |
|
| 235 |
< |
|
| 236 |
< |
Vector3d pos = atom->getPos(); |
| 237 |
< |
|
| 238 |
< |
|
| 239 |
< |
|
| 240 |
< |
G_t[0] += pos.x()*(particleEnergies[thisAtom]-AvgE_a_[thisAtom]); |
| 241 |
< |
G_t[1] += pos.y()*(particleEnergies[thisAtom]-AvgE_a_[thisAtom]); |
| 242 |
< |
G_t[2] += pos.z()*(particleEnergies[thisAtom]-AvgE_a_[thisAtom]); |
| 243 |
< |
|
| 244 |
< |
thisAtom++; |
| 245 |
< |
} |
| 246 |
< |
} |
| 247 |
< |
|
| 248 |
< |
|
| 249 |
< |
G_t_.push_back(G_t); |
| 250 |
< |
frame++; |
| 251 |
< |
std::cerr <<"Frame: " << frame <<"\t" << G_t << std::endl; |
| 252 |
< |
} |
| 253 |
< |
|
| 254 |
< |
|
| 255 |
< |
|
| 256 |
< |
bsMan_->unloadBlock(i); |
| 257 |
< |
} |
| 203 |
> |
int nframes = bsMan_->getNFrames(); |
| 204 |
> |
for (int i = 0; i < AvgE_a_.size(); i++){ |
| 205 |
> |
AvgE_a_[i] /= nframes; |
| 206 |
> |
} |
| 207 |
|
|
| 208 |
+ |
int frame = 0; |
| 209 |
|
|
| 210 |
< |
|
| 210 |
> |
// Do it again to compute G^(kappa)(t) for x,y,z |
| 211 |
> |
for (int i = 0; i < nblocks; ++i) { |
| 212 |
> |
bsMan_->loadBlock(i); |
| 213 |
> |
assert(bsMan_->isBlockActive(i)); |
| 214 |
> |
SnapshotBlock block1 = bsMan_->getSnapshotBlock(i); |
| 215 |
> |
for (int j = block1.first; j < block1.second; ++j) { |
| 216 |
|
|
| 217 |
+ |
// go snapshot-by-snapshot through this block: |
| 218 |
+ |
Snapshot* snap = bsMan_->getSnapshot(j); |
| 219 |
+ |
|
| 220 |
+ |
// update the positions and velocities of the atoms belonging |
| 221 |
+ |
// to rigid bodies: |
| 222 |
+ |
|
| 223 |
+ |
updateFrame(j); |
| 224 |
+ |
|
| 225 |
+ |
// this needs to be updated to the frame value: |
| 226 |
+ |
particleEnergies = E_a_[j]; |
| 227 |
+ |
|
| 228 |
+ |
int thisAtom = 0; |
| 229 |
+ |
Vector3d G_t; |
| 230 |
+ |
|
| 231 |
+ |
for (mol = info_->beginMolecule(mi); mol != NULL; |
| 232 |
+ |
mol = info_->nextMolecule(mi)) { |
| 233 |
+ |
for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) { |
| 234 |
+ |
|
| 235 |
+ |
Vector3d pos = atom->getPos(); |
| 236 |
+ |
|
| 237 |
+ |
G_t[0] += pos.x()*(particleEnergies[thisAtom]-AvgE_a_[thisAtom]); |
| 238 |
+ |
G_t[1] += pos.y()*(particleEnergies[thisAtom]-AvgE_a_[thisAtom]); |
| 239 |
+ |
G_t[2] += pos.z()*(particleEnergies[thisAtom]-AvgE_a_[thisAtom]); |
| 240 |
+ |
|
| 241 |
+ |
thisAtom++; |
| 242 |
+ |
} |
| 243 |
+ |
} |
| 244 |
+ |
|
| 245 |
+ |
G_t_.push_back(G_t); |
| 246 |
+ |
std::cerr <<"Frame: " << j <<"\t" << G_t << std::endl; |
| 247 |
+ |
} |
| 248 |
+ |
bsMan_->unloadBlock(i); |
| 249 |
+ |
} |
| 250 |
|
} |
| 251 |
|
|
| 252 |
|
|
