From 791e94bcb5279291c89b2165aac829142b9b0aaf Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Sun, 4 Jan 2026 20:16:31 +0000 Subject: [PATCH 1/5] Initial plan From bdc5bcdfb76c177bc046b171e5ac3e377d6b1440 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Sun, 4 Jan 2026 20:20:32 +0000 Subject: [PATCH 2/5] Implement correct Birch-Murnaghan equation and add Murnaghan EOS Co-authored-by: jochym <5993422+jochym@users.noreply.github.com> --- docs/source/lib-usage.rst | 12 +++++----- elastic/__init__.py | 2 +- elastic/elastic.py | 46 +++++++++++++++++++++++++++++++++++---- 3 files changed, 49 insertions(+), 11 deletions(-) diff --git a/docs/source/lib-usage.rst b/docs/source/lib-usage.rst index 7f55305..d9d889f 100644 --- a/docs/source/lib-usage.rst +++ b/docs/source/lib-usage.rst @@ -127,16 +127,16 @@ Birch-Murnaghan Equation of State Let us now use the tools provided by the modules to calculate equation of state for the crystal and verify it by plotting the data points -against fitted EOS curve. The EOS used by the module is a well -established Birch-Murnaghan formula (P - pressure, V - volume, B - -parameters): +against fitted EOS curve. The EOS used by the module is the well +established third-order Birch-Murnaghan formula (P - pressure, V - volume, +B - parameters): .. math:: - P(V)= \frac{B_0}{B'_0}\left[ - \left({\frac{V}{V_0}}\right)^{-B'_0} - 1 - \right] + P(V) = \frac{3B_0}{2} \left[\left(\frac{V_0}{V}\right)^{7/3} - + \left(\frac{V_0}{V}\right)^{5/3}\right] + \left\{1 + \frac{3}{4}(B'_0 - 4)\left[\left(\frac{V_0}{V}\right)^{2/3} - 1\right]\right\} Now we repeat the setup and optimization procedure from the example 1 above but using a new Crystal class (see above we skip this part for diff --git a/elastic/__init__.py b/elastic/__init__.py index 63ff5fc..8f168cc 100644 --- a/elastic/__init__.py +++ b/elastic/__init__.py @@ -65,7 +65,7 @@ from __future__ import print_function, division, absolute_import from .elastic import get_BM_EOS, get_elastic_tensor from .elastic import get_elementary_deformations, scan_volumes -from .elastic import get_pressure, get_strain, BMEOS +from .elastic import get_pressure, get_strain, BMEOS, MurnaghanEOS # To reach "elastic.__version__" attribute in other programs from ._version import __version__ diff --git a/elastic/elastic.py b/elastic/elastic.py index 50087a2..6a40236 100644 --- a/elastic/elastic.py +++ b/elastic/elastic.py @@ -80,6 +80,44 @@ def BMEOS(v, v0, b0, b0p): + ''' + Third-order Birch-Murnaghan Equation of State. + + Returns pressure as a function of volume for the Birch-Murnaghan EOS. + + .. math:: + P(V) = \\frac{3B_0}{2} \\left[\\left(\\frac{V_0}{V}\\right)^{7/3} - + \\left(\\frac{V_0}{V}\\right)^{5/3}\\right] + \\left\\{1 + \\frac{3}{4}(B'_0 - 4)\\left[\\left(\\frac{V_0}{V}\\right)^{2/3} - 1\\right]\\right\\} + + :param v: volume + :param v0: equilibrium volume V_0 + :param b0: bulk modulus B_0 + :param b0p: pressure derivative of bulk modulus B'_0 + + :returns: pressure P(V) + ''' + x = (v0/v)**(1.0/3.0) + return (3.0*b0/2.0) * (x**7 - x**5) * (1.0 + (3.0/4.0)*(b0p - 4.0)*(x**2 - 1.0)) + + +def MurnaghanEOS(v, v0, b0, b0p): + ''' + Murnaghan Equation of State. + + Returns pressure as a function of volume for the Murnaghan EOS. + This was the equation previously (incorrectly) used for BMEOS. + + .. math:: + P(V) = \\frac{B_0}{B'_0}\\left[\\left(\\frac{V_0}{V}\\right)^{B'_0} - 1\\right] + + :param v: volume + :param v0: equilibrium volume V_0 + :param b0: bulk modulus B_0 + :param b0p: pressure derivative of bulk modulus B'_0 + + :returns: pressure P(V) + ''' return (b0/b0p)*(pow(v0/v, b0p) - 1) @@ -385,12 +423,12 @@ def get_pressure(s): def get_BM_EOS(cryst, systems): """Calculate Birch-Murnaghan Equation of State for the crystal. - The B-M equation of state is defined by: + The third-order Birch-Murnaghan equation of state is defined by: .. math:: - P(V)= \\frac{B_0}{B'_0}\\left[ - \\left({\\frac{V}{V_0}}\\right)^{-B'_0} - 1 - \\right] + P(V) = \\frac{3B_0}{2} \\left[\\left(\\frac{V_0}{V}\\right)^{7/3} - + \\left(\\frac{V_0}{V}\\right)^{5/3}\\right] + \\left\\{1 + \\frac{3}{4}(B'_0 - 4)\\left[\\left(\\frac{V_0}{V}\\right)^{2/3} - 1\\right]\\right\\} It's coefficients are estimated using n single-point structures ganerated from the crystal (cryst) by the scan_volumes function between two relative From ce435677c2ddfd270d86fd9baf246aafbfa90ddf Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Sun, 4 Jan 2026 20:21:46 +0000 Subject: [PATCH 3/5] Update scan_pressures docstring to clarify Murnaghan approximation Co-authored-by: jochym <5993422+jochym@users.noreply.github.com> --- elastic/elastic.py | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/elastic/elastic.py b/elastic/elastic.py index 6a40236..666a60f 100644 --- a/elastic/elastic.py +++ b/elastic/elastic.py @@ -602,10 +602,16 @@ def get_elastic_tensor(cryst, systems): def scan_pressures(cryst, lo, hi, n=5, eos=None): ''' Scan the pressure axis from lo to hi (inclusive) - using B-M EOS as the volume predictor. + using inverse Murnaghan EOS as the volume predictor. + + Note: This function uses the inverse Murnaghan equation as an approximation + for computational efficiency. The inverse of the full 3rd-order Birch-Murnaghan + equation requires numerical solution and is more complex. For moderate pressure + ranges, the Murnaghan approximation provides reasonable volume estimates. + Pressure (lo, hi) in GPa ''' - # Inverse B-M EOS to get volumes from pressures + # Inverse Murnaghan EOS to get volumes from pressures # This will work only in limited pressure range p>-B/B'. # Warning! Relative, the V0 prefactor is removed. def invbmeos(b, bp, x): From e190ecc4b8b53927b54a7e66853210d83e44cf1d Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Sun, 4 Jan 2026 20:22:46 +0000 Subject: [PATCH 4/5] Fix typo: ganerated -> generated Co-authored-by: jochym <5993422+jochym@users.noreply.github.com> --- elastic/elastic.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/elastic/elastic.py b/elastic/elastic.py index 666a60f..a7c695a 100644 --- a/elastic/elastic.py +++ b/elastic/elastic.py @@ -430,7 +430,7 @@ def get_BM_EOS(cryst, systems): \\left(\\frac{V_0}{V}\\right)^{5/3}\\right] \\left\\{1 + \\frac{3}{4}(B'_0 - 4)\\left[\\left(\\frac{V_0}{V}\\right)^{2/3} - 1\\right]\\right\\} - It's coefficients are estimated using n single-point structures ganerated + It's coefficients are estimated using n single-point structures generated from the crystal (cryst) by the scan_volumes function between two relative volumes. The BM EOS is fitted to the computed points by least squares method. The returned value is a list of fitted From fd4faaa7f423e86ba0958c8074d2afb1fbd507dd Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Mon, 5 Jan 2026 04:25:20 +0000 Subject: [PATCH 5/5] Rename invbmeos to invmurneos for consistent naming Co-authored-by: jochym <5993422+jochym@users.noreply.github.com> --- elastic/elastic.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/elastic/elastic.py b/elastic/elastic.py index a7c695a..2102cc4 100644 --- a/elastic/elastic.py +++ b/elastic/elastic.py @@ -614,18 +614,18 @@ def scan_pressures(cryst, lo, hi, n=5, eos=None): # Inverse Murnaghan EOS to get volumes from pressures # This will work only in limited pressure range p>-B/B'. # Warning! Relative, the V0 prefactor is removed. - def invbmeos(b, bp, x): + def invmurneos(b, bp, x): return array([pow(b/(bp*xv+b), 1/(3*bp)) for xv in x]) if eos is None: raise RuntimeError('Required EOS data missing') # Limit negative pressures to 90% of the singularity value. - # Beyond this B-M EOS is bound to be wrong anyway. + # Beyond this Murnaghan EOS is bound to be wrong anyway. lo = max(lo, -0.9*eos[1]/eos[2]) - scale = (eos[0]/cryst.get_volume())*invbmeos(eos[1], eos[2], - linspace(lo, hi, num=n)) + scale = (eos[0]/cryst.get_volume())*invmurneos(eos[1], eos[2], + linspace(lo, hi, num=n)) # print(scale) uc = cryst.get_cell() systems = [Atoms(cryst) for s in scale]