<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Hello, I am Richard Otis from Prof. Zi-Kui Liu's research group.
Tonight I present pycalphad, a free and open-source library for CALPHAD computations. It works on Mac, Windows, and Linux.
Phase Equilibria Computation
Al-Fe (M.Seiersten et al., 1991)
Step1:  Energy Surfaces of Al-Ni (N. Dupin et al., 2001)
Step2:  Al-Zn (S. Mey et al., 1993)
Next, we show values from Mey et al.'s 1993 assessment of the Al-Zn system. We have not done any fitting here; we are just showing how arbitrary values for the parameters may be specified and easily used to plot the phase diagram. Note the optional parameters keyword argument specified for each call to Model. This feature is what makes it easy to perform arbitrary variable substitutions in our thermodynamic models.
Step4:  Prototyping New Models (Advanced)
What if we have a new kind of functional form for the Gibbs energy we would like to experiment with?
<ASSISTANT_TASK:>
Python Code:
%matplotlib inline
import matplotlib.pyplot as plt
from pycalphad import Database, binplot
db_alfe = Database('alfe_sei.TDB')
my_phases_alfe = ['LIQUID', 'B2_BCC', 'FCC_A1', 'HCP_A3', 'AL5FE2', 'AL2FE', 'AL13FE4', 'AL5FE4']
fig = plt.figure(figsize=(9,6))
pdens = [{'B2_BCC': 20000}, 2000]
%time binplot(db_alfe, ['AL', 'FE', 'VA'] , my_phases_alfe, 'X(AL)', 300, 2000, steps=100, \
              pdens=pdens, ax=fig.gca(), mode='numexpr')
plt.show()
Explanation: Hello, I am Richard Otis from Prof. Zi-Kui Liu's research group.
Tonight I present pycalphad, a free and open-source library for CALPHAD computations. It works on Mac, Windows, and Linux.
Phase Equilibria Computation
Al-Fe (M.Seiersten et al., 1991)
End of explanation
%matplotlib inline
from pycalphad import energy_surf, Database
from pycalphad.plot.utils import phase_legend
import numpy as np
import matplotlib.pyplot as plt
db_alni = Database('NI_AL_DUPIN_2001.TDB')
my_phases_alni = ['LIQUID', 'FCC_L12', 'BCC_B2', 'AL3NI5', 'AL3NI2', 'AL3NI1']
fig = plt.figure(figsize=(9,6))
ax = fig.gca()
%time nrg_df = energy_surf(db_alni, ['AL', 'NI', 'VA'], my_phases_alni, T=1000.0)
legend_handles, colorlist = phase_legend(my_phases_alni)
group_df = nrg_df.groupby('Phase')
for phase, phase_df in group_df:
    ax = phase_df.plot(kind='scatter', x='X(AL)', y='GM', ax=ax, color=colorlist[phase.upper()], s=1, xlim=(0,1), ylim=(-11e4, -3e4))
ax = ax.legend(handles=legend_handles, loc='center left', bbox_to_anchor=(1, 0.6))
plt.title('Al-Ni Gibbs Energy at 1000 K', fontsize=20)
plt.show()
nrg_df
nrg_df.to_excel('alni_energies.xlsx')
Explanation: Energy Surfaces of Al-Ni (N. Dupin et al., 2001)
End of explanation
from pycalphad import Model, Database, binplot
db_alzn = Database('alzn_noparams.tdb')
assessed_parameters = {
    'LIQALZN0H': 10465.3,
    'LIQALZN0S': 3.39259,
    'FCCALZN0H': 7297.5,
    'FCCALZN0S': -0.47512,
    'FCCALZN1H': 6612.9,
    'FCCALZN1S': 4.5911,
    'FCCALZN2H': -3097.2,
    'FCCALZN2S': -3.306635,
    'HCPALZN0H': 18821.0,
    'HCPALZN0S': 8.95255,
    'HCPALZN3H': -702.8,
    'HCPALZN3S': 0.0
    }
# Initialize phase models
assessed_models = {
    'LIQUID': Model(db_alzn, ['AL', 'ZN'], 'LIQUID', parameters=assessed_parameters),
    'FCC_A1': Model(db_alzn, ['AL', 'ZN', 'VA'], 'FCC_A1', parameters=assessed_parameters),
    'HCP_A3': Model(db_alzn, ['AL', 'ZN', 'VA'], 'HCP_A3', parameters=assessed_parameters)
    }
fig = plt.figure(figsize=(9,6))
ax = binplot(db_alzn, ['AL', 'ZN', 'VA'] , ['LIQUID', 'FCC_A1', 'HCP_A3'], \
                   'X(ZN)', 300, 1000, ax=plt.gca(), model=assessed_models)
ax = ax.set_title('Al-Zn', fontsize=20)
blank_parameters = {
    'LIQALZN0H': 0.0,
    'LIQALZN0S': 0.0,
    'FCCALZN0H': 0.0,
    'FCCALZN0S': 0.0,
    'FCCALZN1H': 0.0,
    'FCCALZN1S': 0.0,
    'FCCALZN2H': 0.0,
    'FCCALZN2S': 0.0,
    'HCPALZN0H': 0.0,
    'HCPALZN0S': 0.0,
    'HCPALZN3H': 0.0,
    'HCPALZN3S': 0.0
    }
