> Cool -- feel free to take a look at what I've got and let me know if
> you can spot any possible improvements:
>
> https://github.com/dlonie/avogadro/tree/ENH_intercell_bonds
Hi Dave,
Your algorithm for finding the shortest bond vector doesn't work right for unit
cells with angles other than 90deg. I changed it to something that should
yield correct results but I can't compile it because I get:
/home/tuukka/code/dlonie/libavogadro/src/engines/bsdyengine.cpp:35:34: fatal
error: avogadro/obeigenconv.h: No such file or directory
However I made a patch of my changes, you can see if it compiles.
- Tuukka
diff --git a/libavogadro/src/engines/bsdyengine.cpp b/libavogadro/src/engines/bsdyengine.cpp
index 14ca09c..7daf91d 100644
--- a/libavogadro/src/engines/bsdyengine.cpp
+++ b/libavogadro/src/engines/bsdyengine.cpp
@@ -701,28 +701,30 @@ namespace Avogadro
const Vector3d &begPos = *bond->beginPos();
const Vector3d &endPos = *bond->endPos();
- // Find the shortest bond vector by adjusting the fractional bond
- // vector so that all components are between (-0.5, 0.5]
+ // Find the shortest bond vector
const Vector3d origBondVector (endPos - begPos);
- OpenBabel::vector3 fractionalBondVector =
- cell->CartesianToFractional(Eigen2OB(origBondVector));
- while (fractionalBondVector.x() <= -0.5)
- ++fractionalBondVector.x();
- while (fractionalBondVector.x() > 0.5)
- --fractionalBondVector.x();
- while (fractionalBondVector.y() <= -0.5)
- ++fractionalBondVector.y();
- while (fractionalBondVector.y() > 0.5)
- --fractionalBondVector.y();
- while (fractionalBondVector.z() <= -0.5)
- ++fractionalBondVector.z();
- while (fractionalBondVector.z() > 0.5)
- --fractionalBondVector.z();
- shortestBondVector =
- OB2Eigen(cell->FractionalToCartesian(fractionalBondVector));
- const double shortestBondLength = shortestBondVector.norm();
- const Vector3d normalizedShortestBondVector =
- shortestBondVector / shortestBondLength;
+
+ const OpenBabel::matrix3x3 cellMatrixT = cell->GetCellMatrix().transpose();
+ Vector3d shift(0,0,0);
+ double sqNormDifference = 0;
+
+ // Loop through all neighboring unit cells
+ OpenBabel::vector3 q;
+ for (q.x()=-1; q.x()<2; ++q.x()) {
+ for (q.y()=-1; q.y()<2; ++q.y()) {
+ for (q.z()=-1; q.z()<2; ++q.z()) {
+ // Calculate diff = |d + q1*a + q2*b + q3*c|^2 - |d|^2
+ Vector3d s( (cellMatrixT * q).AsArray() );
+ double diff = 2*s.dot(d) + s.dot(s);
+ if (diff < sqNormDifference) {
+ sqNormDifference = diff;
+ shift = s;
+ }
+ }
+ }
+ }
+
+ shortestBondVector = origBondVector + shift;
// If we're in the original cell, render intracell.
if (shortestBondVector.isApprox(origBondVector, 1e-4)) {
------------------------------------------------------------------------------
For Developers, A Lot Can Happen In A Second.
Boundary is the first to Know...and Tell You.
Monitor Your Applications in Ultra-Fine Resolution. Try it FREE!
http://p.sf.net/sfu/Boundary-d2dvs2
_______________________________________________
Avogadro-devel mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/avogadro-devel
