From f32b6695566e96056cf5ac858b9be45ef3a530b0 Mon Sep 17 00:00:00 2001
From: Ian Bell <ian.bell@nist.gov>
Date: Sun, 21 Jul 2024 10:25:09 -0400
Subject: [PATCH] Fix 2CLJ virials and test

---
 include/teqp/models/model_potentials/2center_ljf.hpp |  6 +++---
 src/tests/catch_test_modelset.cxx                    | 10 ++++++++++
 2 files changed, 13 insertions(+), 3 deletions(-)

diff --git a/include/teqp/models/model_potentials/2center_ljf.hpp b/include/teqp/models/model_potentials/2center_ljf.hpp
index 212ba08..4610b9d 100644
--- a/include/teqp/models/model_potentials/2center_ljf.hpp
+++ b/include/teqp/models/model_potentials/2center_ljf.hpp
@@ -275,7 +275,7 @@ namespace teqp {
             auto alphar(const TauType& tau, const DeltaType& delta, const double& alpha) const {
 
                 auto eta = forceeval((delta / (a + (1.0 - a) * pow(tau, g))));
-                auto r = (pow(alpha, 2) - 1.0) * log(1.0 - eta) + ((pow(alpha, 2) + 3 * alpha) * eta - 3 * alpha * pow(eta, 2)) / (pow(1.0 - eta, 2));
+                auto r = (pow(alpha, 2) - 1.0) * log(1.0 - eta) + ((pow(alpha, 2) + 3 * alpha) * eta - 3 * alpha * powi(eta, 2)) / (pow(1.0 - eta, 2));
                 return forceeval(r);
             }
         };
@@ -337,7 +337,7 @@ namespace teqp {
                 auto delta_eta = forceeval(rho_dimer_star * eta_red);
                 auto alphar_1 = Hard.alphar(tau, delta_eta, alpha);
                 auto alphar_2 = Attr.alphar(tau, delta, alpha);
-                auto val = alphar_1 + alphar_2;
+                auto val = forceeval(alphar_1 + alphar_2);
                 if (mu_sq != 0.0){
                     val += Pole.alphar(tau, delta, mu_sq);
                 }
@@ -406,7 +406,7 @@ namespace teqp {
             return eos;
         }
 
-        // build the 2-center Lennard-Jones model with dipole
+        // build the 2-center Lennard-Jones model without dipole
         inline auto build_two_center_model(const std::string& model_version, const double& L = 0.0) {
 
             // Get reducing for temperature and density
diff --git a/src/tests/catch_test_modelset.cxx b/src/tests/catch_test_modelset.cxx
index 389b268..2c6bcb7 100644
--- a/src/tests/catch_test_modelset.cxx
+++ b/src/tests/catch_test_modelset.cxx
@@ -30,6 +30,15 @@ auto LKPmethane = [](){
     };
     return spec;
 };
+auto TwoLJF = [](){
+    return R"({
+      "kind": "2CLJF",
+      "model": {
+          "author": "2CLJF_Mecke",
+          "L^*": 0.5
+      }
+    })"_json;
+};
 auto TwoLJFDip = [](){
     return R"({
       "kind": "2CLJF-Dipole",
@@ -106,6 +115,7 @@ std::map<std::string, std::pair<nlohmann::json, std::vector<evalpoint>>> PureFlu
     {"SW_EspindolaHeredia2009",{ { {"kind", "SW_EspindolaHeredia2009"}, {"model", {{"lambda", 1.3}}} }, {VirialReference{1.3, 1e-10} } }},
     {"EXP6_Kataoka1992",{ { {"kind", "EXP6_Kataoka1992"}, {"model", {{"alpha", 12}}} }, {VirialReference{1.3, 1e-10} } }},
     
+    {"2CLJF",{ TwoLJF(), {VirialReference{1.3, 1e-10}} }}, 
     {"2CLJF_Dipole",{ TwoLJFDip(), {} }}, // Virials are not valid because non-integer density exponents
     {"2CLJF_Quad",{ TwoLJFQuad(), {} }}, // Virials are not valid because non-integer density exponents