Newer
Older
/***
\brief This file contains the contributions that can be composed together to form SAFT models
*/
#pragma once
#include "nlohmann/json.hpp"
#include "teqp/types.hpp"
#include "teqp/exceptions.hpp"
#include "teqp/constants.hpp"
#include "teqp/math/quadrature.hpp"
#include "teqp/models/saft/polar_terms.hpp"
namespace teqp {
namespace SAFTVRMie {
/// Coefficients for one fluid
struct SAFTVRMieCoeffs {
std::string name; ///< Name of fluid
double m = -1, ///< number of segments
epsilon_over_k = -1, ///< [K] depth of pair potential divided by Boltzman constant
lambda_a = -1, ///< The attractive exponent (the 6 in LJ 12-6 potential)
lambda_r = -1, ///< The repulsive exponent (the 12 in LJ 12-6 potential)
mustar2 = 0, ///< nondimensional, the reduced dipole moment squared
nmu = 0, ///< number of dipolar segments
Qstar2 = 0, ///< nondimensional, the reduced quadrupole squared
nQ = 0; ///< number of quadrupolar segments
std::string BibTeXKey; ///< The BibTeXKey for the reference for these coefficients
};
/// Manager class for SAFT-VR-Mie coefficients
class SAFTVRMieLibrary {
std::map<std::string, SAFTVRMieCoeffs> coeffs;
public:
SAFTVRMieLibrary() {
insert_normal_fluid("Methane", 1.0000, 3.7412e-10, 153.36, 12.650, 6, "Lafitte-JCP-2001");
insert_normal_fluid("Ethane", 1.4373, 3.7257e-10, 206.12, 12.400, 6, "Lafitte-JCP-2001");
insert_normal_fluid("Propane", 1.6845, 3.9056e-10, 239.89, 13.006, 6, "Lafitte-JCP-2001");
void insert_normal_fluid(const std::string& name, double m, const double sigma_m, const double epsilon_over_k, const double lambda_r, const double lambda_a, const std::string& BibTeXKey) {
SAFTVRMieCoeffs coeff;
coeff.name = name;
coeff.m = m;
coeff.epsilon_over_k = epsilon_over_k;
coeff.lambda_r = lambda_r;
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
coeff.BibTeXKey = BibTeXKey;
coeffs.insert(std::pair<std::string, SAFTVRMieCoeffs>(name, coeff));
}
const auto& get_normal_fluid(const std::string& name) {
auto it = coeffs.find(name);
if (it != coeffs.end()) {
return it->second;
}
else {
throw std::invalid_argument("Bad name:" + name);
}
}
auto get_coeffs(const std::vector<std::string>& names){
std::vector<SAFTVRMieCoeffs> c;
for (auto n : names){
c.push_back(get_normal_fluid(n));
}
return c;
}
};
/// Things that only depend on the components themselves, but not on composition, temperature, or density
struct SAFTVRMieChainContributionTerms{
private:
/// The matrix of coefficients needed to evaluate f_k
const Eigen::Matrix<double, 7, 7> phi{(Eigen::Matrix<double, 7, 7>() <<
7.5365557, -359.44, 1550.9, -1.19932, -1911.28, 9236.9, 10,
-37.60463, 1825.6, -5070.1, 9.063632, 21390.175, -129430, 10,
71.745953, -3168.0, 6534.6, -17.9482, -51320.7, 357230, 0.57,
-46.83552, 1884.2, -3288.7, 11.34027, 37064.54, -315530, -6.7,
-2.467982, -0.82376, -2.7171, 20.52142, 1103.742, 1390.2, -8,
-0.50272, -3.1935, 2.0883, -56.6377, -3264.61, -4518.2, 0,
8.0956883, 3.7090, 0, 40.53683, 2556.181, 4241.6, 0 ).finished()};
/// The matrix used to obtain the parameters c_1, c_2, c_3, and c_4 in Eq. A18
const Eigen::Matrix<double, 4, 4> A{(Eigen::Matrix<double, 4, 4>() <<
0.81096, 1.7888, -37.578, 92.284,
1.0205, -19.341, 151.26, -463.50,
-1.9057, 22.845, -228.14, 973.92,
1.0885, -6.1962, 106.98, -677.64).finished()};
// Eq. A48
auto get_lambda_k_ij(const Eigen::ArrayXd& lambda_k) const{
Eigen::ArrayXXd mat(N,N);
for (auto i = 0; i < lambda_k.size(); ++i){
for (auto j = i; j < lambda_k.size(); ++j){
mat(i,j) = 3 + sqrt((lambda_k(i)-3)*(lambda_k(j)-3));
mat(j,i) = mat(i,j);
}
}
return mat;
}
/// Eq. A3
Eigen::ArrayXXd C(N,N);
for (auto i = 0; i < N; ++i){
for (auto j = i; j < N; ++j){
C(i,j) = lambda_r_ij(i,j)/(lambda_r_ij(i,j)-lambda_a_ij(i,j))*pow(lambda_r_ij(i,j)/lambda_a_ij(i,j), lambda_a_ij(i,j)/(lambda_r_ij(i,j)-lambda_a_ij(i,j)));
C(j,i) = C(i,j); // symmetric
}
}
return C;
}
// Eq. A26
std::vector<Eigen::ArrayXXd> f_(8); // 0-th element is present, but not initialized
for (auto k = 1; k < 8; ++k){
f_[k].resize(N,N);
};
for (auto k = 1; k < 8; ++k){
auto phik = phi.col(k-1); // phi matrix is indexed to start at 1, but our matrix starts at 0
Eigen::ArrayXXd num(N,N), den(N,N); num.setZero(), den.setZero();
for (auto n = 0; n < 4; ++n){
num += phik[n]*pow(alpha_ij, n);
}
for (auto n = 4; n < 7; ++n){
den += phik[n]*pow(alpha_ij, n-3);
}
f_[k] = num/(1 + den);
}
return f_;
}
/// Eq. A45
Eigen::ArrayXXd sigma(N,N);
for (auto i = 0; i < N; ++i){
for (auto j = i; j < N; ++j){
sigma(j,i) = sigma(i,j); // symmetric
}
}
return sigma;
}
Eigen::ArrayXXd eps_(N,N);
for (auto i = 0; i < N; ++i){
for (auto j = i; j < N; ++j){
eps_(i,j) = (1.0-kmat(i,j))*sqrt(pow(sigma_ij(i,i),3)*pow(sigma_ij(j,j),3)*epsilon_over_k(i)*epsilon_over_k(j))/pow(sigma_ij(i,j), 3);
eps_(j,i) = eps_(i,j); // symmetric
}
}
return eps_;
}
std::size_t get_N(){
auto sizes = (Eigen::ArrayXd(5) << m.size(), epsilon_over_k.size(), sigma_A.size(), lambda_a.size(), lambda_r.size()).finished();
if (sizes.mean() != sizes.minCoeff()){
throw teqp::InvalidArgument("sizes of pure component arrays are not all the same");
}
return sizes[0];
}
/// Eq. A18 for the attractive exponents
auto get_cij(const Eigen::ArrayXXd& lambdaij) const{
std::vector<Eigen::ArrayXXd> cij(4);
for (auto n = 0; n < 4; ++n){
cij[n].resize(N,N);
};
for (auto i = 0; i < N; ++i){
for (auto j = i; j < N; ++j){
using CV = Eigen::Vector<double, 4>;
const CV b{(CV() << 1, 1.0/lambdaij(i,j), 1.0/pow(lambdaij(i,j),2), 1.0/pow(lambdaij(i,j),3)).finished()};
auto c1234 = (A*b).eval();
cij[0](i,j) = c1234(0);
cij[1](i,j) = c1234(1);
cij[2](i,j) = c1234(2);
cij[3](i,j) = c1234(3);
}
}
return cij;
}
/// Eq. A18 for the attractive exponents
return get_cij(lambda_a_ij);
}
/// Eq. A18 for 2x the attractive exponents
return get_cij(2.0*lambda_a_ij);
}
/// Eq. A18 for the repulsive exponents
return get_cij(lambda_r_ij);
}
/// Eq. A18 for the 2x the repulsive exponents
return get_cij(2.0*lambda_r_ij);
}
/// Eq. A18 for the 2x the repulsive exponents
return get_cij(lambda_r_ij + lambda_a_ij);
}
public:
// One entry per component
const Eigen::ArrayXd m, epsilon_over_k, sigma_A, lambda_a, lambda_r;
const Eigen::ArrayXXd kmat;
const std::size_t N;
// Calculated matrices for the ij pair
const Eigen::ArrayXXd lambda_r_ij, lambda_a_ij, C_ij, alpha_ij, sigma_ij, epsilon_ij; // Matrices of parameters
const std::vector<Eigen::ArrayXXd> crnij, canij, c2rnij, c2anij, carnij;
const std::vector<Eigen::ArrayXXd> fkij; // Matrices of parameters
SAFTVRMieChainContributionTerms(
const Eigen::ArrayXd& m,
const Eigen::ArrayXd& epsilon_over_k,
const Eigen::ArrayXd& sigma_m,
const Eigen::ArrayXd& lambda_r,
const Eigen::ArrayXd& lambda_a,
: m(m), epsilon_over_k(epsilon_over_k), sigma_A(sigma_m*1e10), lambda_a(lambda_a), lambda_r(lambda_r), kmat(kmat),
N(get_N()),
lambda_r_ij(get_lambda_k_ij(lambda_r)), lambda_a_ij(get_lambda_k_ij(lambda_a)),
C_ij(get_C_ij()), alpha_ij(C_ij*(1/(lambda_a_ij-3) - 1/(lambda_r_ij-3))),
sigma_ij(get_sigma_ij()), epsilon_ij(get_epsilon_ij()),
crnij(get_crnij()), canij(get_canij()),
c2rnij(get_c2rnij()), c2anij(get_c2anij()), carnij(get_carnij()),
fkij(get_fkij())
{}
template<typename RType>
auto get_uii_over_kB(std::size_t i, const RType& r) const {
auto rstarinv = forceeval(sigma_A[i]/r);
return forceeval(C_ij(i,i)*epsilon_over_k[i]*(pow(rstarinv, lambda_r[i]) - pow(rstarinv, lambda_a[i])));
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
/// Solve for the value of \f$j=\sigma/r\f$ for which the integrand in \f$d_{ii}\f$ becomes equal to 1 to numerical precision
template <typename TType>
auto get_j_cutoff_dii(std::size_t i, const TType &T) const {
auto lambda_a_ = lambda_a(i), lambda_r_ = lambda_r(i);
auto EPS = std::numeric_limits<decltype(getbaseval(T))>::epsilon();
auto K = forceeval(log(-log(EPS)*T/(C_ij(i,i)*epsilon_ij(i,i))));
auto j0 = forceeval(exp(K/lambda_r_)); // this was proposed by longemen3000 (Andrés Riedemann)
auto kappa = C_ij(i,i)*epsilon_ij(i,i);
// Function to return residual and its derivatives w.r.t.
auto fgh = [&kappa, &lambda_r_, &lambda_a_, &T, &EPS](auto j){
auto jlr = pow(j, lambda_r_), jla = pow(j, lambda_a_);
auto u = kappa*(jlr - jla);
auto uprime = kappa*(lambda_r_*jlr - lambda_a_*jla)/j;
auto uprime2 = kappa*(lambda_r_*(lambda_r_-1.0)*jlr - lambda_a_*(lambda_a_-1.0)*jla)/(j*j);
return std::make_tuple(forceeval(-u/T-log(EPS)), forceeval(-uprime/T), forceeval(-uprime2/T));
};
TType j = j0;
for (auto counter = 0; counter <= 3; ++counter){
// Halley's method steps
auto [R, Rprime, Rprime2] = fgh(j);
auto denominator = 2.0*Rprime*Rprime-R*Rprime2;
if (getbaseval(denominator) < EPS){
break;
}
j -= 2.0*R*Rprime/denominator;
}
double jbase = getbaseval(j);
if (jbase < 1.0){
throw teqp::IterationFailure("Cannot obtain a value of j");
}
return j;
}
/**
\note Eq. A9 from Lafitte
The calculation of the diameter is based upon
\f[
d_{ii} = \int_0^{\sigma_{ii}}(1-\exp(-\beta u_{ii}^{\rm Mie}(r)){\rm d}r
\f]
which is broken up into two parts:
\f[
d = \int_0^{r_{\rm cut}} 1 {\rm d} r + \int_{r_{\rm cut}}^{\sigma_{ii}} [1-\exp(-\beta u_{ii}^{\rm Mie}(r))] {\rm d}r
\f]
but the integrand is basically constant (to numerical precision) from 0 to some cutoff value of \f$r\f$, which we'll call \f$r_{\rm cut}\f$. So first we need to find the value of \f$r_{\rm cut}\f$ that makes the integrand take its constant value, which is explained well in the paper from Aasen (https://github.com/ClapeyronThermo/Clapeyron.jl/issues/152#issuecomment-1480324192). Finding the cutoff value is obtained when
\f[
\exp(-\beta u_{ii}^{\rm Mie}(r)) = EPS
\f]
where EPS is the numerical precision of the floating point type. Taking the logs of both sides,
\f[
-\beta u_{ii}^{\rm Mie} = \ln(EPS)
\f]
To get a starting value, it is first assumed that only the repulsive contribution contributes to the potential, yielding \f$u^{\rm rep} = C\epsilon(\sigma/r)^{\lambda_r}\f$ (with \f$C\f$ the same as the full potential with attraction) which yields
\f[
-\beta C\epsilon(\sigma/r)^{\lambda_r} = \ln(EPS)
\f]
and
\f[
(\sigma/r)_{\rm guess} = (-\ln(EPS)/(\beta C \epsilon))^{1/\lambda_r}
\f]
Then we solve for the residual \f$R(r)=0\f$, where \f$R_0=\exp(-u/T)-EPS\f$. Equivalently we can write the residual in logarithmic terms as \f$R=-u/T-\ln(EPS)\f$. This simplifies the rootfinding as you need \f$R\f$, \f$R'\f$ and \f$R''\f$ to apply Halley's method, which are themselves quite straightforward to obtain because \f$R'=-u'/T\f$, \f$R''=-u''/T\f$, where the primes are derivatives taken with respect to \f$\sigma/r\f$.
*/
TType get_dii(std::size_t i, const TType &T) const{
std::function<TType(TType)> integrand = [this, i, &T](const TType& r){
return forceeval(1.0-exp(-this->get_uii_over_kB(i, r)/T));
// Sum of the two integrals, one is constant, the other is from integration
auto rcut = forceeval(sigma_A[i]/get_j_cutoff_dii(i, T));
auto integral_contribution = quad<10, TType, TType>(integrand, rcut, sigma_A[i]);
auto d = forceeval(rcut + integral_contribution);
if (getbaseval(d) > sigma_A[i]){
throw teqp::IterationFailure("Value of d is larger than sigma; this is impossible");
}
return d;
}
template <typename TType>
auto get_dmat(const TType &T) const{
Eigen::Array<TType, Eigen::Dynamic, Eigen::Dynamic> d(N,N);
// For the pure components, by integration
for (auto i = 0; i < N; ++i){
d(i,i) = get_dii(i, T);
}
// The cross terms, using the linear mixing rule
for (auto i = 0; i < N; ++i){
for (auto j = i+1; j < N; ++j){
d(i,j) = (d(i,i) + d(j,j))/2.0;
d(j,i) = d(i,j);
}
}
return d;
}
// Calculate core parameters that depend on temperature, volume, and composition
template <typename TType, typename RhoType, typename VecType>
auto get_core_calcs(const TType& T, const RhoType& rhomolar, const VecType& molefracs) const{
if (molefracs.size() != N){
throw teqp::InvalidArgument("Length of molefracs of "+std::to_string(molefracs.size()) + " does not match the model size of"+std::to_string(N));
}
using FracType = std::decay_t<decltype(molefracs[0])>;
using NumType = std::common_type_t<TType, RhoType, FracType>;
// Things that are easy to calculate
// ....
auto dmat = get_dmat(T); // Matrix of diameters of pure and cross terms
auto rhoN = forceeval(rhomolar*N_A); // Number density, in molecules/m^3
auto mbar = forceeval((molefracs*m).sum()); // Mean number of segments, dimensionless
auto rhos = forceeval(rhoN*mbar/1e30); // Mean segment number density, in segments/A^3
auto xs = forceeval((m*molefracs/mbar).eval()); // Segment fractions
constexpr double MY_PI = static_cast<double>(EIGEN_PI);
auto pi6 = MY_PI/6;
using TRHOType = std::common_type_t<std::decay_t<TType>, std::decay_t<RhoType>, std::decay_t<decltype(molefracs[0])>, std::decay_t<decltype(m[0])>>;
Eigen::Array<TRHOType, 4, 1> zeta;
for (auto l = 0; l < 4; ++l){
NumType summer_zeta_x = 0.0;
TRHOType summer_zeta_x_bar = 0.0;
for (auto i = 0; i < N; ++i){
for (auto j = 0; j < N; ++j){
summer_zeta_x += xs(i)*xs(j)*powi(dmat(i,j), 3)*rhos;
summer_zeta_x_bar += xs(i)*xs(j)*powi(sigma_ij(i,j), 3);
auto zeta_x = forceeval(pi6*summer_zeta_x); // Eq. A13
auto zeta_x_bar = forceeval(pi6*rhos*summer_zeta_x_bar); // Eq. A23
auto zeta_x_bar5 = forceeval(POW2(zeta_x_bar)*POW3(zeta_x_bar)); // (zeta_x_bar)^5
auto zeta_x_bar8 = forceeval(zeta_x_bar5*POW3(zeta_x_bar)); // (zeta_x_bar)^8
// Coefficients in the gdHSij term, do not depend on component,
// so calculate them here
auto X = forceeval(POW3(1.0 - zeta_x)), X3 = X;
auto X2 = forceeval(POW2(1.0 - zeta_x));
auto k0 = forceeval(-log(1.0-zeta_x) + (42.0*zeta_x - 39.0*POW2(zeta_x) + 9.0*POW3(zeta_x) - 2.0*POW4(zeta_x))/(6.0*X3)); // Eq. A30
auto k1 = forceeval((POW4(zeta_x) + 6.0*POW2(zeta_x) - 12.0*zeta_x)/(2.0*X3));
auto k2 = forceeval(-3.0*POW2(zeta_x)/(8.0*X2));
auto k3 = forceeval((-POW4(zeta_x) + 3.0*POW2(zeta_x) + 3.0*zeta_x)/(6.0*X3));
// Pre-calculate the cubes of the diameters
auto dmat3 = dmat.array().pow(3).eval();
NumType a1kB = 0.0;
NumType a2kB2 = 0.0;
NumType a3kB3 = 0.0;
NumType K_HS = get_KHS(zeta_x);
NumType rho_dK_HS_drho = get_rhos_dK_HS_drhos(zeta_x);
for (auto i = 0; i < N; ++i){
for (auto j = i; j < N; ++j){
NumType x_0_ij = sigma_ij(i,j)/dmat(i, j);
// -----------------------
// Calculations for a_1/kB
// -----------------------
auto I = [&x_0_ij](double lambda_ij){
return forceeval(-(pow(x_0_ij, 3-lambda_ij)-1.0)/(lambda_ij-3.0)); // Eq. A14
};
auto J = [&x_0_ij](double lambda_ij){
return forceeval(-(pow(x_0_ij, 4-lambda_ij)*(lambda_ij-3.0)-pow(x_0_ij, 3.0-lambda_ij)*(lambda_ij-4.0)-1.0)/((lambda_ij-3.0)*(lambda_ij-4.0))); // Eq. A15
auto Bhatij_a = this->get_Bhatij(zeta_x, X, I(lambda_a_ij(i,j)), J(lambda_a_ij(i,j)));
auto Bhatij_2a = this->get_Bhatij(zeta_x, X, I(2*lambda_a_ij(i,j)), J(2*lambda_a_ij(i,j)));
auto Bhatij_r = this->get_Bhatij(zeta_x, X, I(lambda_r_ij(i,j)), J(lambda_r_ij(i,j)));
auto Bhatij_2r = this->get_Bhatij(zeta_x, X, I(2*lambda_r_ij(i,j)), J(2*lambda_r_ij(i,j)));
auto Bhatij_ar = this->get_Bhatij(zeta_x, X, I(lambda_a_ij(i,j)+lambda_r_ij(i,j)), J(lambda_a_ij(i,j)+lambda_r_ij(i,j)));
auto one_term = [this, &x_0_ij, &I, &J, &zeta_x, &X](double lambda_ij, const NumType& zeta_x_eff){
return forceeval(
pow(x_0_ij, lambda_ij)*(
this->get_Bhatij(zeta_x, X, I(lambda_ij), J(lambda_ij))
+ this->get_a1Shatij(zeta_x_eff, lambda_ij)
NumType zeta_x_eff_r = crnij[0](i,j)*zeta_x + crnij[1](i,j)*POW2(zeta_x) + crnij[2](i,j)*POW3(zeta_x) + crnij[3](i,j)*POW4(zeta_x);
NumType zeta_x_eff_a = canij[0](i,j)*zeta_x + canij[1](i,j)*POW2(zeta_x) + canij[2](i,j)*POW3(zeta_x) + canij[3](i,j)*POW4(zeta_x);
NumType dzeta_x_eff_dzetax_r = crnij[0](i,j) + crnij[1](i,j)*2*zeta_x + crnij[2](i,j)*3*POW2(zeta_x) + crnij[3](i,j)*4*POW3(zeta_x);
NumType dzeta_x_eff_dzetax_a = canij[0](i,j) + canij[1](i,j)*2*zeta_x + canij[2](i,j)*3*POW2(zeta_x) + canij[3](i,j)*4*POW3(zeta_x);
NumType a1ij = 2.0*MY_PI*rhos*dmat3(i,j)*epsilon_ij(i,j)*C_ij(i,j)*(
one_term(lambda_a_ij(i,j), zeta_x_eff_a) - one_term(lambda_r_ij(i,j), zeta_x_eff_r)
NumType contribution = xs(i)*xs(j)*a1ij;
double factor = (i == j) ? 1.0 : 2.0; // Off-diagonal terms contribute twice
a1kB += contribution*factor;
// --------------------------
// Calculations for a_2/k_B^2
// --------------------------
NumType zeta_x_eff_2r = c2rnij[0](i,j)*zeta_x + c2rnij[1](i,j)*POW2(zeta_x) + c2rnij[2](i,j)*POW3(zeta_x) + c2rnij[3](i,j)*POW4(zeta_x);
NumType zeta_x_eff_2a = c2anij[0](i,j)*zeta_x + c2anij[1](i,j)*POW2(zeta_x) + c2anij[2](i,j)*POW3(zeta_x) + c2anij[3](i,j)*POW4(zeta_x);
NumType zeta_x_eff_ar = carnij[0](i,j)*zeta_x + carnij[1](i,j)*POW2(zeta_x) + carnij[2](i,j)*POW3(zeta_x) + carnij[3](i,j)*POW4(zeta_x);
NumType dzeta_x_eff_dzetax_2r = c2rnij[0](i,j) + c2rnij[1](i,j)*2*zeta_x + c2rnij[2](i,j)*3*POW2(zeta_x) + c2rnij[3](i,j)*4*POW3(zeta_x);
NumType dzeta_x_eff_dzetax_ar = carnij[0](i,j) + carnij[1](i,j)*2*zeta_x + carnij[2](i,j)*3*POW2(zeta_x) + carnij[3](i,j)*4*POW3(zeta_x);
NumType dzeta_x_eff_dzetax_2a = c2anij[0](i,j) + c2anij[1](i,j)*2*zeta_x + c2anij[2](i,j)*3*POW2(zeta_x) + c2anij[3](i,j)*4*POW3(zeta_x);
NumType chi_ij = fkij[1](i,j)*zeta_x_bar + fkij[2](i,j)*zeta_x_bar5 + fkij[3](i,j)*zeta_x_bar8;
auto a2ij = 0.5*K_HS*(1.0+chi_ij)*epsilon_ij(i,j)*POW2(C_ij(i,j))*(2*MY_PI*rhos*dmat3(i,j)*epsilon_ij(i,j))*(
one_term(2.0*lambda_a_ij(i,j), zeta_x_eff_2a)
-2.0*one_term(lambda_a_ij(i,j)+lambda_r_ij(i,j), zeta_x_eff_ar)
+one_term(2.0*lambda_r_ij(i,j), zeta_x_eff_2r)
); // divided by k_B^2
NumType contributiona2 = xs(i)*xs(j)*a2ij; // Eq. A19
a2kB2 += contributiona2*factor;
// --------------------------
// Calculations for a_3/k_B^3
// --------------------------
auto a3ij = -POW3(epsilon_ij(i,j))*fkij[4](i,j)*zeta_x_bar*exp(
fkij[5](i,j)*zeta_x_bar + fkij[6](i,j)*POW2(zeta_x_bar)
); // divided by k_B^3
NumType contributiona3 = xs(i)*xs(j)*a3ij; // Eq. A25
a3kB3 += contributiona3*factor;
if (i == j){
// ------------------
// Chain contribution
// ------------------
// Eq. A29
auto gdHSii = exp(k0 + k1*x_0_ij + k2*POW2(x_0_ij) + k3*POW3(x_0_ij));
// The g1 terms
// ....
// This is the function for the second part (not the partial) that goes in g_{1,ii},
// divided by 2*PI*d_ij^3*epsilon*rhos
auto g1_term = [&one_term](double lambda_ij, const NumType& zeta_x_eff){
return forceeval(lambda_ij*one_term(lambda_ij, zeta_x_eff));
auto g1_noderivterm = -C_ij(i,i)*(g1_term(lambda_a_ij(i,i), zeta_x_eff_a)-g1_term(lambda_r_ij(i,i), zeta_x_eff_r));
// Bhat = B*rho*kappa; diff(Bhat, rho) = Bhat + rho*dBhat/drho; kappa = 2*pi*eps*d^3
// This is the function for the partial derivative rhos*(da1ij/drhos),
// divided by 2*PI*d_ij^3*epsilon*rhos
auto rhosda1iidrhos_term = [this, &x_0_ij, &I, &J, &zeta_x, &X](double lambda_ij, const NumType& zeta_x_eff, const NumType& dzetaxeff_dzetax, const NumType& Bhatij){
auto I_ = I(lambda_ij);
auto J_ = J(lambda_ij);
auto rhosda1Sdrhos = this->get_rhoda1Shatijdrho(zeta_x, zeta_x_eff, dzetaxeff_dzetax, lambda_ij);
auto rhosdBdrhos = this->get_rhodBijdrho(zeta_x, X, I_, J_, Bhatij);
return forceeval(pow(x_0_ij, lambda_ij)*(rhosda1Sdrhos + rhosdBdrhos));
};
// This is rhos*d(a_1ij)/drhos/(2*pi*d^3*eps*rhos)
auto da1iidrhos_term = C_ij(i,j)*(
rhosda1iidrhos_term(lambda_a_ij(i,i), zeta_x_eff_a, dzeta_x_eff_dzetax_a, Bhatij_a)
-rhosda1iidrhos_term(lambda_r_ij(i,i), zeta_x_eff_r, dzeta_x_eff_dzetax_r, Bhatij_r)
);
auto g1ii = 3.0*da1iidrhos_term + g1_noderivterm;
// The g2 terms
// ....
// This is the second part (not the partial deriv.) that goes in g_{2,ii},
// divided by 2*PI*d_ij^3*epsilon*rhos
auto g2_noderivterm = -POW2(C_ij(i,i))*K_HS*(
lambda_a_ij(i,j)*one_term(2*lambda_a_ij(i,j), zeta_x_eff_2a)
-(lambda_a_ij(i,j)+lambda_r_ij(i,j))*one_term(lambda_a_ij(i,j)+lambda_r_ij(i,j), zeta_x_eff_ar)
+lambda_r_ij(i,j)*one_term(2*lambda_r_ij(i,j), zeta_x_eff_2r)
);
// This is [rhos*d(a_2ij/(1+chi_ij))/drhos]/(2*pi*d^3*eps*rhos)
auto da2iidrhos_term = 0.5*POW2(C_ij(i,j))*(
rho_dK_HS_drho*(
one_term(2.0*lambda_a_ij(i,j), zeta_x_eff_2a)
-2.0*one_term(lambda_a_ij(i,j)+lambda_r_ij(i,j), zeta_x_eff_ar)
+one_term(2.0*lambda_r_ij(i,j), zeta_x_eff_2r))
rhosda1iidrhos_term(2.0*lambda_a_ij(i,i), zeta_x_eff_2a, dzeta_x_eff_dzetax_2a, Bhatij_2a)
-2.0*rhosda1iidrhos_term(lambda_a_ij(i,i)+lambda_r_ij(i,i), zeta_x_eff_ar, dzeta_x_eff_dzetax_ar, Bhatij_ar)
+rhosda1iidrhos_term(2.0*lambda_r_ij(i,i), zeta_x_eff_2r, dzeta_x_eff_dzetax_2r, Bhatij_2r)
)
);
auto g2MCAij = 3.0*da2iidrhos_term + g2_noderivterm;
auto betaepsilon = epsilon_ij(i,i)/T; // (1/(kB*T))/epsilon
auto theta = exp(betaepsilon)-1.0;
auto phi7 = phi.col(6);
auto gamma_cij = phi7(0)*(-tanh(phi7(1)*(phi7(2)-alpha_ij(i,j)))+1.0)*zeta_x_bar*theta*exp(phi7(3)*zeta_x_bar + phi7(4)*POW2(zeta_x_bar)); // Eq. A37
NumType giiMie = gdHSii*exp((betaepsilon*g1ii + POW2(betaepsilon)*g2ii)/gdHSii);
alphar_chain -= molefracs[i]*(m[i]-1.0)*log(giiMie);
auto ahs = get_a_HS(rhos, zeta);
auto alphar_mono = forceeval(mbar*(ahs + a1kB/T + a2kB2/(T*T) + a3kB3/(T*T*T)));
using dmat_t = decltype(dmat);
using rhos_t = decltype(rhos);
using rhoN_t = decltype(rhoN);
using mbar_t = decltype(mbar);
using xs_t = decltype(xs);
using zeta_t = decltype(zeta);
using zeta_x_t = decltype(zeta_x);
using zeta_x_bar_t = decltype(zeta_x_bar);
using alphar_mono_t = decltype(alphar_mono);
using a1kB_t = decltype(a1kB);
using a2kB2_t = decltype(a2kB2);
using a3kB3_t = decltype(a3kB3);
using alphar_chain_t = decltype(alphar_chain);
dmat_t dmat;
rhos_t rhos;
rhoN_t rhoN;
mbar_t mbar;
xs_t xs;
zeta_t zeta;
zeta_x_t zeta_x;
zeta_x_bar_t zeta_x_bar;
alphar_mono_t alphar_mono;
a1kB_t a1kB;
a2kB2_t a2kB2;
a3kB3_t a3kB3;
alphar_chain_t alphar_chain;
return vals{dmat, rhos, rhoN, mbar, xs, zeta, zeta_x, zeta_x_bar, alphar_mono, a1kB, a2kB2, a3kB3, alphar_chain};
template<typename RhoType>
auto get_KHS(const RhoType& pf) const {
return forceeval(pow(1.0-pf,4)/(1.0 + 4.0*pf + 4.0*pf*pf - 4.0*pf*pf*pf + pf*pf*pf*pf));
\rho_s\frac{\partial K_{HS}}{\partial \rho_s} = \zeta\frac{\partial K_{HS}}{\partial \zeta}
\f]
*/
template<typename RhoType>
auto get_rhos_dK_HS_drhos(const RhoType& zeta_x) const {
auto num = -4.0*POW3(zeta_x - 1.0)*(POW2(zeta_x) - 5.0*zeta_x - 2.0);
auto den = POW2(POW4(zeta_x) - 4.0*POW3(zeta_x) + 4.0*POW2(zeta_x) + 4.0*zeta_x + 1.0);
return forceeval(num/den*zeta_x);
/// Eq. A6 from Lafitte, accounting for the case of rho_s=0, for which the limit is zero
template<typename RhoType, typename ZetaType>
auto get_a_HS(const RhoType& rhos, const Eigen::Array<ZetaType, 4, 1>& zeta) const{
constexpr double MY_PI = static_cast<double>(EIGEN_PI);
// The way in which the function goes to zero is subtle, and the factor of 4 accounts for the contributions from each term
return forceeval(4.0*zeta[3]);
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])));
\hat B_{ij} \equiv \frac{B_{ij}}{2\pi\epsilon_{ij}d^3_{ij}\rho_s} = \frac{1-\zeta_x/2}{(1-\zeta_x)^3}I-\frac{9\zeta_x(1+\zeta_x)}{2(1-\zeta_x)^3}J
template<typename ZetaType, typename IJ>
auto get_Bhatij(const ZetaType& zeta_x, const ZetaType& one_minus_zeta_x3, const IJ& I, const IJ& J) const{
(1.0-zeta_x/2.0)/one_minus_zeta_x3*I - 9.0*zeta_x*(1.0+zeta_x)/(2.0*one_minus_zeta_x3)*J
\f[
B = \hat B_{ij}\kappa \rho_s
\f]
\f[
\left(\frac{\partial B_{ij}}{\partial \rho_s}\right)_{T,\vec{z}} = \kappa\left(\hat B + \zeta_x \frac{\partial \hat B}{\partial \zeta_x}\right)
\rho_s \left(\frac{\partial B_{ij}}{\partial \rho_s}\right)_{T,\vec{z}} = \hat B + \zeta_x \frac{\partial \hat B}{\partial \zeta_x}
template<typename ZetaType, typename IJ>
auto get_rhodBijdrho(const ZetaType& zeta_x, const ZetaType& one_minus_zeta_x3, const IJ& I, const IJ& J, const ZetaType& Bhatij) const{
auto dBhatdzetax = (-3.0*I*(zeta_x - 2.0) - 27.0*J*zeta_x*(zeta_x + 1.0) + (zeta_x - 1.0)*(I + 9.0*J*zeta_x + 9.0*J*(zeta_x + 1.0)))/(2.0*POW4(1.0-zeta_x));
return forceeval(Bhatij + dBhatdzetax*zeta_x);
/**
\note Starting from Eq. A16 from Lafitte
\f[
\hat a^S_{1,ii} = \frac{a^S_{1,ii}}{2\pi\epsilon_{ij}d^3_{ij}\rho_s}
\f]
so
\f[
a^S_{1,ii} = \kappa\rho_s\hat a^S_{1,ii}
\f]
*/
template<typename ZetaType>
auto get_a1Shatij(const ZetaType& zeta_x_eff, double lambda_ij) const{
-1.0/(lambda_ij-3.0)*(1.0-zeta_x_eff/2.0)/POW3(forceeval(1.0-zeta_x_eff))
/**
\f[
\left(\frac{\partial a^S_{1,ii}}{\partial \rho_s}\right)_{T,\vec{z}} = \kappa\left(\hat a^S_{1,ii} + \rho_s\frac{\partial \hat a^S_{1,ii}}{\partial \rho_s} \right)
\f]
\f[
\left(\frac{\partial a^S_{1,ii}}{\partial \rho_s}\right)_{T,\vec{z}} = \kappa\left(\hat a^S_{1,ii} + \rho_s\frac{\partial \hat a^S_{1,ii}}{\partial \zeta_{x,eff}}\frac{\partial \zeta_{x,eff}}{\partial \zeta_x}\frac{\partial \zeta_x}{\partial \rho_s} \right)
since \f$\rho_s\frac{\partial \zeta_x}{\partial \rho_s} = \zeta_x\f$
template<typename ZetaType>
auto get_rhoda1Shatijdrho(const ZetaType& zeta_x, const ZetaType& zeta_x_eff, const ZetaType& dzetaxeffdzetax, double lambda_ij) const{
auto zetaxda1Shatdzetax = ((2.0*zeta_x_eff - 5.0)*dzetaxeffdzetax)/(2.0*(lambda_ij-3)*POW4(zeta_x_eff-1.0))*zeta_x;
return forceeval(get_a1Shatij(zeta_x_eff, lambda_ij) + zetaxda1Shatdzetax);
/**
\brief A class used to evaluate mixtures using the SAFT-VR-Mie model
*/
class SAFTVRMieMixture {
private:
std::vector<std::string> names;
const SAFTVRMieChainContributionTerms terms;
const std::optional<SAFTpolar::multipolar_contributions_variant> polar; // Can be present or not
static void check_kmat(const Eigen::ArrayXXd& kmat, std::size_t N) {
if (kmat.size() == 0){
return;
}
if (kmat.cols() != kmat.rows()) {
throw teqp::InvalidArgument("kmat rows and columns are not identical");
}
if (kmat.cols() != N) {
throw teqp::InvalidArgument("kmat needs to be a square matrix the same size as the number of components");
}
};
static auto get_coeffs_from_names(const std::vector<std::string> &names){
SAFTVRMieLibrary library;
return library.get_coeffs(names);
}
public:
static auto build_chain(const std::vector<SAFTVRMieCoeffs> &coeffs, const std::optional<Eigen::ArrayXXd>& kmat){
if (kmat){
check_kmat(kmat.value(), coeffs.size());
}
const std::size_t N = coeffs.size();
Eigen::ArrayXd m(N), epsilon_over_k(N), sigma_m(N), lambda_r(N), lambda_a(N);
auto i = 0;
for (const auto &coeff : coeffs) {
m[i] = coeff.m;
epsilon_over_k[i] = coeff.epsilon_over_k;
sigma_m[i] = coeff.sigma_m;
lambda_r[i] = coeff.lambda_r;
lambda_a[i] = coeff.lambda_a;
i++;
}
if (kmat){
return SAFTVRMieChainContributionTerms(m, epsilon_over_k, sigma_m, lambda_r, lambda_a, std::move(kmat.value()));
}
else{
auto mat = Eigen::ArrayXXd::Zero(N,N);
return SAFTVRMieChainContributionTerms(m, epsilon_over_k, sigma_m, lambda_r, lambda_a, std::move(mat));
auto build_polar(const std::vector<SAFTVRMieCoeffs> &coeffs) -> decltype(this->polar){
Eigen::ArrayXd mustar2(coeffs.size()), nmu(coeffs.size()), Qstar2(coeffs.size()), nQ(coeffs.size());
auto i = 0;
for (const auto &coeff : coeffs) {
mustar2[i] = coeff.mustar2;
nmu[i] = coeff.nmu;
Qstar2[i] = coeff.Qstar2;
nQ[i] = coeff.nQ;
i++;
}
bool has_dipolar = ((mustar2*nmu).cwiseAbs().sum() == 0);
bool has_quadrupolar = ((Qstar2*nQ).cwiseAbs().sum() == 0);
if (!has_dipolar && !has_quadrupolar){
return std::nullopt; // No dipolar or quadrupolar contribution is present
}
else{
// The dispersive and hard chain initialization has already happened at this point
return SAFTpolar::MultipolarContributionGrossVrabec(terms.m, terms.sigma_A, terms.epsilon_over_k, mustar2, nmu, Qstar2, nQ);
SAFTVRMieMixture(const std::vector<std::string> &names, const std::optional<Eigen::ArrayXXd>& kmat = std::nullopt) : SAFTVRMieMixture(get_coeffs_from_names(names), kmat){};
SAFTVRMieMixture(const std::vector<SAFTVRMieCoeffs> &coeffs, const std::optional<Eigen::ArrayXXd> &kmat = std::nullopt) : terms(build_chain(coeffs, kmat)), polar(build_polar(coeffs)) {};
SAFTVRMieMixture(SAFTVRMieChainContributionTerms&& terms, std::optional<SAFTpolar::multipolar_contributions_variant> &&polar = std::nullopt) : terms(std::move(terms)), polar(std::move(polar)) {};
// PCSAFTMixture( const PCSAFTMixture& ) = delete; // non construction-copyable
SAFTVRMieMixture& operator=( const SAFTVRMieMixture& ) = delete; // non copyable
auto chain_factory(const std::vector<SAFTVRMieCoeffs> &coeffs, const std::optional<Eigen::ArrayXXd>& kmat){
SAFTVRMieMixture::build_chain(coeffs, kmat);
}
const auto& get_polar(){ return polar; }
// Checker for whether a polar term is present
bool has_polar() const{ return polar.has_value(); }
const auto& get_terms() const { return terms; }
auto get_core_calcs(double T, double rhomolar, const Eigen::ArrayXd& mole_fractions) const {
auto val = terms.get_core_calcs(T, rhomolar, mole_fractions);
auto fromArrayX = [](const Eigen::ArrayXd &x){std::valarray<double>n(x.size()); for (auto i =0; i < n.size(); ++i){ n[i] = x[i];} return n;};
auto fromArrayXX = [](const Eigen::ArrayXXd &x){
std::size_t N = x.rows();
for (auto i = 0; i < N; ++i){
n[i].resize(N);
for (auto j = 0; j < N; ++j){
n[i][j] = x(i,j);
}
}
return n;
};
return nlohmann::json{
{"dmat", fromArrayXX(val.dmat)},
{"rhos", val.rhos},
{"rhoN", val.rhoN},
{"mbar", val.mbar},
{"xs", fromArrayX(val.xs)},
{"zeta", fromArrayX(val.zeta)},
{"zeta_x", val.zeta_x},
{"zeta_x_bar", val.zeta_x_bar},
{"alphar_mono", val.alphar_mono},
{"a1kB", val.a1kB},
{"a2kB2", val.a2kB2},
{"a3kB3", val.a3kB3},
{"alphar_chain", val.alphar_chain}
};
}
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
auto get_m() const { return terms.m; }
auto get_sigma_Angstrom() const { return (terms.sigma_A).eval(); }
auto get_sigma_m() const { return terms.sigma_A/1e10; }
auto get_epsilon_over_k_K() const { return terms.epsilon_over_k; }
auto get_kmat() const { return terms.kmat; }
auto get_lambda_r() const { return terms.lambda_r; }
auto get_lambda_a() const { return terms.lambda_a; }
// template<typename VecType>
// double max_rhoN(const double T, const VecType& mole_fractions) const {
// auto N = mole_fractions.size();
// Eigen::ArrayX<double> d(N);
// for (auto i = 0; i < N; ++i) {
// d[i] = sigma_Angstrom[i] * (1.0 - 0.12 * exp(-3.0 * epsilon_over_k[i] / T));
// }
// return 6 * 0.74 / EIGEN_PI / (mole_fractions*m*powvec(d, 3)).sum()*1e30; // particles/m^3
// }
template<class VecType>
auto R(const VecType& molefrac) const {
return get_R_gas<decltype(molefrac[0])>();
}
template<typename TTYPE, typename RhoType, typename VecType>
auto alphar(const TTYPE& T, const RhoType& rhomolar, const VecType& mole_fractions) const {
// First values for the Mie chain with dispersion (always included)
error_if_expr(T); error_if_expr(rhomolar);
auto vals = terms.get_core_calcs(T, rhomolar, mole_fractions);
using type = std::common_type_t<TTYPE, RhoType, decltype(mole_fractions[0])>;
type alphar = vals.alphar_mono + vals.alphar_chain;
type packing_fraction = vals.zeta[3];
if (polar){ // polar term is present
using mas = SAFTpolar::multipolar_argument_spec;
auto visitor = [&T, &rhomolar, &mole_fractions, &packing_fraction](const auto& contrib) -> type {
constexpr mas arg_spec = std::decay_t<decltype(contrib)>::arg_spec;
if constexpr(arg_spec == mas::TK_rhoNA3_packingfraction_molefractions){
RhoType rho_A3 = rhomolar*N_A*1e-30;
type alpha = contrib.eval(T, rho_A3, packing_fraction, mole_fractions).alpha;
Ian Bell
committed
else if constexpr(arg_spec == mas::TK_rhoNm3_rhostar_molefractions){
Ian Bell
committed
auto rhostar = contrib.get_rhostar(rhoN_m3, packing_fraction, mole_fractions);
return contrib.eval(T, rhoN_m3, rhostar, mole_fractions).alpha;
}
else{
throw teqp::InvalidArgument("Don't know how to handle this kind of arguments in polar term");
}
};
alphar += std::visit(visitor, polar.value());
}
inline auto SAFTVRMiefactory(const nlohmann::json & spec){
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
std::optional<Eigen::ArrayXXd> kmat;
if (spec.contains("kmat") && spec.at("kmat").is_array() && spec.at("kmat").size() > 0){
kmat = build_square_matrix(spec["kmat"]);
}
if (spec.contains("names")){
std::vector<std::string> names = spec["names"];
if (kmat && kmat.value().rows() != names.size()){
throw teqp::InvalidArgument("Provided length of names of " + std::to_string(names.size()) + " does not match the dimension of the kmat of " + std::to_string(kmat.value().rows()));
}
return SAFTVRMieMixture(names, kmat);
}
else if (spec.contains("coeffs")){
bool something_polar = false;
std::vector<SAFTVRMieCoeffs> coeffs;
for (auto j : spec["coeffs"]) {
SAFTVRMieCoeffs c;
c.name = j.at("name");
c.m = j.at("m");
c.sigma_m = (j.contains("sigma_m")) ? j.at("sigma_m").get<double>() : j.at("sigma_Angstrom").get<double>()/1e10;
c.epsilon_over_k = j.at("epsilon_over_k");
c.lambda_r = j.at("lambda_r");
c.lambda_a = j.at("lambda_a");
c.BibTeXKey = j.at("BibTeXKey");
// These are legacy definitions of the polar moments
if (j.contains("(mu^*)^2") && j.contains("nmu")){
c.mustar2 = j.at("(mu^*)^2");
c.nmu = j.at("nmu");
something_polar = true;
}
if (j.contains("(Q^*)^2") && j.contains("nQ")){
c.Qstar2 = j.at("(Q^*)^2");
c.nQ = j.at("nQ");
something_polar = true;
}
if (j.contains("Q_Cm2") || j.contains("Q_DA") || j.contains("mu_Cm") || j.contains("mu_D")){
something_polar = true;
}
coeffs.push_back(c);
}
if (kmat && kmat.value().rows() != coeffs.size()){
throw teqp::InvalidArgument("Provided length of coeffs of " + std::to_string(coeffs.size()) + " does not match the dimension of the kmat of " + std::to_string(kmat.value().rows()));
}
if (!something_polar){
// Nonpolar, just m, epsilon, sigma and possibly a kmat matrix with kij coefficients
return SAFTVRMieMixture(SAFTVRMieMixture::build_chain(coeffs, kmat));
}
else{
// Polar term is also provided, along with the chain terms
std::string polar_model = "GrossVrabec"; // This is the default, as it was the first one implemented
if (spec.contains("polar_model")){
polar_model = spec["polar_model"];
}
// Go back and extract the dipolar and quadrupolar terms from
// the JSON, in base SI units
const double D_to_Cm = 3.33564e-30; // C m/D
const double mustar2factor = 1.0/(4*static_cast<double>(EIGEN_PI)*8.8541878128e-12*1.380649e-23); // factor=1/(4*pi*epsilon_0*k_B), such that (mu^*)^2 := factor*mu[Cm]^2/((epsilon/kB)[K]*sigma[m]^3)
const double Qstar2factor = 1.0/(4*static_cast<double>(EIGEN_PI)*8.8541878128e-12*1.380649e-23); // same as mustar2factor
auto N = coeffs.size();
Eigen::ArrayXd ms(N), epsks(N), sigma_ms(N), mu_Cm(N), Q_Cm2(N), nQ(N), nmu(N);
Eigen::Index i = 0;
for (auto j : spec["coeffs"]) {
double m = j.at("m");
double sigma_m = (j.contains("sigma_m")) ? j.at("sigma_m").get<double>() : j.at("sigma_Angstrom").get<double>()/1e10;
double epsilon_over_k = j.at("epsilon_over_k");
auto get_dipole_Cm = [&]() -> double {
if (j.contains("(mu^*)^2") && j.contains("nmu")){
// Terms defined like in Gross&Vrabec; backwards-compatibility
double mustar2 = j.at("(mu^*)^2");
return sqrt(mustar2*(m*epsilon_over_k*pow(sigma_m, 3))/mustar2factor);
}
else if (j.contains("mu_Cm")){
return j.at("mu_Cm");
}
else if (j.contains("mu_D")){
return j.at("mu_D").get<double>()*D_to_Cm;
}
else{
return 0.0;
}
};
auto get_quadrupole_Cm2 = [&]() -> double{
if (j.contains("(Q^*)^2") && j.contains("nQ")){
// Terms defined like in Gross&Vrabec; backwards-compatibility
double Qstar2 = j.at("(Q^*)^2");
return sqrt(Qstar2*(m*epsilon_over_k*pow(sigma_m, 5))/Qstar2factor);
}
else if (j.contains("Q_Cm2")){
return j.at("Q_Cm2");
}
else if (j.contains("Q_DA")){
return j.at("Q_DA").get<double>()*D_to_Cm*1e10;
}
else{
return 0.0;
}
};
ms(i) = m; sigma_ms(i) = sigma_m; epsks(i) = epsilon_over_k; mu_Cm(i) = get_dipole_Cm(); Q_Cm2(i) = get_quadrupole_Cm2();
nmu(i) = (j.contains("nmu") ? j["nmu"].get<double>() : 0.0);
nQ(i) = (j.contains("nQ") ? j["nQ"].get<double>() : 0.0);
i++;
};
using namespace SAFTpolar;
if (polar_model == "GrossVrabec"){
auto mustar2 = (mustar2factor*mu_Cm.pow(2)/(ms*epsks*sigma_ms.pow(3))).eval();
auto Qstar2 = (Qstar2factor*Q_Cm2.pow(2)/(ms*epsks*sigma_ms.pow(5))).eval();
auto polar = MultipolarContributionGrossVrabec(ms, sigma_ms*1e10, epsks, mustar2, nmu, Qstar2, nQ);
return SAFTVRMieMixture(SAFTVRMieMixture::build_chain(coeffs, kmat), std::move(polar));
}
if (polar_model == "GubbinsTwu+Luckas"){
using MCGTL = MultipolarContributionGubbinsTwu<LuckasJIntegral, LuckasKIntegral>;
auto mubar2 = (mustar2factor*mu_Cm.pow(2)/(epsks*sigma_ms.pow(3))).eval();
auto Qbar2 = (Qstar2factor*Q_Cm2.pow(2)/(epsks*sigma_ms.pow(5))).eval();
Ian Bell
committed
auto polar = MCGTL(sigma_ms, epsks, mubar2, Qbar2, multipolar_rhostar_approach::use_packing_fraction);
return SAFTVRMieMixture(SAFTVRMieMixture::build_chain(coeffs, kmat), std::move(polar));
}
if (polar_model == "GubbinsTwu+GubbinsTwu"){
using MCGG = MultipolarContributionGubbinsTwu<GubbinsTwuJIntegral, GubbinsTwuKIntegral>;
auto mubar2 = (mustar2factor*mu_Cm.pow(2)/(epsks*sigma_ms.pow(3))).eval();
auto Qbar2 = (Qstar2factor*Q_Cm2.pow(2)/(epsks*sigma_ms.pow(5))).eval();
Ian Bell
committed
auto polar = MCGG(sigma_ms, epsks, mubar2, Qbar2, multipolar_rhostar_approach::use_packing_fraction);
return SAFTVRMieMixture(SAFTVRMieMixture::build_chain(coeffs, kmat), std::move(polar));
}
if (polar_model == "GubbinsTwu+Luckas+GubbinsTwuRhostar"){