From a2712a838599dfacf161f2c6128d7e9737cbe870 Mon Sep 17 00:00:00 2001 From: Philipp Rehner Date: Wed, 22 Mar 2023 10:56:11 +0100 Subject: [PATCH 1/7] More generic association models --- src/pcsaft/parameters.rs | 12 +++++++++++- 1 file changed, 11 insertions(+), 1 deletion(-) diff --git a/src/pcsaft/parameters.rs b/src/pcsaft/parameters.rs index 7053aa6a6..04fa3a067 100644 --- a/src/pcsaft/parameters.rs +++ b/src/pcsaft/parameters.rs @@ -151,13 +151,23 @@ impl FromSegments for PcSaftRecord { impl FromSegments for PcSaftRecord { fn from_segments(segments: &[(Self, usize)]) -> Result { // We do not allow more than a single segment for q, mu, kappa_ab, epsilon_k_ab + let polar_segments: usize = segments + .iter() + .filter_map(|(s, n)| { + if s.q.is_some() || s.mu.is_some() || s.association_record.is_some() { + Some(n) + } else { + None + } + }) + .sum(); let quadpole_segments: usize = segments.iter().filter_map(|(s, n)| s.q.map(|_| n)).sum(); let dipole_segments: usize = segments.iter().filter_map(|(s, n)| s.mu.map(|_| n)).sum(); let assoc_segments: usize = segments .iter() .filter_map(|(s, n)| s.association_record.map(|_| n)) .sum(); - if quadpole_segments + dipole_segments + assoc_segments > 1 { + if polar_segments > 1 { return Err(ParameterError::IncompatibleParameters(format!( "Too many polar/associating segments (dipolar: {dipole_segments}, quadrupolar {quadpole_segments}, associating: {assoc_segments})." ))); From c853858ab1b95846991ef1cdedf1f08ebc4f6118 Mon Sep 17 00:00:00 2001 From: Philipp Rehner Date: Fri, 24 Mar 2023 16:45:48 +0100 Subject: [PATCH 2/7] Rewrite association contribution to be more generic --- docs/theory/models/FeOs_Association.png | Bin 0 -> 29101 bytes docs/theory/models/association.md | 83 ++++ docs/theory/models/hard_spheres.md | 3 + docs/theory/models/index.md | 3 + parameters/pcsaft/sauer2014_hetero.json | 8 +- .../pcsaft/sauer2014_hetero_joback.json | 8 +- parameters/pcsaft/sauer2014_homo.json | 28 +- parameters/pcsaft/sauer2014_homo_joback.json | 8 +- src/association/dft.rs | 193 ++++++--- src/association/mod.rs | 394 +++++++++++------- src/association/python.rs | 19 +- src/gc_pcsaft/dft/mod.rs | 2 +- src/gc_pcsaft/dft/parameter.rs | 8 +- src/gc_pcsaft/eos/mod.rs | 2 +- src/gc_pcsaft/eos/parameter.rs | 50 ++- src/gc_pcsaft/python/mod.rs | 2 +- src/pcsaft/dft/mod.rs | 2 +- src/pcsaft/dft/pure_saft_functional.rs | 35 +- src/pcsaft/eos/mod.rs | 2 +- src/pcsaft/eos/qspr.rs | 9 +- src/pcsaft/parameters.rs | 61 ++- src/pcsaft/python.rs | 11 +- tests/gc_pcsaft/dft.rs | 56 +-- tests/pcsaft/test_parameters.json | 4 +- 24 files changed, 648 insertions(+), 343 deletions(-) create mode 100644 docs/theory/models/FeOs_Association.png create mode 100644 docs/theory/models/association.md create mode 100644 docs/theory/models/hard_spheres.md diff --git a/docs/theory/models/FeOs_Association.png b/docs/theory/models/FeOs_Association.png new file mode 100644 index 0000000000000000000000000000000000000000..dd43fd2f9092ee84f49c454e5858ae9bed022c62 GIT binary patch literal 29101 zcmd?Rc|4R~_&=;qMI}^9MW|#cgG$I=DUujj8@n?0CF@`eeNsubGTbt@>{Q4-7>uP7 znwd}tW6f^L*crR$9?|FfdtT39&+GZ;_x$Gd>Q!^keXetz>s;qt@AJN{Q>d#?xxLNNdA*$F)nQIpya{k6kWPvaU(e&Y!m^ktXbRqd-REQL|*tJb?& zSpHGay=`y?TI#=eV zu9*ufqDsfb1X#{p zwb+&OELW~vxn_G#O^EZlwu-v^1$AfGO<5LAcgb6#T(=CaYojjc zT#?mPV$n6#)D4i<6F+G1NXAHA$@m`JI7II5H5P=_AyeoRu3&!qs=0%qMYNLDO%^LF zD{JY)_a2zsi&H`BYao-Z*vs+Q>#*3PEFIF+oO5otAa1(`A|EIoePF=yAWy^F>W=pl zCm)q#z9uZ_@CN}%Bg{44AS;%jprC)TKEY0APt^HBY*|8Z0bw52;kWoBQvZ!|X2JT} z#^?yf96ZB3d1UaYpvb-*dPJ7~>Zge`U@>O;q%NrA+H({uEpOXq6vw#-D zU~~XTNhdsKArJ^<_e9E?bIV^ns(!&zP0jy+6#JAMQv0*09$@+6aF&!A-sB?Ql*sb^ zZFIXktUZl|3dZ$%OZMfk^p_@$_(_fCvy4@x(J>b$i&!QnC#Pz2809RBl)Q}(EE^je z5;b$?kn8OB($mpr>7&oBXix`VLl5>oxQXy&VY&E$`DcfijY$^^i|cRA>sJkYE$BTP zgJ%peqbrqiG&JeO6LZ(aw6cGQl%!hWho}3`iKoT8BSa1@9XfL{-^1OSgY)r)%dE!K zBSII(PVe1Ob4-#w`pDG-#~!nT8(nmU{Nx*zY>OY3I^N!lrLRhYs$*W>hNuf2@ZXf& zH$nff@SWfB@Zg#Va>scE@;_DqF}pu3g5av)#sn9t!Dw(bIHyH5X}wFHaX@cWXG13t zt*Te8yr#GmQ3S&Vt!%6{@#Bk22PIl*1GMTSwaHW7_RA=6VU| z&KzP|d17H{us8a^1_JuNd;7a7z}tMU55D!j9UDP7@3JB|(X7xo4@pC1(vubuk%FIQyqlM;(L@j%VfG*1I@k zz&xwADsD^e2!IT6C@jc?6{o?TgE>n`s=3Ly>2}kr2p@oFMQ^LAi2-1qSzCqOYRUs3 zHK3*$TTN?`RBe-SzU`)(2%l)4m22BgjR5nE+N!9nrT~CsgPNvqH62D$15L)ewyfg# z1awxGSU|>P8pYhW!&7_Zjlc`(%#8SZ8qKf|dt)^+9^I3_7 zBTm3DHCVcD!GSpasNjkG*|v-#n44ucG5!J7(8oJrgekXt$lCd~dRx{kN1g6fd) zAhjYbRD1MjsJBgeX+^NEn_-xB6V*@Lqaog2gvHjRn{vtLOt2Eb@X|_SL|@srxRUoh zU?;9y)=HZ-4tdpB1iOOC%Dek5mSpp#5#f#~dBL-?SWKRIaJQz8w>qs@dEobfq{B{n zfx8$D?1OQVUh}d(p+=_InFb&gxpccXmN5xOsjK{ycByUfGj8k+dbvSLq)Q(;D@7S; zxtCIijvVyx*VPOHz` z=E$T{k>7xD&+nyjIJ;!`Cv^&T-TZ2KJ2sDWJd;mCjXYE9wJ~UtVlO%CfLi_?L7a8o zP?zq+?LS9(3TauXj%dcewMA8rDh6=y5>2rMi@Zev%PdY49GK79yC-eC_NmV0^2H;QOnNWhcpw%`x7C&ax5lEm;x&UC|1-;(U_fzj9Bm9vvIDt?va)>V~nG2xyUfGOtnX+ zv{_VLV{kfRW>B`KJ4HaIv3;mR!Rg|A7@u|+W@-xjjDIWlcq4o@A}{Q>ZVm2pLupq{ zkotNZ%TBU~{%YgNz7xLf*{n-&jF9VW!NtCLhl%|u&bQh#N?{r;LVQUnvxUB^!aWky z46Pc%jUKgPcWsdVDwacnA_BRGVkyl1O2@F&(~ak=NiJbd9)ZeE9;a80=bl#^x?Vx# zNA?}}QvV+GSv27-rjs?(X{G}1Ni)rK(lZ0Z2N}EUP$=1f{)kiQEvz@2_9gw`zlv|f z?U$smLY7u6f(81_hsJ#hx*@t+DG|m8=W54`xINi;ZZSh3$Hoo!*m(ZjpVD@d%grXl zO7^rlc3gL+Rl!Nm1ZetP5UMiCZF)hmtITU`WXgQPHl@SFPWx#b($k?3&Tj%-zB%!R z(A+mWW%*rU0B06>kKuHPjbt8h&lss$&`y(#cRhR9buY(@d`%SbjPfh9a*-#Hv6XH8 zWy^qsrxG~j_3QnLpLKfPjzx|7x#=nJ+SHT+U^lvyt1H*_QyqO=Onfqej zoI`|H#o9;|KEuf)1@4qlGb|gBDBQDC-USsW1IKi>e@gU;gP*8ZKvE4%#--T+Vb184 z&!*JLr5iz>Zw2;PVhndP^`nYaPKAv52Dm+(0W(6h?FNS)jqjNl2_a`XUE~uzx-|6K z^8qRs1iM=}U6Xh(?Q?-?G|L?KdF_*5Oa@bnIy5-NDeUgp;U}vMrB&lgX2fZO@3`?d zFj}ob(gl0!m*R%Lcf#gaV+@YsTKG-AlxToJ;C9%L6y;g?bi`U%oL{JG@kdf`n~ZaZ z#0n>y7i&F7t)Z)g0%kc{RO$ZO11HNKSYKV#pS$bp_G|#G5U9RT`Z-&sYN(}Q1kRD^ z6di>{PgW8!2zflHL9%`$#p zvdiYy4UY2keM`-NF|zArCFesG-7i&su$oUQJCbJ(aAee~xZXqvu#1CI`~T^3L9DL3 zSgTnl_RZ?fwhi^4w}ytDk?lFFc9s?aYg&wlcfa=1PM&u7ZlZx8lcOV$)aD1Jfz-Re z_ja2#Fg%54rDzA!Zx~OHDoRj7OQ+{z7tgKAcLpM&8%mwE1;X{Ah@>T}BiSYCC*QM^ zmi&hLaA|PJLI*C-br%~AzQbr`e$N-g6rqk zap8A>LTG^-TJWe6dF2p3CV=yxLpCzLTBISvBNM5qkuBVQU<0OD2R<=R4bYC(%E!x; zy$>+9fAbnf?dkFt32c`PH0cQ)jsr_j=E6ZC{kG*;o(TwI9_cKeZCPl89dY}gzY!H1o`zdH>F}OcJ_j_o_Cr=%)L7gpiz&C;D70x zQmx`mszici?`9cEQP!pt(hRv8oFDhGigZGX8US}7MGaP~SABivR~H-yrY?5QOv?AN zWTb6CW^vj9v7Q8yI%sgs2KFxmGa_hXS?i#kSQ z)>%{b^(80XjwlGGuP0;dJDs(oMfZ;TYDYto!mn5(cvMfjf8f5%YFZ&8h_KT(lG5^d zh3wX1VdtNU-2(|qk41hi>ZpeV4N`xE*~wkmTOh=Dzgn5C2GuBwfPC7gJ%f0hfCX&` z`9=g6DG5F^><>z`K=)-zHoY#Z-KGKzkKZPsHk0l=gY6HC@i1cX{OVt?uO69@dUGf24c#%G(m|lL{i^r=(#N&h&OAC`s z-+z6XM-qAZPZZ?B&G!WLj3D>GsWgcSRO~NVLN8qPpfNWT@OWH4X{vrtOl)S5SLbh9 zXm6epzk4sPVrkxH=BEsz3>GN2YHx}vrx#Dl9LK}%XV$8uzk^MWJC>88*Plx=G}kQT z16LP)Q8qJd)j$C$bg*D4MLjLST^tzB9MlYlHYSTm%cB{LhXKm+RN@@l+QMv;YF`=c zYMGzYjIGB-e-chfbuw{|j4Znl2)|^%pY**Zge%H@+NwmvBY`~3**6&U%B~d*kKtK) zmIf|Te|V``*IItRbj#&c-P)_fjM-XgKsq^{;=F6WnRI#7XQ$ zseqNp(p2fP)t~g)e2SzvCGopMi=fxV32%}_w*aB*%@CaXxYY!%BB{GrnD0!t{|YID zQzpv3a%e(SMpZm+xO4H;_ZKx;gigA!v0?`Pt`7_0q#>g(z-l8fD1(GB?!LT{B{>Kn z6N@@T569LD5_qMEo^NA^7l*AQzmxUKMWj>?7+cx)1rHWw+DS{Ab7R6k+((_+N&9NJ zW-K2YD_prlYdr&;`xqM5$f!C&aBY*MtjcJ@dA%q6yK zdz82sDCU@d82!CZ@!Cfr6{I=R$s+`w*ge$|CZ#m9H|WCJ)u;y@^>WE0vrB0U=iYh> zTwd>-4(W=B^YnM^HZ=g`21j@GuszDiem2sRb8&pn?jBRmIFlj~n!2XS$x@8#c(y>$ z(vNyOxQ(KmWq&N?;s)s%~`I3EqMW4SPVL!=)lGYu;Jl)?tRC}xBJu1fPj zW~G-qk4UAZ>m5Y1X|Xp!p;@{KV|s=jsr%j)!IMXV5P>jVTh>nIE9aMYWPTZkG+5y82pz$n)!EMVMy9SGjbU@ebu8!Mw2Vt z*69%CBA?K23Y??5{_0SI`_UENOC$M^l*(rwMk3XA25C7k+MVYUAHRDgWR8lFoi)Vn z5@2C%XI#vxL<(VcMRo~x+4=ZcL@>E|1wAg;b#R{#0-)-L+D09w@!JzZ!FLHfrar^W zbvsM!5_W%qX_e3_7di|@8}4cSu8l#r#VpNl`T(`b&nIr@nsQvIscHj9iU|#8BN}@j zvC!ay8Q@RP^4g^X+2gxI&&QKAH{~*St(b{DyMF`74%}z(441MK7s)2X-&j&>8ExTr zK$9uo9QeU|loEXc>yYBwj8-d;h(R_Kb)-Stqe`|oe~myUl7y1P@ux1Jl=Kh*FA|;I zGO91Bf%l8+YQZx{aCFY4$6sy9 zpoNJu$&fN&XN;vxn3(_p-;6Kol~4YAw9j9*!l6zqj4L{!2mks4gK?twW>7g@<%>UZ zy5Ta_*Te|0`Aw@mZ+)vk>fBi^1%h}C%8{&M#d8o=9F{?N@V!~EntWfGDqKC-h@uV;aiA7{Z%jmM0Ub!rE(>X5q!<#- zrLQ4IT$LnDu#rBPHr_J1Ankyr3wKwe#?d}hpK$TBTIC7f*GMG|ucAW^sbKOqR=ZBW zLjZ>DL*2q~Wyf|=#b%-h!gfPWHGXUle(~$4rp32)N!b{M^^fqzaD8B|{8SI$bEknm zML0V5+T^RB5wS=}mwS*d3_&loQ!C1-z#F2t^ve4{ZF|P(yT(T{GH8*TgyZ2YF!B0lm&hSxcVZ@J&BRfp&AX_A^tS z)lLGoi1#862p(6SBzxX=X=4a-K?}rbt=Ur*k1xNM0m4kMURdNu+%9_Fv6;Mgu@n!O zA1pOICQMHDRgnDm93-MtBsobqj$0<>+s3;k4gmSU#5on}@jhkQ^#`z?Ji7K=9Trxg z*~bxUgus4?ih_07nS9zvU2Z!g5vFRJXNe~p(*z7P5ffo4hrDE2n|~i{!Izs0#X_BC zd`I|u>#f*~1>{WpLTYjrc`;5!)`dME32H1wYPMakOaGj@p(QEm=$6q{-BLsX_iKkq z;Pnu=^c{4uw(HC2WQO%%tgy;ranQ4_b54g6C#5;6SLO|*yYt<&o5MnRFrg{M$!0%L zlvj_y3^utLPsyh2CQhQwN6cxBb5I0bLX?Xfl_cIr;xm6xClHk|lFm((ZiV3M$~wz8 z?v&_vnvD39feL=>xLw+~-o9!F{C9Bw?Jx=O7Q(_I@HsW_A2Wq2OmAfJ@L;!s0SF3fFpC(HH(V&f;C%$idkXBmU_ybc00p1|(W;Pi~ z3+|Q^(eecS*AU=yJvS60iW^Lwg?|Y*1-!IpOc`>0cwZKC{VeM1?rvMkMVLg%=W|F0 zc(YGE4mWc1>iJQEf(fEC#kp^{D%D$i_HBQ_*)Bl3~y1HQGFl#5z7dCN?@Our;f=r7A~>EO_{3t5zl zXC!1dX23HcLi6EamNhM7Wyc>LhvH_6XWzlj)qlpyhzd>1O@Y(YpTf~HH2Fr*OS`i= ziZle2-Cs76pb+uI_*d*)+QUbW;}U?x@Z;$)@q7|^A2wMz5n?m<)k_o6T={|&D<`ND zu=qiSbVG@K;;5Vt0$d?Kky)B6CJnQ(Pnw}VOUb^fKF23Td05mjP9e`Y`=3BsGDa7# zFXxwK+Tm?8z>04Z4{y_dt54*PyK?Bu&E=C!3#eQ^j{T%`VGA=9bd0#ocv{t5odsdy z_zB@B7|GTIHJ@-@V3dxk{`zDM^-Qz!_r`VQ(WvKSk4yaP=Xk<(_$pybog28+6xFI`^r`$)LPM1?fRuUcbh%Q zADd-p>%w-(__|z%5&A8JGH-lRIvVJ?q<^_hN^4~=JG9p{e8g38x7l;4Z!OpglS$1u?$vCe97^i|;g@x;We6^8r2L^@8XRT*C04X8dG{38LAvtTMj}O{QI@eU^=T#Q}xqW z@{FUInK@?cTk49ya{s5UgGlX|+P0zRt5tfG_M!aJU+gtt5_mYP*47+sxEOh%1W`A2E&P*-_3LOFLyLchxU9UVPp@c&bnSz88gZZ$jxrzL?4bU&oKE*PNu07f}3Fg;v@VQ6YYD z?`nO>dH_$ttid^oxOB_{R!qP46r7<4`tJnU^&?TEY~`H*^rSu)2>e8Q(Z?=Zjvu5>w;|BgfbB1s03 zZHD;N(Y4wr;S+k_KMc^%f>|HYX_Bm0@i*F#C?S{tM_TUT6Hg=evLP0bN5>wAoO2r% z3>t}i&rADKe!ligrltY zqwM+QeV__D0}-Pv?~G2KE-x zKL#b(Wf49`z}KZp6rw{1QDviQbAZ|VEh;jP94#-9kwR7CB{%6I3uxOxh)2cuyrWh z*zW#0xjCzLJ$OrR$8%a(9Q9S;jv|K_yu;Hpz!T(5*FHmJYeT$fL=F0j}8 zW|G{RqF>6!+DZSuxk($J#Nn%|>(!DY{TD`7rY8ab^-Z4JM|>%v#6nMl!!V2Vd`JCbsYJ#H2!B#K?~>GQZqF=ssm; zz{S|f|37HVTMUEqC;GK!b67cD+nTdIRp`*vNi_ATQKndLb4uEBv@J0l*Pb^SHi#V5ub~w9lIWcbmjM2&VRd@Rx zBjD(aHV@eY+)TZHCwqiRq?C4Kb6JILrTJU$BIRAfvc)Kxn^gyO=d@0oL7^9$9)gJ( z%sFz*06$;#Fs?kjF5WRZx&Oa*n^Xw=3^IPvbZ1BN_Aqo1P=e%bZyHQgl{w+OnJ|@` z6S?nT1&ay8&JNp3X2Kx0UcwED^dt4b?6kI3ImTpHq;vuO9h>931JUWO5Lp_F#f|F$^CrGW zI6BtdOD=7#y*vbu;F)&Ul*h&gk@E6PFK)#lLxYEK+Q2&)*J6gpe{UK?Z!(@$e+P1S0)1d_shKa!^6TI}SCuvmHCbQQ61dUJLyfQdmP;z+D8#53Zp z&?i!WH?ZU_+01=J-@?%6$?|Fnd$y{$iB`p1ocIs>B6VO^l|8FQUK1!|$ZWZ?Ps z3EJ(Y#87Ia?v?Eh?n0JA=`h@pY&AnUie|z+2lt6lEH+KElQIxIG&G?zONoX$+1AsV zIYc~B^wLP2@MC1(PN_>Pcl7_4AN@NQpDwCTy@Q9r)a4#0WiMEGcM*gWU(WRsa8D3vcBlH>)@W_C|QnmA6B6Ai@GKwQr=SvAK7< zz=s|uT(%drez+E>ys@>-no3f{qcY;IPZAlvP10dvT0b_Y84WFspC1?YXbmI=vk~h} zoiA=JVV44!2P&?OUqz!P>&_o=u^qWO*AI3NJ8N&b@&6Hl>wn*1t#65_dXO~RE^Rq0 zXJWb5)(Y>PB4{v?$JO{>z_xDGegvFFD57V&E?dr?{hhuk)*QLy=YM8c5OZ$; zyWVC~GC7M#xmfl=^j*{@7V?g*Ezwj7NM4h&eKYzPj`n33tDI7fn~Q`am`iXs^36Gl z=R$R?V9>X9&s&t4$ZF}OSMXQMPvF0L(7a3t_2{Bu?C|}T*rkqw1aLJEp0s!Pg5|QYx1{ySKoJp+4eMWJ6RJxcWu3sVKw=OIZId>Yj*9SC^m|E*y7_ zhyt151lkjhoslz-TY~5|4&lLpz;Ucd%G_RPRE}RkC2PU&r8T0Qn@GXX0C<1nqKmpd zFu6DvKf*HMx851>FfPUuPkij53-a4lKHF#7jRFds#=) zKBO(F{$`mrZy}7UOSA|ADQ}X&$5%Z{5~K+*`reVDB&n4Q5+mT#9%v{9u!2Olly3@S z*gp_T))Ww2p>t*dE_zTC%VJd-px{0t_)hRQcE2v2z zfCP6=_N6t2)U28SD4g0k6`hGB!Qy`2$+5gtm}3#xeUBn4$)N4*V8T5<2`zt&b(8i% zy%ObtIPos7u8!~B28tfTN?pWTvLPyRBMpcwXgqIWQXK=iVQUea=Ok?an_CkhvfjqRge^Z2b(G8-l42@4>7@5~GmZf{~=iKAnaQ-1*U{wb}}>NX@u2p?=MAVu?1d%Oi7(DPJE- zdr!+zs7y$iAg&d-+V-`Je;=P59+Q#EGt6}2J`GiZXBz;2>jMapXQ@v5!LXYc~#v^p_L7-RtUm|=tWW#H#SmlR>PjSR>`Fso}h z{GxFaXJiuOv{)&JA?zE@#E>ma0N3(uzYP;dxj-5GMLL@uAvfvQ#u_1bfe6l$=$Fg9 zSj^GN136rp88FAkppuhTU33QWc%c3!WOKZaTza#HP6>G1(|W#=sAgS2!!A6wZ3je% zg)Y$4ROXcR*nU?u&Hoq7IY@?zpn9BwGGL zS6xPeXoIji27cgpB$zW#ftt|SvmZ5sN(X!SKQZxKrvYNUbG^7r0tfCaQyxUtoGcqZ zlv)C#zF%8-sO@1mJMQrAgdO#LF~#F_vO6cUmPUX>|LG=5mNf)QwWNDANy)YE2qhd4 zu=&c{oNS>((CKvdpgFvLRgf!<`E(h=$p2yHha9oO0@cb*B~qRDPwA6fM9jb33^bf2 ze0rzvfbdG}Hb7sutq?8XC8xDXP{c<&JQ^N3AOu3Otbckm0DMCd5UeB{#CPn@%l&T4Ao*wKE(BK@zl$O zU5y5IC^1)WujIq6Dx*Tw`lVR-t*5IS6&#fMx(dT#CJtx9ekWN32&J5{`4^DyG+pL{Oa}3OXHc=# zxP1jeWDk$dn}~twu7z;8KV278oFmX(n~m)>tb*p#Z6>&usi}h5DKCgt3y)uXhsJD# zSap9(dR2P)RB4ZBES`A^lDz|1)(#TrZdJ->6lJ>j$yv)qmeqsYyw=>KeCyXRGa5pI z_x3{5o6W?pR2xmo_BrK;2xPP8WZy0c$p8dseHG8v{isW8E6 z?VZs*8JX)nr=Sx^TqfagXa$T7^K0m-MxbKbh0HJD~0h1(dR2o9cvVQLS-bM4$~K!jkFALAn|~Sm4g&M#)}C4M}MF2+x|1 z`5|ypQ`~L(fg4X^w4p`HTL}wXm8}k-GG3&;ED$^)2puGtry5DsS6We(i-~t&lv?_3 z^}ss+UtKE5Iws*2VF8I>s9Aqwt=gtOPgH}X!@}G=6I*(I+H+;!Mn>%pbt`JWUhRFl zLe~GtLv5(&o|g!#@@~eod7v3``YRI$kawV?Q^E2)vrZBO|rA3(YREw1< zD##Kwwdk9Y#W&OG%Bh8=Fq$`udqwwQB-OM|QmcY|zIrV~Q&J~C5cLeq**tTB%ma1r zfv%h+Y-MUCR8X_%t-72zJSQSiJjn^gjvUCBbmS1d<4G?>0&v`mvQI@iY!?ant5 zkzMbwq^jh8Je12-;ThVU;xMx}<;8o_=IgaleYR0t9#n0kU%h``hfa4*W7K!B@&xQ+ zNmDDs$}izI!coSRIlYbi`1J8caqmrKU&Eq06ZJ8}H>E^O=B(+*x~mkMs9t2fw|6K5 z^JAFvs|8)ARGlA{iq{baeh3Svv~Q>r!pd^31nn04r#~&e7B`|Ba*(0*Zlm4*3iqC9 zH{5#Dm2f=5VhVdAE=i$2d*1$6Q32y?Q6W9hW|IjFZ?N(1A-N;&d#qqlV@@x273~&_H>~R59bjq zrZK7LHSf4e4M?AR2(Tv!IkQ4&eFEmBfl5}hXC&cWGY0ghd$Iz~ADmwKR0qBEkqr8zf})GF2!sgUVyZ=%2!loYnfFF0SRX6EF;*9DT=y zykjbv25~se>os=GJkw3o2YCmWA(kEewSmn0bCEE(XuY`OwCC z!~X>hHo1TjL!KZDCqch~c^`+?ElxkG%i;p8x)j{Gf^piz)XuGg}NZ zrrMRtsE~QGw@ASYl$pnv>w}yrmmf-Or(=YG;m@@T2RWd##wI&>kTUZgb1iI}GOhqC zJgP9)e!qqbeE2J8n(R@4`LG`;uer&M+9UlOQvR3kWwij@Y@+3V@Zayi3ahYe$ptyq z|2Xk?Oh^u0_JD5G{9&xsHyP{cu>R^tnDs`MlFM6bGeV!qvVPX;uRa~*^dC=Lw|QXr zVG1#<8cgOcJq98lm0!jV%LL&Q{}QR`P=$Fk@mC%pdJu+- z2gFlVp3i^z;(7ASUQ>GycOvN;LML!@!}@=Ko12;>fr+KvKp_+GJ+@!GfAzmnX2aG? ze|Q}WN<0K5E&iy%G;y3Ol#`@rMJfAEoj^0V*{Ww;1a}FIW40XCB%Rr8Tm|N?_d05~ ze4tx5#z0FXNNF+L-tmI4wbuo|pnT!B^M6FXtk2t*aseu9{SnqbjGsQ%V9@mbC=)hpUcU@x z9!Mkxq}43|Ph5J5%3<4EQNLkr$Gpj4M3!jiQ{btorgdZTSJhkn+|2 zEw}ym$f2cx;S)hJ!*+lOmM`E_I1{&Rt}N)r%~cPu@+BaQ{;3Po^3bhSO?B%I`}_M} zz8i=UE(=By6_2^pcq|Cg7vs);-M(cl#S~7!9e)l)0g~@m) zx`Gwn8ooj%s~=JqE>rHSSZ^ZDP0xDG)Un5gM4pL(|3(Ej{lJR;9s4HdJ@omTfeG*g znu3G*T8Y0%SSGIB)ch1vbN?|JDGC8&Qea#6)(iV&D~PeiWXiv7F^=q7Q^t?g2+aE} zOY)_`44I@@&J4JCd~}^nE_?Lpa*0D4H`mIMP`SA$+f5-Zd}3;@r!+S=Qj@iNL}0qC zktSbvZ+x2>J|{^LlC^}YVxc0Z7Nn}nc*t}V?9>!AZjZG^ri^bup|1cGQ%^i z^6UArGE1ZfBfDVEWp>ws9LQ$H{SrQe966ylNL(aMh@+-nuIKUmO>;aa zLliu=REd8-!=R~0jN#oR=si-VMEvKq_2gDGHMoPw;HSNe@hJ^_GkT0ZJ;)VU+eg=o z`kU!~lL3>W!g4XRULF5WWj}&E*9jK5p6hWJO3|Ci(O7|x{1gKH_w6Q`e`MpgGuMUy z)#1!U6G&iEm#z8anKt|a=3nQc{~!PY(d#Jb(?*o|?~ys8v;O&)3ub@gKT%u@jhHhI zxJD7#N?5l}kefM)SgGQ_5F`)UgylLS;TGQXOcsx|m;VEI2#2RD_F3+te^EHz@5a5I zUwexrxsD2^fE}nXmxNo}*hYdvUw3>{85?FHZ==8cDiNtEYl0IGIo4BMzc~SmVf^N%G7R2s^jj<@^eVVl_(f8O4hE?rTef9%A>$S7d=lbj(J4N*+g5Ka zAJY3EpF-$$#V1u&b$6p&ERp0_l8bJFMML7#1q|d?-$d<{9^X_qLBiV$qZC zrZ(|l)PCm++wW$>CVOlsm-4IVskvR8<^u(M(EnNR7N5DDYs~|dpr0+~LnkW(`;M!e z^`~YReC64$MS*Go^5N2Mwjoxov$^8Jn|rlCTMu}%89Pjnnhd|*1&8llDl(;AqWNTCYZc{q zlV?DQf&2xBh`DZQqKNIOpX?kBgg2 z*_-T9EHR*92jGQDGuwqbFK+J?(9PxaBHCqnY?5~)*My`UW&2>qi}dRS9>^^f)nx~i zTkoIHe8SAEwuxyelvLw?{#KC}tI`aU*uK;66f-Pq`bl}Pr8_PFq$NjLGy>aZlYz`8 z(LV!Sex!{_-zwlugL!Ypt$w8651t20J|v2R=PfB-99X+;H(x|E2i7@;-s1FhJBmGx z8&)y>AiU-~o1oHQ_R2E5^t5U8U%wFM7ab);px2hNTkaFEQn8UXp?@8O*9<_ptUt*Y zDEfF2$-Z4U@e|A#99j-=?7%+_{|nNCDbdi?^pL<2JcqU5cFhs=FcPbObp4wh-zSxa zTV1lo2`d)|hCh3Pic6%O^G5EdxXZSG^O^u&1W`dkba%~6xE#hL zMQ-qd6&d`q`F{W_5>yThnWr4fM9poK!}TWv+g1Ms$1iL`nfEZ-<+4VOyUicV?1sL_ z?^V2ZwTb5N7me(s%wQ6_yaktECSw-ZpDP4th>evU>ts@Lo8I_=X`*tMn<)$3aygls zbp!&N+tp4=bfW2m%7HVArCpvr_9e`^1m(^8Bk1-Azi2z=iN1a8Xqy6Bv_15|Dn??6qV1cvef8(idE;>?rTr$6<8$}py*Gk;-^Gau_N z3mtwGl@3E6Vd8y`v(Rl!VZrK++skdj?2_JDu(iPH?n9G{;fkB4Q6onu8>`f0--yb# zmM3har9K1aR_-op#ZTTQ>oIB57-4>O^q1^-xfW|ltTaz=RSG!gFS?opo4SzVj@m?+ zR}jsqYb5@=$xQ-Lm*vwH8-+Cx{>n$!b$AwWO{I#=sfF&%qTE1vKEY51T{UKcH9kQh z-Hj-wRWSLu8cPbE@b^$Ks+ zN5HgXeX0k4sS%Xtd>OL35S647S>Z20Pl>f+TrD^m_*}ZvmhakU6 znHeHFLsS_`OH`4RU&tlLMM}IQYolQ>Qw~(YNzL<`oISfYA>a*k56%g>drO1jbnvS{ zhVRqf(`1p90c5UE;t3>84biGYp7HX(Pw@=ua0%1y2}Ak~n(E~91&H>dR0(5V9_V)R z^18sL8&Qycv)@wH=%;|TCoQ?8X`Y0lqm<^}T{=eVf@=|xv-vAVo@Le{eZ?LZmz+{` zxi)q6K!8L~yJ$O!PKxqfz*3gr*c{7+a{0w5te*F~IGrD>ugcgv&rI=)kuTy(SR3df z+-?CM`+j{#S?d-L4vA?AZt{t-9k}sIC@tM3g8dz!cW$!_hg#Itrljk$P)y4~m>qlRZx{{Gl1iZMy+ z*%wdm`^s7M#kzC#-;>^d${;UBqQ=?U7hl6~7opWi3%7}}LK-H&_(e4kP!8kfAVDTt6ZnKKnl)Y(D&7NUr{klAUl@rQS0MPr=T1Qv?faRCGg*#CP~xR z8$vLT9O_HG+-_1M!W^uV76za{pFmFdZdmSXzSAX_uvH#gB1nAVky6&FT z2j{XY-CHWaj;Qa6)QPoH;ty%V|>IaBxB%X}&VEdYp5DY1*2 zpVJ-wTs3JO>)-mn*Y*X z1qCvtqIZ>XK}|iSzc*{_Vs6y#SZFGi!zTXf-G~#R)3;3CXKm{DE&dRN?|^13tgc zh-<0nG!|#4EGPdmllikE%%|(lmZ)MmoW3WRhX_?H*>~BBi&jRJL}C*hCrcX;qgZC` zP^Ng>snwrdrdDuLgx;t1+f}4yKv?l-*X_gztd_u@Xc>F>pN*;c;3S$q?4`wYcMt#b z;30(&OWd;zxk)#+5s`Ywn!-M-N%P6uf2Ki5pgKZB)wani<<~f8ZC<3#u}<(zleiNX z{%d5;YDLenW#_5}{Lw!KCoN5P`IT9ww>O4*KOde7*bh^P?ihN+>H1!!o%%*8=n@fc_$u_A*xmHW4moe9KL|gxi`h005nj) zT|lAt{Sj$UQ6VpH)MXLv{0qI%B>BfinoyeD;cBzmIka2FxX9MC#fH9Sq6Hz!Qr4Pl z)L@fDjuAcxB^%|8Ns(=!;-x=bBkRPYyY)TJVICL#^YM8p#16_Gbj`d_?ngt&KuX*)I+H;+EvJE z4vW*KWk>XVnr``oh-on+o)+W)%F>m-^0x|%dOe$_jzI}!X2+xTZgxk7N}AC|w&`O! zP)!aq0wnf57p~U)zskGvsHU!N-`1+NR;W{PC{RUYFp7eJfC*ZPKv6+aGz_9bM3689 zgaQH%C|b&(KoA+CGC45?2$K*i$dCYGQVFxnND>4>0%YKwo8$)iz4!g`)_QBb_160{ z_ny1&KKq<|&e^}ce>-+fC_KGD_t*i9cvtY6T>Cu80-cQf+tHPVyyNd9ZoECD?ZrV} zW!2#YRxA{Ih4#61b3us#Ws|A1=X@wo)D-IWD*h~%)7EC~OQTppv0>3t(+9;tr}Y(X zd1cN+aXqMuzmQVB_g&mo<~;Pg0z^?ETfF0URh&#o@1VH1;9kH2_XUG!_s53&3YYcX zL%~Dn)kiTFXnWTQ1DZe~_sn>&G1f>oWIAlqUMs5SO- zsiU7Zl1$*bG}ErRTFwG2*%dM~$RueANdn`E(u#aa45GV!%rc=5|1n}J^m?jX-g zbb$DvMjasC9sf|NvXVPjaGEal(T}5Btc;RV3q0lruMZ$pMCS;QYI3G?d$!ikdkk*c z2G~E9y8l;z`~Q1lWQ-OaJ;Lu};YdG70NSh*_z~=8+36htC z+G&R;kJ-3wgYlA7L=k;mwpEnIoKq`Ou{>!voF@$2eLLH?Lt!UV*w2UmYI%7_)x29W z;O7dLcL86E?W%t?oBi~On(be(s~eAD9Vl^biAeYwERbUC`V+=%B3EEU(3u(54SH0v zmwepIn6yKpFbFr$4wU#2LiCZo1|=@--iQCysYhut5`KwIcZaO44N@)+l%Nzh<3g21 zIhP`RJki8S-o+;GJ6kwk!w^5zQ15D86LOp5QZGNn2DQP-?|q6Lal&GzBQQ@4b?QjzT`TaA|1MumKr!_c_l`{Lm-jE~)~3kT zwr^N`AKTMx1k|PYyPMCjL{M)ty%Okn3*6Qfk`#e0*7zPTTttCBdf#hAen-7bol{XH zaTRR*SBSC4aYm@D``zcZk@xPyE9-Guj;Qi}OoZ{4f8YU|ZWHbLHIDMmZl@5Pp{1%# z3FUnJcNRT^|*k>#EC_J_{`BYwcXB7n9;$k*noKeu;b+D?3oA2KV6^g8AwG z^9kPm_#@6^HHBIuHLh;o!;@Bw5U6nT*;(dzAdCMpdWSAD&Cxd5L zKDRKM{w5SX2wq$&ov9V1TtTWy8eEEBZw-tKIY^UEGKt z!X{qHGcsgEawi6i(C*_~&7D2q=*{a5dvaUM7uJlkY6cKV(UTpSi#j^NU{gMRLL z&h?oi_uOP)DZWr5idOJ&&P~sIj`IyVFjo)+HXs7712OZ~eg)Ywld)qOOl^rxz&)=W z=27fD3mlkV1?n`;#mlGbKFZ?HB5z_NIz>m6iU{Ny#hcW)}{9 z9sq$GddW7*&Hq2_$$yjC5`o5-$vA}cy-DYB9lZPAILjzwvWf{lQDpVS zJ70)D=PuuZlTl@;CJ$+RPk#!p=f#MIF~Xs!%ur-Hm2bSwHdgZWt%x7eOzch8Ht0RL zf@dnR4gm8Fu$&`Wzq;?>+=nF;AQW@y-ElC?sn3{ zcBLZV1J5f2AZ1)whLy{eg_9pHfJk!k-qCduV;J7BqM$Ntk>1uy zmQUHkAULaFjcT~zq|%m<{a_ImvFJe6?trhjJ_PH*aNl(H%^adE9Ao(!V@A;Qxo@${ zbC8m|14B=R4uD`cq<_~n<%Ro5P>4@M1d3cpP$pI(F>i7+q4-6v(gs9vbaM|ZM4Nij z*R|F!^ybze8`;1_iMc4$DdWVv6Z=VIx+riV@>hp2b=fAbsVdp$#Ssw5!voAiTzo_O zS^oo2-7ydr5ncRnu(7`tPw@@?hC(zE?Y;i!x3c6jP13?q#s36(E4nJrT2RlGWL*a{HRkRP*fXeBTQaxoy{pX+^QE zO2waki0+d8Tq`G#jm3b1-XN9_tYcxI`}gH{$H1<%9>hOHsOO5rU5b3dPjW^P#B zw1x_xT*`~wEmiY*#`~!Gv^!BYC3FXM3h^VPcw02_=LBp)5}xp_tIXf^iib?l9DT?s z1XIB$b;yRWVgCiRr=cp`9Av{Z$Da318}hV`Sgy;?gr%!f^U5CFD&~100~qasFb}uH zyf2P`G%CUdDk(TduHDMm`Pa1ffEjvKbh^())UfI%?nQ@JVgq6A@8E~po5yKH<$xsL z+q9N2xNR|_m6R4{-DUi-tG(Pw71j}CKn;Jnu`HfNC}q59yRdmY2Y%xfd^%v2jzbgX zA4J6$4~TpLK`^Em#CJHhlRP~c6~ea=_Iq4&b&vBuau{OWKU1enoe4Q_nABFY@OQ6k zmL|GgS8ueC+Q)Mvg%7z&hl%@0FnR));ia8hN6uCFv{cSTPW zjuZKnVGLczU6W@%&GZU_uO(PJ#`#zOCVGv0{16E41#4P`u&oGFtq<>~Z2@A{#bQZ_ zl)G*lg~LT|Z(l)|v$63*+j(0+sfp1G9-??&Az-K_2+n5&s~{WJ%e<(6kw-Z=(1Alo zl0dd>&tYA-Su2qr)Xsi5uSNKg@0{2IDR6`|M673d(@So-Q^eniYbM6!&NIHXj^+^A z*1QAiLM_P#uTec>VB*|;a6)H~BrB*4LG7i}-tqG+M#1|WJ5QonhMr8J0H;{Woghyg z`lG|U5lXrSNI^Vi&Wcy^aDlNt0ElZRox%q)S_AFJymz>X$);t9a9MM7v28vOlnPZv zXp^B%aXDA^k;g@kjDkl}!$wM|Bc_sFt+0~Fzb&w-zd7N6NomH;e8>7zl4H&%xi7*?g0F0(DvL!YfVTfHZWFgQc;PVtcn-o*d0~45t#-Jf?NP}6^Crq`V z^0$U~H@gr{4{w$iRcg>VgJ&NqhB3z9yE*dA?*$kONp84uV z9lYdoqYVB1vS+;{KEhl5XB-Ap@7N}jm8>mdd~rOxmA)rKCGrwopB{qpUTJ$bh?ppQaFJibp6liBUt zc7&`GZW1JZx%Amc*f^*W^y`m!mzQ9UDv77o3d)o&XV=4{Z)0kN>O#E--P;J>r9m-n z(V=0n4bZnEZMQC}S9o}K+M?4|_ZWaE8LE_!FRQG)CehlU=XRo=Zy_TT){{eX%Z)j) zEE$1O;O9)F3e#_#AR0tb?T%-Y)O8BOs5$vwzdBaNFJ(px1Bx9gaY40_sr({CkKbcY zK_8;DsX-H?v`T6FalbHOjYX^D>w5)I2X~4(7y(Lf_Rs4IHOrbwESY115#N{{ac7|E zVkDZx&D3?W*}=6A+6E@z`TAZ#jZZt(X^I!ZcQ&j&l5-TSU}r_csqG6fEC@U*`%a7% zn4`_T3wvQ=CbxaR3`Sir26?i(--%;3^94#G(b1b zaJ7})0?QI37K#HxvGZMd^cZBuFiYlrsAzzd;+0brh4vvE({)`Edg`wIhsd(j`X(gNTW^V1JGj9My6RMyR0hlwvJQX2sXSo?wp&h`r(E9Qa zV#(>VE6V+2o`*;tcf=1FaSRtK*783Ew@1tafE6ma0>~f0RE&#N_=%myGY}FBeSB3! z5L6j^cQ=bOS)VZvzv|3S1*O_Wp{qfH`;HN)@)sJ^;)g4rq@6A z0yX$mRtNkt8#Mp(V@{pCp5Zq6_Yf8gZ9~nVlA`r+SBIK}iv%q&fDjwZ_jTV-Oq}2* zqh5R`HmW-Kmh2lGW`B(nck?V{JByIhsf=K(jHn4n5cl=V91@X$AVZjfla8KUHVf{8 z6`nEHu7%n*&jX}C-vx&o2luVE77rJ+eJ`n21H6ft)uPux+wFg+x=l5L zY9E$k=-JmI;N|Ff_GyKY!GCGw!l{PO@+i2$>RmLVg(-u~`Jg0p-)#F(a)7H=qh=}Q zrUxbkwT-A=SPh@|H(?Sw=OWZnXtqQ@kPtD<080<#WN|6 zK!;9;<2yg{?3VaF`3SDVrSOtg!`7o%T}KvMZOOeHgkO0t9(Sf|!YWGM0B4ir(-FO_ zco8LLv9Es%n2u4^E^>2JO)37FWrbJDA5iOdEG(+N_i=XVofG&}Ir{c&?{%Z5|FsJ< zs$hAKXU3g!;XPE z-~wgV*`vG;a-PF}R2sqr`=Pr17CM9B2&-?3mg#&IaiQoux)vo(!-`cdd(IgN>c6|R z9jEh^dwa6GTs{r{=GQU~szdc68cO%{iIA~x34b$eqUrB-WulWZ4rx;3as9miDvBm9fZ?qXVL|C<}; z&MH|WU7fw+;}&oTK;_wVT}`J4!jmHnQInnPpUN0Kz8;kSy`|sbMyxMOYL?b&ISgK3(n)cv8END^ zc@^$4y5K^%!fGgC!uTjawH^JGS9)oby+6KBQFv`CcjUHWD!~oC>ACuDj-I+)ip@R&r7{<8lqBvhkKe3ss8OH_ zrmrZv1>!BOB)i^o#IAsajUP2$H*_UW{z845Ybi-)-GWIdUHGuR;MsHlwyWrWzdDYiQ67pt;Y!7UYamw|n~pGv32eH`wkKbf=HTxFD+H>{p|`vk z3&ZNzRD?;V$@Ui&o!m~TtRrZ5QdCN+hiSHyVS|lL_djc@NKkf^O4|y+yFXCiinLyu zt=>)B+W)IiFu-RP^evZxdU`wu%S@1S&*Y~pY`*`V;~-Q!xqt%|^T?!smf(Yc;EUjo k7r*+duWi~wShpa1U*%%(CH", "model_record": { @@ -57,7 +56,7 @@ { "identifier": "=C<", "model_record": { - "m": 0.86367, + "m": 0.86367, "sigma": 3.1815, "epsilon_k": 156.31 }, @@ -66,7 +65,7 @@ { "identifier": "C≡CH", "model_record": { - "m": 1.3279, + "m": 1.3279, "sigma": 2.9421, "epsilon_k": 223.05 }, @@ -75,7 +74,7 @@ { "identifier": "CH2_hex", "model_record": { - "m": 0.39496, + "m": 0.39496, "sigma": 3.9126, "epsilon_k": 289.03 }, @@ -84,7 +83,7 @@ { "identifier": "CH_hex", "model_record": { - "m": 0.02880, + "m": 0.02880, "sigma": 8.9779, "epsilon_k": 1306.7 }, @@ -93,7 +92,7 @@ { "identifier": "CH2_pent", "model_record": { - "m": 0.46742, + "m": 0.46742, "sigma": 3.7272, "epsilon_k": 267.16 }, @@ -102,7 +101,7 @@ { "identifier": "CH_pent", "model_record": { - "m": 0.03314, + "m": 0.03314, "sigma": 7.7190, "epsilon_k": 1297.7 }, @@ -111,7 +110,7 @@ { "identifier": "CH_arom", "model_record": { - "m": 0.42335, + "m": 0.42335, "sigma": 3.7270, "epsilon_k": 274.41 }, @@ -169,7 +168,7 @@ { "identifier": "HCOO", "model_record": { - "m": 1.7525, + "m": 1.7525, "sigma": 2.9043, "epsilon_k": 229.63, "mu": 2.7916 @@ -193,7 +192,9 @@ "sigma": 3.2859, "epsilon_k": 488.66, "epsilon_k_ab": 2517.0, - "kappa_ab": 0.006825 + "kappa_ab": 0.006825, + "na": 1.0, + "nb": 1.0 }, "molarweight": 17.00734 }, @@ -204,9 +205,10 @@ "sigma": 3.6456, "epsilon_k": 467.59, "epsilon_k_ab": 1064.6, - "kappa_ab": 0.026662 + "kappa_ab": 0.026662, + "na": 1.0, + "nb": 1.0 }, "molarweight": 16.02238 } -] - +] \ No newline at end of file diff --git a/parameters/pcsaft/sauer2014_homo_joback.json b/parameters/pcsaft/sauer2014_homo_joback.json index 534c9ae75..c56e8632b 100644 --- a/parameters/pcsaft/sauer2014_homo_joback.json +++ b/parameters/pcsaft/sauer2014_homo_joback.json @@ -332,7 +332,9 @@ "sigma": 3.2859, "epsilon_k": 488.66, "epsilon_k_ab": 2517.0, - "kappa_ab": 0.006825 + "kappa_ab": 0.006825, + "na": 1.0, + "nb": 1.0 }, "molarweight": 17.00734, "ideal_gas_record": { @@ -350,7 +352,9 @@ "sigma": 3.6456, "epsilon_k": 467.59, "epsilon_k_ab": 1064.6, - "kappa_ab": 0.026662 + "kappa_ab": 0.026662, + "na": 1.0, + "nb": 1.0 }, "molarweight": 16.02238, "ideal_gas_record": { diff --git a/src/association/dft.rs b/src/association/dft.rs index e824efb24..05c92b378 100644 --- a/src/association/dft.rs +++ b/src/association/dft.rs @@ -52,9 +52,6 @@ where // number of segments let n = self.association_parameters.component_index.len(); - // number of associating segments - let nassoc = self.association_parameters.assoc_comp.len(); - // number of dimensions let dim = (weighted_densities.shape()[0] - 1) / n - 1; @@ -70,7 +67,7 @@ where .collect(); let n3 = weighted_densities.index_axis(Axis(0), n * (dim + 1)); - // calculate rho0 (only associating segments) + // calculate rho0 let [_, _, c2, _] = p.geometry_coefficients(temperature); let diameter = p.hs_diameter(temperature); let mut n2i = n0i.to_owned(); @@ -88,9 +85,6 @@ where *rho0 = n0i; } }); - let rho0 = Array2::from_shape_fn((nassoc, n3.len()), |(i, j)| { - rho0[(self.association_parameters.assoc_comp[i], j)] - }); // calculate xi let n2v: Vec<_> = n2vi.iter().map(|n2vi| n2vi.sum_axis(Axis(0))).collect(); @@ -111,56 +105,143 @@ where // auxiliary variables let n3i = n3.mapv(|n3| (-n3 + 1.0).recip()); - // only one associating component - if nassoc == 1 { - // association strength - let k = &n2 * &n3i * diameter[self.association_parameters.assoc_comp[0]] * 0.5; - let deltarho = (((&k / 18.0 + 0.5) * &k * xi + 1.0) * n3i) - * ((temperature.recip() * self.association_parameters.epsilon_k_aibj[(0, 0)]) - .exp_m1() - * self.association_parameters.sigma3_kappa_aibj[(0, 0)]) - * rho0.index_axis(Axis(0), 0); - - let na = self.association_parameters.na[0]; - let nb = self.association_parameters.nb[0]; - let f = |x: N| x.ln() - x * 0.5 + 0.5; - if nb > 0.0 { - // no cross association, two association sites - let xa = deltarho.mapv(|d| Self::assoc_site_frac_ab(d, na, nb)); - let xb = (&xa - 1.0) * (na / nb) + 1.0; - Ok((xa.mapv(f) * na + xb.mapv(f) * nb) * rho0.index_axis(Axis(0), 0)) - } else { - // no cross association, one association site - let xa = deltarho.mapv(|d| Self::assoc_site_frac_a(d, na)); - - Ok(xa.mapv(f) * na * rho0.index_axis(Axis(0), 0)) + self.calculate_helmholtz_energy_density(temperature, &rho0, &n2, &n3i, &xi) + } +} + +impl Association

{ + pub fn calculate_helmholtz_energy_density< + N: DualNum + ScalarOperand, + S: Data, + >( + &self, + temperature: N, + rho0: &Array2, + n2: &ArrayBase, + n3i: &Array1, + xi: &Array1, + ) -> EosResult> { + let a = &self.association_parameters; + + let d = self.parameters.hs_diameter(temperature); + + match ( + a.sites_a.len() * a.sites_b.len(), + a.sites_c.len(), + self.force_cross_association, + ) { + (0, 0, _) => Ok(Array::zeros(n3i.len())), + (1, 0, false) => { + Ok(self.helmholtz_energy_density_ab_analytic(temperature, rho0, &d, n2, n3i, xi)) + } + (0, 1, false) => { + Ok(self.helmholtz_energy_density_cc_analytic(temperature, rho0, &d, n2, n3i, xi)) + } + (1, 1, false) => { + Ok( + self.helmholtz_energy_density_ab_analytic(temperature, rho0, &d, n2, n3i, xi) + + self.helmholtz_energy_density_cc_analytic( + temperature, + rho0, + &d, + n2, + n3i, + xi, + ), + ) + } + _ => { + let mut x: Array1 = Array::from_elem(a.sites_a.len() + a.sites_b.len(), 0.2); + let (assoc_comp_ab, n_ab): (Vec<_>, Vec<_>) = a + .sites_a + .iter() + .chain(a.sites_b.iter()) + .map(|s| (s.assoc_comp, s.n)) + .unzip(); + let rhoab = Array2::from_shape_fn((x.len(), n3i.len()), |(i, j)| { + rho0[(assoc_comp_ab[i], j)] * n_ab[i] + }); + rhoab + .axis_iter(Axis(1)) + .zip(n2.iter()) + .zip(n3i.iter()) + .zip(xi.iter()) + .map(|(((rho, &n2), &n3i), &xi)| { + let [delta_ab, delta_cc] = + self.association_strength(temperature, &d, n2, n3i, xi); + Self::helmholtz_energy_density_cross_association( + &rho, + &delta_ab, + &delta_cc, + self.max_iter, + self.tol, + Some(&mut x), + ) + }) + .collect() } - } else { - let mut x: Array1 = Array::from_elem(2 * nassoc, 0.2); - Ok(rho0 - .view() - .into_shape([nassoc, rho0.len() / nassoc]) - .unwrap() - .axis_iter(Axis(1)) - .zip(n2.iter()) - .zip(n3i.iter()) - .zip(xi.iter()) - .map(|(((rho0, &n2), &n3i), &xi)| { - self.helmholtz_energy_density_cross_association( - temperature, - &rho0, - &diameter, - n2, - n3i, - xi, - self.max_iter, - self.tol, - Some(&mut x), - ) - }) - .collect::, _>>()? - .into_shape(n2.raw_dim()) - .unwrap()) } } + + fn helmholtz_energy_density_ab_analytic + ScalarOperand, S: Data>( + &self, + temperature: N, + rho0: &Array2, + diameter: &Array1, + n2: &ArrayBase, + n3i: &Array1, + xi: &Array1, + ) -> Array1 { + let a = &self.association_parameters; + + // site densities + let i = a.sites_a[0].assoc_comp; + let j = a.sites_b[0].assoc_comp; + let rhoa = &rho0.index_axis(Axis(0), i) * a.sites_a[0].n; + let rhob = &rho0.index_axis(Axis(0), j) * a.sites_b[0].n; + + // association strength + let di = diameter[i]; + let dj = diameter[j]; + let k = n2 * n3i * (di * dj / (di + dj)); + let delta = (((&k / 18.0 + 0.5) * &k * xi + 1.0) * n3i) + * ((temperature.recip() * a.epsilon_k_ab[(0, 0)]).exp_m1() * a.sigma3_kappa_ab[(0, 0)]); + + // no cross association, two association sites + let aux = &delta * (&rhob - &rhoa) + 1.0; + let xa = ((&aux * &aux + &delta * &rhoa * 4.0).map(N::sqrt) + &aux).map(N::recip) * 2.0; + let aux = -aux + 2.0; + let xb = ((&aux * &aux + delta * &rhob * 4.0).map(N::sqrt) + aux).map(N::recip) * 2.0; + + let f = |x: N| x.ln() - x * 0.5 + 0.5; + rhoa * xa.mapv(f) + rhob * xb.mapv(f) + } + + fn helmholtz_energy_density_cc_analytic + ScalarOperand, S: Data>( + &self, + temperature: N, + rho0: &Array2, + diameter: &Array1, + n2: &ArrayBase, + n3i: &Array1, + xi: &Array1, + ) -> Array1 { + let a = &self.association_parameters; + + // site densities + let i = a.sites_c[0].assoc_comp; + let rhoc = &rho0.index_axis(Axis(0), i) * a.sites_c[0].n; + + // association strength + let di = diameter[i]; + let k = n2 * n3i * (di * 0.5); + let delta = (((&k / 18.0 + 0.5) * &k * xi + 1.0) * n3i) + * ((temperature.recip() * a.epsilon_k_cc[(0, 0)]).exp_m1() * a.sigma3_kappa_cc[(0, 0)]); + + // no cross association, two association sites + let xc = ((delta * 4.0 * &rhoc + 1.0).map(N::sqrt) + 1.0).map(N::recip) * 2.0; + + let f = |x: N| x.ln() - x * 0.5 + 0.5; + rhoc * xc.mapv(f) + } } diff --git a/src/association/mod.rs b/src/association/mod.rs index aa1b74bb7..cd59bf0ec 100644 --- a/src/association/mod.rs +++ b/src/association/mod.rs @@ -1,10 +1,11 @@ //! Generic implementation of the SAFT association contribution //! that can be used across models. use crate::hard_sphere::HardSphereProperties; -use feos_core::{EosError, HelmholtzEnergyDual, StateHD}; +use feos_core::{EosError, EosResult, HelmholtzEnergyDual, StateHD}; use ndarray::*; use num_dual::linalg::{norm, LU}; use num_dual::*; +use num_traits::Zero; use serde::{Deserialize, Serialize}; use std::fmt; use std::ops::SubAssign; @@ -17,6 +18,25 @@ mod python; #[cfg(feature = "python")] pub use python::PyAssociationRecord; +#[derive(Clone, Copy, Debug)] +struct AssociationSite { + assoc_comp: usize, + n: f64, + kappa_ab: f64, + epsilon_k_ab: f64, +} + +impl AssociationSite { + fn new(assoc_comp: usize, n: f64, kappa_ab: f64, epsilon_k_ab: f64) -> Self { + Self { + assoc_comp, + n, + kappa_ab, + epsilon_k_ab, + } + } +} + /// Pure component association parameters. #[derive(Serialize, Deserialize, Clone, Copy, Default)] pub struct AssociationRecord { @@ -25,20 +45,27 @@ pub struct AssociationRecord { /// Association energy parameter in units of Kelvin pub epsilon_k_ab: f64, /// \# of association sites of type A - #[serde(skip_serializing_if = "Option::is_none")] - pub na: Option, + #[serde(skip_serializing_if = "f64::is_zero")] + #[serde(default)] + pub na: f64, /// \# of association sites of type B - #[serde(skip_serializing_if = "Option::is_none")] - pub nb: Option, + #[serde(skip_serializing_if = "f64::is_zero")] + #[serde(default)] + pub nb: f64, + /// \# of association sites of type C + #[serde(skip_serializing_if = "f64::is_zero")] + #[serde(default)] + pub nc: f64, } impl AssociationRecord { - pub fn new(kappa_ab: f64, epsilon_k_ab: f64, na: Option, nb: Option) -> Self { + pub fn new(kappa_ab: f64, epsilon_k_ab: f64, na: f64, nb: f64, nc: f64) -> Self { Self { kappa_ab, epsilon_k_ab, na, nb, + nc, } } } @@ -47,8 +74,16 @@ impl fmt::Display for AssociationRecord { fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { write!(f, "AssociationRecord(kappa_ab={}", self.kappa_ab)?; write!(f, ", epsilon_k_ab={}", self.epsilon_k_ab)?; - write!(f, ", na={}", self.na.unwrap_or(1.0))?; - write!(f, ", nb={})", self.nb.unwrap_or(1.0)) + if self.na > 0.0 { + write!(f, ", na={}", self.na)?; + } + if self.nb > 0.0 { + write!(f, ", nb={})", self.nb)?; + } + if self.nc > 0.0 { + write!(f, ", nc={}", self.nc)?; + } + write!(f, ")") } } @@ -57,61 +92,88 @@ impl fmt::Display for AssociationRecord { #[derive(Clone)] pub struct AssociationParameters { component_index: Array1, - pub assoc_comp: Array1, - pub kappa_ab: Array1, - pub epsilon_k_ab: Array1, - pub sigma3_kappa_aibj: Array2, - pub epsilon_k_aibj: Array2, - pub na: Array1, - pub nb: Array1, + sites_a: Array1, + sites_b: Array1, + sites_c: Array1, + pub sigma3_kappa_ab: Array2, + pub sigma3_kappa_cc: Array2, + pub epsilon_k_ab: Array2, + pub epsilon_k_cc: Array2, } impl AssociationParameters { pub fn new( - records: &[Option], + records: &[Vec], sigma: &Array1, component_index: Option<&Array1>, ) -> Self { - let mut assoc_comp = Vec::new(); - let mut sigma_assoc = Vec::new(); - let mut kappa_ab = Vec::new(); - let mut epsilon_k_ab = Vec::new(); - let mut na = Vec::new(); - let mut nb = Vec::new(); + let mut sites_a = Vec::new(); + let mut sites_b = Vec::new(); + let mut sites_c = Vec::new(); for (i, record) in records.iter().enumerate() { - if let Some(record) = record.as_ref() { - if record.kappa_ab > 0.0 && record.epsilon_k_ab > 0.0 { - assoc_comp.push(i); - sigma_assoc.push(sigma[i]); - kappa_ab.push(record.kappa_ab); - epsilon_k_ab.push(record.epsilon_k_ab); - na.push(record.na.unwrap_or(1.0)); - nb.push(record.nb.unwrap_or(1.0)); + for site in record { + if site.kappa_ab > 0.0 && site.epsilon_k_ab > 0.0 { + if site.na > 0.0 { + sites_a.push(AssociationSite::new( + i, + site.na, + site.kappa_ab, + site.epsilon_k_ab, + )); + } + if site.nb > 0.0 { + sites_b.push(AssociationSite::new( + i, + site.nb, + site.kappa_ab, + site.epsilon_k_ab, + )); + } + if site.nc > 0.0 { + sites_c.push(AssociationSite::new( + i, + site.nc, + site.kappa_ab, + site.epsilon_k_ab, + )); + } } } } - let sigma3_kappa_aibj = Array2::from_shape_fn([kappa_ab.len(); 2], |(i, j)| { - (sigma_assoc[i] * sigma_assoc[j]).powf(1.5) * (kappa_ab[i] * kappa_ab[j]).sqrt() + let sigma3_kappa_ab = Array2::from_shape_fn([sites_a.len(), sites_b.len()], |(i, j)| { + (sigma[sites_a[i].assoc_comp] * sigma[sites_b[j].assoc_comp]).powf(1.5) + * (sites_a[i].kappa_ab * sites_b[j].kappa_ab).sqrt() + }); + let sigma3_kappa_cc = Array2::from_shape_fn([sites_c.len(); 2], |(i, j)| { + (sigma[sites_c[i].assoc_comp] * sigma[sites_c[j].assoc_comp]).powf(1.5) + * (sites_c[i].kappa_ab * sites_c[j].kappa_ab).sqrt() }); - let epsilon_k_aibj = Array2::from_shape_fn([epsilon_k_ab.len(); 2], |(i, j)| { - 0.5 * (epsilon_k_ab[i] + epsilon_k_ab[j]) + let epsilon_k_ab = Array2::from_shape_fn([sites_a.len(), sites_b.len()], |(i, j)| { + 0.5 * (sites_a[i].epsilon_k_ab + sites_b[j].epsilon_k_ab) + }); + let epsilon_k_cc = Array2::from_shape_fn([sites_c.len(); 2], |(i, j)| { + 0.5 * (sites_c[i].epsilon_k_ab + sites_c[j].epsilon_k_ab) }); Self { component_index: component_index .cloned() .unwrap_or_else(|| Array1::from_shape_fn(records.len(), |i| i)), - assoc_comp: Array1::from_vec(assoc_comp), - kappa_ab: Array1::from_vec(kappa_ab), - epsilon_k_ab: Array1::from_vec(epsilon_k_ab), - sigma3_kappa_aibj, - epsilon_k_aibj, - na: Array1::from_vec(na), - nb: Array1::from_vec(nb), + sites_a: Array1::from_vec(sites_a), + sites_b: Array1::from_vec(sites_b), + sites_c: Array1::from_vec(sites_c), + sigma3_kappa_ab, + sigma3_kappa_cc, + epsilon_k_ab, + epsilon_k_cc, } } + + pub fn is_empty(&self) -> bool { + (self.sites_a.is_empty() | self.sites_b.is_empty()) & self.sites_c.is_empty() + } } /// Implementation of the SAFT association Helmholtz energy @@ -158,17 +220,25 @@ impl Association

{ n2: D, n3i: D, xi: D, - ) -> Array2 { - // Calculate association strength - let ac = &self.association_parameters.assoc_comp; - Array2::from_shape_fn([ac.len(); 2], |(i, j)| { - let k = diameter[ac[i]] * diameter[ac[j]] / (diameter[ac[i]] + diameter[ac[j]]) - * (n2 * n3i); + ) -> [Array2; 2] { + let p = &self.association_parameters; + let delta_ab = Array2::from_shape_fn([p.sites_a.len(), p.sites_b.len()], |(i, j)| { + let di = diameter[p.sites_a[i].assoc_comp]; + let dj = diameter[p.sites_b[j].assoc_comp]; + let k = di * dj / (di + dj) * (n2 * n3i); n3i * (k * xi * (k / 18.0 + 0.5) + 1.0) - * self.association_parameters.sigma3_kappa_aibj[(i, j)] - * (temperature.recip() * self.association_parameters.epsilon_k_aibj[(i, j)]) - .exp_m1() - }) + * p.sigma3_kappa_ab[(i, j)] + * (temperature.recip() * p.epsilon_k_ab[(i, j)]).exp_m1() + }); + let delta_cc = Array2::from_shape_fn([p.sites_c.len(); 2], |(i, j)| { + let di = diameter[p.sites_c[i].assoc_comp]; + let dj = diameter[p.sites_c[j].assoc_comp]; + let k = di * dj / (di + dj) * (n2 * n3i); + n3i * (k * xi * (k / 18.0 + 0.5) + 1.0) + * p.sigma3_kappa_cc[(i, j)] + * (temperature.recip() * p.epsilon_k_cc[(i, j)]).exp_m1() + }); + [delta_ab, delta_cc] } } @@ -177,6 +247,7 @@ impl + ScalarOperand, P: HardSphereProperties> HelmholtzEnergyDu { fn helmholtz_energy(&self, state: &StateHD) -> D { let p: &P = &self.parameters; + let a = &self.association_parameters; // temperature dependent segment diameter let diameter = p.hs_diameter(state.temperature); @@ -186,48 +257,43 @@ impl + ScalarOperand, P: HardSphereProperties> HelmholtzEnergyDu let n2 = zeta2 * 6.0; let n3i = (-n3 + 1.0).recip(); - if self.association_parameters.assoc_comp.len() > 1 || self.force_cross_association { - // extract densities of associating segments - let rho_assoc = self - .association_parameters - .assoc_comp - .mapv(|a| state.partial_density[self.association_parameters.component_index[a]]); - - // Helmholtz energy - self.helmholtz_energy_density_cross_association( - state.temperature, - &rho_assoc, - &diameter, - n2, - n3i, - D::one(), - self.max_iter, - self.tol, - None, - ) - .unwrap_or_else(|_| D::from(std::f64::NAN)) - * state.volume - } else { - // association strength - let c = self.association_parameters.component_index - [self.association_parameters.assoc_comp[0]]; - let deltarho = - self.association_strength(state.temperature, &diameter, n2, n3i, D::one())[(0, 0)] - * state.partial_density[c]; - - let na = self.association_parameters.na[0]; - let nb = self.association_parameters.nb[0]; - if nb > 0.0 { - // no cross association, two association sites - let xa = Self::assoc_site_frac_ab(deltarho, na, nb); - let xb = (xa - 1.0) * (na / nb) + 1.0; - - state.moles[c] * ((xa.ln() - xa * 0.5 + 0.5) * na + (xb.ln() - xb * 0.5 + 0.5) * nb) - } else { - // no cross association, one association site - let xa = Self::assoc_site_frac_a(deltarho, na); - - state.moles[c] * (xa.ln() - xa * 0.5 + 0.5) * na + // association strength + let [delta_ab, delta_cc] = + self.association_strength(state.temperature, &diameter, n2, n3i, D::one()); + + match ( + a.sites_a.len() * a.sites_b.len(), + a.sites_c.len(), + self.force_cross_association, + ) { + (0, 0, _) => D::zero(), + (1, 0, false) => self.helmholtz_energy_ab_analytic(state, delta_ab[(0, 0)]), + (0, 1, false) => self.helmholtz_energy_cc_analytic(state, delta_cc[(0, 0)]), + (1, 1, false) => { + self.helmholtz_energy_ab_analytic(state, delta_ab[(0, 0)]) + + self.helmholtz_energy_cc_analytic(state, delta_cc[(0, 0)]) + } + _ => { + // extract site densities of associating segments + let rho: Array1<_> = a + .sites_a + .iter() + .chain(a.sites_b.iter()) + .chain(a.sites_c.iter()) + .map(|s| state.partial_density[a.component_index[s.assoc_comp]] * s.n) + .collect(); + + // Helmholtz energy + Self::helmholtz_energy_density_cross_association( + &rho, + &delta_ab, + &delta_cc, + self.max_iter, + self.tol, + None, + ) + .unwrap_or_else(|_| D::from(std::f64::NAN)) + * state.volume } } } @@ -240,66 +306,72 @@ impl

fmt::Display for Association

{ } impl Association

{ - pub fn assoc_site_frac_ab>(deltarho: D, na: f64, nb: f64) -> D { - (((deltarho * (na - nb) + 1.0).powi(2) + deltarho * nb * 4.0).sqrt() - + (deltarho * (nb - na) + 1.0)) - .recip() - * 2.0 + fn helmholtz_energy_ab_analytic>(&self, state: &StateHD, delta: D) -> D { + let a = &self.association_parameters; + + // site densities + let rhoa = + state.partial_density[a.component_index[a.sites_a[0].assoc_comp]] * a.sites_a[0].n; + let rhob = + state.partial_density[a.component_index[a.sites_b[0].assoc_comp]] * a.sites_b[0].n; + + // fraction of non-bonded association sites + let sqrt = ((delta * (rhoa - rhob) + 1.0).powi(2) + delta * rhob * 4.0).sqrt(); + let xa = (sqrt + (delta * (rhob - rhoa) + 1.0)).recip() * 2.0; + let xb = (sqrt + (delta * (rhoa - rhob) + 1.0)).recip() * 2.0; + + (rhoa * (xa.ln() - xa * 0.5 + 0.5) + rhob * (xb.ln() - xb * 0.5 + 0.5)) * state.volume } - pub fn assoc_site_frac_a>(deltarho: D, na: f64) -> D { - ((deltarho * 4.0 * na + 1.0).sqrt() + 1.0).recip() * 2.0 + fn helmholtz_energy_cc_analytic>(&self, state: &StateHD, delta: D) -> D { + let a = &self.association_parameters; + + // site density + let rhoc = + state.partial_density[a.component_index[a.sites_c[0].assoc_comp]] * a.sites_c[0].n; + + // fraction of non-bonded association sites + let xc = ((delta * 4.0 * rhoc + 1.0).sqrt() + 1.0).recip() * 2.0; + + rhoc * (xc.ln() - xc * 0.5 + 0.5) * state.volume } #[allow(clippy::too_many_arguments)] fn helmholtz_energy_density_cross_association< - S: Data, D: DualNum + ScalarOperand, + S: Data, >( - &self, - temperature: D, - density: &ArrayBase, - diameter: &Array1, - n2: D, - n3i: D, - xi: D, + rho: &ArrayBase, + delta_ab: &Array2, + delta_cc: &Array2, max_iter: usize, tol: f64, x0: Option<&mut Array1>, - ) -> Result { + ) -> EosResult { // check if density is close to 0 - if density.sum().re() < f64::EPSILON { + if rho.sum().re() < f64::EPSILON { if let Some(x0) = x0 { x0.fill(1.0); } return Ok(D::zero()); } - let assoc_comp = &self.association_parameters.assoc_comp; - let nassoc = assoc_comp.len(); - - // association strength - let delta = self.association_strength(temperature, diameter, n2, n3i, xi); - - // extract parameters of associating components - let na = &self.association_parameters.na; - let nb = &self.association_parameters.nb; - // cross-association according to Michelsen2006 // initialize monomer fraction let mut x = match &x0 { Some(x0) => (*x0).clone(), - None => Array::from_elem(2 * nassoc, 0.2), + None => Array::from_elem(rho.len(), 0.2), }; + let delta_ab_re = delta_ab.map(D::re); + let delta_cc_re = delta_cc.map(D::re); + let rho_re = rho.map(D::re); for k in 0..max_iter { - if Self::newton_step_cross_association::<_, f64>( - nassoc, + if Self::newton_step_cross_association( &mut x, - &delta.map(D::re), - na, - nb, - &density.map(D::re), + &delta_ab_re, + &delta_cc_re, + &rho_re, tol, )? { break; @@ -312,7 +384,7 @@ impl Association

{ // calculate derivatives let mut x_dual = x.mapv(D::from); for _ in 0..D::NDERIV { - Self::newton_step_cross_association(nassoc, &mut x_dual, &delta, na, nb, density, tol)?; + Self::newton_step_cross_association(&mut x_dual, delta_ab, delta_cc, rho, tol)?; } // save monomer fraction @@ -321,44 +393,60 @@ impl Association

{ } // Helmholtz energy density - let xa = x_dual.slice(s![..nassoc]); - let xb = x_dual.slice(s![nassoc..]); let f = |x: D| x.ln() - x * 0.5 + 0.5; - Ok((density * (xa.mapv(f) * na + xb.mapv(f) * nb)).sum()) + Ok((rho * x_dual.mapv(f)).sum()) } - fn newton_step_cross_association, D: DualNum + ScalarOperand>( - nassoc: usize, + fn newton_step_cross_association + ScalarOperand, S: Data>( x: &mut Array1, - delta: &Array2, - na: &Array1, - nb: &Array1, + delta_ab: &Array2, + delta_cc: &Array2, rho: &ArrayBase, tol: f64, - ) -> Result { + ) -> EosResult { + let nassoc = x.len(); // gradient let mut g = x.map(D::recip); // Hessian - let mut h: Array2 = Array::zeros((2 * nassoc, 2 * nassoc)); + let mut h: Array2 = Array::zeros([nassoc; 2]); - // split x array - let (xa, xb) = x.view().split_at(Axis(0), nassoc); + // split arrays + let &[a, b] = delta_ab.shape() else { panic!("wrong shape!") }; + let c = delta_cc.shape()[0]; + let (xa, xc) = x.view().split_at(Axis(0), a + b); + let (xa, xb) = xa.split_at(Axis(0), a); + let (rhoa, rhoc) = rho.view().split_at(Axis(0), a + b); + let (rhoa, rhob) = rhoa.split_at(Axis(0), a); - // calculate gradients and approximate Hessian for i in 0..nassoc { - let d = &delta.index_axis(Axis(0), i) * rho; - - let dnx = (&xb * nb * &d).sum() + 1.0; + // calculate gradients + let (d, dnx) = if i < a { + let d = delta_ab.index_axis(Axis(0), i); + (d, (&xb * &rhob * d).sum() + 1.0) + } else if i < a + b { + let d = delta_ab.index_axis(Axis(1), i - a); + (d, (&xa * &rhoa * d).sum() + 1.0) + } else { + let d = delta_cc.index_axis(Axis(0), i - a - b); + (d, (&xc * &rhoc * d).sum() + 1.0) + }; g[i] -= dnx; - for j in 0..nassoc { - h[(i, nassoc + j)] = -d[j] * nb[j]; - h[(nassoc + i, j)] = -d[j] * na[j]; - } - h[(i, i)] = -dnx / xa[i]; - let dnx = (&xa * na * &d).sum() + 1.0; - g[nassoc + i] -= dnx; - h[(nassoc + i, nassoc + i)] = -dnx / xb[i]; + // approximate hessian + h[(i, i)] = -dnx / x[i]; + if i < a { + for j in 0..b { + h[(i, a + j)] = -d[j] * rhob[j]; + } + } else if i < a + b { + for j in 0..a { + h[(i, j)] = -d[j] * rhoa[j]; + } + } else { + for j in 0..c { + h[(i, a + b + j)] -= d[j] * rhoc[j]; + } + } } // Newton step @@ -407,7 +495,7 @@ mod tests_pcsaft { let mut params = water_parameters(); let mut record = params.pure_records.pop().unwrap(); let mut association_record = record.model_record.association_record.unwrap(); - association_record.na = Some(2.0); + association_record.na = 2.0; record.model_record.association_record = Some(association_record); let params = Arc::new(PcSaftParameters::new_pure(record)); let assoc = Association::new(¶ms, ¶ms.association, 50, 1e-10); diff --git a/src/association/python.rs b/src/association/python.rs index cbc8ed680..ac3211bbc 100644 --- a/src/association/python.rs +++ b/src/association/python.rs @@ -4,19 +4,15 @@ use feos_core::parameter::ParameterError; use pyo3::prelude::*; /// Pure component association parameters -#[pyclass( - name = "AssociationRecord", - text_signature = "(kappa_ab, epsilon_k_ab, na=None, nb=None)" -)] +#[pyclass(name = "AssociationRecord")] #[derive(Clone)] pub struct PyAssociationRecord(pub AssociationRecord); #[pymethods] impl PyAssociationRecord { - #[pyo3(signature = (kappa_ab, epsilon_k_ab, na=None, nb=None))] #[new] - fn new(kappa_ab: f64, epsilon_k_ab: f64, na: Option, nb: Option) -> Self { - Self(AssociationRecord::new(kappa_ab, epsilon_k_ab, na, nb)) + fn new(kappa_ab: f64, epsilon_k_ab: f64, na: f64, nb: f64, nc: f64) -> Self { + Self(AssociationRecord::new(kappa_ab, epsilon_k_ab, na, nb, nc)) } #[getter] @@ -30,15 +26,20 @@ impl PyAssociationRecord { } #[getter] - fn get_na(&self) -> Option { + fn get_na(&self) -> f64 { self.0.na } #[getter] - fn get_nb(&self) -> Option { + fn get_nb(&self) -> f64 { self.0.nb } + #[getter] + fn get_nc(&self) -> f64 { + self.0.nc + } + fn __repr__(&self) -> PyResult { Ok(self.0.to_string()) } diff --git a/src/gc_pcsaft/dft/mod.rs b/src/gc_pcsaft/dft/mod.rs index 438af9ba2..e8e502d65 100644 --- a/src/gc_pcsaft/dft/mod.rs +++ b/src/gc_pcsaft/dft/mod.rs @@ -56,7 +56,7 @@ impl GcPcSaftFunctional { contributions.push(Box::new(att)); // Association - if !parameters.association.assoc_comp.is_empty() { + if !parameters.association.is_empty() { let assoc = Association::new( ¶meters, ¶meters.association, diff --git a/src/gc_pcsaft/dft/parameter.rs b/src/gc_pcsaft/dft/parameter.rs index d6c27d9df..d0931dfac 100644 --- a/src/gc_pcsaft/dft/parameter.rs +++ b/src/gc_pcsaft/dft/parameter.rs @@ -80,7 +80,13 @@ impl ParameterHetero for GcPcSaftFunctionalParameters { sigma.push(segment.model_record.sigma); epsilon_k.push(segment.model_record.epsilon_k); - association_records.push(segment.model_record.association_record); + association_records.push( + segment + .model_record + .association_record + .into_iter() + .collect(), + ); psi_dft.push(segment.model_record.psi_dft.unwrap_or(PSI_GC_DFT)); diff --git a/src/gc_pcsaft/eos/mod.rs b/src/gc_pcsaft/eos/mod.rs index c6e0f9906..eec50fea4 100644 --- a/src/gc_pcsaft/eos/mod.rs +++ b/src/gc_pcsaft/eos/mod.rs @@ -60,7 +60,7 @@ impl GcPcSaft { contributions.push(Box::new(Dispersion { parameters: parameters.clone(), })); - if !parameters.association.assoc_comp.is_empty() { + if !parameters.association.is_empty() { contributions.push(Box::new(Association::new( ¶meters, ¶meters.association, diff --git a/src/gc_pcsaft/eos/parameter.rs b/src/gc_pcsaft/eos/parameter.rs index c25add9ad..795ea2792 100644 --- a/src/gc_pcsaft/eos/parameter.rs +++ b/src/gc_pcsaft/eos/parameter.rs @@ -140,10 +140,10 @@ impl ParameterHetero for GcPcSaftEosParameters { let mut assoc = segment.model_record.association_record; if let Some(mut assoc) = assoc.as_mut() { - assoc.na = Some(assoc.na.unwrap_or(1.0) * count); - assoc.nb = Some(assoc.nb.unwrap_or(1.0) * count); + assoc.na *= count; + assoc.nb *= count; }; - association_records.push(assoc); + association_records.push(assoc.into_iter().collect()); m_i += segment.model_record.m * count; sigma_i += segment.model_record.m * segment.model_record.sigma.powi(3) * count; @@ -287,11 +287,17 @@ impl HardSphereProperties for GcPcSaftEosParameters { impl GcPcSaftEosParameters { pub fn to_markdown(&self) -> String { + let gorup_dict: HashMap<&String, &GcPcSaftRecord> = self + .segment_records + .iter() + .map(|r| (&r.identifier, &r.model_record)) + .collect(); + let mut output = String::new(); let o = &mut output; write!( o, - "|component|molarweight|dipole moment|group|$m$|$\\sigma$|$\\varepsilon$|$\\kappa_{{AB}}$|$\\varepsilon_{{AB}}$|$N_A$|$N_B$|\n|-|-|-|-|-|-|-|-|-|-|-|" + "|component|molarweight|dipole moment|group|$m$|$\\sigma$|$\\varepsilon$|$\\kappa_{{AB}}$|$\\varepsilon_{{AB}}$|$N_A$|$N_B$|$N_C$|\n|-|-|-|-|-|-|-|-|-|-|-|-|" ) .unwrap(); for i in 0..self.m.len() { @@ -305,37 +311,39 @@ impl GcPcSaftEosParameters { .as_ref() .unwrap_or(&format!("Component {}", self.component_index[i] + 1)), self.molarweight[self.component_index[i]], - if let Some(d) = self.dipole_comp.iter().position(|&d| d == i) { + if let Some(d) = self + .dipole_comp + .iter() + .position(|&d| d == self.component_index[i]) + { format!("{}", self.mu[d]) } else { "".into() } ) }; - let association = - if let Some(a) = self.association.assoc_comp.iter().position(|&a| a == i) { - format!( - "{}|{}|{}|{}", - self.association.kappa_ab[a], - self.association.epsilon_k_ab[a], - self.association.na[a], - self.association.nb[a] - ) - } else { - "|||".to_string() - }; + let record = gorup_dict[&self.identifiers[i]]; + let association = if let Some(a) = record.association_record { + format!( + "{}|{}|{}|{}|{}", + a.kappa_ab, a.epsilon_k_ab, a.na, a.nb, a.nc + ) + } else { + "||||".to_string() + }; write!( o, "\n|{}|{}|{}|{}|{}|{}|", component, self.identifiers[i], - self.m[i], - self.sigma[i], - self.epsilon_k[i], + record.m, + record.sigma, + record.epsilon_k, association ) .unwrap(); } + write!(o, "\n\n|component|group 1|group 2|bonds|\n|-|-|-|-|").unwrap(); let mut last_component = None; @@ -426,7 +434,7 @@ pub mod test { 2.7702, 334.29, None, - Some(AssociationRecord::new(0.009583, 2575.9, None, None)), + Some(AssociationRecord::new(0.009583, 2575.9, 1.0, 1.0, 0.0)), None, ), None, diff --git a/src/gc_pcsaft/python/mod.rs b/src/gc_pcsaft/python/mod.rs index 6ad10432a..171254a1e 100644 --- a/src/gc_pcsaft/python/mod.rs +++ b/src/gc_pcsaft/python/mod.rs @@ -67,7 +67,7 @@ impl PyGcPcSaftRecord { #[getter] fn get_association_record(&self) -> Option { - self.0.association_record.clone().map(PyAssociationRecord) + self.0.association_record.map(PyAssociationRecord) } fn __repr__(&self) -> PyResult { diff --git a/src/pcsaft/dft/mod.rs b/src/pcsaft/dft/mod.rs index fab80baad..deefc121d 100644 --- a/src/pcsaft/dft/mod.rs +++ b/src/pcsaft/dft/mod.rs @@ -76,7 +76,7 @@ impl PcSaftFunctional { contributions.push(Box::new(att)); // Association - if !parameters.association.assoc_comp.is_empty() { + if !parameters.association.is_empty() { let assoc = Association::new( ¶meters, ¶meters.association, diff --git a/src/pcsaft/dft/pure_saft_functional.rs b/src/pcsaft/dft/pure_saft_functional.rs index 66279ef3d..801a079a8 100644 --- a/src/pcsaft/dft/pure_saft_functional.rs +++ b/src/pcsaft/dft/pure_saft_functional.rs @@ -18,16 +18,18 @@ const PI36M1: f64 = 1.0 / (36.0 * PI); const N3_CUTOFF: f64 = 1e-5; const N0_CUTOFF: f64 = 1e-9; -#[derive(Clone)] pub struct PureFMTAssocFunctional { parameters: Arc, + association: Association, version: FMTVersion, } impl PureFMTAssocFunctional { pub fn new(parameters: Arc, version: FMTVersion) -> Self { + let association = Association::new(¶meters, ¶meters.association, 50, 1e-10); Self { parameters, + association, version, } } @@ -113,33 +115,22 @@ impl + ScalarOperand> FunctionalContributionDual for PureFMTA // association let a = &p.association; - if a.assoc_comp.len() == 1 { + if !a.is_empty() { let mut xi = -(&n2v * &n2v).sum_axis(Axis(0)) / (&n2 * &n2) + 1.0; xi.iter_mut().zip(&n2).for_each(|(xi, &n2)| { if n2.re() < N0_CUTOFF * 4.0 * PI * p.m[0] * r.re().powi(2) { *xi = N::one(); } }); - - let k = &n2 * &n3m1rec * r; - let deltarho = (((&k / 18.0 + 0.5) * &k * &xi + 1.0) * n3m1rec) - * ((temperature.recip() * a.epsilon_k_aibj[(0, 0)]).exp_m1() - * a.sigma3_kappa_aibj[(0, 0)]) - * (&n0 / p.m[0] * &xi); - - let f = |x: N| x.ln() - x * 0.5 + 0.5; - phi = phi - + if a.nb[0] > 0.0 { - let xa = deltarho.mapv(|d| { - Association::::assoc_site_frac_ab(d, a.na[0], a.nb[0]) - }); - let xb = (xa.clone() - 1.0) * a.na[0] / a.nb[0] + 1.0; - (n0 / p.m[0] * xi) * (xa.mapv(f) * a.na[0] + xb.mapv(f) * a.nb[0]) - } else { - let xa = deltarho - .mapv(|d| Association::::assoc_site_frac_a(d, a.na[0])); - n0 / p.m[0] * xi * (xa.mapv(f) * a.na[0]) - }; + let rho0 = (&n0 / p.m[0] * &xi).insert_axis(Axis(0)); + + phi += &(self.association.calculate_helmholtz_energy_density( + temperature, + &rho0, + &n2, + &n3m1rec, + &xi, + ))?; } Ok(phi) diff --git a/src/pcsaft/eos/mod.rs b/src/pcsaft/eos/mod.rs index 48f0b05fa..c11a7714a 100644 --- a/src/pcsaft/eos/mod.rs +++ b/src/pcsaft/eos/mod.rs @@ -86,7 +86,7 @@ impl PcSaft { variant: options.dq_variant, })); }; - if !parameters.association.assoc_comp.is_empty() { + if !parameters.association.is_empty() { contributions.push(Box::new(Association::new( ¶meters, ¶meters.association, diff --git a/src/pcsaft/eos/qspr.rs b/src/pcsaft/eos/qspr.rs index baa2e6436..0f2ce9b4b 100644 --- a/src/pcsaft/eos/qspr.rs +++ b/src/pcsaft/eos/qspr.rs @@ -70,12 +70,13 @@ pub struct QSPR { impl> IdealGasContributionDual for QSPR { fn de_broglie_wavelength(&self, temperature: D, components: usize) -> Array1 { - let (c_300, c_400) = match self.parameters.association.assoc_comp.len() { - 0 => match self.parameters.ndipole + self.parameters.nquadpole { + let (c_300, c_400) = if self.parameters.association.is_empty() { + match self.parameters.ndipole + self.parameters.nquadpole { 0 => (NA_NP_300, NA_NP_400), _ => (NA_P_300, NA_P_400), - }, - _ => (AP_300, AP_400), + } + } else { + (AP_300, AP_400) }; Array1::from_shape_fn(components, |i| { diff --git a/src/pcsaft/parameters.rs b/src/pcsaft/parameters.rs index 04fa3a067..2b4881f93 100644 --- a/src/pcsaft/parameters.rs +++ b/src/pcsaft/parameters.rs @@ -70,14 +70,23 @@ impl FromSegments for PcSaftRecord { [ record.kappa_ab * n, record.epsilon_k_ab * n, - record.na.unwrap_or(1.0) * n, - record.nb.unwrap_or(1.0) * n, + record.na * n, + record.nb * n, + record.nc * n, ] }) }) - .reduce(|a, b| [a[0] + b[0], a[1] + b[1], a[2] + b[2], a[3] + b[3]]) - .map(|[kappa_ab, epsilon_k_ab, na, nb]| { - AssociationRecord::new(kappa_ab, epsilon_k_ab, Some(na), Some(nb)) + .reduce(|a, b| { + [ + a[0] + b[0], + a[1] + b[1], + a[2] + b[2], + a[3] + b[3], + a[4] + b[4], + ] + }) + .map(|[kappa_ab, epsilon_k_ab, na, nb, nc]| { + AssociationRecord::new(kappa_ab, epsilon_k_ab, na, nb, nc) }); // entropy scaling @@ -219,14 +228,19 @@ impl PcSaftRecord { epsilon_k_ab: Option, na: Option, nb: Option, + nc: Option, viscosity: Option<[f64; 4]>, diffusion: Option<[f64; 5]>, thermal_conductivity: Option<[f64; 4]>, ) -> PcSaftRecord { let association_record = match (kappa_ab, epsilon_k_ab) { - (Some(kappa_ab), Some(epsilon_k_ab)) => { - Some(AssociationRecord::new(kappa_ab, epsilon_k_ab, na, nb)) - } + (Some(kappa_ab), Some(epsilon_k_ab)) => Some(AssociationRecord::new( + kappa_ab, + epsilon_k_ab, + na.unwrap_or(0.0), + nb.unwrap_or(0.0), + nc.unwrap_or(0.0), + )), (None, None) => None, _ => { panic!("To model association, both kappa_ab and epsilon_k_ab need to be specified.") @@ -338,7 +352,7 @@ impl Parameter for PcSaftParameters { epsilon_k[i] = r.epsilon_k; mu[i] = r.mu.unwrap_or(0.0); q[i] = r.q.unwrap_or(0.0); - association_records.push(r.association_record); + association_records.push(r.association_record.into_iter().collect()); viscosity.push(r.viscosity); diffusion.push(r.diffusion); thermal_conductivity.push(r.thermal_conductivity); @@ -470,7 +484,7 @@ impl PcSaftParameters { let o = &mut output; write!( o, - "|component|molarweight|$m$|$\\sigma$|$\\varepsilon$|$\\mu$|$Q$|$\\kappa_{{AB}}$|$\\varepsilon_{{AB}}$|$N_A$|$N_B$|\n|-|-|-|-|-|-|-|-|-|-|-|" + "|component|molarweight|$m$|$\\sigma$|$\\varepsilon$|$\\mu$|$Q$|$\\kappa_{{AB}}$|$\\varepsilon_{{AB}}$|$N_A$|$N_B$|$N_C$|\n|-|-|-|-|-|-|-|-|-|-|-|-|" ) .unwrap(); for (i, record) in self.pure_records.iter().enumerate() { @@ -479,10 +493,10 @@ impl PcSaftParameters { let association = record .model_record .association_record - .unwrap_or_else(|| AssociationRecord::new(0.0, 0.0, None, None)); + .unwrap_or_else(|| AssociationRecord::new(0.0, 0.0, 0.0, 0.0, 0.0)); write!( o, - "\n|{}|{}|{}|{}|{}|{}|{}|{}|{}|{}|{}|", + "\n|{}|{}|{}|{}|{}|{}|{}|{}|{}|{}|{}|{}|", component, record.molarweight, record.model_record.m, @@ -492,8 +506,9 @@ impl PcSaftParameters { record.model_record.q.unwrap_or(0.0), association.kappa_ab, association.epsilon_k_ab, - association.na.unwrap_or(1.0), - association.nb.unwrap_or(1.0) + association.na, + association.nb, + association.nc ) .unwrap(); } @@ -515,13 +530,13 @@ impl std::fmt::Display for PcSaftParameters { if !self.quadpole_comp.is_empty() { write!(f, "\n\tq={}", self.q)?; } - if !self.association.assoc_comp.is_empty() { - write!(f, "\n\tassociating={}", self.association.assoc_comp)?; - write!(f, "\n\tkappa_ab={}", self.association.kappa_ab)?; - write!(f, "\n\tepsilon_k_ab={}", self.association.epsilon_k_ab)?; - write!(f, "\n\tna={}", self.association.na)?; - write!(f, "\n\tnb={}", self.association.nb)?; - } + // if !self.association.is_empty() { + // write!(f, "\n\tassociating={}", self.association.assoc_comp)?; + // write!(f, "\n\tkappa_ab={}", self.association.kappa_ab)?; + // write!(f, "\n\tepsilon_k_ab={}", self.association.epsilon_k_ab)?; + // write!(f, "\n\tna={}", self.association.na)?; + // write!(f, "\n\tnb={}", self.association.nb)?; + // } if !self.k_ij.iter().all(|k| k.is_zero()) { write!(f, "\n\tk_ij=\n{}", self.k_ij)?; } @@ -649,7 +664,9 @@ pub mod utils { "sigma": 3.000683, "epsilon_k": 366.5121, "kappa_ab": 0.034867983, - "epsilon_k_ab": 2500.6706 + "epsilon_k_ab": 2500.6706, + "na": 1.0, + "nb": 1.0 }, "molarweight": 18.0152 }"#; diff --git a/src/pcsaft/python.rs b/src/pcsaft/python.rs index 658722422..4da44bc6b 100644 --- a/src/pcsaft/python.rs +++ b/src/pcsaft/python.rs @@ -36,6 +36,7 @@ impl PyPcSaftRecord { epsilon_k_ab: Option, na: Option, nb: Option, + nc: Option, viscosity: Option<[f64; 4]>, diffusion: Option<[f64; 5]>, thermal_conductivity: Option<[f64; 4]>, @@ -50,6 +51,7 @@ impl PyPcSaftRecord { epsilon_k_ab, na, nb, + nc, viscosity, diffusion, thermal_conductivity, @@ -93,12 +95,17 @@ impl PyPcSaftRecord { #[getter] fn get_na(&self) -> Option { - self.0.association_record.and_then(|a| a.na) + self.0.association_record.map(|a| a.na) } #[getter] fn get_nb(&self) -> Option { - self.0.association_record.and_then(|a| a.nb) + self.0.association_record.map(|a| a.nb) + } + + #[getter] + fn get_nc(&self) -> Option { + self.0.association_record.map(|a| a.nc) } #[getter] diff --git a/tests/gc_pcsaft/dft.rs b/tests/gc_pcsaft/dft.rs index b777405b5..836643501 100644 --- a/tests/gc_pcsaft/dft.rs +++ b/tests/gc_pcsaft/dft.rs @@ -120,13 +120,13 @@ fn test_bulk_association() -> Result<(), Box> { segment_records.clone(), None, )?); - let eos = Arc::new(GcPcSaft::new(eos_parameters.clone())); + let eos = Arc::new(GcPcSaft::new(eos_parameters)); let func_parameters = Arc::new(GcPcSaftFunctionalParameters::from_segments( vec![ethylene_glycol], segment_records, None, )?); - let func = Arc::new(GcPcSaftFunctional::new(func_parameters.clone())); + let func = Arc::new(GcPcSaftFunctional::new(func_parameters)); let t = 200.0 * KELVIN; let v = 0.002 * METER.powi(3); @@ -135,35 +135,35 @@ fn test_bulk_association() -> Result<(), Box> { let state_func = State::new_nvt(&func, t, v, &n)?; let p_eos = state_eos.pressure_contributions(); let p_func = state_func.pressure_contributions(); - println!( - "Equation of state: - \tcomps: {} - \tkappa_ab: {} - \tepsilon_k_ab: {} - \tna: {} - \tnb: {}", - eos_parameters.association.assoc_comp, - eos_parameters.association.kappa_ab, - eos_parameters.association.epsilon_k_ab, - eos_parameters.association.na, - eos_parameters.association.nb, - ); + // println!( + // "Equation of state: + // \tcomps: {} + // \tkappa_ab: {} + // \tepsilon_k_ab: {} + // \tna: {} + // \tnb: {}", + // eos_parameters.association.assoc_comp, + // eos_parameters.association.kappa_ab, + // eos_parameters.association.epsilon_k_ab, + // eos_parameters.association.na, + // eos_parameters.association.nb, + // ); for (s, x) in &p_eos { println!("{s:18}: {x:21.16}"); } - println!( - "\nHelmholtz energy functional: - \tcomps: {} - \tkappa_ab: {} - \tepsilon_k_ab: {} - \tna: {} - \tnb: {}", - func_parameters.association.assoc_comp, - func_parameters.association.kappa_ab, - func_parameters.association.epsilon_k_ab, - func_parameters.association.na, - func_parameters.association.nb, - ); + // println!( + // "\nHelmholtz energy functional: + // \tcomps: {} + // \tkappa_ab: {} + // \tepsilon_k_ab: {} + // \tna: {} + // \tnb: {}", + // func_parameters.association.assoc_comp, + // func_parameters.association.kappa_ab, + // func_parameters.association.epsilon_k_ab, + // func_parameters.association.na, + // func_parameters.association.nb, + // ); for (s, x) in &p_func { println!("{s:26}: {x:21.16}"); } diff --git a/tests/pcsaft/test_parameters.json b/tests/pcsaft/test_parameters.json index 988fcb3b8..c49a5de31 100644 --- a/tests/pcsaft/test_parameters.json +++ b/tests/pcsaft/test_parameters.json @@ -93,7 +93,9 @@ "sigma": 3.000683, "epsilon_k": 366.5121, "kappa_ab": 0.034867983, - "epsilon_k_ab": 2500.6706 + "epsilon_k_ab": 2500.6706, + "na": 1.0, + "nb": 1.0 }, "molarweight": 18.0152 }, From 8b3e47411edb5c1357e38f04336f89d6a788eec3 Mon Sep 17 00:00:00 2001 From: Philipp Rehner Date: Sat, 25 Mar 2023 11:43:04 +0100 Subject: [PATCH 3/7] Michael L. Michelsen to the rescue once again --- docs/theory/models/association.md | 2 +- src/association/mod.rs | 11 +++++++++-- 2 files changed, 10 insertions(+), 3 deletions(-) diff --git a/docs/theory/models/association.md b/docs/theory/models/association.md index d8f28342e..7cfc38193 100644 --- a/docs/theory/models/association.md +++ b/docs/theory/models/association.md @@ -34,7 +34,7 @@ with the dirac delta function $\delta_{\alpha\beta}$. However, [Michelsen 2006]( $$\hat{H}_{\alpha\beta}=-\frac{\delta_{\alpha\beta}}{X_\alpha}\left(1+\sum_\beta\rho_\beta X_\beta\Delta^{\alpha\beta}\right)-\rho_\beta\Delta^{\alpha\beta}$$ -instead. $X_\alpha$ can then be solved robustly using a standard undamped Newton algorithm. With the split in $AB$ and $CC$ association the two kinds different versions could be solved separately from each other. This is currently not implemented, because use cases are rare to nonexistent and the benefit is small. +instead. $X_\alpha$ can then be solved robustly using a Newton algorithm with a check that ensures that the values of $X_\alpha$ remain positive. With the split in $AB$ and $CC$ association the two kinds different versions could be solved separately from each other. This is currently not implemented, because use cases are rare to nonexistent and the benefit is small. A drastic improvement in performance, however, can be achieved by solving eq. {eq}`eqn:x_assoc` analytically for simple cases. If there is only one $A$ and one $B$ site the corresponding fractions of non-bonded association sites can be calculated from diff --git a/src/association/mod.rs b/src/association/mod.rs index cd59bf0ec..627cc0ebd 100644 --- a/src/association/mod.rs +++ b/src/association/mod.rs @@ -8,7 +8,6 @@ use num_dual::*; use num_traits::Zero; use serde::{Deserialize, Serialize}; use std::fmt; -use std::ops::SubAssign; use std::sync::Arc; #[cfg(feature = "dft")] @@ -450,7 +449,15 @@ impl Association

{ } // Newton step - x.sub_assign(&LU::new(h)?.solve(&g)); + // avoid stepping to negative values for x (see Michelsen 2006) + let delta_x = LU::new(h)?.solve(&g); + Zip::from(x).and(&delta_x).for_each(|x, &delta_x| { + if delta_x.re() < x.re() * 0.8 { + *x -= delta_x + } else { + *x *= 0.2 + } + }); // check convergence Ok(norm(&g.map(D::re)) < tol) From bbd54eb8ce09c43ce5bec77ed47c7c16bfbf1c6e Mon Sep 17 00:00:00 2001 From: Philipp Rehner Date: Mon, 3 Apr 2023 10:19:45 +0200 Subject: [PATCH 4/7] Add hard-sphere theory guide --- docs/theory/models/association.md | 2 +- docs/theory/models/hard_spheres.md | 56 +++++++++++++++++++++++++++++- 2 files changed, 56 insertions(+), 2 deletions(-) diff --git a/docs/theory/models/association.md b/docs/theory/models/association.md index 7cfc38193..a0c359329 100644 --- a/docs/theory/models/association.md +++ b/docs/theory/models/association.md @@ -55,7 +55,7 @@ $$\Delta^{\alpha\beta}=\left(\frac{1}{1-\zeta_3}+\frac{\frac{3}{2}d_{ij}\zeta_2} with -$$d_{ij}=\frac{2d_id_j}{d_i+d_j},~~~~d_k=\sigma_k\left(1-0.12e^{\frac{-3\varepsilon_i}{k_\mathrm{B}T}}\right)$$ +$$d_{ij}=\frac{2d_id_j}{d_i+d_j},~~~~d_k=\sigma_k\left(1-0.12e^{\frac{-3\varepsilon_k}{k_\mathrm{B}T}}\right)$$ and diff --git a/docs/theory/models/hard_spheres.md b/docs/theory/models/hard_spheres.md index 5962c582a..26c2a278f 100644 --- a/docs/theory/models/hard_spheres.md +++ b/docs/theory/models/hard_spheres.md @@ -1,3 +1,57 @@ # Hard spheres -under construction \ No newline at end of file +$\text{FeO}_\text{s}$ provides an implementation of the Boublík-Mansoori-Carnahan-Starling-Leland (BMCSL) equation of state ([Boublík, 1970](https://doi.org/10.1063/1.1673824), [Mansoori et al., 1971](https://doi.org/10.1063/1.1675048)) for hard-sphere mixtures which is often used as reference contribution in SAFT equations of state. The implementation is generalized to allow the description of non-sperical or fused-sphere reference fluids. + +The reduced Helmholtz energy density is calculated according to + +$$\frac{\beta A}{V}=\frac{6}{\pi}\left(\frac{3\zeta_1\zeta_2}{1-\zeta_3}+\frac{\zeta_2^3}{\zeta_3\left(1-\zeta_3\right)^2}+\left(\frac{\zeta_2^3}{\zeta_3^2}-\zeta_0\right)\ln\left(1-\zeta_3\right)\right)$$ + +with the packing fractions + +$$\zeta_k=\frac{\pi}{6}\sum_\alpha C_{k,\alpha}\rho_\alpha d_\alpha^k,~~~~~~~~k=0\ldots 3$$ + +The geometry coefficients $C_{k,\alpha}$ and the segment diameters $d_\alpha$ depend on the context in which the model is used. The following table shows how the expression can be reused in various models. For details on the fused-sphere chain model, check the [repository](https://github.com/feos-org/feos-fused-chains) or the [publication](https://doi.org/10.1103/PhysRevE.105.034110). + +||Hard spheres|PC-SAFT|Fused-sphere chains| +|-|:-:|:-:|:-:| +|$d_\alpha$|$\sigma_\alpha$|$\sigma_\alpha\left(1-0.12e^{\frac{-3\varepsilon_\alpha}{k_\mathrm{B}T}}\right)$|$\sigma_\alpha$| +|$C_{0,\alpha}$|$1$|$m_\alpha$|$1$| +|$C_{1,\alpha}$|$1$|$m_\alpha$|$A_\alpha^*$| +|$C_{1,\alpha}$|$1$|$m_\alpha$|$A_\alpha^*$| +|$C_{1,\alpha}$|$1$|$m_\alpha$|$V_\alpha^*$| + +## Fundamental measure theory + +An model for inhomogeneous mixtues of hard spheres is provided by fundamental measure theory (FMT, [Rosenfeld, 1989](https://doi.org/10.1103/PhysRevLett.63.980)). Different variants have been proposed of which only those that are consistent with the BMCSL equation of state in the homogeneous limit are currently considered in $\text{FeO}_\text{s}$ (exluding, e.g., the original Rosenfeld and White-Bear II variants). + +The Helmholtz energy density is calculated according to + +$$\beta f=-n_0\ln\left(1-n_3\right)+\frac{n_{12}}{1-n_3}+\frac{1}{36\pi}n_2n_{22}f_3(n_3)$$ + +The expressions for $n_{12}$ and $n_{22}$ depend on the FMT version. + +|version|$n_{12}$|$n_{22}$|references| +|-|:-:|:-:|:-:| +|WhiteBear|$n_1n_2-\vec n_1\cdot\vec n_2$|$n_2^2-3\vec n_2\cdot\vec n_2$|[Roth et al., 2002](https://doi.org/10.1088/0953-8984/14/46/313), [Yu and Wu, 2002](https://doi.org/10.1063/1.1520530)| +|KierlikRosinberg|$n_1n_2$|$n_2^2$|[Kierlik and Rosinberg, 1990](https://doi.org/10.1103/PhysRevA.42.3382)| +|AntiSymWhiteBear|$n_1n_2-\vec n_1\cdot\vec n_2$|$n_2^2\left(1-\frac{\vec n_2\cdot\vec n_2}{n_2^2}\right)^3$|[Rosenfeld et al., 1997](https://doi.org/10.1103/PhysRevE.55.4245), [Kessler et al., 2021](https://doi.org/10.1016/j.micromeso.2021.111263)| + +For small $n_3$, the value of $f(n_3)$ numerically diverges. Therefore, it is approximated with a Taylor expansion. + +$$f_3=\begin{cases}\frac{n_3+\left(1-n_3\right)^2\ln\left(1-n_3\right)}{n_3^2\left(1-n_3\right)^2}&\text{if }n_3>10^{-5}\\ +\frac{3}{2}+\frac{8}{3}n_3+\frac{15}{4}n_3^2+\frac{24}{5}n_3^3+\frac{35}{6}n_3^4&\text{else}\end{cases}$$ + +The weighted densities $n_k(\mathbf{r})$ are calculated by convolving the density profiles $\rho_\alpha(\mathbf{r})$ with weight functions $\omega_k^\alpha(\mathbf{r})$ + +$$n_k(\mathbf{r})=\sum_\alpha n_k^\alpha(\mathbf{r})=\sum_\alpha\int\rho_\alpha(\mathbf{r}')\omega_k^\alpha(\mathbf{r}-\mathbf{r}')\mathrm{d}\mathbf{r}'$$ + +which differ between the different FMT versions. + +||WhiteBear/AntiSymWhiteBear|KierlikRosinberg| +|-|:-:|:-:| +|$\omega_0^\alpha(\mathbf{r})$|$\frac{C_{0,\alpha}}{\pi\sigma_\alpha^2}\,\delta\!\left(\frac{d_\alpha}{2}-\|\mathbf{r}\|\right)$|$C_{0,\alpha}\left(-\frac{1}{8\pi}\,\delta''\!\left(\frac{d_\alpha}{2}-\|\mathbf{r}\|\right)+\frac{1}{2\pi\|\mathbf{r}\|}\,\delta'\!\left(\frac{d_\alpha}{2}-\|\mathbf{r}\|\right)\right)$| +|$\omega_1^\alpha(\mathbf{r})$|$\frac{C_{1,\alpha}}{2\pi\sigma_\alpha}\,\delta\!\left(\frac{d_\alpha}{2}-\|\mathbf{r}\|\right)$|$\frac{C_{1,\alpha}}{8\pi}\,\delta'\!\left(\frac{d_\alpha}{2}-\|\mathbf{r}\|\right)$| +|$\omega_2^\alpha(\mathbf{r})$|$C_{2,\alpha}\,\delta\!\left(\frac{d_\alpha}{2}-\|\mathbf{r}\|\right)$|$C_{2,\alpha}\,\delta\!\left(\frac{d_\alpha}{2}-\|\mathbf{r}\|\right)$| +|$\omega_3^\alpha(\mathbf{r})$|$C_{3,\alpha}\,\Theta\!\left(\frac{d_\alpha}{2}-\|\mathbf{r}\|\right)$|$C_{3,\alpha}\,\Theta\!\left(\frac{d_\alpha}{2}-\|\mathbf{r}\|\right)$| +|$\vec\omega_1^\alpha(\mathbf{r})$|$C_{3,\alpha}\frac{\mathbf{r}}{2\pi\sigma_\alpha\|\mathbf{r}\|}\,\delta\!\left(\frac{d_\alpha}{2}-\|\mathbf{r}\|\right)$|-| +|$\vec\omega_2^\alpha(\mathbf{r})$|$C_{3,\alpha}\frac{\mathbf{r}}{\|\mathbf{r}\|}\,\delta\!\left(\frac{d_\alpha}{2}-\|\mathbf{r}\|\right)$|-| \ No newline at end of file From acb5b35a4721f799d723cb5d0eca055f2413c50d Mon Sep 17 00:00:00 2001 From: Philipp Rehner Date: Mon, 3 Apr 2023 13:27:36 +0200 Subject: [PATCH 5/7] update parameter files --- parameters/pcsaft/gross2002.json | 72 ++++++++++---- parameters/pcsaft/loetgeringlin2015_homo.json | 4 + parameters/pcsaft/loetgeringlin2018.json | 72 ++++++++++++++ parameters/pcsaft/rehner2020.json | 94 +++++++++++-------- 4 files changed, 186 insertions(+), 56 deletions(-) diff --git a/parameters/pcsaft/gross2002.json b/parameters/pcsaft/gross2002.json index b888201d8..3b8737154 100644 --- a/parameters/pcsaft/gross2002.json +++ b/parameters/pcsaft/gross2002.json @@ -13,7 +13,9 @@ "sigma": 3.23, "epsilon_k": 188.9, "kappa_ab": 0.035176, - "epsilon_k_ab": 2899.5 + "epsilon_k_ab": 2899.5, + "na": 1.0, + "nb": 1.0 }, "molarweight": 32.042 }, @@ -31,7 +33,9 @@ "sigma": 3.1771, "epsilon_k": 198.24, "kappa_ab": 0.032384, - "epsilon_k_ab": 2653.4 + "epsilon_k_ab": 2653.4, + "na": 1.0, + "nb": 1.0 }, "molarweight": 46.069 }, @@ -49,7 +53,9 @@ "sigma": 3.2522, "epsilon_k": 233.4, "kappa_ab": 0.015268, - "epsilon_k_ab": 2276.8 + "epsilon_k_ab": 2276.8, + "na": 1.0, + "nb": 1.0 }, "molarweight": 60.096 }, @@ -67,7 +73,9 @@ "sigma": 3.6139, "epsilon_k": 259.59, "kappa_ab": 0.006692, - "epsilon_k_ab": 2544.6 + "epsilon_k_ab": 2544.6, + "na": 1.0, + "nb": 1.0 }, "molarweight": 74.123 }, @@ -85,7 +93,9 @@ "sigma": 3.4508, "epsilon_k": 247.28, "kappa_ab": 0.010319, - "epsilon_k_ab": 2252.1 + "epsilon_k_ab": 2252.1, + "na": 1.0, + "nb": 1.0 }, "molarweight": 88.15 }, @@ -103,7 +113,9 @@ "sigma": 3.6735, "epsilon_k": 262.32, "kappa_ab": 0.005747, - "epsilon_k_ab": 2538.9 + "epsilon_k_ab": 2538.9, + "na": 1.0, + "nb": 1.0 }, "molarweight": 102.177 }, @@ -121,7 +133,9 @@ "sigma": 3.545, "epsilon_k": 253.46, "kappa_ab": 0.001155, - "epsilon_k_ab": 2878.5 + "epsilon_k_ab": 2878.5, + "na": 1.0, + "nb": 1.0 }, "molarweight": 116.203 }, @@ -139,7 +153,9 @@ "sigma": 3.7145, "epsilon_k": 262.74, "kappa_ab": 0.002197, - "epsilon_k_ab": 2754.8 + "epsilon_k_ab": 2754.8, + "na": 1.0, + "nb": 1.0 }, "molarweight": 130.23 }, @@ -157,7 +173,9 @@ "sigma": 3.7292, "epsilon_k": 263.64, "kappa_ab": 0.001427, - "epsilon_k_ab": 2941.9 + "epsilon_k_ab": 2941.9, + "na": 1.0, + "nb": 1.0 }, "molarweight": 144.257 }, @@ -175,7 +193,9 @@ "sigma": 3.2085, "epsilon_k": 208.42, "kappa_ab": 0.024675, - "epsilon_k_ab": 2253.9 + "epsilon_k_ab": 2253.9, + "na": 1.0, + "nb": 1.0 }, "molarweight": 60.096 }, @@ -193,7 +213,9 @@ "sigma": 3.9053, "epsilon_k": 266.01, "kappa_ab": 0.001863, - "epsilon_k_ab": 2618.8 + "epsilon_k_ab": 2618.8, + "na": 1.0, + "nb": 1.0 }, "molarweight": 88.15 }, @@ -211,7 +233,9 @@ "sigma": 3.0007, "epsilon_k": 366.51, "kappa_ab": 0.034868, - "epsilon_k_ab": 2500.7 + "epsilon_k_ab": 2500.7, + "na": 1.0, + "nb": 1.0 }, "molarweight": 18.015 }, @@ -229,7 +253,9 @@ "sigma": 2.8906, "epsilon_k": 214.94, "kappa_ab": 0.095103, - "epsilon_k_ab": 684.3 + "epsilon_k_ab": 684.3, + "na": 1.0, + "nb": 1.0 }, "molarweight": 31.06 }, @@ -247,7 +273,9 @@ "sigma": 3.1343, "epsilon_k": 221.53, "kappa_ab": 0.017275, - "epsilon_k_ab": 854.7 + "epsilon_k_ab": 854.7, + "na": 1.0, + "nb": 1.0 }, "molarweight": 45.09 }, @@ -265,7 +293,9 @@ "sigma": 3.5347, "epsilon_k": 250.52, "kappa_ab": 0.022674, - "epsilon_k_ab": 1028.1 + "epsilon_k_ab": 1028.1, + "na": 1.0, + "nb": 1.0 }, "molarweight": 59.11 }, @@ -283,7 +313,9 @@ "sigma": 3.4777, "epsilon_k": 231.8, "kappa_ab": 0.02134, - "epsilon_k_ab": 932.2 + "epsilon_k_ab": 932.2, + "na": 1.0, + "nb": 1.0 }, "molarweight": 59.11 }, @@ -301,7 +333,9 @@ "sigma": 3.7021, "epsilon_k": 335.47, "kappa_ab": 0.074883, - "epsilon_k_ab": 1351.6 + "epsilon_k_ab": 1351.6, + "na": 1.0, + "nb": 1.0 }, "molarweight": 93.13 }, @@ -319,7 +353,9 @@ "sigma": 3.8582, "epsilon_k": 211.59, "kappa_ab": 0.07555, - "epsilon_k_ab": 3044.4 + "epsilon_k_ab": 3044.4, + "na": 1.0, + "nb": 1.0 }, "molarweight": 60.053 } diff --git a/parameters/pcsaft/loetgeringlin2015_homo.json b/parameters/pcsaft/loetgeringlin2015_homo.json index ff9df6b16..030332969 100644 --- a/parameters/pcsaft/loetgeringlin2015_homo.json +++ b/parameters/pcsaft/loetgeringlin2015_homo.json @@ -277,6 +277,8 @@ "epsilon_k": 488.66, "epsilon_k_ab": 2517.0, "kappa_ab": 0.006825, + "na": 1.0, + "nb": 1.0, "viscosity": [ -15.7583e-3, -2.5654e-1, @@ -294,6 +296,8 @@ "epsilon_k": 467.59, "epsilon_k_ab": 1064.6, "kappa_ab": 0.026662, + "na": 1.0, + "nb": 1.0, "viscosity": [ -4.4048e-3, -0.6089e-1, diff --git a/parameters/pcsaft/loetgeringlin2018.json b/parameters/pcsaft/loetgeringlin2018.json index 989e3439c..153e40ed3 100644 --- a/parameters/pcsaft/loetgeringlin2018.json +++ b/parameters/pcsaft/loetgeringlin2018.json @@ -896,6 +896,8 @@ "mu": 0.9204, "kappa_ab": 0.03, "epsilon_k_ab": 1547.4422, + "na": 1.0, + "nb": 1.0, "viscosity": [ -0.6824, -2.129, @@ -921,6 +923,8 @@ "mu": 1.0313, "kappa_ab": 0.03, "epsilon_k_ab": 959.5053, + "na": 1.0, + "nb": 1.0, "viscosity": [ -0.3807, 2.1585, @@ -946,6 +950,8 @@ "mu": 1.07, "kappa_ab": 0.03, "epsilon_k_ab": 39.2881, + "na": 1.0, + "nb": 1.0, "viscosity": [ -0.8987, -2.1083, @@ -1040,6 +1046,8 @@ "mu": 1.391, "kappa_ab": 0.03, "epsilon_k_ab": 1497.5173, + "na": 1.0, + "nb": 1.0, "viscosity": [ -0.8691, -2.3216, @@ -1065,6 +1073,8 @@ "mu": 1.6698, "kappa_ab": 0.03, "epsilon_k_ab": 1846.2029, + "na": 1.0, + "nb": 1.0, "viscosity": [ -1.0065, -1.88906, @@ -1090,6 +1100,8 @@ "mu": 1.6608, "kappa_ab": 0.03, "epsilon_k_ab": 1718.0488, + "na": 1.0, + "nb": 1.0, "viscosity": [ -0.9864, -1.3818, @@ -1363,6 +1375,8 @@ "mu": 1.409, "kappa_ab": 0.03, "epsilon_k_ab": 2644.5966, + "na": 1.0, + "nb": 1.0, "viscosity": [ -1.3992, -3.2065, @@ -1388,6 +1402,8 @@ "mu": 1.6189, "kappa_ab": 0.03, "epsilon_k_ab": 1967.8312, + "na": 1.0, + "nb": 1.0, "viscosity": [ -1.5371, -2.97397, @@ -1526,6 +1542,8 @@ "mu": 1.6908, "kappa_ab": 0.03, "epsilon_k_ab": 2002.9613, + "na": 1.0, + "nb": 1.0, "viscosity": [ -1.6684, -3.9471, @@ -1639,6 +1657,8 @@ "mu": 2.4103, "kappa_ab": 0.03, "epsilon_k_ab": 2711.6591, + "na": 1.0, + "nb": 1.0, "viscosity": [ -0.7689, -0.9258, @@ -1664,6 +1684,8 @@ "mu": 1.6908, "kappa_ab": 0.03, "epsilon_k_ab": 2514.0609, + "na": 1.0, + "nb": 1.0, "viscosity": [ -1.08153, -1.1998, @@ -1870,6 +1892,8 @@ "mu": 1.6608, "kappa_ab": 0.03, "epsilon_k_ab": 2220.2939, + "na": 1.0, + "nb": 1.0, "viscosity": [ -1.1698, -3.4917, @@ -1895,6 +1919,8 @@ "mu": 1.7388, "kappa_ab": 0.03, "epsilon_k_ab": 1956.8424, + "na": 1.0, + "nb": 1.0, "viscosity": [ -1.3075, -2.24414, @@ -1920,6 +1946,8 @@ "mu": 1.6488, "kappa_ab": 0.03, "epsilon_k_ab": 1878.5922, + "na": 1.0, + "nb": 1.0, "viscosity": [ -1.3001, -1.7492, @@ -2080,6 +2108,8 @@ "mu": 1.6698, "kappa_ab": 0.03, "epsilon_k_ab": 2140.2545, + "na": 1.0, + "nb": 1.0, "viscosity": [ -1.8993, -3.0227, @@ -2127,6 +2157,8 @@ "mu": 1.5889, "kappa_ab": 0.03, "epsilon_k_ab": 82.5108, + "na": 1.0, + "nb": 1.0, "viscosity": [ -1.0803, -2.3348, @@ -2152,6 +2184,8 @@ "mu": 1.6488, "kappa_ab": 0.03, "epsilon_k_ab": 1930.6231, + "na": 1.0, + "nb": 1.0, "viscosity": [ -1.218, -2.0634, @@ -2177,6 +2211,8 @@ "mu": 1.6578, "kappa_ab": 0.03, "epsilon_k_ab": 1793.3615, + "na": 1.0, + "nb": 1.0, "viscosity": [ -1.2075, -1.7111, @@ -2291,6 +2327,8 @@ "mu": 1.3101, "kappa_ab": 0.03, "epsilon_k_ab": 1103.6088, + "na": 1.0, + "nb": 1.0, "viscosity": [ -0.369, -1.2753, @@ -2338,6 +2376,8 @@ "mu": 1.6998, "kappa_ab": 0.03, "epsilon_k_ab": 2519.7116, + "na": 1.0, + "nb": 1.0, "viscosity": [ -0.577, -0.44059, @@ -2454,6 +2494,8 @@ "mu": 0.8994, "kappa_ab": 0.03, "epsilon_k_ab": 2731.5672, + "na": 1.0, + "nb": 1.0, "viscosity": [ -1.3279, -3.1517, @@ -2479,6 +2521,8 @@ "mu": 1.6099, "kappa_ab": 0.03, "epsilon_k_ab": 2073.3555, + "na": 1.0, + "nb": 1.0, "viscosity": [ -1.4657, -2.71167, @@ -2592,6 +2636,8 @@ "mu": 1.6608, "kappa_ab": 0.03, "epsilon_k_ab": 2150.7482, + "na": 1.0, + "nb": 1.0, "viscosity": [ -2.0036, -3.7363, @@ -2639,6 +2685,8 @@ "mu": 1.421, "kappa_ab": 0.03, "epsilon_k_ab": 2476.9879, + "na": 1.0, + "nb": 1.0, "viscosity": [ -1.2518, -2.9995, @@ -2664,6 +2712,8 @@ "mu": 1.6488, "kappa_ab": 0.03, "epsilon_k_ab": 1869.8446, + "na": 1.0, + "nb": 1.0, "viscosity": [ -1.3896, -2.58867, @@ -2689,6 +2739,8 @@ "mu": 1.6488, "kappa_ab": 0.03, "epsilon_k_ab": 1721.3869, + "na": 1.0, + "nb": 1.0, "viscosity": [ -1.3845, -1.1782, @@ -2781,6 +2833,8 @@ "mu": 1.5499, "kappa_ab": 0.03, "epsilon_k_ab": 1923.3059, + "na": 1.0, + "nb": 1.0, "viscosity": [ -0.9812, -2.6533, @@ -2806,6 +2860,8 @@ "mu": 1.7, "kappa_ab": 0.03, "epsilon_k_ab": 1835.0315, + "na": 1.0, + "nb": 1.0, "viscosity": [ -1.1188, -2.22062, @@ -2831,6 +2887,8 @@ "mu": 1.6668, "kappa_ab": 0.03, "epsilon_k_ab": 1713.2429, + "na": 1.0, + "nb": 1.0, "viscosity": [ -1.1043, -1.0389, @@ -2856,6 +2914,8 @@ "mu": 1.6399, "kappa_ab": 0.03, "epsilon_k_ab": 1806.0011, + "na": 1.0, + "nb": 1.0, "viscosity": [ -1.1043, -0.9823, @@ -2971,6 +3031,8 @@ "mu": 1.1692, "kappa_ab": 0.03, "epsilon_k_ab": 86.0572, + "na": 1.0, + "nb": 1.0, "viscosity": [ -0.7383, -1.7723, @@ -2996,6 +3058,8 @@ "mu": 1.6788, "kappa_ab": 0.03, "epsilon_k_ab": 2044.5298, + "na": 1.0, + "nb": 1.0, "viscosity": [ -0.8756, -1.47721, @@ -3021,6 +3085,8 @@ "mu": 1.6608, "kappa_ab": 0.03, "epsilon_k_ab": 1871.8788, + "na": 1.0, + "nb": 1.0, "viscosity": [ -0.847, -1.0735, @@ -3181,6 +3247,8 @@ "mu": 1.55, "kappa_ab": 0.03, "epsilon_k_ab": 2238.1838, + "na": 1.0, + "nb": 1.0, "viscosity": [ -1.7882, -2.8133, @@ -3294,6 +3362,8 @@ "mu": 1.6488, "kappa_ab": 0.03, "epsilon_k_ab": 2212.4364, + "na": 1.0, + "nb": 1.0, "viscosity": [ -1.7295, -2.8003, @@ -3341,6 +3411,8 @@ "mu": 1.6698, "kappa_ab": 0.03, "epsilon_k_ab": 1888.8681, + "na": 1.0, + "nb": 1.0, "viscosity": [ -1.6044, -3.44231, diff --git a/parameters/pcsaft/rehner2020.json b/parameters/pcsaft/rehner2020.json index 87d9d42e6..df289c35b 100644 --- a/parameters/pcsaft/rehner2020.json +++ b/parameters/pcsaft/rehner2020.json @@ -14,7 +14,9 @@ "sigma": 2.937523956051823, "epsilon_k": 272.02757407828676, "kappa_ab": 0.04448012165716923, - "epsilon_k_ab": 3125.3202766200056 + "epsilon_k_ab": 3125.3202766200056, + "na": 1.0, + "nb": 1.0 } }, { @@ -33,7 +35,8 @@ "epsilon_k": 238.31794591948915, "kappa_ab": 0.037807044304069184, "epsilon_k_ab": 2749.004567423258, - "na": 2 + "na": 2.0, + "nb": 1.0 } }, { @@ -52,8 +55,8 @@ "epsilon_k": 169.7751872692369, "kappa_ab": 0.13373757842938191, "epsilon_k_ab": 1772.0393059972052, - "na": 2, - "nb": 2 + "na": 2.0, + "nb": 2.0 } }, { @@ -72,7 +75,9 @@ "epsilon_k": 166.6147951235982, "mu": 1.6152087869692175, "kappa_ab": 0.09819448826630345, - "epsilon_k_ab": 2667.2518268470913 + "epsilon_k_ab": 2667.2518268470913, + "na": 1.0, + "nb": 1.0 } }, { @@ -92,7 +97,8 @@ "mu": 1.9373715367486852, "kappa_ab": 0.038235772849856756, "epsilon_k_ab": 2377.871413812318, - "na": 2 + "na": 2.0, + "nb": 1.0 } }, { @@ -112,8 +118,8 @@ "mu": 1.5050185176397637, "kappa_ab": 0.082906656369032, "epsilon_k_ab": 1784.1459137848506, - "na": 2, - "nb": 2 + "na": 2.0, + "nb": 2.0 } }, { @@ -132,8 +138,8 @@ "epsilon_k": 101.0845669390466, "kappa_ab": 0.11953491759977186, "epsilon_k_ab": 1834.845540500216, - "na": 2, - "nb": 2 + "na": 2.0, + "nb": 2.0 } }, { @@ -152,8 +158,8 @@ "epsilon_k": 124.58467252505255, "kappa_ab": 0.10067627754508424, "epsilon_k_ab": 1810.4338396651838, - "na": 2, - "nb": 2 + "na": 2.0, + "nb": 2.0 } }, { @@ -172,8 +178,8 @@ "epsilon_k": 144.37733584205245, "kappa_ab": 0.05810650285845608, "epsilon_k_ab": 1959.8294504365563, - "na": 2, - "nb": 2 + "na": 2.0, + "nb": 2.0 } }, { @@ -192,8 +198,8 @@ "epsilon_k": 163.6648110273199, "kappa_ab": 0.07216895211562263, "epsilon_k_ab": 1862.0561323915151, - "na": 2, - "nb": 2 + "na": 2.0, + "nb": 2.0 } }, { @@ -212,8 +218,8 @@ "epsilon_k": 179.754303082289, "kappa_ab": 0.0762520371571717, "epsilon_k_ab": 1824.6154575116475, - "na": 2, - "nb": 2 + "na": 2.0, + "nb": 2.0 } }, { @@ -231,7 +237,9 @@ "sigma": 3.7482827652225086, "epsilon_k": 270.86189939448803, "kappa_ab": 0.002566318648210027, - "epsilon_k_ab": 2778.9261790720657 + "epsilon_k_ab": 2778.9261790720657, + "na": 1.0, + "nb": 1.0 } }, { @@ -249,7 +257,9 @@ "sigma": 3.620930465025628, "epsilon_k": 252.7630330258869, "kappa_ab": 0.003060843493411822, - "epsilon_k_ab": 2700.937487295063 + "epsilon_k_ab": 2700.937487295063, + "na": 1.0, + "nb": 1.0 } }, { @@ -267,7 +277,9 @@ "sigma": 3.8281719743325215, "epsilon_k": 268.72917896983984, "kappa_ab": 0.0031104598923094376, - "epsilon_k_ab": 2804.2762044751535 + "epsilon_k_ab": 2804.2762044751535, + "na": 1.0, + "nb": 1.0 } }, { @@ -285,7 +297,9 @@ "sigma": 3.942380472821992, "epsilon_k": 278.68277806898686, "kappa_ab": 0.0014907242131781597, - "epsilon_k_ab": 3103.6530564736395 + "epsilon_k_ab": 3103.6530564736395, + "na": 1.0, + "nb": 1.0 } }, { @@ -303,7 +317,9 @@ "sigma": 4.042094734564122, "epsilon_k": 281.2912655243266, "kappa_ab": 0.002493036173908781, - "epsilon_k_ab": 3023.609567972825 + "epsilon_k_ab": 3023.609567972825, + "na": 1.0, + "nb": 1.0 } }, { @@ -321,7 +337,9 @@ "sigma": 3.977088533521953, "epsilon_k": 274.5543622578429, "kappa_ab": 0.0012672291257527345, - "epsilon_k_ab": 3223.361822388517 + "epsilon_k_ab": 3223.361822388517, + "na": 1.0, + "nb": 1.0 } }, { @@ -340,8 +358,8 @@ "epsilon_k": 136.16495641117777, "kappa_ab": 0.09762722816926367, "epsilon_k_ab": 1718.8359393009746, - "na": 2, - "nb": 2 + "na": 2.0, + "nb": 2.0 } }, { @@ -360,8 +378,8 @@ "epsilon_k": 149.97423748959227, "kappa_ab": 0.16261464675774262, "epsilon_k_ab": 1469.6277485887485, - "na": 2, - "nb": 2 + "na": 2.0, + "nb": 2.0 } }, { @@ -380,8 +398,8 @@ "epsilon_k": 159.24901703474268, "kappa_ab": 0.0787091926637352, "epsilon_k_ab": 1771.398524849797, - "na": 2, - "nb": 2 + "na": 2.0, + "nb": 2.0 } }, { @@ -400,8 +418,8 @@ "epsilon_k": 162.73433804276743, "kappa_ab": 0.12199296128478557, "epsilon_k_ab": 1516.0966359637514, - "na": 2, - "nb": 2 + "na": 2.0, + "nb": 2.0 } }, { @@ -420,8 +438,8 @@ "epsilon_k": 211.90069259634961, "kappa_ab": 0.041207808637705026, "epsilon_k_ab": 2550.2519514335218, - "na": 2, - "nb": 2 + "na": 2.0, + "nb": 2.0 } }, { @@ -440,8 +458,8 @@ "epsilon_k": 183.09717953416344, "kappa_ab": 0.08499573407848357, "epsilon_k_ab": 2268.8453113675146, - "na": 2, - "nb": 2 + "na": 2.0, + "nb": 2.0 } }, { @@ -460,8 +478,8 @@ "epsilon_k": 299.7826857254255, "kappa_ab": 0.019561736415810844, "epsilon_k_ab": 3069.7332849899476, - "na": 2, - "nb": 2 + "na": 2.0, + "nb": 2.0 } } ] \ No newline at end of file From fba05c232e482e552f143ef195fe04e49ac528ef Mon Sep 17 00:00:00 2001 From: Philipp Rehner Date: Mon, 3 Apr 2023 15:25:20 +0200 Subject: [PATCH 6/7] cleanup --- src/pcsaft/parameters.rs | 28 ---------------------------- src/pcsaft/python.rs | 4 ---- tests/gc_pcsaft/dft.rs | 26 -------------------------- 3 files changed, 58 deletions(-) diff --git a/src/pcsaft/parameters.rs b/src/pcsaft/parameters.rs index 2b4881f93..5cfef1564 100644 --- a/src/pcsaft/parameters.rs +++ b/src/pcsaft/parameters.rs @@ -7,7 +7,6 @@ use feos_core::parameter::{ }; use ndarray::{Array, Array1, Array2}; use num_dual::DualNum; -use num_traits::Zero; use quantity::si::{JOULE, KB, KELVIN}; use serde::{Deserialize, Serialize}; use std::collections::HashMap; @@ -517,33 +516,6 @@ impl PcSaftParameters { } } -impl std::fmt::Display for PcSaftParameters { - fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { - write!(f, "PcSaftParameters(")?; - write!(f, "\n\tmolarweight={}", self.molarweight)?; - write!(f, "\n\tm={}", self.m)?; - write!(f, "\n\tsigma={}", self.sigma)?; - write!(f, "\n\tepsilon_k={}", self.epsilon_k)?; - if !self.dipole_comp.is_empty() { - write!(f, "\n\tmu={}", self.mu)?; - } - if !self.quadpole_comp.is_empty() { - write!(f, "\n\tq={}", self.q)?; - } - // if !self.association.is_empty() { - // write!(f, "\n\tassociating={}", self.association.assoc_comp)?; - // write!(f, "\n\tkappa_ab={}", self.association.kappa_ab)?; - // write!(f, "\n\tepsilon_k_ab={}", self.association.epsilon_k_ab)?; - // write!(f, "\n\tna={}", self.association.na)?; - // write!(f, "\n\tnb={}", self.association.nb)?; - // } - if !self.k_ij.iter().all(|k| k.is_zero()) { - write!(f, "\n\tk_ij=\n{}", self.k_ij)?; - } - write!(f, "\n)") - } -} - #[cfg(test)] pub mod utils { use super::*; diff --git a/src/pcsaft/python.rs b/src/pcsaft/python.rs index 4da44bc6b..077e4d7e6 100644 --- a/src/pcsaft/python.rs +++ b/src/pcsaft/python.rs @@ -179,10 +179,6 @@ impl PyPcSaftParameters { fn _repr_markdown_(&self) -> String { self.0.to_markdown() } - - fn __repr__(&self) -> PyResult { - Ok(self.0.to_string()) - } } #[pymodule] diff --git a/tests/gc_pcsaft/dft.rs b/tests/gc_pcsaft/dft.rs index 836643501..cc5a76827 100644 --- a/tests/gc_pcsaft/dft.rs +++ b/tests/gc_pcsaft/dft.rs @@ -135,35 +135,9 @@ fn test_bulk_association() -> Result<(), Box> { let state_func = State::new_nvt(&func, t, v, &n)?; let p_eos = state_eos.pressure_contributions(); let p_func = state_func.pressure_contributions(); - // println!( - // "Equation of state: - // \tcomps: {} - // \tkappa_ab: {} - // \tepsilon_k_ab: {} - // \tna: {} - // \tnb: {}", - // eos_parameters.association.assoc_comp, - // eos_parameters.association.kappa_ab, - // eos_parameters.association.epsilon_k_ab, - // eos_parameters.association.na, - // eos_parameters.association.nb, - // ); for (s, x) in &p_eos { println!("{s:18}: {x:21.16}"); } - // println!( - // "\nHelmholtz energy functional: - // \tcomps: {} - // \tkappa_ab: {} - // \tepsilon_k_ab: {} - // \tna: {} - // \tnb: {}", - // func_parameters.association.assoc_comp, - // func_parameters.association.kappa_ab, - // func_parameters.association.epsilon_k_ab, - // func_parameters.association.na, - // func_parameters.association.nb, - // ); for (s, x) in &p_func { println!("{s:26}: {x:21.16}"); } From b44014403b06b05c433b4d4f6055c981fc1c4440 Mon Sep 17 00:00:00 2001 From: Philipp Rehner Date: Fri, 14 Apr 2023 13:18:24 +0200 Subject: [PATCH 7/7] add changelog --- CHANGELOG.md | 3 +++ docs/theory/models/association.md | 2 +- src/association/mod.rs | 2 +- src/association/python.rs | 1 + 4 files changed, 6 insertions(+), 2 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index d8740a33f..87df11d7f 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -5,6 +5,9 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/), and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html). ## [Unreleased] +### Changed +- Changed the internal implementation of the association contribution to accomodate more general association schemes. [#150](https://github.com/feos-org/feos/pull/150) +- To comply with the new association implementation, the default values of `na` and `nb` are now `0` rather than `1`. Parameter files have been adapted accordingly. [#150](https://github.com/feos-org/feos/pull/150) ## [0.4.2] - 2023-03-20 - Python only: Release the changes introduced in `feos-core` 0.4.1 and `feos-dft` 0.4.1. diff --git a/docs/theory/models/association.md b/docs/theory/models/association.md index a0c359329..9313bbe46 100644 --- a/docs/theory/models/association.md +++ b/docs/theory/models/association.md @@ -30,7 +30,7 @@ The analytic Hessian is $$H_{\alpha\beta}=-\frac{\delta_{\alpha\beta}}{X_\alpha^2}-\rho_\beta\Delta^{\alpha\beta}$$ -with the dirac delta function $\delta_{\alpha\beta}$. However, [Michelsen 2006](https://pubs.acs.org/doi/full/10.1021/ie060029x) shows that it is beneficial to use +with the Kronecker delta $\delta_{\alpha\beta}$. However, [Michelsen 2006](https://pubs.acs.org/doi/full/10.1021/ie060029x) shows that it is beneficial to use $$\hat{H}_{\alpha\beta}=-\frac{\delta_{\alpha\beta}}{X_\alpha}\left(1+\sum_\beta\rho_\beta X_\beta\Delta^{\alpha\beta}\right)-\rho_\beta\Delta^{\alpha\beta}$$ diff --git a/src/association/mod.rs b/src/association/mod.rs index 627cc0ebd..8dfcbc9de 100644 --- a/src/association/mod.rs +++ b/src/association/mod.rs @@ -77,7 +77,7 @@ impl fmt::Display for AssociationRecord { write!(f, ", na={}", self.na)?; } if self.nb > 0.0 { - write!(f, ", nb={})", self.nb)?; + write!(f, ", nb={}", self.nb)?; } if self.nc > 0.0 { write!(f, ", nc={}", self.nc)?; diff --git a/src/association/python.rs b/src/association/python.rs index ac3211bbc..9744f5552 100644 --- a/src/association/python.rs +++ b/src/association/python.rs @@ -11,6 +11,7 @@ pub struct PyAssociationRecord(pub AssociationRecord); #[pymethods] impl PyAssociationRecord { #[new] + #[pyo3(signature = (kappa_ab, epsilon_k_ab, na=0.0, nb=0.0, nc=0.0))] fn new(kappa_ab: f64, epsilon_k_ab: f64, na: f64, nb: f64, nc: f64) -> Self { Self(AssociationRecord::new(kappa_ab, epsilon_k_ab, na, nb, nc)) }