From 2929bc747437a95551f9a1af8ac03d3a2ce5868f Mon Sep 17 00:00:00 2001
From: "Daniel R. Roe" <daniel.r.roe@gmail.com>
Date: Wed, 10 Aug 2022 10:37:45 -0400
Subject: [PATCH 01/13] Make diffusion correction exactly like unwrap

---
 src/Action_Diffusion.cpp | 55 +++++++++++++++++++++++++++++++++++++---
 1 file changed, 51 insertions(+), 4 deletions(-)

diff --git a/src/Action_Diffusion.cpp b/src/Action_Diffusion.cpp
index e9cf89186b..d8ede4c2f1 100644
--- a/src/Action_Diffusion.cpp
+++ b/src/Action_Diffusion.cpp
@@ -346,20 +346,65 @@ Action::RetType Action_Diffusion::DoAction(int frameNum, ActionFrame& frm) {
       // If the particle moved more than half the box, assume it was imaged
       // and adjust the distance of the total movement with respect to the
       // original frame.
-      if      (delx >  boxcenter_[0]) fixedXYZ[0] -= frm.Frm().BoxCrd().Param(Box::X);
+      mprintf("DEBUG: %4i %6s ortho prevdist %8.3f %8.3f %8.3f", frameNum, avg_x_->legend(), delx, dely, delz);
+/*      if      (delx >  boxcenter_[0]) fixedXYZ[0] -= frm.Frm().BoxCrd().Param(Box::X);
       else if (delx < -boxcenter_[0]) fixedXYZ[0] += frm.Frm().BoxCrd().Param(Box::X);
       if      (dely >  boxcenter_[1]) fixedXYZ[1] -= frm.Frm().BoxCrd().Param(Box::Y);
       else if (dely < -boxcenter_[1]) fixedXYZ[1] += frm.Frm().BoxCrd().Param(Box::Y);
       if      (delz >  boxcenter_[2]) fixedXYZ[2] -= frm.Frm().BoxCrd().Param(Box::Z);
-      else if (delz < -boxcenter_[2]) fixedXYZ[2] += frm.Frm().BoxCrd().Param(Box::Z);
+      else if (delz < -boxcenter_[2]) fixedXYZ[2] += frm.Frm().BoxCrd().Param(Box::Z);*/
+      double t0 = -floor( delx / frm.Frm().BoxCrd().Param(Box::X) + 0.5 ) * frm.Frm().BoxCrd().Param(Box::X);
+      double t1 = -floor( dely / frm.Frm().BoxCrd().Param(Box::Y) + 0.5 ) * frm.Frm().BoxCrd().Param(Box::Y);
+      double t2 = -floor( delz / frm.Frm().BoxCrd().Param(Box::Z) + 0.5 ) * frm.Frm().BoxCrd().Param(Box::Z);
+      fixedXYZ[0] += t0;
+      fixedXYZ[1] += t1;
+      fixedXYZ[2] += t2;
       // Calculate the distance between this "fixed" coordinate
       // and the reference (initial) frame.
       delx = fixedXYZ[0] - iXYZ[0];
       dely = fixedXYZ[1] - iXYZ[1];
       delz = fixedXYZ[2] - iXYZ[2];
+      mprintf("  fixeddist %8.3f %8.3f %8.3f\n", delx, dely, delz);
     } else if ( imageOpt_.ImagingType() == ImageOption::NONORTHO ) {
       // Non-orthorhombic imaging
-      // Calculate distance to previous frames coordinates.
+      Vec3 boxTrans(0.0);
+      Vec3 vtgt( XYZ[0], XYZ[1], XYZ[2] );
+      Vec3 vref( previous_[idx], previous_[idx+1], previous_[idx+2] );
+      // Calculate original distance from the ref (previous) position. 
+      Vec3 vd = vtgt - vref; // dx dy dz
+      double minDistanceSquare = vd.Magnitude2();
+      // Reciprocal coordinates
+      vd = frm.Frm().BoxCrd().FracCell() * vd ; // recip * dxyz
+      double cx = floor(vd[0]);
+      double cy = floor(vd[1]);
+      double cz = floor(vd[2]);
+      // Loop over all possible translations 
+      for (int ix = -1; ix < 2; ++ix) {
+        for (int iy = -1; iy < 2; ++iy) {
+          for (int iz = -1; iz < 2; ++iz) {
+            // Calculate the translation.
+            Vec3 vcc = frm.Frm().BoxCrd().UnitCell().TransposeMult( Vec3( cx+(double)ix,
+                                                                          cy+(double)iy,
+                                                                          cz+(double)iz ) ); // ucell^T * ccxyz
+            // Calc. the potential new coordinate for tgt
+            Vec3 vnew = vtgt - vcc;
+            // Calc. the new distance from the ref (previous) position
+            Vec3 vr = vref - vnew;
+            double distanceSquare = vr.Magnitude2();
+            // If the orig. distance is greater than the new distance, unwrap. 
+            if ( minDistanceSquare > distanceSquare ) {
+                minDistanceSquare = distanceSquare;
+                boxTrans = vcc;
+            }
+          }
+        }
+      } // END loop over ix
+      // Fix position
+      fixedXYZ[0] -= boxTrans[0];
+      fixedXYZ[1] -= boxTrans[1];
+      fixedXYZ[2] -= boxTrans[2];
+      
+/*      // Calculate distance to previous frames coordinates.
       delx = XYZ[0] - previous_[idx  ];
       dely = XYZ[1] - previous_[idx+1];
       delz = XYZ[2] - previous_[idx+2];
@@ -401,17 +446,19 @@ Action::RetType Action_Diffusion::DoAction(int frameNum, ActionFrame& frm) {
         fixedXYZ[0] = minCurr[0];
         fixedXYZ[1] = minCurr[1];
         fixedXYZ[2] = minCurr[2];
-      }
+      }*/
       // Calculate the distance between this "fixed" coordinate
       // and the reference (initial) frame.
       delx = fixedXYZ[0] - iXYZ[0];
       dely = fixedXYZ[1] - iXYZ[1];
       delz = fixedXYZ[2] - iXYZ[2];
     } else {
+      mprintf("DEBUG: %4i %6s noimg prevdist %8.3f %8.3f %8.3f", frameNum, avg_x_->legend(), XYZ[0] - previous_[idx], XYZ[1] - previous_[idx+1], XYZ[2] - previous_[idx+2]);
       // No imaging. Calculate distance from current position to initial position.
       delx = XYZ[0] - iXYZ[0];
       dely = XYZ[1] - iXYZ[1];
       delz = XYZ[2] - iXYZ[2];
+      mprintf("  fixeddist %8.3f %8.3f %8.3f\n", frameNum, avg_x_->legend(), " ", delx, dely, delz);
     }
     // Calc distances for this atom
     double distx = delx * delx;

From 432734e34a6480be4da00077f16496a4b7583e47 Mon Sep 17 00:00:00 2001
From: "Daniel R. Roe" <daniel.r.roe@gmail.com>
Date: Wed, 10 Aug 2022 10:53:49 -0400
Subject: [PATCH 02/13] Add orthogonal unwrap function

---
 src/Action_Diffusion.cpp |  9 +++++++--
 src/ImageRoutines.h      | 17 ++++++++++++++++-
 2 files changed, 23 insertions(+), 3 deletions(-)

diff --git a/src/Action_Diffusion.cpp b/src/Action_Diffusion.cpp
index d8ede4c2f1..52760a1d71 100644
--- a/src/Action_Diffusion.cpp
+++ b/src/Action_Diffusion.cpp
@@ -2,6 +2,7 @@
 #include "Action_Diffusion.h"
 #include "CpptrajStdio.h"
 #include "StringRoutines.h" // validDouble
+#include "ImageRoutines.h"
 #include "DataSet_1D.h" // LinearRegression
 #ifdef TIMER
 # include "Timer.h"
