# Third-order high-temperature Birch-Murnaghan equation of state

Solve for volume, and for integral of V dP, for solid phases at high pressure and temperature.

#pragma version=1.3

// set bracketing values and tolerance for volume
// absolute brackets will be (1 +/- kBracketFactor) * V0
// should be okay for EOS for solids
static constant kBracketFactor=0.1
static constant kTol=1e-5

// V = volume at P, T
// Tref = reference temperature
// V0 = volume at Tref
// thermal expansion coefficient = c1 + c2*T + c3/T + c4/T^2 + c5/T^0.5
// K0 is the bulk modulus at Tref
// Kprime is the pressure deriviative of bulk modulus
// dKdT is temperature derivative of bulk modulus

// For room temperature EOS, set a, b, c, and dKdT to zero.

// 3rd order high-temperature Birch-Murnaghan equation of state, returns P
threadsafe function HTBM(V, T, Tref, V0, c1, c2, c3, c4, c5, K0, Kprime, dKdT)
variable V, T, Tref, V0, c1, c2, c3, c4, c5, K0, Kprime, dKdT

variable KT=K0+dKdT*(T-Tref)
variable V0T=V0*exp(c1*(T-Tref)+0.5*c2*(T^2-Tref^2)+c3*ln(T/Tref)-c4*(T^-1-Tref^-1)+2*c5*(T^0.5-Tref^0.5))
return 3/2*KT*((V0T/V)^(7/3) - (V0T/V)^(5/3))*(1 + 3/4*(Kprime-4)*((V0T/V)^(2/3) - 1))
end

// wrapper function for FindRoots
threadsafe function HTBM_wrapper(wave w, variable V)
// w contains V0, T, Tref, a, b, c, dKdT, K0, Kprime
return HTBM(V, w, w, w, w, w, w, w, w, w, w, w)
end

// solve for V that satisfies P=pressure
threadsafe function solveHTBM(pressure, temperature, Tref, V0, c1, c2, c3, c4, c5, K0, Kprime, dKdT, [limH,limL])
variable pressure, temperature, Tref, V0, c1, c2, c3, c4, c5, K0, Kprime, dKdT, limL, limH
limL = ParamIsDefault(limL) ? (1-kBracketFactor)*V0 : limL
limH = ParamIsDefault(limH) ? (1+kBracketFactor)*V0 : limH

Make /free/D w={temperature, Tref, V0, c1, c2, c3, c4, c5, K0, Kprime, dKdT}
FindRoots /H=(limH)/L=(limL) /Q /T=(kTol)/Z=(pressure) HTBM_wrapper, w
return V_flag==0 ? V_Root : NaN
end

// calculate integral of V dP from P=P0 to P=pressure at T=temperature
threadsafe function integrateVolume(P0, pressure, temperature, Tref, V0, c1, c2, c3, c4, c5, K0, Kprime, dKdT)
variable P0, pressure, temperature, Tref, V0, c1, c2, c3, c4, c5, K0, Kprime, dKdT

variable volume=solveHTBM(pressure, temperature, Tref, V0, c1, c2, c3, c4, c5, K0, Kprime, dKdT)
return (pressure-P0)*volume-integrateP(volume, temperature, Tref, V0, c1, c2, c3, c4, c5, K0, Kprime, dKdT)
end

// isothermal integral of 3rd order Birch-Murnaghan wrt V from V0 to V at T=T
threadsafe function integrateP(V, T, Tref, V0, c1, c2, c3, c4, c5, K0, Kprime, dKdT)
variable V, T, Tref, V0, c1, c2, c3, c4, c5, K0, Kprime, dKdT

variable KT=K0+dKdT*(T-Tref)
variable V0T=V0*exp(c1*(T-Tref)+0.5*c2*(T^2-Tref^2)+c3*ln(T/Tref)-c4*(T^-1-Tref^-1)+2*c5*(T^0.5-Tref^0.5))

variable integral = 0
integral = -(Kprime-4)*V0T^3*V^-2 - (2-3*(Kprime-4))*V0T^(7/3)*V^(-4/3) - (3*(Kprime-4)-4)*V0T^(5/3)*V^(-2/3)
integral += V0T*((Kprime-4) + (2-3*(Kprime-4)) + (3*(Kprime-4)-4))
integral *= 9/16*KT
return integral
end

Tony: Do you have reference values to create an example? I am especially interested from the academic side, i.e. for demonstration in either a physical chemistry or a materials science as a companion to the isotherms that are created for gases (https://www.wavemetrics.com/project/vdWIsothermPlots).

Well, yes, I should be able to provide plenty of examples.

I posted this after digging out my EOS-solving code to make a calculation that involved iron at high pressure. Of course iron tends to sit in the centre of planets where the VdP contribution to enthalpy is important...

I used values from here:

https://www.sciencedirect.com/science/article/abs/pii/S0012821X1300294X

And a more general collection of reference values can be found here:

https://agupubs.onlinelibrary.wiley.com/doi/book/10.1029/RF002

Incidentally, if you're interested in EOS for real fluids, I have code for plenty of EOS for water, including IAPWS, and some other formulations that are optimised for high pressure. I'd be happy to share. Forum Support Gallery