Fix imaging correction in diffusion calculation. OpenMP parallelize diffusion and unwrap.

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
@@ -1177,6 +1181,10 @@ spam
 surf
 \end_layout
 
+\begin_layout LyX-Code
+unwrap
+\end_layout
+
 \begin_layout LyX-Code
 velocityautocorr
 \end_layout
@@ -25698,7 +25706,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

src/Action_Diffusion.cpp

@@ -2,10 +2,16 @@
 #include "Action_Diffusion.h"
 #include "CpptrajStdio.h"
 #include "StringRoutines.h" // validDouble
+#include "Unwrap.h"
 #include "DataSet_1D.h" // LinearRegression
 #ifdef TIMER
 # include "Timer.h"
 #endif
+#ifdef _OPENMP
+# include <omp.h>
+#endif
+
+using namespace Cpptraj;
 
 // CONSTRUCTOR
 Action_Diffusion::Action_Diffusion() :
@@ -21,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)
 {}
 
@@ -203,6 +208,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;
 }
 
@@ -317,91 +331,55 @@ Action::RetType Action_Diffusion::DoAction(int frameNum, ActionFrame& frm) {
     }
   }
   // Diffusion calculation
+  Vec3 boxLengths(0.0), limxyz(0.0);
   if (imageOpt_.ImagingEnabled()) {
     imageOpt_.SetImageType( frm.Frm().BoxCrd().Is_X_Aligned_Ortho() );
-    boxcenter_ = frm.Frm().BoxCrd().Center();
+    boxLengths = frm.Frm().BoxCrd().Lengths();
+    limxyz = boxLengths / 2.0;
   }
   // 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;
+# ifdef _OPENMP
+# pragma omp parallel private(imask) reduction(+ : average2, avgx, avgy, avgz)
+  {
+# pragma omp for
+# endif
+  for (imask = 0; imask < mask_.Nselected(); imask++)
+  {
+    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 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      (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);
+      // Calculate vector needed to correct XYZ for imaging.
+      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];
       // 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 if ( imageOpt_.ImagingType() == ImageOption::NONORTHO ) {
       // Non-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.
-      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 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(), limxyz);
+      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];
@@ -426,22 +404,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;

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_;

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;
 }
@@ -125,11 +137,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;

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,72 +250,55 @@ 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 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)
   {
-    vtgt = AtomPairs.GetCoord(idx, tgtIn);
-    vref = AtomPairs.GetCoord(idx, refIn);
-
-    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();
+# 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_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);
 }
 
 // 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();
+  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);
-
-    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];
+# 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);
 }