diff --git a/bench/vdw_new.cpp b/bench/vdw_new.cpp
new file mode 100644
index 0000000..d260404
--- /dev/null
+++ b/bench/vdw_new.cpp
@@ -0,0 +1,56 @@
+#include <iostream>
+#include <Eigen/Core>
+
+using namespace Eigen;
+
+#ifndef SCALAR
+#define SCALAR float
+#endif
+
+#ifndef SIZE
+#define SIZE 10000
+#endif
+
+#ifndef REPEAT
+#define REPEAT 10000
+#endif
+
+typedef Matrix<SCALAR, Eigen::Dynamic, 1> Vec;
+
+using namespace std;
+
+SCALAR E_VDW(const Vec &interactions1, const Vec &interactions2)
+{
+  return (interactions2.cwise()/interactions1)
+         .cwise().cube()
+         .cwise().square()
+         .cwise().square()
+         .sum();
+}
+
+int main() 
+{
+  //
+  //          1   2   3   4  ... (interactions)
+  // ka       .   .   .   .  ...
+  // rab      .   .   .   .  ...
+  // energy   .   .   .   .  ...
+  // ...     ... ... ... ... ...
+  // (variables
+  //    for
+  // interaction)
+  //
+  Vec interactions1(SIZE), interactions2(SIZE); // SIZE is the number of vdw interactions in our system
+  // SetupCalculations()
+  SCALAR rab = 1.0;  
+  interactions1.setConstant(2.4);
+  interactions2.setConstant(rab);
+  
+  // Energy()
+  SCALAR energy = 0.0;
+  for (unsigned int i = 0; i<REPEAT; ++i) {
+    energy += E_VDW(interactions1, interactions2);
+    energy *= 1 + 1e-20 * i; // prevent compiler from optimizing the loop
+  }
+  cout << "energy = " << energy << endl;
+}