Improve the speed of imaged distance calculations. Check for duplicate atom names in residues.

doc/DocumentChecksums.txt

@@ -1,3 +1,3 @@
 f6f8cb1a79951d80a9d2656fd9c30f55  CpptrajDevelopmentGuide.lyx
-cbbe066a99553d484a9113510c740f1a  cpptraj.lyx
+676d7874e4942cf94f91fee003cf3c28  cpptraj.lyx
 5d9b5b5ed47a3ded57b6464df99b3585  CpptrajManual.lyx

src/Action_AutoImage.cpp

@@ -473,7 +473,7 @@ Action::RetType Action_AutoImage::autoimage_by_vector(int frameNum, ActionFrame&
         Vec3 translated_coord_frac = mobile_unit_center_frac + image_vec;
         // Closest distance of mobile unit center to anchor center
         int ixyz[3];
-        Cpptraj::Dist2_Imaged_Frac( translated_coord_frac, anchor_center_frac, box.UnitCell(), box.FracCell(), ixyz );
+        Cpptraj::Dist2_Imaged_Frac( translated_coord_frac, anchor_center_frac, box.UnitCell(), ixyz );
         if (ixyz[0] != 0 || ixyz[1] != 0 || ixyz[2] != 0) {
           // The reflection is closer to the center, so move the mobile unit.
           translate_unit(fracCoords_, Vec3(ixyz[0], ixyz[1], ixyz[2]), *it);

src/Action_Distance.cpp

@@ -2,6 +2,7 @@
 #include "CpptrajStdio.h"
 #include "DistRoutines.h"
 #include "AssociatedData_NOE.h"
+#include <cmath>
 
 // CONSTRUCTOR
 Action_Distance::Action_Distance() :
@@ -168,7 +169,8 @@ Action::RetType Action_Distance::DoAction(int frameNum, ActionFrame& frm) {
       a1 = frm.Frm().VGeometricCenter( Mask1_ );
   }
 
-  double Dist = DIST(imageOpt_.ImagingType(), a1, a2_, frm.Frm().BoxCrd());
+  double Dist2 = DIST2(imageOpt_.ImagingType(), a1, a2_, frm.Frm().BoxCrd());
+  double Dist = sqrt(Dist2);
 
   dist_->Add(frameNum, &Dist);
 

src/Action_NMRrst.cpp

@@ -551,7 +551,8 @@ Action::RetType Action_NMRrst::DoAction(int frameNum, ActionFrame& frm) {
         a2 = frm.Frm().VGeometricCenter( noe->dMask2_ );
       }
 
-      double Dist = DIST(imageOpt_.ImagingType(), a1, a2, frm.Frm().BoxCrd());
+      double Dist2 = DIST2(imageOpt_.ImagingType(), a1, a2, frm.Frm().BoxCrd());
+      double Dist = sqrt(Dist2);
 
       noe->dist_->Add(frameNum, &Dist);
     }

src/Action_PairDist.cpp

@@ -143,7 +143,8 @@ Action::RetType Action_PairDist::Setup(ActionSetup& setup) {
 void Action_PairDist::BinHist( std::vector<double>& tmp,
                                const double* XYZ1, const double* XYZ2, Box const& box)
 {
-  double Dist = DIST(imageOpt_.ImagingType(), XYZ1, XYZ2, box);
+  double Dist2 = DIST2(imageOpt_.ImagingType(), XYZ1, XYZ2, box);
+  double Dist = sqrt(Dist2);
 
   unsigned long bin = (unsigned long) (Dist / delta_);
   // Resize on the fly if necessary

src/Action_Radial.cpp

@@ -3,7 +3,7 @@
 #include "Action_Radial.h"
 #include "CpptrajStdio.h"
 #include "Constants.h" // FOURTHIRDSPI
-#include "DistRoutines.h"
+#include "Dcalc.h"
 #ifdef _OPENMP
 #  include <omp.h>
 #endif
@@ -483,45 +483,23 @@ void Action_Radial::calcRDF_singleMask(Frame const& frmIn) {
   int outer_max = OuterMask_.Nselected();
   int idx1;
 
-  double boxX = frmIn.BoxCrd().Param(Box::X);
-  double boxY = frmIn.BoxCrd().Param(Box::Y);
-  double boxZ = frmIn.BoxCrd().Param(Box::Z);
-  bool isOrtho = (imageOpt_.ImagingType() == ImageOption::ORTHO &&
-                  boxX != 0.0 && boxY != 0.0 && boxZ != 0.0);
-
   long unsigned* myRDF = &RDF_[0];
 # ifdef _OPENMP
 # pragma omp parallel private(idx1, myRDF)
   {
-    //mprintf("OPENMP: %i threads\n",omp_get_num_threads());
-    myRDF = &(rdf_thread_[omp_get_thread_num()][0]);
+  //mprintf("OPENMP: %i threads\n",omp_get_num_threads());
+  myRDF = &(rdf_thread_[omp_get_thread_num()][0]);
+  Cpptraj::Dcalc dcalc(frmIn, imageOpt_.ImagingType());
 # pragma omp for schedule(dynamic)
+# else
+  Cpptraj::Dcalc dcalc(frmIn, imageOpt_.ImagingType());
 # endif
   for (idx1 = 0; idx1 < outer_max; idx1++) {
     const double* xyz1 = frmIn.XYZ( OuterMask_[idx1] );
+    dcalc.SetXYZ1( xyz1 );
     for (int idx2 = idx1 + 1; idx2 < outer_max; idx2++) {
       const double* xyz2 = frmIn.XYZ( OuterMask_[idx2] );
-      double D2;
-      if (isOrtho) {
-        double x = xyz1[0] - xyz2[0];
-        double y = xyz1[1] - xyz2[1];
-        double z = xyz1[2] - xyz2[2];
-        x = std::fabs(x);
-        y = std::fabs(y);
-        z = std::fabs(z);
-        while (x > boxX) x -= boxX;
-        while (y > boxY) y -= boxY;
-        while (z > boxZ) z -= boxZ;
-        double D = boxX - x;
-        if (D < x) x = D;
-        D = boxY - y;
-        if (D < y) y = D;
-        D = boxZ - z;
-        if (D < z) z = D;
-        D2 = x*x + y*y + z*z;
-      }
-      else
-        D2 = DIST2(imageOpt_.ImagingType(), xyz1, xyz2, frmIn.BoxCrd());
+      double D2 = dcalc.Dist2(xyz2);
       if (D2 <= maximum2_) {
       // NOTE: Can we modify the histogram to store D^2?
         double dist = sqrt(D2);
@@ -542,42 +520,24 @@ void Action_Radial::calcRDF_twoMask(Frame const& frmIn) {
   int outer_max = OuterMask_.Nselected();
   int inner_max = InnerMask_.Nselected();
   int idx1;
-  double boxX = frmIn.BoxCrd().Param(Box::X);
-  double boxY = frmIn.BoxCrd().Param(Box::Y);
-  double boxZ = frmIn.BoxCrd().Param(Box::Z);
-  bool isOrtho = (imageOpt_.ImagingType() == ImageOption::ORTHO &&
-                boxX != 0.0 && boxY != 0.0 && boxZ != 0.0);
 
   long unsigned* myRDF = &RDF_[0];
 # ifdef _OPENMP
 # pragma omp parallel private(idx1, myRDF)
   {
-    //mprintf("OPENMP: %i threads\n",omp_get_num_threads());
-    myRDF = &(rdf_thread_[omp_get_thread_num()][0]);
+  //mprintf("OPENMP: %i threads\n",omp_get_num_threads());
+  myRDF = &(rdf_thread_[omp_get_thread_num()][0]);
+  Cpptraj::Dcalc dcalc(frmIn, imageOpt_.ImagingType());
 # pragma omp for
+# else
+  Cpptraj::Dcalc dcalc(frmIn, imageOpt_.ImagingType());
 # endif
   for (idx1 = 0; idx1 < outer_max; idx1++) {
     const double* xyz1 = frmIn.XYZ( OuterMask_[idx1] );
+    dcalc.SetXYZ1( xyz1 );
     for (int idx2 = 0; idx2 < inner_max; idx2++) {
       const double* xyz2 = frmIn.XYZ( InnerMask_[idx2] );
-      double D2;
-      if (isOrtho) {
-        double x = std::fabs(xyz1[0] - xyz2[0]);
-        double y = std::fabs(xyz1[1] - xyz2[1]);
-        double z = std::fabs(xyz1[2] - xyz2[2]);
-        while (x > boxX) x -= boxX;
-        while (y > boxY) y -= boxY;
-        while (z > boxZ) z -= boxZ;
-        double D = boxX - x;
-        if (D < x) x = D;
-        D = boxY - y;
-        if (D < y) y = D;
-        D = boxZ - z;
-        if (D < z) z = D;
-        D2 = x*x + y*y + z*z;
-      }
-      else
-        D2 = DIST2(imageOpt_.ImagingType(), xyz1, xyz2, frmIn.BoxCrd());
+      double D2 = dcalc.Dist2(xyz2);
       if (D2 <= maximum2_) {
       // NOTE: Can we modify the histogram to store D^2?
         double dist = sqrt(D2);
@@ -658,14 +618,19 @@ Action::RetType Action_Radial::DoAction(int frameNum, ActionFrame& frm) {
     {
     //mprintf("OPENMP: %i threads\n",omp_get_num_threads());
     mythread = omp_get_thread_num();
+    Cpptraj::Dcalc dcalc(frm.Frm(), imageOpt_.ImagingType());
 #   pragma omp for
+#   else
+    Cpptraj::Dcalc dcalc(frm.Frm(), imageOpt_.ImagingType());
 #   endif
     for (nmask1 = 0; nmask1 < outer_max; nmask1++) {
       atom1 = OuterMask_[nmask1];
+      dcalc.SetXYZ1( frm.Frm().XYZ(atom1) );
       for (nmask2 = 0; nmask2 < inner_max; nmask2++) {
         atom2 = InnerMask_[nmask2];
         if ( (*currentParm_)[atom1].MolNum() != (*currentParm_)[atom2].MolNum() ) {
-          D = DIST2( imageOpt_.ImagingType(), frm.Frm().XYZ(atom1), frm.Frm().XYZ(atom2), frm.Frm().BoxCrd());
+          //D = DIST2( imageOpt_.ImagingType(), frm.Frm().XYZ(atom1), frm.Frm().XYZ(atom2), frm.Frm().BoxCrd());
+          D = dcalc.Dist2( frm.Frm().XYZ(atom2) );
           if (D <= maximum2_) {
             // NOTE: Can we modify the histogram to store D^2?
             D = sqrt(D);
@@ -692,17 +657,22 @@ Action::RetType Action_Radial::DoAction(int frameNum, ActionFrame& frm) {
 #   pragma omp parallel private(nmask1,D,idx,mythread)
     {
     mythread = omp_get_thread_num();
+    Cpptraj::Dcalc dcalc(frm.Frm(), imageOpt_.ImagingType());
 #   pragma omp for
+#   else
+    Cpptraj::Dcalc dcalc(frm.Frm(), imageOpt_.ImagingType());
 #   endif
     for (nmask1 = 0; nmask1 < mask1_max; nmask1++)
     {
       AtomMask const& site1 = Sites1_[nmask1];
       Vec3 com1 = getCenter( frm.Frm(), site1 );
+      dcalc.SetXYZ1( com1.Dptr() );
       for (Marray::const_iterator site2 = Sites2_.begin(); site2 != Sites2_.end(); ++site2)
       {
         if (site1 != *site2) {
           Vec3 com2 = getCenter( frm.Frm(), *site2 );
-          D = DIST2(imageOpt_.ImagingType(), com1.Dptr(), com2.Dptr(), frm.Frm().BoxCrd());
+          //D = DIST2(imageOpt_.ImagingType(), com1.Dptr(), com2.Dptr(), frm.Frm().BoxCrd());
+          D = dcalc.Dist2( com2.Dptr() );
           if (D <= maximum2_) {
             D = sqrt(D);
             //mprintf("MASKLOOP: %10i %10i %10.4f\n",atom1,atom2,D);
@@ -734,11 +704,17 @@ Action::RetType Action_Radial::DoAction(int frameNum, ActionFrame& frm) {
 #   pragma omp parallel private(nmask2,atom2,D,idx,mythread)
     {
     mythread = omp_get_thread_num();
+    Cpptraj::Dcalc dcalc(frm.Frm(), imageOpt_.ImagingType());
+    dcalc.SetXYZ1( coord_center.Dptr() );
 #   pragma omp for
+#   else
+    Cpptraj::Dcalc dcalc(frm.Frm(), imageOpt_.ImagingType());
+    dcalc.SetXYZ1( coord_center.Dptr() );
 #   endif
     for (nmask2 = 0; nmask2 < mask2_max; nmask2++) {
       atom2 = InnerMask_[nmask2];
-      D = DIST2(imageOpt_.ImagingType(), coord_center.Dptr(), frm.Frm().XYZ(atom2), frm.Frm().BoxCrd());
+      D = dcalc.Dist2( frm.Frm().XYZ(atom2) );
+      //D = DIST2(imageOpt_.ImagingType(), coord_center.Dptr(), frm.Frm().XYZ(atom2), frm.Frm().BoxCrd());
       if (D <= maximum2_) {
         // NOTE: Can we modify the histogram to store D^2?
         D = sqrt(D);

src/Dcalc.h

@@ -0,0 +1,144 @@
+#ifndef INC_CPPTRAJ_DCALC_H
+#define INC_CPPTRAJ_DCALC_H
+#include "Frame.h" // will bring in Vec3, Matrix_3x3, Box
+#include <limits> // std::numeric_limits
+#include <cmath> // std::fabs
+namespace Cpptraj {
+/** Class for implementing multiple distance calcs for a single frame. */
+class Dcalc {
+  public:
+    /// CONSTRUCTOR
+    Dcalc( Frame const& frmIn, ImageOption::Type itype) :
+      ucell_( frmIn.BoxCrd().UnitCell() ),
+      fcell_( frmIn.BoxCrd().FracCell() ),
+      boxX_( frmIn.BoxCrd().Param(Box::X) ),
+      boxY_( frmIn.BoxCrd().Param(Box::Y) ),
+      boxZ_( frmIn.BoxCrd().Param(Box::Z) ),
+      xyz1_(0),
+      itype_(ImageOption::NO_IMAGE)
+    {
+      if (boxX_ != 0.0 && boxY_ != 0.0 && boxZ_ != 0.0) {
+        itype_ = itype;
+      }
+    }
+    /// Set xyz1 - typically the coordinate from an outer loop
+    void SetXYZ1(const double* xyz1) {
+      xyz1_ = xyz1;
+      if (itype_ == ImageOption::NONORTHO) {
+        f1_ = fcell_ * Vec3(xyz1_);
+        wrap_frac(f1_);
+      }
+    }
+  private:
+    // ---------------------------------------------------------------
+    /// \return Dist^2 for xyz2 from previously set xyz1, orthogonal imaging
+    double Dist2_ortho(const double* xyz2) const
+    {
+      double x = xyz1_[0] - xyz2[0];
+      double y = xyz1_[1] - xyz2[1];
+      double z = xyz1_[2] - xyz2[2];
+      x = std::fabs(x);
+      y = std::fabs(y);
+      z = std::fabs(z);
+      while (x > boxX_) x -= boxX_;
+      while (y > boxY_) y -= boxY_;
+      while (z > boxZ_) z -= boxZ_;
+      double D = boxX_ - x;
+      if (D < x) x = D;
+      D = boxY_ - y;
+      if (D < y) y = D;
+      D = boxZ_ - z;
+      if (D < z) z = D;
+      return (x*x + y*y + z*z);
+    }
+    // ---------------------------------------------------------------
+    /// Bring each fractional coordinate back into the main unit cell (0 to 1)
+    static inline void wrap_frac(Vec3& frac) {
+      for (int idx = 0; idx != 3; idx++) {
+        while (frac[idx] < 0.0)
+          frac[idx] += 1.0;
+        while (frac[idx] > 1.0)
+          frac[idx] -= 1.0;
+      }
+    }
+    /// \return Dist^2 for xyz2 from previously set xyz1, non-orthogonal imaging
+    double Dist2_nonortho(const double* xyz2) const
+    {
+      // Wrap each coordinate in fractional coordinates back to the main unit cell
+      // (so each frac coord is between 0 and 1).
+      Vec3 f2 = fcell_ * Vec3(xyz2);
+      wrap_frac(f2);
+      // Calculate f2 to f1 in fractional space
+      Vec3 fdelta = f2 - f1_;
+      //fdelta.Print("fdelta");
+      // Set the min and max cell indices we need to look for based
+      // on the f1->f2 vector.
+      int imin[3];
+      int imax[3];
+      for (int idx = 0; idx != 3; idx++) {
+        if (fdelta[idx] < 0) {
+          imin[idx] = -1;
+          imax[idx] = 1;
+        } else {
+          imin[idx] = 0;
+          imax[idx] = 2;
+        }
+      }
+      // Calculate wrapped f2 back in Cartesian space
+      Vec3 c2 = ucell_.TransposeMult( f2 );
+      // Searching only images we are close to.
+      // Check each image of f1 to determine which distance is the smallest.
+      double minD2 = std::numeric_limits<double>::max();
+      for (int ix = imin[0]; ix != imax[0]; ix++) {
+        for (int iy = imin[1]; iy != imax[1]; iy++) {
+          for (int iz = imin[2]; iz != imax[2]; iz++) {
+            // Image in frac space
+            Vec3 f1img = f1_ + Vec3(ix, iy, iz);
+            // Image in Cartesian space
+            Vec3 c1img = ucell_.TransposeMult( f1img );
+            // Calculate distance in Cartesian space
+            double dx = (c1img[0] - c2[0]);
+            double dy = (c1img[1] - c2[1]);
+            double dz = (c1img[2] - c2[2]);
+            double D2 = (dx*dx) + (dy*dy) + (dz*dz);
+            if (D2 < minD2) {
+              minD2 = D2;
+            }
+          } // END loop over Z image indices
+        } // END loop over Y image indices
+      } // END loop over X image indices
+      //mprintf("minD %12.4f\n", sqrt(minD2));
+      return minD2;
+    }
+    // ---------------------------------------------------------------
+    /// \return Dist^2 for xyz2 from previously set xyz1, no imaging
+    double Dist2_noimage(const double* xyz2) const
+    {
+      double x = xyz1_[0] - xyz2[0];
+      double y = xyz1_[1] - xyz2[1];
+      double z = xyz1_[2] - xyz2[2];
+      //double D = x*x + y*y + z*z;
+      //fprintf(stdout,"Mask1=%8.3f %8.3f %8.3f Mask2=%8.3f %8.3f %8.3f D=%8.3f\n",
+      //        a1[0],a1[1],a1[2],a2[0],a2[1],a2[2],D);
+      return (x*x + y*y + z*z);
+     }
+  public:
+    /// \return Dist^2 for xyz2 from previously set xyz1
+    double Dist2(const double* xyz2) const {
+      switch (itype_) {
+        case ImageOption::ORTHO : return Dist2_ortho(xyz2);
+        case ImageOption::NONORTHO : return Dist2_nonortho(xyz2);
+        case ImageOption::NO_IMAGE : return Dist2_noimage(xyz2);
+      }
+      return -1.0;
+    }
+  private:
+    Matrix_3x3 const& ucell_;
+    Matrix_3x3 const& fcell_;
+    Vec3 f1_;
+    double boxX_, boxY_, boxZ_;
+    const double* xyz1_;
+    ImageOption::Type itype_;
+};
+}
+#endif