Help Needed Translating ENKI MELTS Examples into Python Code

[liushanke]

Hello everyone,

I am currently learning to use alphaMELTS by following the examples provided on the ENKI website. Two specific examples demonstrate modeling systematic variation for the same initial conditions:

Fractionation example: 1. https://enki-portal.gitlab.io/ThermoEngine/MELTS-v1.0.2-fractionation.html
Equilibrium example: 2. https://enki-portal.gitlab.io/ThermoEngine/MELTS-v1.0.2-equilibrium.html

To deepen my understanding, I tried translating these examples into Python scripts. However, I encountered issues running the Python code. Below is one of my attempts:
~~~~~~~~
from thermoengine import equilibrate
import matplotlib.pyplot as plt
import numpy as np

# Initialize MELTS v1.0.2 instance
melts = equilibrate.MELTSmodel()

# Set initial bulk composition (wt% or grams)
feasible = melts.set_bulk_composition({
    'SiO2': 77.5, 'TiO2': 0.08, 'Al2O3': 12.5, 'Fe2O3': 0.207,
    'Cr2O3': 0.0, 'FeO': 0.473, 'MnO': 0.0, 'MgO': 0.03,
    'NiO': 0.0, 'CoO': 0.0, 'CaO': 0.43, 'Na2O': 3.98,
    'K2O': 4.88, 'P2O5': 0.0, 'H2O': 5.5
})

# Optional: Suppress unwanted phases
melts.set_phase_inclusion_status({'Nepheline': False, 'OrthoOxide': False})

# Compute equilibrium at specific T (°C) and P (MPa)
output = melts.equilibrate_tp(760.0, 175.0, initialize=True)
(status, t, p, xmlout) = output

print(status, t, p)

# Display summary output
melts.output_summary(xmlout)

# List thermodynamic properties for the system
props = melts.get_list_of_properties()
for prop in props:
    print(f"{prop:<20s} {melts.get_property_of_phase(xmlout, 'System', prop):13.6e} {melts.get_units_of_property(prop):<10s}")

# List chemical affinities and undersaturated phase compositions
affinity_dict = melts.get_dictionary_of_affinities(xmlout, sort=True)
for phase in affinity_dict:
    affinity, formulae = affinity_dict[phase]
    print(f"{phase:<20s} {affinity:10.2f} {formulae:<60s}")

# Set up for sequential calculations and plotting
number_of_steps = 20
t_increment = -1.0
p_increment = 0.0

plotPhases = ['Liquid', 'Sanidine', 'Plagioclase', 'Quartz', 'Water']
plotColors = ['ro', 'bo', 'go', 'co', 'mo']

# Initialize data arrays for plotting
n = len(plotPhases)
xPlot = np.zeros(number_of_steps + 1)
yPlot = np.zeros((n, number_of_steps + 1))

xPlot

• = t

for i in range(n):
    yPlot

• = melts.get_property_of_phase(xmlout, plotPhases)

# Create a plot for the phase data
plt.ion()
fig, ax = plt.subplots()
ax.set_xlim([min(t, t + t_increment * number_of_steps), max(t, t + t_increment * number_of_steps)])
ax.set_ylim([0.0, 100.0])

graphs = [ax.plot(xPlot, yPlot, plotColors)

• for i in range(n)]

ax.legend(graphs, plotPhases, loc='upper left')

# Perform sequential calculations and update the plot
for i in range(1, number_of_steps + 1):
    output = melts.equilibrate_tp(t + i * t_increment, p + i * p_increment)
    status, t, p, xmlout = output

    print(f"{status:<30s} {t:8.2f} {p:8.2f}")

    xPlot = t
    for j in range(n):
        yPlot[j] = melts.get_property_of_phase(xmlout, plotPhases[j])

    for j, graph in enumerate(graphs):
        graph.set_xdata(xPlot)
        graph.set_ydata(yPlot[j])

    fig.canvas.draw()

# Show the plot (important for standard Python scripts)
plt.show()
~~~~~~~~~~~

I know these examples are designed to run on the ENKI server. However, I encountered difficulties translating them into a Python script, particularly because I couldn't find the thermoengine  module referenced in the original examples. This has been a blocking issue, preventing me from running the code successfully on my local machine.

Could someone kindly help me verify the code or guide me in correctly setting it up and running it in Python, so I can run it smoothly, just like the example from this web link (https://github.com/magmasource/alphaMELTS/blob/main/examples/py/tutorial.py), using the "alphamelts-py-2.3.1-win64" package on my personal computer?

I would greatly appreciate any advice or insights.

Thank you in advance for your help!

Best regards,
Shanke

[Paula]

Hi Shanke,

If you want to run the ENKI examples you need to install Thermoengine locally or run it in Docker. If you are on Windows then Docker is your only option. There's more information here: https://gitlab.com/ENKI-portal/ThermoEngine and you should be able to get help on the ENKI Slack channel (enki-portal.slack.com). It's the free version of Slack so you can't see old posts, unfortunately, but queries usually get answered by community members quickly.

However, this will not help with running alphaMELTS for Python. The underlying thermodynamic models are the same but it's a different implementation (forked from Mark Ghiorso's xMELTS code) and a different interface, although both use Python. There are alphaMELTS examples from past workshops, including fractionation crystallization and degassing. For example, see here: https://magmasource.caltech.edu/gitlist/Workshops.git/tree/main/workshop_files/matlab_python_rc12/

There isn't an isobaric equilibrium (batch) crystallization example but you can make one by replacing

Code Select Expand

ptpath.engine.setSystemProperties(["Log fO2 Path: -1FMQ", "Mode: Fractionate Solids"])with

Code Select Expand

ptpath.engine.setSystemProperties(["Log fO2 Path: -1FMQ"])
Or have a look at PetThermoTools, which builds on alphaMELTS for Python, and includes documentation and examples for common calculations like modeling liquid lines of descent. See https://github.com/gleesonm1/PetThermoTools

Best wishes,
Paula

[liushanke]

Hi Paul,

Thank you for your kind response—it has been very helpful. However, I was unable to access the link: "https://magmasource.caltech.edu/gitlist/Workshops.git/tree/main/workshop_files/matlab_python_rc12/".

I realize that fully understanding the different ways to run alphaMELTS is still a long journey ahead.

Thank you again for your help!

Sincerely,
Shanke

[Paula]

It should be ok now.

Paula

[liushanke]

Yes，thanks for your response.