In [1]:

    

from pycalphad import Database, Model, variables as v
from pycalphad import equilibrium, calculate
import numpy as np
from pycalphad.tests.datasets import *
dbf = Database('alzn_mey.tdb')


    

In [2]:

    

from pycalphad.core.solver import InteriorPointSolver
class ProblemSaver(InteriorPointSolver):
    saved_problem = [None]
    def solve(self, prob):
        self.saved_problem[0] = prob
        self.verbose = True
        return super(ProblemSaver, self).solve(prob)
v.T.default_value = 930
eq2 = equilibrium(dbf, ['AL', 'ZN', 'VA'], ['FCC_A1', 'LIQUID', 'HCP_A3'],
                 {v.X('ZN'): 1e-4, v.NP('LIQUID'): 1e-4,
                  v.P: 1e5}, verbose=False, solver=ProblemSaver())
print(eq2)


    

    




('Composition Sets', [CompositionSet(FCC_A1, [0.11022044 0.88977956], NP=0.3333333333333333, GM=-49634.65671841146), CompositionSet(HCP_A3, [0.03006012 0.96993988], NP=0.3333333333333333, GM=-50490.08145164182), CompositionSet(LIQUID, [0.05811623 0.94188377], NP=0.0001 [fixed], GM=-53231.39019434171)])
('Chemical Potentials', array([0., 0.]))
Trying to improve poor solution
Calculation Failed:  OrderedDict([('N', array(1.)), ('NP_LIQUID', array(0.0001)), ('P', array(100000.)), ('X_ZN', array(0.0001))]) b"Restoration phase failed, algorithm doesn't know how to proceed."
Chemical Potentials [-5625.11667848 -9260.85608864]
[1.00000000e+00 1.00000000e+05 1.00000000e-06 1.79170456e-10
 1.46201273e-10 1.16415322e-10 1.00000001e+06 2.41173741e+03
 4.43595388e+03 5.49002284e-04]
[1.00000000e+00 1.00000000e+05 3.00000000e+02 9.97257972e-01
 2.74202829e-03 9.98058935e-01 1.94106472e-03 2.62233591e-03
 9.97377664e-01 5.86390262e-01 4.13608748e-01 1.00000000e-06]
Status: -2 b"Restoration phase failed, algorithm doesn't know how to proceed."
<xarray.Dataset>
Dimensions:    (N: 1, NP_LIQUID: 1, P: 1, X_ZN: 1, component: 2, internal_dof: 2, vertex: 3)
Coordinates:
  * N          (N) float64 1.0
  * NP_LIQUID  (NP_LIQUID) float64 0.0001
  * P          (P) float64 1e+05
  * X_ZN       (X_ZN) float64 0.0001
  * vertex     (vertex) int64 0 1 2
  * component  (component) <U2 'AL' 'ZN'
Dimensions without coordinates: internal_dof
Data variables:
    NP         (N, P, NP_LIQUID, X_ZN, vertex) float64 nan nan nan
    GM         (N, P, NP_LIQUID, X_ZN) float64 nan
    MU         (N, P, NP_LIQUID, X_ZN, component) float64 nan nan
    X          (N, P, NP_LIQUID, X_ZN, vertex, component) float64 nan ... nan
    Y          (N, P, NP_LIQUID, X_ZN, vertex, internal_dof) float64 nan ... nan
    Phase      (N, P, NP_LIQUID, X_ZN, vertex) <U6 '' '' ''
    T          (N, P, NP_LIQUID, X_ZN) float64 nan
Attributes:
    engine:   pycalphad 0.7.1.post2+33.gf8d132d
    created:  2019-05-22T16:10:24.949663

In [3]:

    

eq2.GM - ((eq2.NP.isel(vertex=0) * eq2.MU * eq2.X.isel(vertex=0)).sum() + 
          (eq2.NP.isel(vertex=1) * (eq2.MU * eq2.X.isel(vertex=1)).sum()) + 
          (eq2.NP.isel(vertex=2) * (eq2.MU * eq2.X.isel(vertex=2)).sum()))


    

    
Out[3]:





<xarray.DataArray (N: 1, P: 1, NP_LIQUID: 1, X_ZN: 1)>
array([[[[nan]]]])
Coordinates:
  * N          (N) float64 1.0
  * NP_LIQUID  (NP_LIQUID) float64 0.0001
  * P          (P) float64 1e+05
  * X_ZN       (X_ZN) float64 0.0001
    vertex     int64 2

