From 71a0fd853f4a32d40d94f06bd7ca14876dc0070c Mon Sep 17 00:00:00 2001
From: Gareth Tribello <gareth.tribello@gmail.com>
Date: Tue, 6 Jun 2017 09:55:21 +0100
Subject: [PATCH] Added option to caclulate sum and made small change to
 Steinhardt to avoid if statements

---
 src/crystallization/Steinhardt.cpp | 8 ++++++--
 src/multicolvar/ActionVolume.cpp   | 2 +-
 2 files changed, 7 insertions(+), 3 deletions(-)

diff --git a/src/crystallization/Steinhardt.cpp b/src/crystallization/Steinhardt.cpp
index e6f7eb28e..261c07996 100644
--- a/src/crystallization/Steinhardt.cpp
+++ b/src/crystallization/Steinhardt.cpp
@@ -99,14 +99,16 @@ void Steinhardt::calculateVector( multicolvar::AtomValuePack& myatoms ) const {
 
       // The complex number of which we have to take powers
       std::complex<double> com1( distance[0]/dlen,distance[1]/dlen );
+      powered = std::complex<double>(1.0,0.0);
 
       // Do stuff for all other m values
       for(unsigned m=1; m<=tmom; ++m) {
         // Calculate Legendre Polynomial
         poly_ass=deriv_poly( m, distance[2]/dlen, dpoly_ass );
         // Calculate power of complex number
-        if(std::abs(com1)>epsilon) powered=pow(com1,m-1);
-        else powered = std::complex<double>(0.,0.);
+        // if(std::abs(com1)>epsilon) powered=pow(com1,m-1);
+        // else if(m==1) powered=std::complex<double>(1.,0);
+        // else powered = std::complex<double>(0.,0.);
         // Real and imaginary parts of z
         real_z = real(com1*powered); imag_z = imag(com1*powered );
 
@@ -140,6 +142,8 @@ void Steinhardt::calculateVector( multicolvar::AtomValuePack& myatoms ) const {
         accumulateSymmetryFunction( 2+tmom-m, i, pref*sw*tq6, pref*myrealvec, pref*Tensor( -myrealvec,distance ), myatoms );
         // Imaginary part
         accumulateSymmetryFunction( 2+ncomp+tmom-m, i, -pref*sw*itq6, -pref*myimagvec, pref*Tensor( myimagvec,distance ), myatoms );
+        // Calculate next power of complex number
+        powered *= com1;
       }
     }
   }
diff --git a/src/multicolvar/ActionVolume.cpp b/src/multicolvar/ActionVolume.cpp
index d57f55b68..b54290464 100644
--- a/src/multicolvar/ActionVolume.cpp
+++ b/src/multicolvar/ActionVolume.cpp
@@ -28,7 +28,7 @@ void ActionVolume::registerKeywords( Keywords& keys ) {
   VolumeGradientBase::registerKeywords( keys );
   if( keys.reserved("VMEAN") ) keys.use("VMEAN");
   keys.use("MEAN"); keys.use("LESS_THAN"); keys.use("MORE_THAN");
-  keys.use("BETWEEN"); keys.use("HISTOGRAM");
+  keys.use("BETWEEN"); keys.use("HISTOGRAM"); keys.use("SUM");
   keys.add("compulsory","SIGMA","the width of the function to be used for kernel density estimation");
   keys.add("compulsory","KERNEL","gaussian","the type of kernel function to be used");
   keys.addFlag("OUTSIDE",false,"calculate quantities for colvars that are on atoms outside the region of interest");
-- 
GitLab