In [1]:
%matplotlib inline
import matplotlib.pyplot as plt
import pycalphad.io.tdb_keywords
pycalphad.io.tdb_keywords.TDB_PARAM_TYPES.append('TH')
from pycalphad import Database, binplot, Model
import pycalphad.variables as v
from dask.distributed import Client

In [5]:
TDB = """
 ELEMENT /-   ELECTRON_GAS              0.0000E+00  0.0000E+00  0.0000E+00!
 ELEMENT VA   VACUUM                    0.0000E+00  0.0000E+00  0.0000E+00!
 ELEMENT AL   FCC_A1                    2.6982E+01  4.5773E+03  2.8322E+01!
 ELEMENT SB   RHOMBOHEDRAL_A7           1.2175E+02  5.8702E+03  4.5522E+01!


TYPE_DEFINITION % SEQ * !
DEFINE_SYSTEM_DEFAULT SPECIE 2 !
DEFAULT_COMMAND DEF_SYS_ELEMENT VA !

FUN GHSERSB   298.15   -9242.858+156.154689*T-30.5130752*T*LN(T)
                          +.007748768*T**2-3.003415E-06*T**3
                          +100625*T**(-1);
              903.9 Y -11738.671+169.485713*T-31.38*T*LN(T)
                          +1.610442E+27*T**(-9);
             2000.00 N 91DIN !

 FUNCTION GHSERAL    2.98150E+02  -7976.15+137.093038*T-24.3671976*T*LN(T)
     -.001884662*T**2-8.77664E-07*T**3+74092*T**(-1);  7.00000E+02  Y
      -11276.24+223.048446*T-38.5844296*T*LN(T)+.018531982*T**2
     -5.764227E-06*T**3+74092*T**(-1);  9.33470E+02  Y
      -11278.378+188.684153*T-31.748192*T*LN(T)-1.230524E+28*T**(-9);  
     2.90000E+03  N !


SPECIES AL1SB1 AL1SB1 !

PHASE LIQUID:L %  1  1.0 !
   CONSTITUENT LIQUID:L  : AL,SB,AL1SB1 : !

PHASE RHOMBOHEDRAL_A7 %  1  1.0 !
   CONSTITUENT RHOMBOHEDRAL_A7  : AL,SB : !

PHASE ALSB % 2 1.0 1.0 !
   CONSTITUENT ALSB : AL : SB : !


PARA G(RHOMBOHEDRAL_A7,SB;0)   298.15
 -9242.858+156.154689*T-30.5130752*T*LN(T)+7.748768E-3*T**2-3.003415E-6*T**3
 +100625*T**(-1);   HICKEL Y
 -11738.83+169.485872*T-31.38*T*LN(T)+1616.849E24*T**(-9);  2000.00 N !

PARA G(RHOMBOHEDRAL_A7,AL;0) 298.15 20000+GHSERAL#; 2000 N !

PARA TH(RHOMBOHEDRAL_A7,AL;0) 1 11111; 10000 N !
PARA TH(RHOMBOHEDRAL_A7,SB;0) 1 222222; 10000 N !

PARA G(LIQUID,AL;0) 2.98140E+02  +11005-11.8419*T+7.934E-20*T**7+GHSERAL#;  
                 933  Y
      +10482.3-11.254*T+1.231E+28*T**(-9)+GHSERAL#;  6.00000E+03  N !

PARA G(LIQUID,SB;0)  298.15   +19822.329-21.923164*T-1.7485E-20*T**7
                                 +GHSERSB#;
                     903.78 Y +19914.189-22.029886*T-1.6168E+27*T**(-9)
                                 +GHSERSB#;
                    2000.00 N 91Din !
PARA G(LIQUID,AL1SB1;0) 298.15 -24797 + 6.6259*T + 0.82937*T*(1-LN(T)) + 2.6817E-3*T**2
                                 +GHSERSB# + GHSERAL#; 2000 N 90Cou !

PARA G(ALSB,AL:SB;0) 298.15 -24486 + 6.7006*T + 0.82783*T*(1-LN(T)) + 2.6797E-3*T**2
                                 +GHSERSB# + GHSERAL#; 2000 N 90Cou !

"""

In [7]:
from sympy import Symbol

dbf = Database(TDB)

class CustomModel(Model):
    def __init__(self, dbe, comps, phase_name, parameters=None):
        super(self, CustomModel).__init__(dbe, comps, phase_name, parameters=parameters)
        phase = dbe.phases[self.phase_name]
        param_search = dbe.search
        th_param_query = (
            (where('phase_name') == phase.name) & \
            (where('parameter_type') == 'TH') & \
            (where('constituent_array').test(self._array_validity))
        )

        th = self.redlich_kister_sum(phase, param_search, th_param_query)
        #for name, value in self.models.items():
        #    self.models[name] = self.models[name].xreplace({Symbol('HICKEL'): th})

mod = Model(dbf, ['AL'], 'RHOMBOHEDRAL_A7')
mod.GM


Out[7]:
1.0*Piecewise((Piecewise((-8.77664e-7*T**3 - 0.001884662*T**2 - 24.3671976*T*log(T) + 137.093038*T - 7976.15 + 74092/T, (T >= 298.15) & (T < 700.0)), (-5.764227e-6*T**3 + 0.018531982*T**2 - 38.5844296*T*log(T) + 223.048446*T - 11276.24 + 74092/T, (T >= 700.0) & (T < 933.47)), (-31.748192*T*log(T) + 188.684153*T - 11278.378 - 1.230524e+28/T**9, (T >= 933.47) & (T < 2900.0)), (0, True)) + 20000, (T >= 298.15) & (T < 2000.0)), (0, True)) + 8.3145*T*Piecewise((RHOMBOHEDRAL_A70AL*log(RHOMBOHEDRAL_A70AL), RHOMBOHEDRAL_A70AL > 1.0e-13), (0, True))/RHOMBOHEDRAL_A70AL + 1.0*(24.9435*T*log(1 - exp(-exp(RHOMBOHEDRAL_A70AL*Piecewise((11111, (T >= 1.0) & (T < 10000.0)), (0, True)))/T)) + 12.47175*exp(RHOMBOHEDRAL_A70AL*Piecewise((11111, (T >= 1.0) & (T < 10000.0)), (0, True))))/RHOMBOHEDRAL_A70AL

In [ ]:
fig = plt.figure(figsize=(9,6))
scheduler = Client()

%time binplot(dbf, ['AL', 'SB', 'AL1SB1'] , sorted(dbf.phases.keys()), {v.X('SB'):(0,1,0.02), v.T: (300, 1300, 50), v.P:101325},  ax=fig.gca(), eq_kwargs={'scheduler': 'distributed'})

plt.show()

In [ ]:
dbf.phases

In [ ]: