From 3a871b3e19ad7ee03bca35ec3d7a650fea3ef98b Mon Sep 17 00:00:00 2001 From: Ian Bell Date: Wed, 28 Aug 2024 12:16:30 -0400 Subject: [PATCH] Add COSMO-SAC-2010 as an option Other COSMO-SAC variants are possible, would require a bit of refactoring --- include/teqp/models/activity/COSMOSAC.hpp | 360 ++++++++++++++++++ .../teqp/models/activity/activity_models.hpp | 21 +- 2 files changed, 380 insertions(+), 1 deletion(-) create mode 100644 include/teqp/models/activity/COSMOSAC.hpp diff --git a/include/teqp/models/activity/COSMOSAC.hpp b/include/teqp/models/activity/COSMOSAC.hpp new file mode 100644 index 0000000..cd05e50 --- /dev/null +++ b/include/teqp/models/activity/COSMOSAC.hpp @@ -0,0 +1,360 @@ +/** + This header is derived from the code in https://github.com/usnistgov/COSMOSAC, and the methodology is described in the paper: + + Ian H. Bell, Erik Mickoleit, Chieh-Ming Hsieh, Shiang-Tai Lin, Jadran Vrabec, Cornelia Breitkopf, and Andreas Jäger + Journal of Chemical Theory and Computation 2020 16 (4), 2635-2646 + DOI: 10.1021/acs.jctc.9b01016 + */ + +namespace teqp::activity::activity_models::COSMOSAC{ + +/** + A single sigma profile. In this data structure (for consistency with the 2005 Virginia Tech + database of COSMO-SAC parameters), the first column is the electron density in [e/A^2], and the second column + is the probability of finding a segment with this electron density, multiplied by the segment area, in A^2 + */ +class SigmaProfile { +public: + Eigen::ArrayXd m_sigma, m_psigmaA; + /// Default constructor + SigmaProfile() {}; + /// Copy-by-reference constructor + SigmaProfile(const Eigen::ArrayXd &sigma, const Eigen::ArrayXd &psigmaA) : m_sigma(sigma), m_psigmaA(psigmaA) {}; + /// Copy-by-reference constructor with std::vector + SigmaProfile(const std::vector &sigma, const std::vector &psigmaA) : m_sigma(Eigen::Map(&(sigma[0]), sigma.size())), m_psigmaA(Eigen::Map(&(psigmaA[0]), psigmaA.size())) {}; + /// Move constructor + SigmaProfile(Eigen::ArrayXd &&sigma, Eigen::ArrayXd &&psigmaA) : m_sigma(sigma), m_psigmaA(psigmaA) {}; + const Eigen::ArrayXd &psigmaA() const { return m_psigmaA; } + const Eigen::ArrayXd &sigma() const { return m_sigma; } + const Eigen::ArrayXd psigma(double A_i) const { return m_psigmaA/A_i; } +}; + +// A holder class for the sigma profiles for a given fluid +struct FluidSigmaProfiles { + SigmaProfile nhb, ///< The profile for the non-hydrogen-bonding segments + oh, ///< The profile for the OH-bonding segments + ot; ///< The profile for the "other" segments +}; + +struct COSMO3Constants { + double + AEFFPRIME = 7.25, // [A^2] + c_OH_OH = 4013.78, // [kcal A^4/(mol e^2)] + c_OT_OT = 932.31, // [kcal A^4 /(mol e^2)] + c_OH_OT = 3016.43, // [kcal A^4 /(mol e^2)] + A_ES = 6525.69, // [kcal A^4 /(mol e^2)] + B_ES = 1.4859e8, // [kcal A^4 K^2/(mol e^2)] + N_A = 6.022140758e23, // [mol^{-1}] + k_B = 1.38064903e-23, // [J K^{-1}] + R = k_B*N_A/4184, // [kcal/(mol*K)]; Universal gas constant of new redefinition of 2018, https://doi.org/10.1088/1681-7575/aa99bc + Gamma_rel_tol = 1e-8; // relative tolerance for Gamma in iterative loop + bool fast_Gamma = false; + std::string to_string() { + return "NOT IMPLEMENTED YET"; + //return "c_hb: " + std::to_string(c_hb) + " kcal A^4 /(mol*e^2) \nsigma_hb: " + std::to_string(sigma_hb) + " e/A^2\nalpha_prime: " + std::to_string(alpha_prime) + " kcal A^4 /(mol*e^2)\nAEFFPRIME: " + std::to_string(AEFFPRIME) + " A\nR: " + std::to_string(R) + " kcal/(mol*K)"; + } +}; + +enum class profile_type { NHB_PROFILE, OH_PROFILE, OT_PROFILE }; + +class COSMO3 { +private: + std::vector A_COSMO_A2; ///< The area per fluid, in \AA^2 + std::vector profiles; ///< The vector of profiles, one per fluid + COSMO3Constants m_consts; + Eigen::Index ileft, w; +public: + COSMO3(const std::vector& A_COSMO_A2, const std::vector &SigmaProfiles, const COSMO3Constants &constants = COSMO3Constants()) + : A_COSMO_A2(A_COSMO_A2), profiles(SigmaProfiles), m_consts(constants) { + Eigen::Index iL, iR; + std::tie(iL, iR) = get_nonzero_bounds(); + this->ileft = iL; this->w = iR - iL + 1; + }; + + /// Determine the left and right bounds for the psigma that is nonzero to accelerate the iterative calculations + std::tuple get_nonzero_bounds(){ + // Determine the range of entries in p(sigma) that are greater than zero, we + // will only calculate segment activity coefficients for these segments + Eigen::Index min_ileft = 51, max_iright = 0; + for (auto i = 0U; i < profiles.size(); ++i) { + const Eigen::ArrayXd psigma = profiles[i].nhb.psigma(A_COSMO_A2[i]); + Eigen::Index ileft = 0, iright = psigma.size(); + for (auto ii = 0; ii < psigma.size(); ++ii) { if (std::abs(psigma(ii)) > 1e-16) { ileft = ii; break; } } + for (auto ii = psigma.size() - 1; ii > ileft; --ii) { if (std::abs(psigma(ii)) > 1e-16) { iright = ii; break; } } + if (ileft < min_ileft) { min_ileft = ileft; } + if (iright > max_iright) { max_iright = iright; } + } + return std::make_tuple(min_ileft, max_iright); + } + + template + auto get_psigma_mix(const MoleFractions &z, profile_type type = profile_type::NHB_PROFILE) const { + Eigen::ArrayX> psigma_mix(51); psigma_mix.fill(0); + std::decay_t xA = 0.0; + for (auto i = 0; i < z.size(); ++i) { + switch (type) { + case profile_type::NHB_PROFILE: + psigma_mix += z[i] * profiles[i].nhb.psigmaA(); break; + case profile_type::OH_PROFILE: + psigma_mix += z[i] * profiles[i].oh.psigmaA(); break; + case profile_type::OT_PROFILE: + psigma_mix += z[i] * profiles[i].ot.psigmaA(); break; + } + xA += z[i] * A_COSMO_A2[i]; + } + psigma_mix /= xA; + return psigma_mix; + } + + /// Get access to the parameters that are in use + COSMO3Constants &get_mutable_COSMO_constants() { return m_consts; } + + std::tuple get_ileftw() const { return std::make_tuple(ileft, w);} ; + + double get_c_hb(profile_type type1, profile_type type2) const{ + + if (type1 == type2){ + if (type1 == profile_type::OH_PROFILE) { + return m_consts.c_OH_OH; + } + else if (type1 == profile_type::OT_PROFILE) { + return m_consts.c_OT_OT; + } + else { + return 0.0; + } + } + else if ((type1 == profile_type::OH_PROFILE && type2 == profile_type::OT_PROFILE) + || (type1 == profile_type::OT_PROFILE && type2 == profile_type::OH_PROFILE)) { + return m_consts.c_OH_OT; + } + else { + return 0.0; + } + } + template + auto get_DELTAW(const TType& T, profile_type type_t, profile_type type_s) const { + auto delta_sigma = 2*0.025/50; + Eigen::ArrayXX DELTAW(51, 51); + double cc = get_c_hb(type_t, type_s); + for (auto m = 0; m < 51; ++m) { + for (auto n = 0; n < 51; ++n) { + double sigma_m = -0.025 + delta_sigma*m, + sigma_n = -0.025 + delta_sigma*n, + c_hb = (sigma_m*sigma_n >= 0) ? 0 : cc; + auto c_ES = m_consts.A_ES + m_consts.B_ES/(T*T); + DELTAW(m, n) = c_ES*POW2(sigma_m + sigma_n) - c_hb*POW2(sigma_m-sigma_n); + } + } + return DELTAW; + } + template + auto get_DELTAW_fast(TType T, profile_type type_t, profile_type type_s) const { + auto delta_sigma = 2 * 0.025 / 50; + Eigen::ArrayXX DELTAW(51, 51); DELTAW.setZero(); + double cc = get_c_hb(type_t, type_s); + for (auto m = ileft; m < ileft+w+1; ++m) { + for (auto n = ileft; n < ileft+w+1; ++n) { + double sigma_m = -0.025 + delta_sigma*m, + sigma_n = -0.025 + delta_sigma*n, + c_hb = (sigma_m*sigma_n >= 0) ? 0 : cc; + auto c_ES = m_consts.A_ES + m_consts.B_ES / (T*T); + DELTAW(m, n) = c_ES*POW2(sigma_m + sigma_n) - c_hb*POW2(sigma_m - sigma_n); + } + } + return DELTAW; + } + + template + auto get_AA(const TType& T, PSigma psigmas){ + // Build the massive \Delta W matrix that is 153*153 in size + Eigen::ArrayXX DELTAW(51, 51); // Depends on temperature + double R = m_consts.R; + std::vector types = { profile_type::NHB_PROFILE, profile_type::OH_PROFILE, profile_type::OT_PROFILE }; + for (auto i = 0; i < 3; ++i) { + for (auto j = 0; j < 3; ++j) { + DELTAW.matrix().block(51 * i, 51 * j, 51, 51) = get_DELTAW(T, types[i], types[j]); + } + } + return Eigen::exp(-DELTAW/(R*T)).rowwise()*psigmas.transpose(); + } + /** + Obtain the segment activity coefficients \f$\Gamma\f$ for given a set of charge densities + + If fast_Gamma is enabled, then only segments that have some contribution are included in the iteration (can be a significant acceleration for nonpolar+nonpolar mixtures) + + \param T Temperature, in K + \param psigmas Charge densities, in the order of NHB, OH, OT. Length is 153. + */ + template + auto get_Gamma(const TType& T, PSigmaType psigmas) const { + auto startTime = std::chrono::high_resolution_clock::now(); + + using TXType = std::decay_t>; + + double R = m_consts.R; + Eigen::ArrayX Gamma(153), Gammanew(153); Gamma.setOnes(); Gammanew.setOnes(); + + // A convenience function to convert double values to string in scientific format + auto to_scientific = [](double val) { std::ostringstream out; out << std::scientific << val; return out.str(); }; + + auto max_iter = 2000; + if (!m_consts.fast_Gamma){ + // The slow and simple branch + + // Build the massive \Delta W matrix that is 153*153 in size + Eigen::ArrayXX DELTAW(153, 153); + std::vector types = { profile_type::NHB_PROFILE, profile_type::OH_PROFILE, profile_type::OT_PROFILE }; + for (auto i = 0; i < 3; ++i) { + for (auto j = 0; j < 3; ++j) { + DELTAW.matrix().block(51*i, 51*j, 51, 51) = get_DELTAW(T, types[i], types[j]); + } + } + + auto AA = (Eigen::exp(-DELTAW / (R*T)).template cast().rowwise()*psigmas.template cast().transpose()).eval(); + for (auto counter = 0; counter <= max_iter; ++counter) { + Gammanew = 1 / (AA.rowwise()*Gamma.transpose()).rowwise().sum(); + Gamma = (Gamma + Gammanew) / 2; + double maxdiff = getbaseval(((Gamma - Gammanew) / Gamma).cwiseAbs().real().maxCoeff()); + if (maxdiff < m_consts.Gamma_rel_tol) { + break; + } + if (counter == max_iter) { + throw std::invalid_argument("Could not obtain the desired tolerance of " + + to_scientific(m_consts.Gamma_rel_tol) + + " after " + + std::to_string(max_iter) + + " iterations in get_Gamma; current value is " + + to_scientific(maxdiff)); + } + } + return Gamma; + } + else { + // The fast branch! + // ---------------- + + // Build the massive AA matrix that is 153*153 in size + auto midTime = std::chrono::high_resolution_clock::now(); + std::vector types = { profile_type::NHB_PROFILE, profile_type::OH_PROFILE, profile_type::OT_PROFILE }; + std::vector offsets = {0*51, 1*51, 2*51}; + Eigen::ArrayXX AA(153, 153); + for (auto i = 0; i < 3; ++i) { + for (auto j = 0; j < 3; ++j) { + Eigen::Index rowoffset = offsets[i], coloffset = offsets[j]; + AA.matrix().block(rowoffset + ileft, coloffset + ileft, w, w) = Eigen::exp(-get_DELTAW_fast(T, types[i], types[j]).block(ileft, ileft, w, w).array() / (R*T)).template cast().rowwise()*psigmas.template cast().segment(coloffset+ileft,w).transpose(); + } + } + auto midTime2 = std::chrono::high_resolution_clock::now(); + + for (auto counter = 0; counter <= max_iter; ++counter) { + for (Eigen::Index offset : {51*0, 51*1, 51*2}){ + Gammanew.segment(offset + ileft, w) = 1 / ( + AA.matrix().block(offset+ileft,51*0+ileft,w,w).array().rowwise()*Gamma.segment(51*0+ileft, w).transpose() + + AA.matrix().block(offset+ileft,51*1+ileft,w,w).array().rowwise()*Gamma.segment(51*1+ileft, w).transpose() + + AA.matrix().block(offset+ileft,51*2+ileft,w,w).array().rowwise()*Gamma.segment(51*2+ileft, w).transpose() + ).rowwise().sum(); + } + for (Eigen::Index offset : {51 * 0, 51 * 1, 51 * 2}) { + Gamma.segment(offset + ileft, w) = (Gamma.segment(offset + ileft, w) + Gammanew.segment(offset + ileft, w)) / 2; + } + double maxdiff = getbaseval(((Gamma - Gammanew) / Gamma).cwiseAbs().real().maxCoeff()); + if (maxdiff < m_consts.Gamma_rel_tol) { + break; + } + if (counter == max_iter){ + throw std::invalid_argument("Could not obtain the desired tolerance of " + + to_scientific(m_consts.Gamma_rel_tol) + +" after " + +std::to_string(max_iter) + +" iterations in get_Gamma; current value is " + + to_scientific(maxdiff)); + } + } + auto endTime = std::chrono::high_resolution_clock::now(); + //std::cout << std::chrono::duration(midTime - startTime).count() << " s elapsed (DELTAW)\n"; + //std::cout << std::chrono::duration(midTime2 - midTime).count() << " s elapsed (AA)\n"; + //std::cout << std::chrono::duration(endTime - midTime2).count() << " s elapsed (comps)\n"; + //std::cout << std::chrono::duration(endTime - startTime).count() << " s elapsed (total)\n"; + return Gamma; + } + } + + /** + The residual part of ln(γ_i), for the i-th component + */ + template + auto get_lngamma_resid(std::size_t i, TType T, const Array &lnGamma_mix) const + { + double AEFFPRIME = m_consts.AEFFPRIME; + Eigen::ArrayX psigmas(3*51); // For a pure fluid, p(sigma) does not depend on temperature + double A_i = A_COSMO_A2[i]; + psigmas << profiles[i].nhb.psigma(A_i), profiles[i].oh.psigma(A_i), profiles[i].ot.psigma(A_i); + // double check_sum = psigmas.sum(); //// Should sum to 1.0 + auto lnGammai = get_Gamma(T, psigmas).log().eval(); + return A_i/AEFFPRIME*(psigmas*(lnGamma_mix - lnGammai)).sum(); + } + /** + This overload is a convenience overload, less computationally + efficient, but simpler to use and more in keeping with the other + contributions. It does not require the lnGamma_mix to be pre-calculated + and passed into this function + */ + + template + auto get_lngamma_resid(TType T, const MoleFracs& molefracs) const + { + using TXType = std::decay_t>; + + Eigen::Array psigmas; + psigmas << get_psigma_mix(molefracs, profile_type::NHB_PROFILE), get_psigma_mix(molefracs, profile_type::OH_PROFILE), get_psigma_mix(molefracs, profile_type::OT_PROFILE); + + Eigen::ArrayX lngamma(molefracs.size()); + // double check_sum = psigmas.sum(); //// Should sum to 1.0 + Eigen::ArrayX lnGamma_mix = get_Gamma(T, psigmas).log(); + for (Eigen::Index i = 0; i < molefracs.size(); ++i) { + lngamma(i) = get_lngamma_resid(i, T, lnGamma_mix); + } + return lngamma; + } + + // EigenArray get_lngamma_disp(const EigenArray &x) const { + // if (x.size() != 2){ throw std::invalid_argument("Multi-component mixtures not supported for dispersive contribution yet"); } + // double w = 0.27027; // default value + // auto cls0 = m_fluids[0].dispersion_flag, cls1 = m_fluids[1].dispersion_flag; + // using d = FluidProfiles::dispersion_classes; + // if ((cls0 == d::DISP_WATER && cls1 == d::DISP_ONLY_ACCEPTOR) + // | + // (cls1 == d::DISP_WATER && cls0 == d::DISP_ONLY_ACCEPTOR)){ + // w = -0.27027; + // } + // else if ((cls0 == d::DISP_WATER && cls1 == d::DISP_COOH) + // | + // (cls1 == d::DISP_WATER && cls0 == d::DISP_COOH)) { + // w = -0.27027; + // } + // else if ((cls0 == d::DISP_COOH && (cls1 == d::DISP_NHB || cls1 == d::DISP_DONOR_ACCEPTOR)) + // | + // (cls1 == d::DISP_COOH && (cls0 == d::DISP_NHB || cls0 == d::DISP_DONOR_ACCEPTOR))) { + // w = -0.27027; + // } + // + // double ekB0 = m_fluids[0].dispersion_eoverkB, ekB1 = m_fluids[1].dispersion_eoverkB; + // double A = w*(0.5*(ekB0+ekB1) - sqrt(ekB0*ekB1)); + // EigenArray lngamma_dsp(2); + // lngamma_dsp(0) = A*x[1]*x[1]; + // lngamma_dsp(1) = A*x[0]*x[0]; + // return lngamma_dsp; + // } + // EigenArray get_lngamma(double T, const EigenArray &x) const override { + // return get_lngamma_comb(T, x) + get_lngamma_resid(T, x) + get_lngamma_disp(x); + // } + + // Returns ares/RT + template + auto operator () (const TType& T, const MoleFractions& molefracs) const { + return contiguous_dotproduct(molefracs, get_lngamma_resid(T, molefracs)); + } +}; + +}; diff --git a/include/teqp/models/activity/activity_models.hpp b/include/teqp/models/activity/activity_models.hpp index e725afa..c5ea73d 100644 --- a/include/teqp/models/activity/activity_models.hpp +++ b/include/teqp/models/activity/activity_models.hpp @@ -1,5 +1,6 @@ #pragma once +#include "teqp/models/activity/COSMOSAC.hpp" namespace teqp::activity::activity_models{ @@ -140,7 +141,7 @@ class BinaryInvariantResidualHelmholtzOverRT { } }; -using ResidualHelmholtzOverRTVariant = std::variant, WilsonResidualHelmholtzOverRT, BinaryInvariantResidualHelmholtzOverRT>; +using ResidualHelmholtzOverRTVariant = std::variant, WilsonResidualHelmholtzOverRT, BinaryInvariantResidualHelmholtzOverRT, COSMOSAC::COSMO3>; inline ResidualHelmholtzOverRTVariant ares_model_factory(const nlohmann::json& armodel) { @@ -159,6 +160,24 @@ inline ResidualHelmholtzOverRTVariant ares_model_factory(const nlohmann::json& a std::vector c = armodel.at("c"); return BinaryInvariantResidualHelmholtzOverRT(c); } + else if (type == "COSMO-SAC-2010"){ + std::vector A_COSMOSAC_A2 = armodel.at("A_COSMOSAC / A^2"); + std::vector profiles; + for (auto& el : armodel.at("profiles")){ + COSMOSAC::FluidSigmaProfiles prof; + auto get_ = [](const auto& j){ + return COSMOSAC::SigmaProfile{ + j.at("sigma / e/A^2").template get>(), + j.at("p(sigma)*A / A^2").template get>() + }; + }; + prof.nhb = get_(el.at("nhb")); + prof.oh = get_(el.at("oh")); + prof.ot = get_(el.at("ot")); + profiles.push_back(prof); + } + return COSMOSAC::COSMO3(A_COSMOSAC_A2, profiles); + } else{ throw teqp::InvalidArgument("bad type of ares model: " + type); }