Fix bug when polar term is present but no polar fluids are actually b…

…eing evaluated
usnistgov · Oct 23, 2023 · 3928724 · 3928724
1 parent 1caa2dc
commit 3928724
Show file tree

Hide file tree

Showing 2 changed files with 26 additions and 0 deletions.
diff --git a/include/teqp/models/saft/polar_terms.hpp b/include/teqp/models/saft/polar_terms.hpp
@@ -1143,6 +1143,11 @@ class MultipolarContributionGrayGubbins {
                 alpha2 = get_alpha2(T, rhoN, rhostar, mole_fractions, mu);
                 alpha3 = get_alpha3(T, rhoN, rhostar, mole_fractions, mu);
                 alpha = forceeval(alpha2/(1.0-alpha3/alpha2));
+                // Handle the case where the polar term is present but the contribution is zero
+                if (getbaseval(alpha2) == 0){
+                    alpha2 = 0;
+                    alpha = 0;
+                }
             }
             return MultipolarContributionGubbinsTwuTermsGT<type>{alpha2, alpha3, alpha};
         }

diff --git a/src/tests/catch_test_SAFTpolar.cxx b/src/tests/catch_test_SAFTpolar.cxx
@@ -333,5 +333,26 @@ TEST_CASE("Check Stockmayer critical points with polarity terms", "[SAFTVRMiepol
 //            CHECK(rhostar_guess == Approx(0.221054));
         }
     }
+}
+
+TEST_CASE("Critical points with Stockmayer-LJ model of Mooji", "[polar]"){
+    auto j = R"({"kind": "SAFT-VR-Mie", "model": {"coeffs": [{"m": 1.0, "lambda_r": 12, "lambda_a": 6, "sigma_m": 3.405e-10, "epsilon_over_k": 119.8, "mu_Cm": 5.390822314007849e-30, "Q_Cm2": 0.0, "name": "Stockmayer", "BibTeXKey": "?", "nQ": 1.0, "nmu": 1.0}, {"m": 1.0, "lambda_r": 12, "lambda_a": 6, "sigma_m": 3.405e-10, "epsilon_over_k": 119.8, "mu_Cm": 0.0, "Q_Cm2": 0.0, "name": "LJ", "BibTeXKey": "?", "nQ": 1.0, "nmu": 1.0}], "polar_model": "GrayGubbins+GubbinsTwu"}} )"_json;
+    auto model = teqp::cppinterface::make_model(j);
+    SECTION("Lennard-Jones"){
+        double Tstar_guess = 1.3, rhostar_guess = 0.3;
+        double ek = 119.8, sigma_m = 3.405e-10;
+        nlohmann::json flags = {{"alternative_pure_index", 1}, {"alternative_length", 2}};
+        auto pure = model->solve_pure_critical(Tstar_guess*ek, rhostar_guess/(N_A*pow(sigma_m, 3)), flags);
+        double Tstar = std::get<0>(pure)/ek;
+        int t = 0;
+    }
+    SECTION("Stockmayer"){
+        double Tstar_guess = 2.3, rhostar_guess = 0.3;
+        double ek = 119.8, sigma_m = 3.405e-10;
+        nlohmann::json flags = {{"alternative_pure_index", 0}, {"alternative_length", 2}};
+        auto pure = model->solve_pure_critical(Tstar_guess*ek, rhostar_guess/(N_A*pow(sigma_m, 3)), flags);
+        double Tstar = std::get<0>(pure)/ek;
+        int t = 0;
+    }
 
 }