vdw update

Awallace3 · Jul 9, 2024 · 7a2397a · 7a2397a
1 parent 6d9bbd2
commit 7a2397a
Show file tree

Hide file tree

Showing 2 changed files with 46 additions and 15 deletions.
diff --git a/src/constants.hpp b/src/constants.hpp
@@ -1,20 +1,51 @@
 #include <vector>
 
 namespace constants {
+// std::vector<double> vdw_ls = {
+//     // indexing starts at 1 for Hydrogen
+//     /* 0.0,   3.1,    2.65,  4.16,  4.17,   3.89, 3.59, 3.34, 3.19, 3.04, */
+//     0.0,   2.63,   2.65,  4.16,  4.17,   3.89, 3.34, 3.18, 3.09, 3.04,
+//     2.91,  3.73,   4.27,  4.33,  4.2,    4.01, 3.86, 3.71, 3.55, 3.71,
+//     4.65,  4.59,   4.51,  4.44,  3.99,   3.97, 4.23, 4.18, 3.82, 3.76,
+//     4.02,  4.19,   4.2,   4.11,  4.04,   3.93, 3.82, 3.72, 4.54, 4.8151,
+//     4.53,  4.2365, 4.099, 4.076, 3.9953, 3.95, 3.66, 3.82, 3.99, 4.23198,
+//     4.303, 4.276,  4.22,  4.17,  4.08,   3.78, 4.77, 3.14, 3.26, 3.28,
+//     3.3,   3.27,   3.32,  3.4,   3.62,   3.42, 3.26, 3.24, 3.3,  3.26,
+//     3.22,  3.2,    4.21,  4.15,  4.08,   4.02, 3.84, 4.0,  3.92, 3.86,
+//     3.98,  3.91,   4.31,  4.32,  4.097,  4.07, 4.23, 3.9,  4.98, 2.75,
+//     2.85,  2.71,   3.0,   3.28,  3.45,   3.51, 3.47, 3.56, 3.55, 3.76,
+//     3.89,  3.93,   3.78,
+// };
 std::vector<double> vdw_ls = {
     // indexing starts at 1 for Hydrogen
-    /* 0.0,   3.1,    2.65,  4.16,  4.17,   3.89, 3.59, 3.34, 3.19, 3.04, */
-    0.0,   2.63,   2.65,  4.16,  4.17,   3.89, 3.34, 3.18, 3.09, 3.04,
-    2.91,  3.73,   4.27,  4.33,  4.2,    4.01, 3.86, 3.71, 3.55, 3.71,
-    4.65,  4.59,   4.51,  4.44,  3.99,   3.97, 4.23, 4.18, 3.82, 3.76,
-    4.02,  4.19,   4.2,   4.11,  4.04,   3.93, 3.82, 3.72, 4.54, 4.8151,
-    4.53,  4.2365, 4.099, 4.076, 3.9953, 3.95, 3.66, 3.82, 3.99, 4.23198,
-    4.303, 4.276,  4.22,  4.17,  4.08,   3.78, 4.77, 3.14, 3.26, 3.28,
-    3.3,   3.27,   3.32,  3.4,   3.62,   3.42, 3.26, 3.24, 3.3,  3.26,
-    3.22,  3.2,    4.21,  4.15,  4.08,   4.02, 3.84, 4.0,  3.92, 3.86,
-    3.98,  3.91,   4.31,  4.32,  4.097,  4.07, 4.23, 3.9,  4.98, 2.75,
-    2.85,  2.71,   3.0,   3.28,  3.45,   3.51, 3.47, 3.56, 3.55, 3.76,
-    3.89,  3.93,   3.78,
+    -1.0, 2.0786987380036113, 2.6456165756409593,
+    3.420404287078669, 2.891280971950477,
+    3.62827416087903, 3.212534413278308,
+    2.9290754944596338, 2.8723837106958987,
+    2.7778974044230074, 2.9101782332050554,
+    4.28967830478927, 3.269226197042043,
+    3.477096070842404, 3.9684248634614394,
+    3.401507025824091, 3.401507025824091,
+    3.307020719551199, 3.552685115860717,
+    5.196746845009027, 4.365267349807583, 0.0, 0.0, 0.0,
+    0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 3.533787854606139,
+    3.9873221247160173, 3.4959933320969823,
+    3.590479638369873, 3.4581988095878256,
+    3.817246773424813, 5.725870160137219,
+    4.705418052389993, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
+    0.0, 0.0, 0.0, 3.6471714221336082, 4.100705692243487,
+    3.892835818443126, 3.892835818443126,
+    3.7416577284064996, 4.081808430988909,
+    6.481760610320351, 5.06446601622698, 0.0, 0.0, 0.0,
+    0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
+    0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
+    0.0, 3.7038632058973433, 3.817246773424813,
+    3.911733079697704, 3.722760467151921,
+    3.817246773424813, 4.1573974760072225,
+    6.576246916593242, 5.347924935045654, 0.0, 0.0, 0.0,
+    0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
+    0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
+    0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
 };
 std::vector<double> r4r2_ls = {
     // indexing starts at 0 for Hydrogen

diff --git a/src/disp.cpp b/src/disp.cpp
@@ -1276,10 +1276,10 @@ double disp_2B_TT(Ref<VectorXi> pos, py::EigenDRef<MatrixXd> carts,
                     (carts(i, 1) - carts(j, 1)) * (carts(i, 1) - carts(j, 1)) +
                     (carts(i, 2) - carts(j, 2)) * (carts(i, 2) - carts(j, 2)); 
 
-        f6_ij = f_n_TT(b_ij, dis, 6);
-        f8_ij = f_n_TT(b_ij, dis, 8);
+        f6_ij = f_n_TT(b_ij, pow(dis, 0.5), 6);
+        // f8_ij = f_n_TT(b_ij, pow(dis, 0.5), 8);
         C8 = -C6s(i, j) * 3 * pow(Q_A * Q_B, 0.5);
-        de = s6 * f6_ij * C6s(i, j) / pow(dis, 3) + s8 * f8_ij * C8 / pow(dis, 4);
+        de = s6 * f6_ij * C6s(i, j) / pow(dis, 3) + s8 * f6_ij * C8 / pow(dis, 4);
         energy -= de;
         // printf("f6_ij: %f, C6: %f, e6: %f, f8_ij: %f, C8: %f, e8: %f, de: %f, e: %f\n", f6_ij, C6s(i, j), e6, f8_ij, C8, e8, de, energy);
       }