models = {
    'LIQUID': Model(db_alzn, ['AL', 'ZN'], 'LIQUID', parameters=blank_parameters),
    'FCC_A1': Model(db_alzn, ['AL', 'ZN', 'VA'], 'FCC_A1', parameters=blank_parameters),
    'HCP_A3': Model(db_alzn, ['AL', 'ZN', 'VA'], 'HCP_A3', parameters=blank_parameters)
    }
fig = plt.figure(figsize=(9,6))
ax = binplot(db_alzn, ['AL', 'ZN', 'VA'] , ['LIQUID', 'FCC_A1', 'HCP_A3'], \
             'X(ZN)', 300, 1000, ax=plt.gca(), model=models)
ax = ax.set_title('Al-Zn with endmembers only', fontsize=20)
Explanation: Al-Zn (S. Mey et al., 1993)
Next, we show values from Mey et al.'s 1993 assessment of the Al-Zn system. We have not done any fitting here; we are just showing how arbitrary values for the parameters may be specified and easily used to plot the phase diagram. Note the optional parameters keyword argument specified for each call to Model. This feature is what makes it easy to perform arbitrary variable substitutions in our thermodynamic models.
End of explanation
from sympy import Piecewise, log
from tinydb import where
import pycalphad.variables as v
class MyMagneticModel(Model):
    def __init__(self, db, comps, phase, **kwargs):
        super().__init__(db, comps, phase, **kwargs)
        self.models['mag'] = self.my_magnetic_energy(db.phases[phase], db.search)
    def my_magnetic_energy(self, phase, param_search):
        
        Return the energy from magnetic ordering in symbolic form.
        The implemented model is the Inden-Hillert-Jarl formulation.
        The approach follows from the background section of W. Xiong, 2011.
        
        Note that this is a demonstration model; use pycalphad's builtin
        implementation as a reference.
        
        site_ratio_normalization = self._site_ratio_normalization(phase)
        # define basic variables
        afm_factor = phase.model_hints['ihj_magnetic_afm_factor']
        bm_param_query = (
            (where('phase_name') == phase.name) & \
            (where('parameter_type') == 'BMAGN') & \
            (where('constituent_array').test(self._array_validity))
        )
        tc_param_query = (
            (where('phase_name') == phase.name) & \
            (where('parameter_type') == 'TC') & \
            (where('constituent_array').test(self._array_validity))
        )
        mean_magnetic_moment = \
            self.redlich_kister_sum(phase, param_search, bm_param_query)
        beta = Piecewise(
            (mean_magnetic_moment, mean_magnetic_moment > 0),
            (mean_magnetic_moment/afm_factor, mean_magnetic_moment <= 0)
            )
        curie_temp = \
            self.redlich_kister_sum(phase, param_search, tc_param_query)
        tc = Piecewise(
            (curie_temp, curie_temp > 0),
            (curie_temp/afm_factor, curie_temp <= 0)
            )
        tau = Piecewise(
            (v.T / tc, tc > 0),
            (10000., tc == 0) # tau is 'infinity'
            )
        # define model parameters
        p = phase.model_hints['ihj_magnetic_structure_factor']
        A = 518/1125 + (11692/15975)*(1/p - 1)
        # factor when tau < 1
        sub_tau = 1 - (1/A) * ((79/(140*p))*(tau**(-1)) + (474/497)*(1/p - 1) \
            * ((tau**3)/6 + (tau**9)/135 + (tau**15)/600)
                              )
        # factor when tau >= 1
        super_tau = -(1/A) * ((tau**-5)/10 + (tau**-15)/315 + (tau**-25)/1500)
        expr_cond_pairs = [(sub_tau, tau < 1),
                           (super_tau, tau >= 1)
                          ]
        g_term = Piecewise(*expr_cond_pairs)
        return 1e6 * v.R * v.T * log(beta+1) * \
            g_term / site_ratio_normalization
%matplotlib inline
import matplotlib.pyplot as plt
from pycalphad import Database, binplot
db_alfe = Database('alfe_sei.TDB')
my_phases_alfe = ['LIQUID', 'B2_BCC', 'FCC_A1', 'HCP_A3', 'AL5FE2', 'AL2FE', 'AL13FE4', 'AL5FE4']
models = [{'FCC_A1': MyMagneticModel}, Model]
fig = plt.figure(figsize=(9,6))
pdens = [{'B2_BCC': 20000}, 2000]
%time binplot(db_alfe, ['AL', 'FE', 'VA'] , my_phases_alfe, 'X(AL)', 300, 2000, steps=100, \
              pdens=pdens, ax=fig.gca(), model=models, mode='numexpr')
plt.title('Al-Fe with toy magnetic model', fontsize=20)
plt.show()
Explanation: Prototyping New Models (Advanced)
What if we have a new kind of functional form for the Gibbs energy we would like to experiment with?
End of explanation
<END_TASK>