In [4]:

    

eq2.X


    

    
Out[4]:





<xarray.DataArray 'X' (N: 1, P: 1, NP_LIQUID: 1, X_ZN: 1, vertex: 3, component: 2)>
array([[[[[[nan, nan],
           [nan, nan],
           [nan, nan]]]]]])
Coordinates:
  * N          (N) float64 1.0
  * NP_LIQUID  (NP_LIQUID) float64 0.0001
  * P          (P) float64 1e+05
  * X_ZN       (X_ZN) float64 0.0001
  * vertex     (vertex) int64 0 1 2
  * component  (component) <U2 'AL' 'ZN'

In [5]:

    

%matplotlib inline
import matplotlib.pyplot as plt
print(eq2['X_ZN'].values.flatten())
print(eq2['T'].values.flatten())
plt.scatter(eq2['X_ZN'].values.flatten(), eq2['T'].values.flatten())
plt.xlim((0,1))
plt.xlabel('X(ZN)')
plt.ylabel('T')


    

    




[0.0001]
[nan]






    
Out[5]:





Text(0, 0.5, 'T')






    

In [6]:

    

fcc = calculate(dbf, ['AL', 'ZN', 'VA'], 'FCC_A1',
                 T=float(eq2['T']), P=1e5)
hcp = calculate(dbf, ['AL', 'ZN', 'VA'], 'HCP_A3',
                 T=float(eq2['T']), P=1e5)
liq = calculate(dbf, ['AL', 'ZN', 'VA'], 'LIQUID',
                 T=float(eq2['T']), P=1e5)


    

In [7]:

    

%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
plt.figure(figsize=(12,12))
plt.scatter(fcc.X.sel(component='ZN').values.flatten(),
            fcc.GM.values.flatten(), c='g')
plt.scatter(hcp.X.sel(component='ZN').values.flatten(),
            hcp.GM.values.flatten(), c='y')
plt.scatter(liq.X.sel(component='ZN').values.flatten(),
            liq.GM.values.flatten(), c='b')
plt.plot([0,1], eq2.MU.values.flatten())
plt.scatter([1e-4], [float(eq2.GM)])
plt.grid(True)


    

In [8]:

    

eq2.Y


    

    
Out[8]:





<xarray.DataArray 'Y' (N: 1, P: 1, NP_LIQUID: 1, X_ZN: 1, vertex: 3, internal_dof: 2)>
array([[[[[[nan, nan],
           [nan, nan],
           [nan, nan]]]]]])
Coordinates:
  * N          (N) float64 1.0
  * NP_LIQUID  (NP_LIQUID) float64 0.0001
  * P          (P) float64 1e+05
  * X_ZN       (X_ZN) float64 0.0001
  * vertex     (vertex) int64 0 1 2
Dimensions without coordinates: internal_dof

In [9]:

    

eq2.MU


    

    
Out[9]:





<xarray.DataArray 'MU' (N: 1, P: 1, NP_LIQUID: 1, X_ZN: 1, component: 2)>
array([[[[[nan, nan]]]]])
Coordinates:
  * N          (N) float64 1.0
  * NP_LIQUID  (NP_LIQUID) float64 0.0001
  * P          (P) float64 1e+05
  * X_ZN       (X_ZN) float64 0.0001
  * component  (component) <U2 'AL' 'ZN'

In [10]:

    

eq2.GM


    

    
Out[10]:





<xarray.DataArray 'GM' (N: 1, P: 1, NP_LIQUID: 1, X_ZN: 1)>
array([[[[nan]]]])
Coordinates:
  * N          (N) float64 1.0
  * NP_LIQUID  (NP_LIQUID) float64 0.0001
  * P          (P) float64 1e+05
  * X_ZN       (X_ZN) float64 0.0001

In [11]:

    

import numpy as np
def constraint_calculator(temps):
    # Should be zero at T=933
    x_in = [1.00000000e+00,   1.00000000e+05,   3.00000000e+02,   9.97257972e-01,
   2.74202827e-03,   9.98058935e-01,   1.94106495e-03,   2.62233490e-03,
   9.97377665e-01,   5.86390285e-01,   4.13608725e-01,   1.00000000e-06]
    desired_chempots = [-37840.74636944, -110371.97774282]
    cons = []
    total_chempots = []
    fixed_chempots = []
    free_chempots_1 = []
    free_chempots_2 = []
    for temp in temps:
        x_in[2] = temp
        cons.append(ProblemSaver.saved_problem[0].constraints(x_in))
        total_chempots.append(ProblemSaver.saved_problem[0].chemical_potentials(x_in))
        fixed_chempots.append(ProblemSaver.saved_problem[0].chemical_potentials(x_in, selected_phase=1))
        free_chempots_1.append(ProblemSaver.saved_problem[0].chemical_potentials(x_in, selected_phase=0))
        free_chempots_2.append(ProblemSaver.saved_problem[0].chemical_potentials(x_in, free_only=True))
    cons = np.array(cons - np.array(ProblemSaver.saved_problem[0].cl))
    total_chempots = np.array(total_chempots) - np.array(desired_chempots)
    fixed_chempots = np.array(fixed_chempots) - np.array(desired_chempots)
    free_chempots_1 = np.array(free_chempots_1) - np.array(desired_chempots)
    free_chempots_2 = np.array(free_chempots_2) - np.array(desired_chempots)
    return cons, total_chempots, fixed_chempots, free_chempots_1, free_chempots_2

def constraint_calculator_2(xfccs):
    # Should be zero at T=933
    x_in = [1.00000000e+00,   1.00000000e+05,   3.00000000e+02,   9.97257972e-01,
   2.74202827e-03,   9.98058935e-01,   1.94106495e-03,   2.62233490e-03,
   9.97377665e-01,   5.86390285e-01,   4.13608725e-01,   1.00000000e-06]
    desired_chempots = [-37840.74636944, -110371.97774282]
    cons = []
    total_chempots = []
    fixed_chempots = []
    free_chempots_1 = []
    free_chempots_2 = []
    for xfcc in xfccs:
        x_in[2] = 933
        x_in[3] = xfcc
        x_in[4] - 1-xfcc
        cons.append(ProblemSaver.saved_problem[0].constraints(x_in))
        total_chempots.append(ProblemSaver.saved_problem[0].chemical_potentials(x_in))
        fixed_chempots.append(ProblemSaver.saved_problem[0].chemical_potentials(x_in, selected_phase=1))
        free_chempots_1.append(ProblemSaver.saved_problem[0].chemical_potentials(x_in, selected_phase=0))
        free_chempots_2.append(ProblemSaver.saved_problem[0].chemical_potentials(x_in, free_only=True))
    cons = np.array(cons - np.array(ProblemSaver.saved_problem[0].cl))
    total_chempots = np.array(total_chempots) - np.array(desired_chempots)
    fixed_chempots = np.array(fixed_chempots) - np.array(desired_chempots)
    free_chempots_1 = np.array(free_chempots_1) - np.array(desired_chempots)
    free_chempots_2 = np.array(free_chempots_2) - np.array(desired_chempots)
    return cons, total_chempots, fixed_chempots, free_chempots_1, free_chempots_2

temps = np.arange(300., 2000., step=100)
xfccs = np.linspace(1e-6,1-1e-6, num=100)
results = constraint_calculator(temps)
results_2 = constraint_calculator_2(xfccs)

%matplotlib inline
import matplotlib.pyplot as plt
plt.plot(temps, results[0][:, -2], label='cons 1')
plt.plot(temps, results[0][:, -1], label='cons 2')
plt.plot(temps, results[1][:, 0], label='total')
plt.plot(temps, results[2][:, 0], label='fixed')
plt.plot(temps, results[3][:, 0], label='free 1')
plt.plot(temps, results[4][:, 0], label='free 2')
#plt.ylim((-1000, 1000))
#plt.xlim((700, 900))
#plt.plot(xfccs, results_2[0][:, -2], label='cons 1')
#plt.plot(xfccs, results_2[0][:, -1], label='cons 2')
#plt.xlim((0,0.1))
plt.legend(loc='best')


    

    
Out[11]:





<matplotlib.legend.Legend at 0x7f6a6a401828>






    

In [12]:

    

np.argmin(np.abs(results_2[0][:, -1]))


    

    
Out[12]:





0

In [13]:

    

xfccs[1]


    

    
Out[13]:





0.010101989898989897

In [14]:

    

results[0][:, -2]


    

    
Out[14]:





array([ 4435.95365963,  8862.03732084, 13286.41980792, 17705.7904437 ,
       22112.93527331, 26503.39500828, 30888.03652269, 35271.18282894,
       39653.88934555, 44036.44486113, 48418.94323598, 52801.41818906,
       57183.88286977, 61566.34277611, 76923.97786827, 82272.80919738,
       87621.64057337])

In [ ]: