Among the opened MPDS data there is a considerable number of the diagrams cell parameters vs. temperature VT
, as well as cell parameters vs. pressure pV
. They provide us the link to the thermodynamic properties. Let's explore this link.
We will calculate an isothermal bulk modulus, as well as its pressure derivative, using the pV
-data from the MPDS API and fitting to the equation of state (EoS).
Important! Before you proceed: the notebooks running at the third-party servers are not secure. Using this notebook assumes you authenticate at the MPDS server with your own API key. Please run this notebook only if you have an open-access account (i.e. an access section of your MPDS account reads: Programmatic data access: only open data
).
Please do not run this notebook at the third-party servers if you have an elevated API access to the MPDS, since there's a nonzero probability of key leakage!
Be sure to always invalidate (revoke) your API key at your MPDS account after using the notebooks.
Now let's proceed with the authentication part. First, apply for an MPDS account, if you have none. Then copy your API key, run the next cell, paste the key in the appeared prompt input, and hit Enter.
import os, getpass
os.environ['MPDS_KEY'] = getpass.getpass()
OK, now you may talk to the MPDS server programmatically from this notebook on your behalf.
For fitting we will use a pytheos
pure-Python library written by Prof. Sang-Heon Dan Shim. (Note there's a wide range of the other EoS fitting tools!)
!pip install pytheos
Make sure pytheos
works, using this simple example:
import numpy as np
from pytheos import BM3Model
v0, k0, k0p = 10, 200, 4
exp_bm3 = BM3Model()
v_data = v0 * np.linspace(0.6, 1, 20)
fit_params = exp_bm3.make_params(v0=v0, k0=k0, k0p=k0p)
p_bm3 = exp_bm3.eval(fit_params, v=v_data)
p_data = p_bm3 + np.random.normal(0.0, 2, 20)
fit_params['v0'].vary = False
fitresult_bm3 = exp_bm3.fit(p_data, fit_params, v=v_data, weights=None)
print(fitresult_bm3.fit_report())
If there were no errors at the previous step, let's continue.
We will download the available isothermal bulk moduli and their pressure derivatives from the MPDS.
!pip install mpds_client
import warnings
import pandas as pd
import numpy as np
from numpy.linalg import det
from ase.geometry import cellpar_to_cell
from mpds_client import MPDSDataRetrieval
Note, if your API key isn't valid, the API returns an HTTP error 403
.
client = MPDSDataRetrieval()
dfrm_k0p = client.get_dataframe({"classes": "binary", "elements": "O", "props": "pressure derivative of isothermal bulk modulus"})
dfrm_k0p = dfrm_k0p[np.isfinite(dfrm_k0p['Phase'])] # only data for the existing distinct phases
avg_k0p = dfrm_k0p.groupby('Phase')['Value'].median().to_frame().reset_index().rename(columns={'Value': 'avg_k0p'})
dfrm_k0 = client.get_dataframe({"props": "isothermal bulk modulus"}, phases=set(dfrm_k0p['Phase'].tolist()))
avg_k0 = dfrm_k0.groupby('Phase')['Value'].median().to_frame().reset_index().rename(columns={'Value': 'avg_k0'})
avg_k0 = avg_k0.merge(avg_k0p, how='inner', on='Phase')
Let's fit the isothermal bulk moduli and their pressure derivatives from the MPDS pV
-data and compare with the experimental values:
def get_cell_volume(a, b, c, alpha, beta, gamma):
'''
Calculate V from cell parameters.
NB ab_normal and a_direction are standard.
'''
return abs(
det(
cellpar_to_cell([a, b, c, alpha, beta, gamma])
)
)
pvts = {}
for matrix in client.get_data(
{"props": "cell parameters - pressure diagram"},
phases=set(dfrm_k0p['Phase'].tolist()), # only those phases we have experimental bulk modulus
fields={} # all fields
):
try: decks = matrix['sample']['measurement'][0]['property']['matrix']
except KeyError:
warnings.warn('Error: there is no expected property in the entry %s' % matrix)
continue
p_all, v_all, t_all = [], [], []
for deck in decks:
p, t, v = deck[0], deck[1], get_cell_volume(*deck[2:])
p_all.append(p)
v_all.append(v)
t_all.append(t)
# TODO here in principle we should do something smarter
# than just omitting the data for the same phase
if matrix['sample']['material']['phase_id'] in pvts and len(pvts[ matrix['sample']['material']['phase_id'] ][0]) > len(p_all):
warnings.warn('Skipping entry %s' % matrix['sample']['material']['entry'])
continue
pvts[ matrix['sample']['material']['phase_id'] ] = [matrix['sample']['material']['entry'], p_all, v_all, t_all]
pvts = [[key] + value for key, value in pvts.items()]
pvts = pd.DataFrame(pvts, columns=['Phase', 'Entry', 'P', 'V', 'T'])
print(pvts)
dfrm = avg_k0.merge(pvts, how='inner', on='Phase')
for n, system in dfrm.iterrows():
params = exp_bm3.make_params(v0=system['V'][0], k0=system['avg_k0'], k0p=system['avg_k0p'])
fitresult_bm3 = exp_bm3.fit(system['P'], params, v=system['V'], weights=None)
if abs(fitresult_bm3.params['k0'].value - system['avg_k0']) > 50:
# show the discrepancies, if any
print("*" * 30 + " Distinct phase https://mpds.io/#phase_id/%s " % system['Phase'] + "*" * 30)
print("BM_fit: %4.1f BM_exp: %4.1f" % (fitresult_bm3.params['k0'].value, system['avg_k0']))
print("BM0p_fit: %1.2f BM0p_exp: %1.2f" % (fitresult_bm3.params['k0p'].value, system['avg_k0p']))
Were you able to follow everything? Please, try to answer:
PS don't forget to invalidate (revoke) your API key.