@@ -353,12 +354,16 @@ Action::RetType Action_Diffusion::DoAction(int frameNum, ActionFrame& frm) {
       else if (dely < -boxcenter_[1]) fixedXYZ[1] += frm.Frm().BoxCrd().Param(Box::Y);
       if      (delz >  boxcenter_[2]) fixedXYZ[2] -= frm.Frm().BoxCrd().Param(Box::Z);
       else if (delz < -boxcenter_[2]) fixedXYZ[2] += frm.Frm().BoxCrd().Param(Box::Z);*/
-      double t0 = -floor( delx / frm.Frm().BoxCrd().Param(Box::X) + 0.5 ) * frm.Frm().BoxCrd().Param(Box::X);
+/*      double t0 = -floor( delx / frm.Frm().BoxCrd().Param(Box::X) + 0.5 ) * frm.Frm().BoxCrd().Param(Box::X);
       double t1 = -floor( dely / frm.Frm().BoxCrd().Param(Box::Y) + 0.5 ) * frm.Frm().BoxCrd().Param(Box::Y);
       double t2 = -floor( delz / frm.Frm().BoxCrd().Param(Box::Z) + 0.5 ) * frm.Frm().BoxCrd().Param(Box::Z);
       fixedXYZ[0] += t0;
       fixedXYZ[1] += t1;
-      fixedXYZ[2] += t2;
+      fixedXYZ[2] += t2;*/
+      Vec3 transVec = Image::UnwrapVec_Ortho<double>(Vec3(XYZ), Vec3((&previous_[0])+idx), frm.Frm().BoxCrd().Lengths());
+      fixedXYZ[0] += transVec[0];
+      fixedXYZ[1] += transVec[1];
+      fixedXYZ[2] += transVec[2];
       // Calculate the distance between this "fixed" coordinate
       // and the reference (initial) frame.
       delx = fixedXYZ[0] - iXYZ[0];
diff --git a/src/ImageRoutines.h b/src/ImageRoutines.h
index 78e9115401..f546f87714 100644
--- a/src/ImageRoutines.h
+++ b/src/ImageRoutines.h
@@ -3,11 +3,11 @@
 #include <string>
 #include <vector>
 #include "ImageTypes.h"
+#include "Vec3.h"
 class Topology;
 class Frame;
 class AtomMask;
 class Matrix_3x3;
-class Vec3;
 class Box;
 class Unit;
 
@@ -34,6 +34,21 @@ namespace Image {
   /// Perform orthogonal imaging on given coordinates using given box boundaries
   Vec3 Ortho(Vec3 const&, Vec3 const&, Vec3 const&, Box const&);
 
+  /// \return Vector required for unwrapping in orthogonal box
+  template <typename T> Vec3 UnwrapVec_Ortho(Vec3 const& vtgt, Vec3 const& vref, Vec3 const& boxVec)
+  {
+    T dxyz[3];
+    dxyz[0] = (vtgt[0] - vref[0]);
+    dxyz[1] = (vtgt[1] - vref[1]);
+    dxyz[2] = (vtgt[2] - vref[2]);
+    return Vec3(
+      -floor( dxyz[0] / boxVec[0] + 0.5 ) * boxVec[0],
+      -floor( dxyz[1] / boxVec[1] + 0.5 ) * boxVec[1],
+      -floor( dxyz[2] / boxVec[2] + 0.5 ) * boxVec[2]
+    );
+  }
+    
+
   /// Perform unwrap of non-orthogonal cell using given reference.
   void UnwrapNonortho( Frame&, Frame&, List const&, Unit const&,
                        Matrix_3x3 const&, Matrix_3x3 const& );

From 559830ee032cf9339bc30e011fab16159b7fac77 Mon Sep 17 00:00:00 2001
From: "Daniel R. Roe" <daniel.r.roe@gmail.com>
Date: Wed, 10 Aug 2022 11:07:31 -0400
Subject: [PATCH 03/13] Put the unwrap routines in a separate file

---
 src/Action_Diffusion.cpp | 15 +++++++---
 src/ImageRoutines.h      | 17 +----------
 src/Unwrap.h             | 61 ++++++++++++++++++++++++++++++++++++++++
 3 files changed, 73 insertions(+), 20 deletions(-)
 create mode 100644 src/Unwrap.h

diff --git a/src/Action_Diffusion.cpp b/src/Action_Diffusion.cpp
index 52760a1d71..dbd6750a1d 100644
--- a/src/Action_Diffusion.cpp
+++ b/src/Action_Diffusion.cpp
@@ -2,12 +2,14 @@
 #include "Action_Diffusion.h"
 #include "CpptrajStdio.h"
 #include "StringRoutines.h" // validDouble
-#include "ImageRoutines.h"
+#include "Unwrap.h"
 #include "DataSet_1D.h" // LinearRegression
 #ifdef TIMER
 # include "Timer.h"
 #endif
 
+using namespace Cpptraj;
+
 // CONSTRUCTOR
 Action_Diffusion::Action_Diffusion() :
   avg_x_(0), avg_y_(0), avg_z_(0), avg_r_(0), avg_a_(0),
@@ -360,7 +362,7 @@ Action::RetType Action_Diffusion::DoAction(int frameNum, ActionFrame& frm) {
       fixedXYZ[0] += t0;
       fixedXYZ[1] += t1;
       fixedXYZ[2] += t2;*/
-      Vec3 transVec = Image::UnwrapVec_Ortho<double>(Vec3(XYZ), Vec3((&previous_[0])+idx), frm.Frm().BoxCrd().Lengths());
+      Vec3 transVec = Unwrap::UnwrapVec_Ortho<double>(Vec3(XYZ), Vec3((&previous_[0])+idx), frm.Frm().BoxCrd().Lengths());
       fixedXYZ[0] += transVec[0];
       fixedXYZ[1] += transVec[1];
       fixedXYZ[2] += transVec[2];
@@ -372,10 +374,14 @@ Action::RetType Action_Diffusion::DoAction(int frameNum, ActionFrame& frm) {
       mprintf("  fixeddist %8.3f %8.3f %8.3f\n", delx, dely, delz);
     } else if ( imageOpt_.ImagingType() == ImageOption::NONORTHO ) {
       // Non-orthorhombic imaging
-      Vec3 boxTrans(0.0);
+      //Vec3 boxTrans(0.0);
       Vec3 vtgt( XYZ[0], XYZ[1], XYZ[2] );
       Vec3 vref( previous_[idx], previous_[idx+1], previous_[idx+2] );
-      // Calculate original distance from the ref (previous) position. 
+      Vec3 transVec = Unwrap::UnwrapVec_Nonortho<double>(vtgt, vref, frm.Frm().BoxCrd().UnitCell(), frm.Frm().BoxCrd().FracCell());
+      fixedXYZ[0] += transVec[0];
+      fixedXYZ[1] += transVec[1];
+      fixedXYZ[2] += transVec[2];
+/*      // Calculate original distance from the ref (previous) position. 
       Vec3 vd = vtgt - vref; // dx dy dz
       double minDistanceSquare = vd.Magnitude2();
       // Reciprocal coordinates
@@ -408,6 +414,7 @@ Action::RetType Action_Diffusion::DoAction(int frameNum, ActionFrame& frm) {
       fixedXYZ[0] -= boxTrans[0];
       fixedXYZ[1] -= boxTrans[1];
       fixedXYZ[2] -= boxTrans[2];
+*/
       
 /*      // Calculate distance to previous frames coordinates.
       delx = XYZ[0] - previous_[idx  ];
diff --git a/src/ImageRoutines.h b/src/ImageRoutines.h
index f546f87714..78e9115401 100644
--- a/src/ImageRoutines.h
+++ b/src/ImageRoutines.h
@@ -3,11 +3,11 @@
 #include <string>
 #include <vector>
 #include "ImageTypes.h"
-#include "Vec3.h"
 class Topology;
 class Frame;
 class AtomMask;
 class Matrix_3x3;
+class Vec3;
 class Box;
 class Unit;
 
@@ -34,21 +34,6 @@ namespace Image {
   /// Perform orthogonal imaging on given coordinates using given box boundaries
   Vec3 Ortho(Vec3 const&, Vec3 const&, Vec3 const&, Box const&);
 
-  /// \return Vector required for unwrapping in orthogonal box
-  template <typename T> Vec3 UnwrapVec_Ortho(Vec3 const& vtgt, Vec3 const& vref, Vec3 const& boxVec)
-  {
-    T dxyz[3];
-    dxyz[0] = (vtgt[0] - vref[0]);
-    dxyz[1] = (vtgt[1] - vref[1]);
-    dxyz[2] = (vtgt[2] - vref[2]);
-    return Vec3(
-      -floor( dxyz[0] / boxVec[0] + 0.5 ) * boxVec[0],
-      -floor( dxyz[1] / boxVec[1] + 0.5 ) * boxVec[1],
-      -floor( dxyz[2] / boxVec[2] + 0.5 ) * boxVec[2]
-    );
-  }
-    
-
   /// Perform unwrap of non-orthogonal cell using given reference.
   void UnwrapNonortho( Frame&, Frame&, List const&, Unit const&,
                        Matrix_3x3 const&, Matrix_3x3 const& );
diff --git a/src/Unwrap.h b/src/Unwrap.h
new file mode 100644
index 0000000000..c657ae7301
--- /dev/null
+++ b/src/Unwrap.h
@@ -0,0 +1,61 @@
+#ifndef INC_UNWRAP_H
+#define INC_UNWRAP_H
+#include "Matrix_3x3.h"
+#include "Vec3.h"
+namespace Cpptraj {
+namespace Unwrap {
+  /// \return Vector required for unwrapping in orthogonal box
+  template <typename T> Vec3 UnwrapVec_Ortho(Vec3 const& vtgt, Vec3 const& vref, Vec3 const& boxVec)
+  {
+    T dxyz[3];
+    dxyz[0] = (vtgt[0] - vref[0]);
+    dxyz[1] = (vtgt[1] - vref[1]);
+    dxyz[2] = (vtgt[2] - vref[2]);
+    return Vec3(
+      -floor( dxyz[0] / boxVec[0] + 0.5 ) * boxVec[0],
+      -floor( dxyz[1] / boxVec[1] + 0.5 ) * boxVec[1],
+      -floor( dxyz[2] / boxVec[2] + 0.5 ) * boxVec[2]
+    );
+  }
+
+  /// \return Vector required for unwrapping in non-orthogonal box
+  template <typename T> Vec3 UnwrapVec_Nonortho(Vec3 const& vtgt, Vec3 const& vref, Matrix_3x3 const& ucell, Matrix_3x3 const& frac)
+  {
+    Vec3 boxTrans(0.0);
+    // Calculate original distance from the ref (previous) position. 
+    Vec3 vd = vtgt - vref; // dx dy dz
+    T minDistanceSquare = vd.Magnitude2();
+    // Reciprocal coordinates
+    vd = frac * vd ; // recip * dxyz
+    T cx = floor(vd[0]);
+    T cy = floor(vd[1]);
+    T cz = floor(vd[2]);
+    // Loop over all possible translations 
+    for (int ix = -1; ix < 2; ++ix) {
+      for (int iy = -1; iy < 2; ++iy) {
+        for (int iz = -1; iz < 2; ++iz) {
+          // Calculate the translation.
+          Vec3 vcc = ucell.TransposeMult( Vec3( cx+(T)ix, 
+                                                cy+(T)iy, 
+                                                cz+(T)iz ) ); // ucell^T * ccxyz
+          // Calc. the potential new coordinate for tgt
+          Vec3 vnew = vtgt - vcc; 
+          // Calc. the new distance from the ref (previous) position
+          Vec3 vr = vref - vnew; 
+          T distanceSquare = vr.Magnitude2();
+          // If the orig. distance is greater than the new distance, unwrap. 
+          if ( minDistanceSquare > distanceSquare ) {
+              minDistanceSquare = distanceSquare;
+              boxTrans = vcc;
+          }
+        }
+      }
+    }
+    // Translate tgt atoms
+    boxTrans.Neg();
+    return boxTrans;
+  }
+
+} // END namespace Unwrap
+} // END namespace Cpptraj
+#endif

From 4110eb0d0ae61bac69409dfa8b535fbcd57d3108 Mon Sep 17 00:00:00 2001
From: "Daniel R. Roe" <daniel.r.roe@gmail.com>
Date: Wed, 10 Aug 2022 11:14:29 -0400
Subject: [PATCH 04/13] Use the unwrap functions

---
 src/ImageRoutines.cpp | 15 ++++++++++-----
 src/cpptrajdepend     |  4 ++--
 2 files changed, 12 insertions(+), 7 deletions(-)

diff --git a/src/ImageRoutines.cpp b/src/ImageRoutines.cpp
index e79b5746a3..20484e3ec9 100644
--- a/src/ImageRoutines.cpp
+++ b/src/ImageRoutines.cpp
@@ -2,6 +2,7 @@
 #include "ImageRoutines.h"
 #include "DistRoutines.h"
 #include "CpptrajStdio.h"
+#include "Unwrap.h"
 #include "Image_List.h"
 #include "Image_List_Pair_CoM.h"
 #include "Image_List_Pair_Geom.h"
@@ -11,6 +12,8 @@
 #include "Image_List_Unit_First.h"
 #include "Image_List_Mask.h"
 
+using namespace Cpptraj;
+
 /** \return Empty list for imaging. */
 Image::List* Image::CreateImageList(Mode modeIn, bool useMass, bool center) {
   Image::List* listOut = 0;
@@ -247,7 +250,7 @@ Vec3 Image::Ortho(Vec3 const& Coord, Vec3 const& bp, Vec3 const& bm, Box const&
 // Image::UnwrapNonortho()
 void Image::UnwrapNonortho(Frame& tgtIn, Frame& refIn, List const& AtomPairs,
                            Unit const& allEntities,
-                           Matrix_3x3 const& ucell, Matrix_3x3 const& recip)
+                           Matrix_3x3 const& ucell, Matrix_3x3 const& frac)
 {
   Vec3 vtgt, vref, boxTrans;
   // Loop over atom pairs
@@ -255,7 +258,8 @@ void Image::UnwrapNonortho(Frame& tgtIn, Frame& refIn, List const& AtomPairs,
   {
     vtgt = AtomPairs.GetCoord(idx, tgtIn);
     vref = AtomPairs.GetCoord(idx, refIn);
-
+    boxTrans = Unwrap::UnwrapVec_Nonortho<double>(vtgt, vref, ucell, frac);
+/*
     boxTrans.Zero();
     // Calculate original distance from the ref (previous) position. 
     Vec3 vd = vtgt - vref; // dx dy dz
@@ -287,7 +291,7 @@ void Image::UnwrapNonortho(Frame& tgtIn, Frame& refIn, List const& AtomPairs,
       }
     }
     // Translate tgt atoms
-    boxTrans.Neg();
+    boxTrans.Neg();*/
     AtomPairs.DoTranslation( tgtIn, idx, boxTrans );
   } // END loop over atom pairs 
   // Save new ref positions
@@ -305,11 +309,12 @@ void Image::UnwrapOrtho(Frame& tgtIn, Frame& refIn, List const& AtomPairs,
   {
     vtgt = AtomPairs.GetCoord(idx, tgtIn);
     vref = AtomPairs.GetCoord(idx, refIn);
-
+    boxTrans = Unwrap::UnwrapVec_Ortho<double>(vtgt, vref, boxVec);
+/*
     Vec3 dxyz = vtgt - vref;
     boxTrans[0] = -floor( dxyz[0] / boxVec[0] + 0.5 ) * boxVec[0];
     boxTrans[1] = -floor( dxyz[1] / boxVec[1] + 0.5 ) * boxVec[1];
-    boxTrans[2] = -floor( dxyz[2] / boxVec[2] + 0.5 ) * boxVec[2];
+    boxTrans[2] = -floor( dxyz[2] / boxVec[2] + 0.5 ) * boxVec[2];*/
     // Translate atoms from first to last
     AtomPairs.DoTranslation( tgtIn, idx, boxTrans );
   } // END loop over atom pairs
diff --git a/src/cpptrajdepend b/src/cpptrajdepend
index 8e704c1e10..7b734af701 100644
--- a/src/cpptrajdepend
+++ b/src/cpptrajdepend
@@ -23,7 +23,7 @@ Action_CreateReservoir.o : Action_CreateReservoir.cpp Action.h ActionState.h Act
 Action_DNAionTracker.o : Action_DNAionTracker.cpp Action.h ActionState.h Action_DNAionTracker.h ArgList.h AssociatedData.h Atom.h AtomMask.h AtomType.h BaseIOtype.h Box.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataFile.h DataFileList.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_REF.h Dimension.h DispatchObject.h DistRoutines.h FileIO.h FileName.h FileTypes.h Frame.h ImageOption.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h Segment.h SymbolExporting.h TextFormat.h Timer.h Topology.h TypeNameHolder.h Unit.h Vec3.h
 Action_DSSP.o : Action_DSSP.cpp Action.h ActionState.h Action_DSSP.h ArgList.h AssociatedData.h Atom.h AtomMask.h AtomType.h BaseIOtype.h Box.h CharMask.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataFile.h DataFileList.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_REF.h Dimension.h DispatchObject.h DistRoutines.h FileIO.h FileName.h FileTypes.h Frame.h ImageOption.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h Segment.h SymbolExporting.h TextFormat.h Timer.h Topology.h TypeNameHolder.h Unit.h Vec3.h
 Action_Density.o : Action_Density.cpp Action.h ActionState.h Action_Density.h ArgList.h AssociatedData.h Atom.h AtomMask.h AtomType.h BaseIOtype.h Box.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataFile.h DataFileList.h DataSet.h DataSetList.h DataSet_1D.h DataSet_Coords.h DataSet_Coords_REF.h Dimension.h DispatchObject.h FileIO.h FileName.h FileTypes.h Frame.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h OnlineVarT.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h Segment.h StringRoutines.h SymbolExporting.h TextFormat.h Timer.h Topology.h TypeNameHolder.h Unit.h Vec3.h
-Action_Diffusion.o : Action_Diffusion.cpp Action.h ActionState.h Action_Diffusion.h ArgList.h AssociatedData.h Atom.h AtomMask.h AtomType.h BaseIOtype.h Box.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataFile.h DataFileList.h DataSet.h DataSetList.h DataSet_1D.h DataSet_Coords.h DataSet_Coords_REF.h Dimension.h DispatchObject.h FileIO.h FileName.h FileTypes.h Frame.h ImageOption.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h Segment.h StringRoutines.h SymbolExporting.h TextFormat.h Timer.h Topology.h TypeNameHolder.h Unit.h Vec3.h
+Action_Diffusion.o : Action_Diffusion.cpp Action.h ActionState.h Action_Diffusion.h ArgList.h AssociatedData.h Atom.h AtomMask.h AtomType.h BaseIOtype.h Box.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataFile.h DataFileList.h DataSet.h DataSetList.h DataSet_1D.h DataSet_Coords.h DataSet_Coords_REF.h Dimension.h DispatchObject.h FileIO.h FileName.h FileTypes.h Frame.h ImageOption.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h Segment.h StringRoutines.h SymbolExporting.h TextFormat.h Timer.h Topology.h TypeNameHolder.h Unit.h Unwrap.h Vec3.h
 Action_Dihedral.o : Action_Dihedral.cpp Action.h ActionState.h Action_Dihedral.h ArgList.h AssociatedData.h Atom.h AtomMask.h AtomType.h BaseIOtype.h Box.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataFile.h DataFileList.h DataSet.h DataSetList.h DataSet_1D.h DataSet_Coords.h DataSet_Coords_REF.h DataSet_double.h Dimension.h DispatchObject.h FileIO.h FileName.h FileTypes.h Frame.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h Segment.h SymbolExporting.h TextFormat.h Timer.h Topology.h TorsionRoutines.h TypeNameHolder.h Unit.h Vec3.h
 Action_DihedralRMS.o : Action_DihedralRMS.cpp Action.h ActionState.h Action_DihedralRMS.h ArgList.h AssociatedData.h Atom.h AtomMask.h AtomType.h BaseIOtype.h Box.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataFile.h DataFileList.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_REF.h DihedralSearch.h Dimension.h DispatchObject.h FileIO.h FileName.h FileTypes.h Frame.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReferenceAction.h ReferenceFrame.h ReplicaDimArray.h Residue.h Segment.h SymbolExporting.h TextFormat.h Timer.h Topology.h TorsionRoutines.h TypeNameHolder.h Unit.h Vec3.h
 Action_Dipole.o : Action_Dipole.cpp Action.h ActionState.h Action_Dipole.h ArgList.h ArrayIterator.h AssociatedData.h Atom.h AtomMask.h AtomType.h BaseIOtype.h Box.h CharMask.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataFile.h DataFileList.h DataSet.h DataSetList.h DataSet_3D.h DataSet_Coords.h DataSet_Coords_REF.h DataSet_GridFlt.h Dimension.h DispatchObject.h FileIO.h FileName.h FileTypes.h Frame.h Grid.h GridAction.h GridBin.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h Segment.h SymbolExporting.h TextFormat.h Timer.h Topology.h TypeNameHolder.h Unit.h Vec3.h
@@ -344,7 +344,7 @@ GridAction.o : GridAction.cpp ArgList.h ArrayIterator.h AssociatedData.h Atom.h
 GridBin.o : GridBin.cpp Box.h CpptrajStdio.h GridBin.h Matrix_3x3.h Parallel.h Vec3.h
 HistBin.o : HistBin.cpp Constants.h CpptrajStdio.h Dimension.h HistBin.h
 Hungarian.o : Hungarian.cpp ArrayIterator.h Constants.h CpptrajStdio.h Hungarian.h Matrix.h
-ImageRoutines.o : ImageRoutines.cpp Atom.h AtomMask.h Box.h CoordinateInfo.h CpptrajStdio.h DistRoutines.h Frame.h ImageOption.h ImageRoutines.h ImageTypes.h Image_List.h Image_List_Mask.h Image_List_Pair.h Image_List_Pair_CoM.h Image_List_Pair_First.h Image_List_Pair_Geom.h Image_List_Unit.h Image_List_Unit_CoM.h Image_List_Unit_First.h Image_List_Unit_Geom.h MaskToken.h Matrix_3x3.h Molecule.h NameType.h Parallel.h ReplicaDimArray.h Residue.h Segment.h SymbolExporting.h Unit.h Vec3.h
+ImageRoutines.o : ImageRoutines.cpp Atom.h AtomMask.h Box.h CoordinateInfo.h CpptrajStdio.h DistRoutines.h Frame.h ImageOption.h ImageRoutines.h ImageTypes.h Image_List.h Image_List_Mask.h Image_List_Pair.h Image_List_Pair_CoM.h Image_List_Pair_First.h Image_List_Pair_Geom.h Image_List_Unit.h Image_List_Unit_CoM.h Image_List_Unit_First.h Image_List_Unit_Geom.h MaskToken.h Matrix_3x3.h Molecule.h NameType.h Parallel.h ReplicaDimArray.h Residue.h Segment.h SymbolExporting.h Unit.h Unwrap.h Vec3.h
 Image_List_Mask.o : Image_List_Mask.cpp Atom.h AtomMask.h AtomType.h Box.h Constants.h CoordinateInfo.h CpptrajStdio.h FileName.h Frame.h Image_List.h Image_List_Mask.h MaskToken.h Matrix_3x3.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReplicaDimArray.h Residue.h Segment.h StringRoutines.h SymbolExporting.h Topology.h TypeNameHolder.h Unit.h Vec3.h
 Image_List_Pair.o : Image_List_Pair.cpp Atom.h AtomMask.h AtomType.h Box.h Constants.h CoordinateInfo.h CpptrajStdio.h FileName.h Frame.h Image_List.h Image_List_Pair.h MaskToken.h Matrix_3x3.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReplicaDimArray.h Residue.h Segment.h SymbolExporting.h Topology.h TypeNameHolder.h Unit.h Vec3.h
 Image_List_Unit.o : Image_List_Unit.cpp Atom.h AtomMask.h AtomType.h Box.h Constants.h CoordinateInfo.h CpptrajStdio.h FileName.h Frame.h Image_List.h Image_List_Unit.h MaskToken.h Matrix_3x3.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReplicaDimArray.h Residue.h Segment.h SymbolExporting.h Topology.h TypeNameHolder.h Unit.h Vec3.h

From ff3e13d4aa9ceac908b470a1658785a9f795bfbd Mon Sep 17 00:00:00 2001
From: "Daniel R. Roe" <daniel.r.roe@gmail.com>
Date: Wed, 10 Aug 2022 11:19:52 -0400
Subject: [PATCH 05/13] Get rid of old code

---
 src/Action_Diffusion.cpp | 108 ++-------------------------------------
 1 file changed, 3 insertions(+), 105 deletions(-)

diff --git a/src/Action_Diffusion.cpp b/src/Action_Diffusion.cpp
index dbd6750a1d..f8d5329ca6 100644
--- a/src/Action_Diffusion.cpp
+++ b/src/Action_Diffusion.cpp
@@ -342,26 +342,7 @@ Action::RetType Action_Diffusion::DoAction(int frameNum, ActionFrame& frm) {
     double delx, dely, delz;
     if ( imageOpt_.ImagingType() == ImageOption::ORTHO ) {
       // Orthorhombic imaging
-      // Calculate distance to previous frames coordinates.
-      delx = XYZ[0] - previous_[idx  ];
-      dely = XYZ[1] - previous_[idx+1];
-      delz = XYZ[2] - previous_[idx+2];
-      // If the particle moved more than half the box, assume it was imaged
-      // and adjust the distance of the total movement with respect to the
-      // original frame.
-      mprintf("DEBUG: %4i %6s ortho prevdist %8.3f %8.3f %8.3f", frameNum, avg_x_->legend(), delx, dely, delz);
-/*      if      (delx >  boxcenter_[0]) fixedXYZ[0] -= frm.Frm().BoxCrd().Param(Box::X);
-      else if (delx < -boxcenter_[0]) fixedXYZ[0] += frm.Frm().BoxCrd().Param(Box::X);
-      if      (dely >  boxcenter_[1]) fixedXYZ[1] -= frm.Frm().BoxCrd().Param(Box::Y);
-      else if (dely < -boxcenter_[1]) fixedXYZ[1] += frm.Frm().BoxCrd().Param(Box::Y);
-      if      (delz >  boxcenter_[2]) fixedXYZ[2] -= frm.Frm().BoxCrd().Param(Box::Z);
-      else if (delz < -boxcenter_[2]) fixedXYZ[2] += frm.Frm().BoxCrd().Param(Box::Z);*/
-/*      double t0 = -floor( delx / frm.Frm().BoxCrd().Param(Box::X) + 0.5 ) * frm.Frm().BoxCrd().Param(Box::X);
-      double t1 = -floor( dely / frm.Frm().BoxCrd().Param(Box::Y) + 0.5 ) * frm.Frm().BoxCrd().Param(Box::Y);
-      double t2 = -floor( delz / frm.Frm().BoxCrd().Param(Box::Z) + 0.5 ) * frm.Frm().BoxCrd().Param(Box::Z);
-      fixedXYZ[0] += t0;
-      fixedXYZ[1] += t1;
-      fixedXYZ[2] += t2;*/
+      // Calculate vector needed to correct XYZ for imaging.
       Vec3 transVec = Unwrap::UnwrapVec_Ortho<double>(Vec3(XYZ), Vec3((&previous_[0])+idx), frm.Frm().BoxCrd().Lengths());
       fixedXYZ[0] += transVec[0];
       fixedXYZ[1] += transVec[1];
@@ -371,106 +352,23 @@ Action::RetType Action_Diffusion::DoAction(int frameNum, ActionFrame& frm) {
       delx = fixedXYZ[0] - iXYZ[0];
       dely = fixedXYZ[1] - iXYZ[1];
       delz = fixedXYZ[2] - iXYZ[2];
-      mprintf("  fixeddist %8.3f %8.3f %8.3f\n", delx, dely, delz);
     } else if ( imageOpt_.ImagingType() == ImageOption::NONORTHO ) {
       // Non-orthorhombic imaging
-      //Vec3 boxTrans(0.0);
-      Vec3 vtgt( XYZ[0], XYZ[1], XYZ[2] );
-      Vec3 vref( previous_[idx], previous_[idx+1], previous_[idx+2] );
-      Vec3 transVec = Unwrap::UnwrapVec_Nonortho<double>(vtgt, vref, frm.Frm().BoxCrd().UnitCell(), frm.Frm().BoxCrd().FracCell());
+      // Calculate vector needed to correct XYZ for imaging.
+      Vec3 transVec = Unwrap::UnwrapVec_Nonortho<double>(Vec3(XYZ), Vec3((&previous_[0])+idx), frm.Frm().BoxCrd().UnitCell(), frm.Frm().BoxCrd().FracCell());
       fixedXYZ[0] += transVec[0];
       fixedXYZ[1] += transVec[1];
       fixedXYZ[2] += transVec[2];
-/*      // Calculate original distance from the ref (previous) position. 
-      Vec3 vd = vtgt - vref; // dx dy dz
-      double minDistanceSquare = vd.Magnitude2();
-      // Reciprocal coordinates
-      vd = frm.Frm().BoxCrd().FracCell() * vd ; // recip * dxyz
-      double cx = floor(vd[0]);
-      double cy = floor(vd[1]);
-      double cz = floor(vd[2]);
-      // Loop over all possible translations 
-      for (int ix = -1; ix < 2; ++ix) {
-        for (int iy = -1; iy < 2; ++iy) {
-          for (int iz = -1; iz < 2; ++iz) {
-            // Calculate the translation.
-            Vec3 vcc = frm.Frm().BoxCrd().UnitCell().TransposeMult( Vec3( cx+(double)ix,
-                                                                          cy+(double)iy,
-                                                                          cz+(double)iz ) ); // ucell^T * ccxyz
-            // Calc. the potential new coordinate for tgt
-            Vec3 vnew = vtgt - vcc;
-            // Calc. the new distance from the ref (previous) position
-            Vec3 vr = vref - vnew;
-            double distanceSquare = vr.Magnitude2();
-            // If the orig. distance is greater than the new distance, unwrap. 
-            if ( minDistanceSquare > distanceSquare ) {
-                minDistanceSquare = distanceSquare;
-                boxTrans = vcc;
-            }
-          }
-        }
-      } // END loop over ix
-      // Fix position
-      fixedXYZ[0] -= boxTrans[0];
-      fixedXYZ[1] -= boxTrans[1];
-      fixedXYZ[2] -= boxTrans[2];
-*/
-      
-/*      // Calculate distance to previous frames coordinates.
-      delx = XYZ[0] - previous_[idx  ];
-      dely = XYZ[1] - previous_[idx+1];
-      delz = XYZ[2] - previous_[idx+2];
-      // If the particle moved more than half the box, assume it was imaged
-      // and adjust the distance of the total movement with respect to the
-      // original frame.
-      if (fabs(delx) > boxcenter_[0] ||
-          fabs(dely) > boxcenter_[1] ||
-          fabs(delz) > boxcenter_[2])
-      {
-        // Previous position in Cartesian space
-        Vec3 pCart( previous_[idx], previous_[idx+1], previous_[idx+2] );
-        // Current position in fractional coords
-        Vec3 cFrac = frm.Frm().BoxCrd().FracCell() * Vec3( XYZ[0], XYZ[1], XYZ[2] );
-        // Look for imaged distance closer to previous than current position
-        double minDist2 = frm.Frm().BoxCrd().Param(Box::X) *
-                          frm.Frm().BoxCrd().Param(Box::Y) *
-                          frm.Frm().BoxCrd().Param(Box::Z);
-        Vec3 minCurr(0.0);
-        for (int ix = -1; ix < 2; ix++) {
-          for (int iy = -1; iy < 2; iy++) {
-            for (int iz = -1; iz < 2; iz++) {
-              if (ix != 0 || iy != 0 || iz != 0) { // Ignore current position
-                Vec3 ixyz(ix, iy, iz);
-                // Current position shifted and back in Cartesian space
-                Vec3 IMG = frm.Frm().BoxCrd().UnitCell().TransposeMult(cFrac + ixyz);
-                // Distance from previous position to imaged current position
-                Vec3 dxyz = IMG - pCart;
-                double dist2 = dxyz.Magnitude2();
-                if (dist2 < minDist2) {
-                  minDist2 = dist2;
-                  minCurr = IMG;
-                }
-              }
-            }
-          }
-        }
-        // minCurr contains the shortest imaged position from previous coords
-        fixedXYZ[0] = minCurr[0];
-        fixedXYZ[1] = minCurr[1];
-        fixedXYZ[2] = minCurr[2];
-      }*/
       // Calculate the distance between this "fixed" coordinate
       // and the reference (initial) frame.
       delx = fixedXYZ[0] - iXYZ[0];
       dely = fixedXYZ[1] - iXYZ[1];
       delz = fixedXYZ[2] - iXYZ[2];
     } else {
-      mprintf("DEBUG: %4i %6s noimg prevdist %8.3f %8.3f %8.3f", frameNum, avg_x_->legend(), XYZ[0] - previous_[idx], XYZ[1] - previous_[idx+1], XYZ[2] - previous_[idx+2]);
       // No imaging. Calculate distance from current position to initial position.
       delx = XYZ[0] - iXYZ[0];
       dely = XYZ[1] - iXYZ[1];
       delz = XYZ[2] - iXYZ[2];
-      mprintf("  fixeddist %8.3f %8.3f %8.3f\n", frameNum, avg_x_->legend(), " ", delx, dely, delz);
     }
     // Calc distances for this atom
     double distx = delx * delx;

From 8f1618596038664d67683ec6e5cba9cfb4900e69 Mon Sep 17 00:00:00 2001
From: "Daniel R. Roe" <daniel.r.roe@gmail.com>
Date: Wed, 10 Aug 2022 11:20:30 -0400
Subject: [PATCH 06/13] Remove old code

---
 src/ImageRoutines.cpp | 38 --------------------------------------
 1 file changed, 38 deletions(-)

diff --git a/src/ImageRoutines.cpp b/src/ImageRoutines.cpp
index 20484e3ec9..aa1ef2a28e 100644
--- a/src/ImageRoutines.cpp
+++ b/src/ImageRoutines.cpp
@@ -259,39 +259,6 @@ void Image::UnwrapNonortho(Frame& tgtIn, Frame& refIn, List const& AtomPairs,
     vtgt = AtomPairs.GetCoord(idx, tgtIn);
     vref = AtomPairs.GetCoord(idx, refIn);
     boxTrans = Unwrap::UnwrapVec_Nonortho<double>(vtgt, vref, ucell, frac);
-/*
-    boxTrans.Zero();
-    // Calculate original distance from the ref (previous) position. 
-    Vec3 vd = vtgt - vref; // dx dy dz
-    double minDistanceSquare = vd.Magnitude2();
-    // Reciprocal coordinates
-    vd = recip * vd ; // recip * dxyz
-    double cx = floor(vd[0]);
-    double cy = floor(vd[1]);
-    double cz = floor(vd[2]);
-    // Loop over all possible translations 
-    for (int ix = -1; ix < 2; ++ix) {
-      for (int iy = -1; iy < 2; ++iy) {
-        for (int iz = -1; iz < 2; ++iz) {
-          // Calculate the translation.
-          Vec3 vcc = ucell.TransposeMult( Vec3( cx+(double)ix, 
-                                                cy+(double)iy, 
-                                                cz+(double)iz ) ); // ucell^T * ccxyz
-          // Calc. the potential new coordinate for tgt
-          Vec3 vnew = vtgt - vcc; 
-          // Calc. the new distance from the ref (previous) position
-          Vec3 vr = vref - vnew; 
-          double distanceSquare = vr.Magnitude2();
-          // If the orig. distance is greater than the new distance, unwrap. 
-          if ( minDistanceSquare > distanceSquare ) {
-              minDistanceSquare = distanceSquare;
-              boxTrans = vcc;
-          }
-        }
-      }
-    }
-    // Translate tgt atoms
-    boxTrans.Neg();*/
     AtomPairs.DoTranslation( tgtIn, idx, boxTrans );
   } // END loop over atom pairs 
   // Save new ref positions
@@ -310,11 +277,6 @@ void Image::UnwrapOrtho(Frame& tgtIn, Frame& refIn, List const& AtomPairs,
     vtgt = AtomPairs.GetCoord(idx, tgtIn);
     vref = AtomPairs.GetCoord(idx, refIn);
     boxTrans = Unwrap::UnwrapVec_Ortho<double>(vtgt, vref, boxVec);
-/*
-    Vec3 dxyz = vtgt - vref;
-    boxTrans[0] = -floor( dxyz[0] / boxVec[0] + 0.5 ) * boxVec[0];
-    boxTrans[1] = -floor( dxyz[1] / boxVec[1] + 0.5 ) * boxVec[1];
-    boxTrans[2] = -floor( dxyz[2] / boxVec[2] + 0.5 ) * boxVec[2];*/
     // Translate atoms from first to last
     AtomPairs.DoTranslation( tgtIn, idx, boxTrans );
   } // END loop over atom pairs

From 92b3332d5bae26b1bb8f4abd9319facca3b2a2a5 Mon Sep 17 00:00:00 2001
From: "Daniel R. Roe" <daniel.r.roe@gmail.com>
Date: Wed, 10 Aug 2022 13:19:36 -0400
Subject: [PATCH 07/13] Add openmp diffusion

---
 src/Action_Diffusion.cpp | 53 +++++++++++++++++++++++++++++++---------
 1 file changed, 41 insertions(+), 12 deletions(-)

diff --git a/src/Action_Diffusion.cpp b/src/Action_Diffusion.cpp
index f8d5329ca6..197b3c85c1 100644
--- a/src/Action_Diffusion.cpp
+++ b/src/Action_Diffusion.cpp
@@ -7,6 +7,9 @@
 #ifdef TIMER
 # include "Timer.h"
 #endif
+#ifdef _OPENMP
+# include <omp.h>
+#endif
 
 using namespace Cpptraj;
 
@@ -206,6 +209,15 @@ Action::RetType Action_Diffusion::Init(ArgList& actionArgs, ActionInit& init, in
     mprintf("\tTo calculate diffusion constant from mean squared displacement plots,\n"
             "\t  calculate the slope of MSD vs time and multiply by 10.0/2*N (where N\n"
             "\t  is # of dimensions); this will give units of 1x10^-5 cm^2/s.\n");
+# ifdef _OPENMP
+# pragma omp parallel
+  {
+# pragma omp master
+  {
+  mprintf("\tParallelizing calculation with %i threads.\n", omp_get_num_threads());
+  }
+  }
+# endif
   return Action::OK;
 }
 
@@ -320,30 +332,44 @@ Action::RetType Action_Diffusion::DoAction(int frameNum, ActionFrame& frm) {
     }
   }
   // Diffusion calculation
+  Vec3 boxLengths(0.0);
   if (imageOpt_.ImagingEnabled()) {
     imageOpt_.SetImageType( frm.Frm().BoxCrd().Is_X_Aligned_Ortho() );
     boxcenter_ = frm.Frm().BoxCrd().Center();
+    if (imageOpt_.ImagingType() == ImageOption::ORTHO)
+      boxLengths = frm.Frm().BoxCrd().Lengths();
   }
   // For averaging over selected atoms
   double average2 = 0.0;
   double avgx = 0.0;
   double avgy = 0.0;
   double avgz = 0.0;
-  unsigned int idx = 0; // Index into previous_
-  double fixedXYZ[3];
-  for (AtomMask::const_iterator at = mask_.begin(); at != mask_.end(); ++at, idx += 3)
-  { // Get current and initial coords for this atom.
-    const double* XYZ = frm.Frm().XYZ(*at);
+  int imask;
+//  unsigned int idx = 0; // Index into previous_
+//  double fixedXYZ[3];
+# ifdef _OPENMP
+# pragma omp parallel private(imask) reduction(+ : average2, avgx, avgy, avgz)
+  {
+# pragma omp for
+# endif
+  for (imask = 0; imask < mask_.Nselected(); imask++)
+  //for (AtomMask::const_iterator at = mask_.begin(); at != mask_.end(); ++at, idx += 3)
+  {
+    int at = mask_[imask];
+    int idx = imask * 3; // Index into previous_
+    // Get current and initial coords for this atom.
+    const double* XYZ = frm.Frm().XYZ(at);
+    double fixedXYZ[3];
     fixedXYZ[0] = XYZ[0];
     fixedXYZ[1] = XYZ[1];
     fixedXYZ[2] = XYZ[2];
-    const double* iXYZ = initial_.XYZ(*at);
+    const double* iXYZ = initial_.XYZ(at);
     // Calculate distance from initial position. 
     double delx, dely, delz;
     if ( imageOpt_.ImagingType() == ImageOption::ORTHO ) {
       // Orthorhombic imaging
       // Calculate vector needed to correct XYZ for imaging.
-      Vec3 transVec = Unwrap::UnwrapVec_Ortho<double>(Vec3(XYZ), Vec3((&previous_[0])+idx), frm.Frm().BoxCrd().Lengths());
+      Vec3 transVec = Unwrap::UnwrapVec_Ortho<double>(Vec3(XYZ), Vec3((&previous_[0])+idx), boxLengths);
       fixedXYZ[0] += transVec[0];
       fixedXYZ[1] += transVec[1];
       fixedXYZ[2] += transVec[2];
@@ -383,22 +409,25 @@ Action::RetType Action_Diffusion::DoAction(int frameNum, ActionFrame& frm) {
     // Store distances for this atom
     if (printIndividual_) {
       float fval = (float)distx;
-      atom_x_[*at]->Add(frameNum, &fval);
+      atom_x_[at]->Add(frameNum, &fval);
       fval = (float)disty;
-      atom_y_[*at]->Add(frameNum, &fval);
+      atom_y_[at]->Add(frameNum, &fval);
       fval = (float)distz;
-      atom_z_[*at]->Add(frameNum, &fval);
+      atom_z_[at]->Add(frameNum, &fval);
       fval = (float)dist2;
-      atom_r_[*at]->Add(frameNum, &fval);
+      atom_r_[at]->Add(frameNum, &fval);
       dist2 = sqrt(dist2);
       fval = (float)dist2;
-      atom_a_[*at]->Add(frameNum, &fval);
+      atom_a_[at]->Add(frameNum, &fval);
     }
     // Update the previous coordinate set to match the current coordinates
     previous_[idx  ] = fixedXYZ[0];
     previous_[idx+1] = fixedXYZ[1];
     previous_[idx+2] = fixedXYZ[2];
   } // END loop over selected atoms
+# ifdef _OPENMP
+  } // END pragma omp parallel
+# endif
   // Calc averages
   double dNselected = 1.0 / (double)mask_.Nselected();
   avgx *= dNselected;

From cee856855da4368416cd35aee82fc629a245f9b1 Mon Sep 17 00:00:00 2001
From: "Daniel R. Roe" <daniel.r.roe@gmail.com>
Date: Wed, 10 Aug 2022 13:21:25 -0400
Subject: [PATCH 08/13] Remove deprecated excludeions keyword from GIST. Add
 diffusion to list of openmp-parallelized commands.

---
 doc/cpptraj.lyx | 6 +++++-
 1 file changed, 5 insertions(+), 1 deletion(-)

diff --git a/doc/cpptraj.lyx b/doc/cpptraj.lyx
index 6ce3580f4a..41b27424fc 100644
--- a/doc/cpptraj.lyx
+++ b/doc/cpptraj.lyx
@@ -1117,6 +1117,10 @@ closest
 cluster (pair-wise distance calculation and sieved frame restore only)
 \end_layout
 
+\begin_layout LyX-Code
+diffusion
+\end_layout
+
 \begin_layout LyX-Code
 dssp/secstruct
 \end_layout
@@ -25698,7 +25702,7 @@ gist
 \end_layout
 
 \begin_layout LyX-Code
-     [noimage] [gridcntr <xval> <yval> <zval>] [excludeions]
+     [noimage] [gridcntr <xval> <yval> <zval>]
 \end_layout
 
 \begin_layout LyX-Code

From c63124d9ffa5e59a23880dbe33be88a37b7792b6 Mon Sep 17 00:00:00 2001
From: "Daniel R. Roe" <daniel.r.roe@gmail.com>
Date: Wed, 10 Aug 2022 13:48:23 -0400
Subject: [PATCH 09/13] Improve the speed of unwrapping by only calculating the
 unwrap vector if object has moved more than half a box length in any
 dimension

---
 src/Action_Diffusion.cpp | 12 +++++-------
 src/Action_Diffusion.h   |  1 -
 src/Action_Unwrap.cpp    |  7 ++++---
 src/ImageRoutines.cpp    |  9 +++++----
 src/ImageRoutines.h      |  4 ++--
 src/Unwrap.h             | 26 +++++++++++++++++++-------
 6 files changed, 35 insertions(+), 24 deletions(-)

diff --git a/src/Action_Diffusion.cpp b/src/Action_Diffusion.cpp
index 197b3c85c1..9eae14846b 100644
--- a/src/Action_Diffusion.cpp
+++ b/src/Action_Diffusion.cpp
@@ -27,7 +27,6 @@ Action_Diffusion::Action_Diffusion() :
   debug_(0),
   outputx_(0), outputy_(0), outputz_(0), outputr_(0), outputa_(0),
   diffout_(0),
-  boxcenter_(0.0),
   masterDSL_(0)
 {}
 
@@ -332,12 +331,11 @@ Action::RetType Action_Diffusion::DoAction(int frameNum, ActionFrame& frm) {
     }
   }
   // Diffusion calculation
-  Vec3 boxLengths(0.0);
+  Vec3 boxLengths(0.0), limxyz(0.0);
   if (imageOpt_.ImagingEnabled()) {
     imageOpt_.SetImageType( frm.Frm().BoxCrd().Is_X_Aligned_Ortho() );
-    boxcenter_ = frm.Frm().BoxCrd().Center();
-    if (imageOpt_.ImagingType() == ImageOption::ORTHO)
-      boxLengths = frm.Frm().BoxCrd().Lengths();
+    boxLengths = frm.Frm().BoxCrd().Lengths();
+    limxyz = boxLengths / 2.0;
   }
   // For averaging over selected atoms
   double average2 = 0.0;
@@ -369,7 +367,7 @@ Action::RetType Action_Diffusion::DoAction(int frameNum, ActionFrame& frm) {
     if ( imageOpt_.ImagingType() == ImageOption::ORTHO ) {
       // Orthorhombic imaging
       // Calculate vector needed to correct XYZ for imaging.
-      Vec3 transVec = Unwrap::UnwrapVec_Ortho<double>(Vec3(XYZ), Vec3((&previous_[0])+idx), boxLengths);
+      Vec3 transVec = Unwrap::UnwrapVec_Ortho<double>(Vec3(XYZ), Vec3((&previous_[0])+idx), boxLengths, limxyz);
       fixedXYZ[0] += transVec[0];
       fixedXYZ[1] += transVec[1];
       fixedXYZ[2] += transVec[2];
@@ -381,7 +379,7 @@ Action::RetType Action_Diffusion::DoAction(int frameNum, ActionFrame& frm) {
     } else if ( imageOpt_.ImagingType() == ImageOption::NONORTHO ) {
       // Non-orthorhombic imaging
       // Calculate vector needed to correct XYZ for imaging.
-      Vec3 transVec = Unwrap::UnwrapVec_Nonortho<double>(Vec3(XYZ), Vec3((&previous_[0])+idx), frm.Frm().BoxCrd().UnitCell(), frm.Frm().BoxCrd().FracCell());
+      Vec3 transVec = Unwrap::UnwrapVec_Nonortho<double>(Vec3(XYZ), Vec3((&previous_[0])+idx), frm.Frm().BoxCrd().UnitCell(), frm.Frm().BoxCrd().FracCell(), limxyz);
       fixedXYZ[0] += transVec[0];
       fixedXYZ[1] += transVec[1];
       fixedXYZ[2] += transVec[2];
diff --git a/src/Action_Diffusion.h b/src/Action_Diffusion.h
index 95430aedff..176591c465 100644
--- a/src/Action_Diffusion.h
+++ b/src/Action_Diffusion.h
@@ -49,7 +49,6 @@ class Action_Diffusion : public Action {
     DataFile* outputr_;
     DataFile* outputa_;
     DataFile* diffout_;
-    Vec3 boxcenter_; ///< Hold center of box each frame
     DataSetList* masterDSL_;
     std::string dsname_;
     Dimension Xdim_;
diff --git a/src/Action_Unwrap.cpp b/src/Action_Unwrap.cpp
index 4b639d0ae6..e7f63e7c57 100644
--- a/src/Action_Unwrap.cpp
+++ b/src/Action_Unwrap.cpp
@@ -125,11 +125,12 @@ Action::RetType Action_Unwrap::DoAction(int frameNum, ActionFrame& frm) {
     RefFrame_ = frm.Frm();
     return Action::OK;
   }
- 
+
+  Vec3 limxyz = frm.Frm().BoxCrd().Lengths() / 2.0; 
   if (frm.Frm().BoxCrd().Is_X_Aligned_Ortho())
-    Image::UnwrapOrtho( frm.ModifyFrm(), RefFrame_, *imageList_, allEntities_ );
+    Image::UnwrapOrtho( frm.ModifyFrm(), RefFrame_, *imageList_, allEntities_, limxyz );
   else {
-    Image::UnwrapNonortho( frm.ModifyFrm(), RefFrame_, *imageList_, allEntities_, frm.Frm().BoxCrd().UnitCell(), frm.Frm().BoxCrd().FracCell() );
+    Image::UnwrapNonortho( frm.ModifyFrm(), RefFrame_, *imageList_, allEntities_, frm.Frm().BoxCrd().UnitCell(), frm.Frm().BoxCrd().FracCell(), limxyz );
   }
 
   return Action::MODIFY_COORDS;
diff --git a/src/ImageRoutines.cpp b/src/ImageRoutines.cpp
index aa1ef2a28e..c3baff05ca 100644
--- a/src/ImageRoutines.cpp
+++ b/src/ImageRoutines.cpp
@@ -250,7 +250,8 @@ Vec3 Image::Ortho(Vec3 const& Coord, Vec3 const& bp, Vec3 const& bm, Box const&
 // Image::UnwrapNonortho()
 void Image::UnwrapNonortho(Frame& tgtIn, Frame& refIn, List const& AtomPairs,
                            Unit const& allEntities,
-                           Matrix_3x3 const& ucell, Matrix_3x3 const& frac)
+                           Matrix_3x3 const& ucell, Matrix_3x3 const& frac,
+                           Vec3 const& limxyz)
 {
   Vec3 vtgt, vref, boxTrans;
   // Loop over atom pairs
@@ -258,7 +259,7 @@ void Image::UnwrapNonortho(Frame& tgtIn, Frame& refIn, List const& AtomPairs,
   {
     vtgt = AtomPairs.GetCoord(idx, tgtIn);
     vref = AtomPairs.GetCoord(idx, refIn);
-    boxTrans = Unwrap::UnwrapVec_Nonortho<double>(vtgt, vref, ucell, frac);
+    boxTrans = Unwrap::UnwrapVec_Nonortho<double>(vtgt, vref, ucell, frac, limxyz);
     AtomPairs.DoTranslation( tgtIn, idx, boxTrans );
   } // END loop over atom pairs 
   // Save new ref positions
@@ -267,7 +268,7 @@ void Image::UnwrapNonortho(Frame& tgtIn, Frame& refIn, List const& AtomPairs,
 
 // Image::UnwrapOrtho()
 void Image::UnwrapOrtho(Frame& tgtIn, Frame& refIn, List const& AtomPairs,
-                        Unit const& allEntities)
+                        Unit const& allEntities, Vec3 const& limxyz)
 {
   Vec3 vtgt, vref, boxTrans;
   Vec3 boxVec = tgtIn.BoxCrd().Lengths();
@@ -276,7 +277,7 @@ void Image::UnwrapOrtho(Frame& tgtIn, Frame& refIn, List const& AtomPairs,
   {
     vtgt = AtomPairs.GetCoord(idx, tgtIn);
     vref = AtomPairs.GetCoord(idx, refIn);
-    boxTrans = Unwrap::UnwrapVec_Ortho<double>(vtgt, vref, boxVec);
+    boxTrans = Unwrap::UnwrapVec_Ortho<double>(vtgt, vref, boxVec, limxyz);
     // Translate atoms from first to last
     AtomPairs.DoTranslation( tgtIn, idx, boxTrans );
   } // END loop over atom pairs
diff --git a/src/ImageRoutines.h b/src/ImageRoutines.h
index 78e9115401..bce75e4c9a 100644
--- a/src/ImageRoutines.h
+++ b/src/ImageRoutines.h
@@ -36,9 +36,9 @@ namespace Image {
 
   /// Perform unwrap of non-orthogonal cell using given reference.
   void UnwrapNonortho( Frame&, Frame&, List const&, Unit const&,
-                       Matrix_3x3 const&, Matrix_3x3 const& );
+                       Matrix_3x3 const&, Matrix_3x3 const&, Vec3 const& );
   /// Perform unwrap of orthogonal cell using given reference.
-  void UnwrapOrtho( Frame&, Frame&, List const&, Unit const& );
+  void UnwrapOrtho( Frame&, Frame&, List const&, Unit const&, Vec3 const& );
 
   /// Wrap selected atom coords from given frame into primary cell, store in result.
   void WrapToCell0(std::vector<double>&, Frame const&, AtomMask const&,
diff --git a/src/Unwrap.h b/src/Unwrap.h
index c657ae7301..4e08bc46d5 100644
--- a/src/Unwrap.h
+++ b/src/Unwrap.h
@@ -4,26 +4,38 @@
 #include "Vec3.h"
 namespace Cpptraj {
 namespace Unwrap {
+  /// \return True if absolute value of any component of dxyz is greater than limxyz
+  template <typename T> bool vGreaterThan(T const* dxyz, T const* limxyz) {
+    return ( fabs(dxyz[0]) > limxyz[0] ||
+             fabs(dxyz[1]) > limxyz[1] ||
+             fabs(dxyz[2]) > limxyz[2] );
+  }
+
   /// \return Vector required for unwrapping in orthogonal box
-  template <typename T> Vec3 UnwrapVec_Ortho(Vec3 const& vtgt, Vec3 const& vref, Vec3 const& boxVec)
+  template <typename T> Vec3 UnwrapVec_Ortho(Vec3 const& vtgt, Vec3 const& vref, Vec3 const& boxVec, Vec3 const& limxyz)
   {
     T dxyz[3];
     dxyz[0] = (vtgt[0] - vref[0]);
     dxyz[1] = (vtgt[1] - vref[1]);
     dxyz[2] = (vtgt[2] - vref[2]);
-    return Vec3(
-      -floor( dxyz[0] / boxVec[0] + 0.5 ) * boxVec[0],
-      -floor( dxyz[1] / boxVec[1] + 0.5 ) * boxVec[1],
-      -floor( dxyz[2] / boxVec[2] + 0.5 ) * boxVec[2]
-    );
+    if (vGreaterThan<T>(dxyz, limxyz.Dptr()))
+      return Vec3(
+        -floor( dxyz[0] / boxVec[0] + 0.5 ) * boxVec[0],
+        -floor( dxyz[1] / boxVec[1] + 0.5 ) * boxVec[1],
+        -floor( dxyz[2] / boxVec[2] + 0.5 ) * boxVec[2]
+      );
+    else
+      return Vec3(0.0);
   }
 
   /// \return Vector required for unwrapping in non-orthogonal box
-  template <typename T> Vec3 UnwrapVec_Nonortho(Vec3 const& vtgt, Vec3 const& vref, Matrix_3x3 const& ucell, Matrix_3x3 const& frac)
+  template <typename T> Vec3 UnwrapVec_Nonortho(Vec3 const& vtgt, Vec3 const& vref, Matrix_3x3 const& ucell, Matrix_3x3 const& frac, Vec3 const& limxyz)
   {
     Vec3 boxTrans(0.0);
     // Calculate original distance from the ref (previous) position. 
     Vec3 vd = vtgt - vref; // dx dy dz
+    if (!vGreaterThan<T>(vd.Dptr(), limxyz.Dptr()))
+      return boxTrans;
     T minDistanceSquare = vd.Magnitude2();
     // Reciprocal coordinates
     vd = frac * vd ; // recip * dxyz

From 3120e918bf9cb4caa4a1c61042f2f2901c7e9b8a Mon Sep 17 00:00:00 2001
From: "Daniel R. Roe" <daniel.r.roe@gmail.com>
Date: Wed, 10 Aug 2022 13:55:01 -0400
Subject: [PATCH 10/13] 6.14.0. Minor version bump for fixes to diffusion calc

---
 src/Version.h | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/Version.h b/src/Version.h
index 8fa4222147..a1fb809ad0 100644
--- a/src/Version.h
+++ b/src/Version.h
@@ -12,7 +12,7 @@
  * Whenever a number that precedes <revision> is incremented, all subsequent
  * numbers should be reset to 0.
  */
-#define CPPTRAJ_INTERNAL_VERSION "V6.13.0"
+#define CPPTRAJ_INTERNAL_VERSION "V6.14.0"
 /// PYTRAJ relies on this
 #define CPPTRAJ_VERSION_STRING CPPTRAJ_INTERNAL_VERSION
 #endif

From 087863a34ab054a346d778b9e21bab78577f360d Mon Sep 17 00:00:00 2001
From: "Daniel R. Roe" <daniel.r.roe@gmail.com>
Date: Wed, 10 Aug 2022 14:23:47 -0400
Subject: [PATCH 11/13] Openmp parallelize unwrap

---
 src/Action_Unwrap.cpp | 12 ++++++++++++
 src/ImageRoutines.cpp | 40 +++++++++++++++++++++++++++++-----------
 2 files changed, 41 insertions(+), 11 deletions(-)

diff --git a/src/Action_Unwrap.cpp b/src/Action_Unwrap.cpp
index e7f63e7c57..8ec969bb89 100644
--- a/src/Action_Unwrap.cpp
+++ b/src/Action_Unwrap.cpp
@@ -2,6 +2,9 @@
 #include "CpptrajStdio.h"
 #include "ImageRoutines.h"
 #include "Image_List.h"
+#ifdef _OPENMP
+# include <omp.h>
+#endif
 
 // CONSTRUCTOR
 Action_Unwrap::Action_Unwrap() :
@@ -74,6 +77,15 @@ Action::RetType Action_Unwrap::Init(ArgList& actionArgs, ActionInit& init, int d
   else
     mprintf("\tReference is first frame.");
   mprintf("\n");
+  # ifdef _OPENMP
+# pragma omp parallel
+  {
+# pragma omp master
+  {
+  mprintf("\tParallelizing calculation with %i threads.\n", omp_get_num_threads());
+  }
+  }
+# endif
 
   return Action::OK;
 }
diff --git a/src/ImageRoutines.cpp b/src/ImageRoutines.cpp
index c3baff05ca..18ea5966a1 100644
--- a/src/ImageRoutines.cpp
+++ b/src/ImageRoutines.cpp
@@ -253,15 +253,24 @@ void Image::UnwrapNonortho(Frame& tgtIn, Frame& refIn, List const& AtomPairs,
                            Matrix_3x3 const& ucell, Matrix_3x3 const& frac,
                            Vec3 const& limxyz)
 {
-  Vec3 vtgt, vref, boxTrans;
+  int idx;
+  int maxidx = (int)AtomPairs.nEntities();
   // Loop over atom pairs
-  for (unsigned int idx = 0; idx != AtomPairs.nEntities(); idx++)
+# ifdef _OPENMP
+# pragma omp parallel private(idx)
+  {
+# pragma omp for
+# endif
+  for (idx = 0; idx < maxidx; idx++)
   {
-    vtgt = AtomPairs.GetCoord(idx, tgtIn);
-    vref = AtomPairs.GetCoord(idx, refIn);
-    boxTrans = Unwrap::UnwrapVec_Nonortho<double>(vtgt, vref, ucell, frac, limxyz);
+    Vec3 vtgt = AtomPairs.GetCoord(idx, tgtIn);
+    Vec3 vref = AtomPairs.GetCoord(idx, refIn);
+    Vec3 boxTrans = Unwrap::UnwrapVec_Nonortho<double>(vtgt, vref, ucell, frac, limxyz);
     AtomPairs.DoTranslation( tgtIn, idx, boxTrans );
-  } // END loop over atom pairs 
+  } // END loop over atom pairs
+# ifdef _OPENMP
+  }
+# endif
   // Save new ref positions
   refIn.CopyFrom(tgtIn, allEntities);
 }
@@ -270,17 +279,26 @@ void Image::UnwrapNonortho(Frame& tgtIn, Frame& refIn, List const& AtomPairs,
 void Image::UnwrapOrtho(Frame& tgtIn, Frame& refIn, List const& AtomPairs,
                         Unit const& allEntities, Vec3 const& limxyz)
 {
-  Vec3 vtgt, vref, boxTrans;
   Vec3 boxVec = tgtIn.BoxCrd().Lengths();
+  int idx;
+  int maxidx = (int)AtomPairs.nEntities();
   // Loop over atom pairs
-  for (unsigned int idx = 0; idx != AtomPairs.nEntities(); idx++)
+# ifdef _OPENMP
+# pragma omp parallel private(idx)
   {
-    vtgt = AtomPairs.GetCoord(idx, tgtIn);
-    vref = AtomPairs.GetCoord(idx, refIn);
-    boxTrans = Unwrap::UnwrapVec_Ortho<double>(vtgt, vref, boxVec, limxyz);
+# pragma omp for
+# endif
+  for (idx = 0; idx < maxidx; idx++)
+  {
+    Vec3 vtgt = AtomPairs.GetCoord(idx, tgtIn);
+    Vec3 vref = AtomPairs.GetCoord(idx, refIn);
+    Vec3 boxTrans = Unwrap::UnwrapVec_Ortho<double>(vtgt, vref, boxVec, limxyz);
     // Translate atoms from first to last
     AtomPairs.DoTranslation( tgtIn, idx, boxTrans );
   } // END loop over atom pairs
+# ifdef _OPENMP
+  }
+# endif
   // Save new ref positions
   refIn.CopyFrom(tgtIn, allEntities);
 }

From 9e791eb618921ebe95b960220e4930e8b38816db Mon Sep 17 00:00:00 2001
From: "Daniel R. Roe" <daniel.r.roe@gmail.com>
Date: Wed, 10 Aug 2022 14:24:29 -0400
Subject: [PATCH 12/13] Note unwrap is openmp parallelized

---
 doc/cpptraj.lyx | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/doc/cpptraj.lyx b/doc/cpptraj.lyx
index 41b27424fc..f71a76cf9e 100644
--- a/doc/cpptraj.lyx
+++ b/doc/cpptraj.lyx
@@ -1181,6 +1181,10 @@ spam
 surf
 \end_layout
 
+\begin_layout LyX-Code
+unwrap
+\end_layout
+
 \begin_layout LyX-Code
 velocityautocorr
 \end_layout

From 0f0236a8e742b8ae7735bb3a5f2e3f90ae2dce3c Mon Sep 17 00:00:00 2001
From: "Daniel R. Roe" <daniel.r.roe@gmail.com>
Date: Wed, 10 Aug 2022 14:51:56 -0400
Subject: [PATCH 13/13] Remove old code

---
 src/Action_Diffusion.cpp | 3 ---
 1 file changed, 3 deletions(-)

diff --git a/src/Action_Diffusion.cpp b/src/Action_Diffusion.cpp
index 9eae14846b..4363f04112 100644
--- a/src/Action_Diffusion.cpp
+++ b/src/Action_Diffusion.cpp
@@ -343,15 +343,12 @@ Action::RetType Action_Diffusion::DoAction(int frameNum, ActionFrame& frm) {
   double avgy = 0.0;
   double avgz = 0.0;
   int imask;
-//  unsigned int idx = 0; // Index into previous_
-//  double fixedXYZ[3];
 # ifdef _OPENMP
 # pragma omp parallel private(imask) reduction(+ : average2, avgx, avgy, avgz)
   {
 # pragma omp for
 # endif
   for (imask = 0; imask < mask_.Nselected(); imask++)
-  //for (AtomMask::const_iterator at = mask_.begin(); at != mask_.end(); ++at, idx += 3)
   {
     int at = mask_[imask];
     int idx = imask * 3; // Index into previous_