| 237 |
|
stat[Stats::PRESSURE_TENSOR_ZX] = tensor(2, 0); |
| 238 |
|
stat[Stats::PRESSURE_TENSOR_ZY] = tensor(2, 1); |
| 239 |
|
stat[Stats::PRESSURE_TENSOR_ZZ] = tensor(2, 2); |
| 240 |
< |
|
| 240 |
> |
Vector3d GKappa_t = getThermalHelfand(); |
| 241 |
> |
stat[Stats::THERMAL_HELFANDMOMENT_X] = GKappa_t.x(); |
| 242 |
> |
stat[Stats::THERMAL_HELFANDMOMENT_Y] = GKappa_t.y(); |
| 243 |
> |
stat[Stats::THERMAL_HELFANDMOMENT_Z] = GKappa_t.z(); |
| 244 |
|
|
| 245 |
|
Globals* simParams = info_->getSimParams(); |
| 246 |
|
|
| 420 |
|
return boxDipole; |
| 421 |
|
} |
| 422 |
|
|
| 423 |
+ |
Vector3d Thermo::getThermalHelfand() { |
| 424 |
+ |
Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot(); |
| 425 |
+ |
SimInfo::MoleculeIterator miter; |
| 426 |
+ |
std::vector<Atom*>::iterator aiter; |
| 427 |
+ |
Molecule* mol; |
| 428 |
+ |
Atom* atom; |
| 429 |
+ |
RealType mass; |
| 430 |
+ |
Vector3d velocity; |
| 431 |
+ |
Vector3d x_a; |
| 432 |
+ |
RealType kinetic; |
| 433 |
+ |
RealType potential; |
| 434 |
+ |
RealType eatom; |
| 435 |
+ |
RealType AvgE_a_ = 0; |
| 436 |
+ |
Vector3d GKappa_t = V3Zero; |
| 437 |
+ |
Vector3d ThermalHelfandMoment; |
| 438 |
+ |
|
| 439 |
+ |
for (mol = info_->beginMolecule(miter); mol != NULL; |
| 440 |
+ |
mol = info_->nextMolecule(miter)) { |
| 441 |
|
|
| 442 |
+ |
for (atom = mol->beginAtom(aiter); atom != NULL; |
| 443 |
+ |
atom = mol->nextAtom(aiter)) { |
| 444 |
+ |
|
| 445 |
+ |
mass = atom->getMass(); |
| 446 |
+ |
velocity = atom->getVel(); |
| 447 |
+ |
kinetic = mass * (velocity[0]*velocity[0] + velocity[1]*velocity[1] + |
| 448 |
+ |
velocity[2]*velocity[2]) / PhysicalConstants::energyConvert; |
| 449 |
+ |
potential = atom->getParticlePot(); |
| 450 |
+ |
eatom += (kinetic + potential)/2.0; |
| 451 |
+ |
} |
| 452 |
+ |
} |
| 453 |
+ |
|
| 454 |
+ |
int natoms = info_->getNGlobalAtoms(); |
| 455 |
+ |
#ifdef IS_MPI |
| 456 |
+ |
|
| 457 |
+ |
MPI_Allreduce(&eatom, &AvgE_a_, 1, MPI_REALTYPE, MPI_SUM, |
| 458 |
+ |
MPI_COMM_WORLD); |
| 459 |
+ |
#else |
| 460 |
+ |
AvgE_a_ = eatom; |
| 461 |
+ |
#endif |
| 462 |
+ |
AvgE_a_ = AvgE_a_/RealType(natoms); |
| 463 |
+ |
|
| 464 |
+ |
for (mol = info_->beginMolecule(miter); mol != NULL; |
| 465 |
+ |
mol = info_->nextMolecule(miter)) { |
| 466 |
+ |
|
| 467 |
+ |
for (atom = mol->beginAtom(aiter); atom != NULL; |
| 468 |
+ |
atom = mol->nextAtom(aiter)) { |
| 469 |
+ |
|
| 470 |
+ |
/* We think that x_a is relative to the total box and should be a wrapped coordinate */ |
| 471 |
+ |
x_a = atom->getPos(); |
| 472 |
+ |
currSnapshot->wrapVector(x_a); |
| 473 |
+ |
potential = atom->getParticlePot(); |
| 474 |
+ |
|
| 475 |
+ |
GKappa_t += x_a*(potential-AvgE_a_); |
| 476 |
+ |
} |
| 477 |
+ |
} |
| 478 |
+ |
#ifdef IS_MPI |
| 479 |
+ |
MPI_Allreduce(GKappa_t.getArrayPointer(), ThermalHelfandMoment.getArrayPointer(), 3, |
| 480 |
+ |
MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD); |
| 481 |
+ |
#else |
| 482 |
+ |
ThermalHelfandMoment = GKappa_t; |
| 483 |
+ |
#endif |
| 484 |
+ |
return ThermalHelfandMoment; |
| 485 |
+ |
|
| 486 |
+ |
} |
| 487 |
+ |
|
| 488 |
+ |
|
| 489 |
+ |
|
| 490 |
|
} //end namespace OpenMD |
