Commit dd8025c1 authored by Alexander Stukowski's avatar Alexander Stukowski

Extended the GSD file reader to parse the PBC image information. Extended the...

Extended the GSD file reader to parse the PBC image information. Extended the Unwrap Trajectories modifier to use this information to unwrap particle coordinates.
parent 7b19e3ab
......@@ -49,15 +49,15 @@
<para>
The modifier uses the <link xlink:href="https://en.wikipedia.org/wiki/Periodic_boundary_conditions#Practical_implementation:_continuity_and_the_minimum_image_convention">minimum image principle</link>
to detect transitions of a particle through a periodic boundary of the simulation cell from one frame to the next, i.e. in an incremental fashion.
After unfolding, the continuous trajectories can span distances larger than one cell size.
After unfolding, the continuous trajectories can span distances larger than the cell dimensions.
</para>
<para>
<informalfigure><screenshot><mediaobject><imageobject>
<imagedata fileref="images/modifiers/unwrap_trajectories_panel.png" format="PNG" scale="50" />
</imageobject></mediaobject></screenshot></informalfigure>
In the current program version, you must press the <guibutton>Update</guibutton> button of the modifier to manually trigger the
analysis of the input particle trajectories. The modifier will step through the entire sequence of input simulation frames
and record all crossings through the periodic cell boundaries. Subsequently, it will use this information to unfold the trajectories
In the current program version, you must press the <guibutton>Update</guibutton> button of the modifier to manually trigger an
analysis of the input particle trajectories. The modifier will step through the entire sequence of simulation frames
and detect and record all crossings of the periodic cell boundaries. Subsequently, it will use this information to unfold the trajectories
and replace particle coordinates with unwrapped versions.
</para>
......@@ -67,6 +67,12 @@
input simulation file if available and can be manually set in the <link linkend="scene_objects.simulation_cell">Simulation cell</link> panel.
</para>
<para>
If the <literal>Periodic Image</literal> particle property is present, which some MD simulation codes store in the trajectory file,
then the modifier will directly use this information to unwrap the particle positions. In this case you don't have to use the <guibutton>Update</guibutton> button,
because a detection of the particles crossing through the cell boundaries is not necessary.
</para>
<simplesect>
<title>See also</title>
<para>
......
......@@ -41,12 +41,15 @@ public:
};
/// Default constructor that creates a status object with status StatusType::Success and an empty info text.
PipelineStatus() : _type(Success) {}
PipelineStatus() = default;
/// Constructs a status object with the given status and optional text string describing the status.
PipelineStatus(StatusType t, const QString& text = QString()) :
PipelineStatus(StatusType t, const QString& text = {}) :
_type(t), _text(text) {}
/// Constructs a status object with success status and a text string describing the status.
PipelineStatus(const QString& text) : _text(text) {}
/// Returns the type of status stores in this object.
StatusType type() const { return _type; }
......@@ -70,7 +73,7 @@ public:
private:
/// The status.
StatusType _type;
StatusType _type = Success;
/// A human-readable string describing the status.
QString _text;
......
......@@ -26,7 +26,7 @@
namespace Ovito { namespace Particles { OVITO_BEGIN_INLINE_NAMESPACE(Import) OVITO_BEGIN_INLINE_NAMESPACE(Formats)
IMPLEMENT_OVITO_CLASS(GSDImporter);
IMPLEMENT_OVITO_CLASS(GSDImporter);
/******************************************************************************
* Checks if the given file has format that can be read by this importer.
......@@ -142,9 +142,10 @@ FileSourceImporter::FrameDataPtr GSDImporter::FrameLoader::loadFile(QFile& file)
else
std::fill(typeProperty->dataInt(), typeProperty->dataInt() + typeProperty->size(), 0);
PropertyStorage* massProperty = readOptionalParticleProperty(gsd, "particles/mass", frameNumber, numParticles, ParticlesObject::MassProperty, frameData);
readOptionalParticleProperty(gsd, "particles/mass", frameNumber, numParticles, ParticlesObject::MassProperty, frameData);
readOptionalParticleProperty(gsd, "particles/charge", frameNumber, numParticles, ParticlesObject::ChargeProperty, frameData);
PropertyStorage* velocityProperty = readOptionalParticleProperty(gsd, "particles/velocity", frameNumber, numParticles, ParticlesObject::VelocityProperty, frameData);
readOptionalParticleProperty(gsd, "particles/velocity", frameNumber, numParticles, ParticlesObject::VelocityProperty, frameData);
readOptionalParticleProperty(gsd, "particles/image", frameNumber, numParticles, ParticlesObject::PeriodicImageProperty, frameData);
PropertyStorage* radiusProperty = readOptionalParticleProperty(gsd, "particles/diameter", frameNumber, numParticles, ParticlesObject::RadiusProperty, frameData);
if(radiusProperty) {
// Convert particle diameter to radius by dividing by 2.
......@@ -154,7 +155,9 @@ FileSourceImporter::FrameDataPtr GSDImporter::FrameLoader::loadFile(QFile& file)
if(orientationProperty) {
// Convert quaternion representation from GSD format to internal format.
// Left-shift all quaternion components by one: (W,X,Y,Z) -> (X,Y,Z,W).
std::for_each(orientationProperty->dataQuaternion(), orientationProperty->dataQuaternion() + orientationProperty->size(), [](Quaternion& q) { std::rotate(q.begin(), q.begin() + 1, q.end()); });
std::for_each(orientationProperty->dataQuaternion(), orientationProperty->dataQuaternion() + orientationProperty->size(), [](Quaternion& q) {
std::rotate(q.begin(), q.begin() + 1, q.end());
});
}
// Parse number of bonds.
......@@ -215,7 +218,12 @@ PropertyStorage* GSDImporter::FrameLoader::readOptionalParticleProperty(GSDFile&
if(gsd.hasChunk(chunkName, frameNumber)) {
PropertyPtr prop = ParticlesObject::OOClass().createStandardStorage(numParticles, propertyType, false);
frameData->addParticleProperty(prop);
gsd.readFloatArray(chunkName, frameNumber, prop->dataFloat(), numParticles, prop->componentCount());
if(prop->dataType() == PropertyStorage::Float)
gsd.readFloatArray(chunkName, frameNumber, prop->dataFloat(), numParticles, prop->componentCount());
else if(prop->dataType() == PropertyStorage::Int)
gsd.readIntArray(chunkName, frameNumber, prop->dataInt(), numParticles, prop->componentCount());
else
throw Exception(tr("Particle property '%1' cannot be read from GSD file, because it has an unsupported data type.").arg(prop->name()));
return prop.get();
}
else return nullptr;
......
......@@ -80,6 +80,51 @@ void UnwrapTrajectoriesModifier::evaluatePreliminary(TimePoint time, ModifierApp
******************************************************************************/
void UnwrapTrajectoriesModifier::unwrapParticleCoordinates(TimePoint time, ModifierApplication* modApp, PipelineFlowState& state)
{
const ParticlesObject* inputParticles = state.expectObject<ParticlesObject>();
inputParticles->verifyIntegrity();
// If the periodic image flags particle property is present, use it to unwrap particle positions.
if(const PropertyObject* particlePeriodicImageProperty = inputParticles->getProperty(ParticlesObject::PeriodicImageProperty)) {
// Get current simulation cell.
const SimulationCellObject* simCellObj = state.expectObject<SimulationCellObject>();
const SimulationCell cell = simCellObj->data();
// Make a modifiable copy of the particles object.
ParticlesObject* outputParticles = state.expectMutableObject<ParticlesObject>();
// Make a modifiable copy of the particle position property.
PropertyPtr posProperty = outputParticles->expectMutableProperty(ParticlesObject::PositionProperty)->modifiableStorage();
const Vector3I* pbcShift = particlePeriodicImageProperty->constDataVector3I();
for(Point3& p : posProperty->point3Range()) {
p += cell.matrix() * Vector3(*pbcShift++);
}
// Unwrap bonds by adjusting their PBC shift vectors.
if(outputParticles->bonds()) {
if(ConstPropertyPtr topologyProperty = outputParticles->bonds()->getPropertyStorage(BondsObject::TopologyProperty)) {
outputParticles->makeBondsMutable();
PropertyObject* periodicImageProperty = outputParticles->bonds()->createProperty(BondsObject::PeriodicImageProperty, true);
for(size_t bondIndex = 0; bondIndex < topologyProperty->size(); bondIndex++) {
size_t particleIndex1 = topologyProperty->getInt64Component(bondIndex, 0);
size_t particleIndex2 = topologyProperty->getInt64Component(bondIndex, 1);
if(particleIndex1 >= particlePeriodicImageProperty->size() || particleIndex2 >= particlePeriodicImageProperty->size())
continue;
const Vector3I& particleShift1 = particlePeriodicImageProperty->getVector3I(particleIndex1);
const Vector3I& particleShift2 = particlePeriodicImageProperty->getVector3I(particleIndex2);
periodicImageProperty->dataVector3I()[bondIndex] += particleShift1 - particleShift2;
}
}
}
// After unwrapping the particles, the PBC image flags are obsolete.
// It's time to remove the particle property.
outputParticles->removeProperty(particlePeriodicImageProperty);
state.setStatus(tr("Unwrapping particle positions using stored image information."));
return;
}
// Obtain the precomputed list of periodic cell crossings from the modifier application that is needed to determine
// the unwrapped particle positions.
UnwrapTrajectoriesModifierApplication* myModApp = dynamic_object_cast<UnwrapTrajectoriesModifierApplication>(modApp);
......@@ -114,6 +159,8 @@ void UnwrapTrajectoriesModifier::unwrapParticleCoordinates(TimePoint time, Modif
state.setStatus(PipelineStatus(PipelineStatus::Success, tr("Detected %1 periodic cell boundary crossing(s) of particle trajectories.").arg(unwrapRecords.size())));
if(unwrapRecords.empty()) return;
state.setStatus(tr("Unwrapping particle positions based on detected boundary traversals."));
// Get current simulation cell.
const SimulationCellObject* simCellObj = state.expectObject<SimulationCellObject>();
const SimulationCell cell = simCellObj->data();
......
from ovito.io import import_file
from ovito.modifiers import UnwrapTrajectoriesModifier, ComputePropertyModifier
import numpy as np
pipeline = import_file("../files/gsd/image_flags.gsd")
pipeline.modifiers.append(ComputePropertyModifier(operate_on = "bonds", output_property = "length1", expressions = ["BondLength"]))
pipeline.modifiers.append(UnwrapTrajectoriesModifier())
pipeline.modifiers.append(ComputePropertyModifier(operate_on = "bonds", output_property = "length2", expressions = ["BondLength"]))
data1 = pipeline.source.compute()
data2 = pipeline.compute()
assert(not np.all(data1.particles.bonds["Periodic Image"] == [0,0,0]))
assert(np.all(data2.particles.bonds["Periodic Image"] == [0,0,0]))
assert(np.allclose(data2.particles.bonds["length1"], data2.particles.bonds["length2"]))
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment