Merge pull request #143 from usnistgov/virialsfixtest

Virialsfixtest

CMakeLists.txt

@@ -178,6 +178,8 @@ if (NOT TEQP_NO_TEQPCPP)
   target_sources(teqpcpp PRIVATE "${CMAKE_CURRENT_SOURCE_DIR}/src/data/FEANN.cpp")
   target_sources(teqpcpp PRIVATE "${CMAKE_CURRENT_SOURCE_DIR}/src/data/Dufal_assoc.cpp")
   target_sources(teqpcpp PRIVATE "${CMAKE_CURRENT_SOURCE_DIR}/src/data/PCSAFT.cpp")
+
+  target_compile_definitions(teqpcpp PUBLIC -DTEQP_MULTIPRECISION_ENABLED)
 
   # target_compile_options(teqpcpp PRIVATE
   #    $<$<OR:$<CXX_COMPILER_ID:Clang>,$<CXX_COMPILER_ID:AppleClang>,$<CXX_COMPILER_ID:GNU>>:

include/teqp/cpp/deriv_adapter.hpp

@@ -4,6 +4,13 @@
 #include "teqp/cpp/teqpcpp.hpp"
 #include "teqp/exceptions.hpp"
 
+#if defined(TEQP_MULTIPRECISION_ENABLED)
+// Imports from boost
+#include <boost/multiprecision/cpp_bin_float.hpp>
+using namespace boost::multiprecision;
+#include "teqp/finite_derivs.hpp"
+#endif
+
 namespace teqp{
 namespace cppinterface{
 namespace adapter{
@@ -128,6 +135,31 @@ class DerivativeAdapter : public teqp::cppinterface::AbstractModel{
     ARN0_args
 #undef X
 
+    virtual double get_Ar01ep(const double T, const double rho, const EArrayd& molefrac) const  override {
+        using namespace boost::multiprecision;
+        using my_float_t = number<cpp_bin_float<100U>>;
+        auto f = [&](const auto& rhoep){
+            return mp.get_cref().alphar(T, rhoep, molefrac);
+        };
+        return rho*static_cast<double>(centered_diff<1,4>(f, static_cast<my_float_t>(rho), 1e-16*static_cast<my_float_t>(rho)));
+    }
+    virtual double get_Ar02ep(const double T, const double rho, const EArrayd& molefrac) const  override {
+        using namespace boost::multiprecision;
+        using my_float_t = number<cpp_bin_float<100U>>;
+        auto f = [&](const auto& rhoep){
+            return mp.get_cref().alphar(T, rhoep, molefrac);
+        };
+        return rho*rho*static_cast<double>(centered_diff<2,4>(f, static_cast<my_float_t>(rho), 1e-16*static_cast<my_float_t>(rho)));
+    }
+    virtual double get_Ar03ep(const double T, const double rho, const EArrayd& molefrac) const  override {
+        using namespace boost::multiprecision;
+        using my_float_t = number<cpp_bin_float<100U>>;
+        auto f = [&](const auto& rhoep){
+            return mp.get_cref().alphar(T, rhoep, molefrac);
+        };
+        return rho*rho*rho*static_cast<double>(centered_diff<3,4>(f, static_cast<my_float_t>(rho), 1e-16*static_cast<my_float_t>(rho)));
+    }
+
     // Virial derivatives
     virtual double get_B2vir(const double T, const EArrayd& z) const override {
         return VirialDerivatives<decltype(mp.get_cref()), double, EArrayd>::get_B2vir(mp.get_cref(), T, z);

include/teqp/cpp/teqpcpp.hpp

@@ -116,6 +116,10 @@ namespace teqp {
                 ARN0_args
             #undef X
 
+            // Extended precision evaluations
+            virtual double get_Ar01ep(const double, const double, const EArrayd&) const = 0;
+            virtual double get_Ar02ep(const double, const double, const EArrayd&) const = 0;
+            virtual double get_Ar03ep(const double, const double, const EArrayd&) const = 0;
 
             // Virial derivatives
             virtual double get_B2vir(const double T, const EArrayd& z) const = 0;

include/teqp/math/quadrature.hpp

@@ -30,7 +30,7 @@ inline auto quad(const std::function<T(Double)>& F, const Double& a, const Doubl
         {
             5,
             {{0, 1.0/3.0*sqrt(5-2*sqrt(10.0/7.0)), -1.0/3.0*sqrt(5-2*sqrt(10.0/7.0)), 1.0/3.0*sqrt(5+2*sqrt(10.0/7.0)), -1.0/3.0*sqrt(5+2*sqrt(10.0/7.0))},
-            {128/225.0, (322.0+13*sqrt(70))/900.0, (322.0+13*sqrt(70))/900.0, (322.0-13*sqrt(70))/900.0, (322.0-13*sqrt(70))/900.0}}
+            {128/225.0, (322.0+13*sqrt(70.0))/900.0, (322.0+13*sqrt(70.0))/900.0, (322.0-13*sqrt(70.0))/900.0, (322.0-13*sqrt(70.0))/900.0}}
         },
         {
             7,

include/teqp/models/CPA.hpp

@@ -131,7 +131,7 @@ class CPACubic {
     CPACubic(cubic_flag flag, const std::valarray<double> &a0, const std::valarray<double> &bi, const std::valarray<double> &c1, const std::valarray<double> &Tc, const double R_gas, const std::optional<std::vector<std::vector<double>>> & kmat = std::nullopt) : a0(a0), bi(bi), c1(c1), Tc(Tc), R_gas(R_gas), kmat(kmat) {
         switch (flag) {
         case cubic_flag::PR:
-        { delta_1 = 1 + sqrt(2); delta_2 = 1 - sqrt(2); break; }
+        { delta_1 = 1 + sqrt(2.0); delta_2 = 1 - sqrt(2.0); break; }
         case cubic_flag::SRK:
         { delta_1 = 1; delta_2 = 0; break; }
         default:

include/teqp/models/association/association.hpp

@@ -374,11 +374,11 @@ class Association{
         using rDDXtype = std::decay_t<std::common_type_t<typename decltype(Delta)::Scalar, decltype(rhomolar), decltype(molefracs[0])>>; // Type promotion, without the const-ness
         Eigen::MatrixX<rDDXtype> rDDX = rhomolar*N_A*(Delta.array()*D.cast<resulttype>().array()).matrix();
         for (auto j = 0; j < rDDX.rows(); ++j){
-            rDDX.row(j).array() = rDDX.row(j).array()*xj.array();
+            rDDX.row(j).array() = rDDX.row(j).array()*xj.array().template cast<rDDXtype>();
         }
 //        rDDX.rowwise() *= xj;
 
-        Eigen::ArrayX<std::decay_t<rDDXtype>> X = X_init, Xnew;
+        Eigen::ArrayX<std::decay_t<rDDXtype>> X = X_init.template cast<rDDXtype>(), Xnew;
 
         for (auto counter = 0; counter < options.max_iters; ++counter){
             // calculate the new array of non-bonded site fractions X

include/teqp/models/mie/lennardjones/johnson.hpp

@@ -111,7 +111,7 @@ class LJ126KolafaNezbeda1994{
     auto get_a(const TTYPE& Tstar, const RHOTYPE& rhostar) const{
         std::common_type_t<TTYPE, RHOTYPE> summer = 0.0;
         for (auto [i, j, Cij] : c_Cij){
-            summer += Cij*pow(Tstar, i/2.0)*pow(rhostar, j);
+            summer += Cij*pow(Tstar, i/2.0)*powi(rhostar, j);
         }
         return forceeval(get_ahs(Tstar, rhostar) + exp(-gamma*POW2(rhostar))*rhostar*Tstar*get_DeltaB2hBH(Tstar)+summer);
     }

include/teqp/models/mie/mie.hpp

@@ -59,9 +59,9 @@ namespace Mie{
         auto alphar(const TTYPE& Tstar, const RHOTYPE& rhostar, const MoleFracType& /*molefrac*/) const {
             auto tau = forceeval(Tc / Tstar); auto delta = forceeval(rhostar / rhoc);
             using _t = std::decay_t<std::common_type_t<TTYPE, RHOTYPE>>;
-            _t s1 = 0; for (auto i = 0; i < n_pol.size(); ++i){ s1 += n_pol[i] * pow(tau, t_pol[i]) * powi(delta, d_pol[i]); }
-            _t s2 = 0; for (auto i = 0; i < n_exp.size(); ++i){ s2 += n_exp[i] * pow(tau, t_exp[i]) * powi(delta, d_exp[i]) * exp(-powi(delta, p[i])); }
-            _t s3 = 0; for (auto i = 0; i < n_gbs.size(); ++i){ s3 += n_gbs[i] * pow(tau, t_gbs[i]) * powi(delta, d_gbs[i]) * exp(-eta[i] * powi(forceeval(delta - eps[i]), 2) - beta[i] * powi(forceeval(tau - gam[i]), 2)); }
+            _t s1 = 0; for (auto i = 0; i < n_pol.size(); ++i){ s1 += n_pol[i] * pow(tau, t_pol[i]) * powi(delta, static_cast<int>(d_pol[i])); }
+            _t s2 = 0; for (auto i = 0; i < n_exp.size(); ++i){ s2 += n_exp[i] * pow(tau, t_exp[i]) * powi(delta, static_cast<int>(d_exp[i])) * exp(-powi(delta, static_cast<int>(p[i]))); }
+            _t s3 = 0; for (auto i = 0; i < n_gbs.size(); ++i){ s3 += n_gbs[i] * pow(tau, t_gbs[i]) * powi(delta, static_cast<int>(d_gbs[i])) * exp(-eta[i] * powi(forceeval(delta - eps[i]), 2) - beta[i] * powi(forceeval(tau - gam[i]), 2)); }
             return forceeval(s1 + s2 + s3);
         }
 

include/teqp/models/model_potentials/2center_ljf.hpp

@@ -258,7 +258,7 @@ namespace teqp {
 
                 result r = 0.0;
                 for (auto i = 0U; i < c.size(); ++i) {
-                    r = r + c[i] * pow(tau, m[i]) * pow(delta, n[i]) * pow(alpha, o[i]) * exp(p[i] * pow(delta, q[i]));
+                    r = r + c[i] * pow(tau, m[i]) * powi(delta, static_cast<int>(n[i])) * pow(alpha, o[i]) * exp(p[i] * powi(delta, static_cast<int>(q[i])));
                 }
                 return forceeval(r);
             }
@@ -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,8 +337,10 @@ 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 alphar_3 = Pole.alphar(tau, delta, mu_sq);
-                auto val = alphar_1 + alphar_2 + alphar_3;
+                auto val = forceeval(alphar_1 + alphar_2);
+                if (mu_sq != 0.0){
+                    val += Pole.alphar(tau, delta, mu_sq);
+                }
                 return forceeval(val);
             }
 
@@ -404,6 +406,25 @@ namespace teqp {
             return eos;
         }
 
+        // 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
+            auto DC_funcs = get_density_reducing(model_version);
+            auto TC_func = get_temperature_reducing(model_version);
+
+            //// Get contributions to EOS ( Attractive and regular part)
+            auto EOS_hard = get_HardSphere_contribution();
+            auto EOS_att = get_Attractive_contribution(model_version);
+            auto EOS_dipolar = get_Dipolar_contribution(model_version);
+            double mu_sq = 0.0;
+
+            // Build the 2-center Lennard-Jones model
+            auto model = Twocenterljf(std::move(DC_funcs), std::move(TC_func), std::move(EOS_hard), std::move(EOS_att), std::move(EOS_dipolar), L, mu_sq);
+
+            return model;
+        }
+
         // build the 2-center Lennard-Jones model with dipole
         inline auto build_two_center_model_dipole(const std::string& model_version, const double& L = 0.0, const double& mu_sq = 0.0) {
 

include/teqp/models/model_potentials/exp6.hpp

@@ -124,8 +124,9 @@ class Kataoka1992{
     {
         std::common_type_t<TType, RhoType> o = 0.0;
         for (auto el : c){
-            auto p = el[0], q = el[1], r = el[2], Apqr = el[3];
-            o = o + Apqr * pow(rhostar, p) * pow(Tstar, -q) * pow(alphastar, r);
+            int p = static_cast<int>(el[0]);
+            auto q = el[1], r = el[2], Apqr = el[3];
+            o = o + Apqr * powi(rhostar, p) * pow(Tstar, -q) * pow(alphastar, r);
         }
         return forceeval(o);
     }

include/teqp/models/model_potentials/squarewell.hpp

@@ -127,7 +127,7 @@ class EspindolaHeredia2009{
         for (auto n = 1; n < 5; ++n){
             num += thetai[n]*pow(lambda_, n);
         }
-        num *= pow(rhostar, 2);
+        num *= powi(rhostar, 2);
         RhoType den = 0;
         for (auto n = 5; n < 8; ++n){
             den += thetai[n]*pow(lambda_, n-4);
@@ -137,25 +137,25 @@ class EspindolaHeredia2009{
     }
 
     template<typename RhoType>
-    auto Chi(const RhoType & rhostar, double lambda_) const { return forceeval(a2i(2, lambda_)*rhostar*(1.0-pow(rhostar,2)/1.5129)); }
+    auto Chi(const RhoType & rhostar, double lambda_) const { return forceeval(a2i(2, lambda_)*rhostar*(1.0-powi(rhostar,2)/1.5129)); }
 
     template<typename RhoType>
     auto aHS(const RhoType & rhostar) const{
-        return forceeval(-3.0*m_pi*rhostar*(m_pi*rhostar-8.0)/pow(m_pi*rhostar-6.0, 2));
+        return forceeval(-3.0*m_pi*rhostar*(m_pi*rhostar-8.0)/powi(forceeval(m_pi*rhostar-6.0), 2));
     }
 
     template<typename RhoType>
     auto get_a1(const RhoType & rhostar, double lambda_) const{
-        RhoType o = a2i(1, lambda_)*pow(rhostar, 2-1) + a31(lambda_)*pow(rhostar, 3-1);
+        RhoType o = a2i(1, lambda_)*powi(rhostar, 2-1) + a31(lambda_)*powi(rhostar, 3-1);
         for (auto i = 1; i < 5; ++i){
-            o = o + gamman(i, lambda_)*pow(rhostar, i+2);
+            o = o + gamman(i, lambda_)*powi(rhostar, i+2);
         }
         return forceeval(o);
     }
 
     template<typename RhoType>
     auto get_a2(const RhoType & rhostar, double lambda_) const{
-        return forceeval(Chi(rhostar, lambda_)*exp(xi2(lambda_)*rhostar + phii(1, lambda_)*pow(rhostar,3) + phii(2,lambda_)*pow(rhostar,4)));
+        return forceeval(Chi(rhostar, lambda_)*exp(xi2(lambda_)*rhostar + phii(1, lambda_)*powi(rhostar,3) + phii(2,lambda_)*powi(rhostar,4)));
     }
 
     template<typename RhoType>

include/teqp/models/pcsaft.hpp

@@ -123,18 +123,23 @@ auto get_alphar_hs(const VecType& zeta, const VecType2& D) {
      for Nderiv in [1, 2, 3, 4, 5]:
          display(simplify((simplify(diff(alpha, rho, Nderiv)).subs(rho,0)*rho**Nderiv/factorial(Nderiv)).subs(D_0, zeta_0/rho).subs(D_1, zeta_1/rho).subs(D_2, zeta_2/rho).subs(D_3, zeta_3/rho)))
     </sympy>
+     
+     The updated approach is simpler; start off with the expression for alphar_hs, and simplify
+     ratios where 0/0 for which the l'hopital limit would be ok until you remove all the terms
+     in the denominator, allowing it to be evaluated. All terms have a zeta_x/zeta_0 and a few
+     have zeta_3*zeta_0 in the denominator
+     
     */
+    auto Upsilon = 1.0 - zeta[3];
     if (getbaseval(zeta[3]) == 0){
         return forceeval(
-             0.0 // 0-th order term, the limit of the function at zero density is zero
-             + zeta[3] + 3.0*D[1]*zeta[2]/D[0] // 1st order term f'(x=0)*x/1!
-             + (zeta[3]*zeta[3] + 6.0*D[1]/D[0]*zeta[2]*zeta[3] + 3.0*zeta[2]*zeta[2]*D[2]/D[0])/2.0 // 2nd order term f''(x=0)*x^2/2!
-             + D[3]/D[0]*(zeta[0]*zeta[3]*zeta[3] + 9.0*zeta[1]*zeta[2]*zeta[3] + 8.0*zeta[2]*zeta[2]*zeta[2])/3.0 // 3rd order term f'''(x=0)*x^3/3!
-             + zeta[3]*D[3]/D[0]*(zeta[0]*zeta[3]*zeta[3] + 12.0*zeta[1]*zeta[2]*zeta[3] + 15.0*zeta[2]*zeta[2]*zeta[2])/4.0 // 4th order term f''''(x=0)*x^4/4!
-             // ... and so on
-         );
+             3.0*D[1]/D[0]*zeta[2]/Upsilon
+             + D[2]*D[2]*zeta[2]/(D[3]*D[0]*Upsilon*Upsilon)
+             - log(Upsilon)
+             + (D[2]*D[2]*D[2])/(D[3]*D[3]*D[0])*log(Upsilon)
+        );
     }
-    auto Upsilon = 1.0 - zeta[3];
+
     return forceeval(1.0 / zeta[0] * (3.0 * zeta[1] * zeta[2] / Upsilon
         + zeta[2] * zeta[2] * zeta[2] / zeta[3] / Upsilon / Upsilon
         + (zeta[2] * zeta[2] * zeta[2] / (zeta[3] * zeta[3]) - zeta[0]) * log(1.0 - zeta[3])

include/teqp/models/saft/genericsaft.hpp

@@ -4,13 +4,15 @@
 #include "teqp/models/saftvrmie.hpp"
 #include "teqp/models/association/association.hpp"
 #include "teqp/models/saft/softsaft.hpp"
+#include "teqp/models/model_potentials/2center_ljf.hpp"
 
 namespace teqp::saft::genericsaft{
 
 struct GenericSAFT{
 
 public:
-    using NonPolarTerms = std::variant<saft::pcsaft::PCSAFTMixture, SAFTVRMie::SAFTVRMieNonpolarMixture, saft::softsaft::SoftSAFT>;
+    using TwoCLJ = twocenterljf::Twocenterljf<twocenterljf::DipolarContribution>;
+    using NonPolarTerms = std::variant<saft::pcsaft::PCSAFTMixture, SAFTVRMie::SAFTVRMieNonpolarMixture, saft::softsaft::SoftSAFT, TwoCLJ>;
 //    using PolarTerms = EOSTermContainer<>;
     using AssociationTerms = std::variant<association::Association>;
 
@@ -26,7 +28,9 @@ struct GenericSAFT{
         else if (kind == "Johnson+Johnson" || kind == "softSAFT"){
             return saft::softsaft::SoftSAFT(j.at("model"));
         }
-        // TODO: 2CLJ
+        else if (kind == "2CLJF"){
+            return twocenterljf::build_two_center_model(j.at("author"), j.at("L^*"));
+        }
         else{
             throw std::invalid_argument("Not valid nonpolar kind:" + kind);
         }

include/teqp/models/saft/pcsaftpure.hpp

@@ -41,7 +41,8 @@ class PCSAFTPureGrossSadowski2001{
 
         auto d = forceeval(sigma_A*(1.0-0.12*exp(-3.0*eps_k/T)));
         Eigen::Array<decltype(d), 4, 1> dpowers; dpowers(0) = 1.0; for (auto i = 1U; i <= 3; ++i){ dpowers(i) = d*dpowers(i-1); }
-        auto zeta = pi/6.0*rhoN_A3*m*dpowers;
+        auto D = pi/6.0*m*dpowers;
+        auto zeta = rhoN_A3*D.template cast<std::common_type_t<TTYPE, RhoType>>();
 
         auto zeta2_to2 = zeta[2]*zeta[2];
         auto zeta2_to3 = zeta2_to2*zeta[2];
@@ -51,9 +52,18 @@ class PCSAFTPureGrossSadowski2001{
         auto onemineta_to3 = onemineta*onemineta_to2;
         auto onemineta_to4 = onemineta*onemineta_to3;
 
-        auto alpha_hs = (3.0*zeta[1]*zeta[2]/onemineta
+        std::decay_t<decltype(zeta[0])> alpha_hs = forceeval((3.0*zeta[1]*zeta[2]/onemineta
                          + zeta2_to3/(zeta[3]*onemineta_to2)
-                         + (zeta2_to3/zeta3_to2-zeta[0])*log(1.0-zeta[3]))/zeta[0];
+                         + (zeta2_to3/zeta3_to2-zeta[0])*log(1.0-zeta[3]))/zeta[0]);
+        if (getbaseval(zeta[0]) == 0){
+            auto Upsilon = 1.0-zeta[3];
+            alpha_hs = forceeval(
+                3.0*D[1]/D[0]*zeta[2]/Upsilon
+                + D[2]*D[2]*zeta[2]/(D[3]*D[0]*Upsilon*Upsilon)
+                - log(Upsilon)
+                + (D[2]*D[2]*D[2])/(D[3]*D[3]*D[0])*log(Upsilon)
+            );
+        }
 
         auto fac_g_hs = d/2.0; // d*d/(2*d)
         auto gii = (1.0/onemineta

include/teqp/models/saftvrmie.hpp

@@ -681,15 +681,20 @@ struct SAFTVRMieChainContributionTerms{
              for Nderiv in [1, 2, 3, 4, 5]:
                  display(simplify((simplify(diff(alpha, rho, Nderiv)).subs(rho,0)*rho**Nderiv/factorial(Nderiv)).subs(D_0, zeta_0/rho).subs(D_1, zeta_1/rho).subs(D_2, zeta_2/rho).subs(D_3, zeta_3/rho)))
             </sympy>
+             
+             The updated approach is simpler; start off with the expression for alphar_hs, and simplify
+             ratios where 0/0 for which the l'hopital limit would be ok until you remove all the terms
+             in the denominator, allowing it to be evaluated. All terms have a zeta_x/zeta_0 and a few
+             have zeta_3*zeta_0 in the denominator
+             
             */
+            auto Upsilon = 1.0 - zeta[3];
             return forceeval(
-                 0.0 // 0-th order term, the limit of the function at zero density is zero
-                 + zeta[3] + 3.0*D[1]*zeta[2]/D[0] // 1st order term f'(x=0)*x/1!
-                 + (zeta[3]*zeta[3] + 6.0*D[1]/D[0]*zeta[2]*zeta[3] + 3.0*zeta[2]*zeta[2]*D[2]/D[0])/2.0 // 2nd order term f''(x=0)*x^2/2!
-                 + D[3]/D[0]*(zeta[0]*zeta[3]*zeta[3] + 9.0*zeta[1]*zeta[2]*zeta[3] + 8.0*zeta[2]*zeta[2]*zeta[2])/3.0 // 3rd order term f'''(x=0)*x^3/3!
-                 + zeta[3]*D[3]/D[0]*(zeta[0]*zeta[3]*zeta[3] + 12.0*zeta[1]*zeta[2]*zeta[3] + 15.0*zeta[2]*zeta[2]*zeta[2])/4.0 // 4th order term f''''(x=0)*x^4/4!
-                 // ... and so on
-             );
+                 3.0*D[1]/D[0]*zeta[2]/Upsilon
+                 + D[2]*D[2]*zeta[2]/(D[3]*D[0]*Upsilon*Upsilon)
+                 - log(Upsilon)
+                 + (D[2]*D[2]*D[2])/(D[3]*D[3]*D[0])*log(Upsilon)
+            );
         }
         else{
             return forceeval(6.0/(MY_PI*rhos)*(3.0*zeta[1]*zeta[2]/(1.0-zeta[3]) + POW3(zeta[2])/(zeta[3]*POW2(1.0-zeta[3])) + (POW3(zeta[2])/POW2(zeta[3])-zeta[0])*log(1.0-zeta[3])));

include/teqp/types.hpp

@@ -179,6 +179,11 @@ namespace teqp {
         //return e.cast<T>().unaryExpr([&x](const auto& e_) {return powi(x, e_); }).eval();
     }
 
+    template<typename T>
+    inline auto powIVd(const T& x, const Eigen::ArrayXd& e) {
+        return e.cast<T>().unaryExpr([&x](const auto& e_) {return forceeval(pow(x, e_)); }).eval();
+    }
+
     //template<typename T>
     //auto powIV(const T& x, const Eigen::ArrayXd& e) {
     //    Eigen::Array<T, Eigen::Dynamic, 1> o = e.cast<T>();

interface/CPP/model_factory_model_potentials.cpp

@@ -22,6 +22,9 @@ namespace teqp{
         std::unique_ptr<teqp::cppinterface::AbstractModel> make_Mie_Chaparro2023(const nlohmann::json &spec){
             return make_owned(FEANN::ChaparroJCP2023(spec.at("lambda_r"), spec.at("lambda_a")));
         }
+        std::unique_ptr<teqp::cppinterface::AbstractModel> make_2CLJF(const nlohmann::json &spec){
+            return make_owned(twocenterljf::build_two_center_model(spec.at("author"), spec.at("L^*")));
+        }
         std::unique_ptr<teqp::cppinterface::AbstractModel> make_2CLJF_Dipole(const nlohmann::json &spec){
             return make_owned(twocenterljf::build_two_center_model_dipole(spec.at("author"), spec.at("L^*"), spec.at("(mu^*)^2")));
         }