In [1]:
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 CR BCC_A2 5.1996E+01 4.0500E+03 2.3560E+01!
ELEMENT V BCC_A2 5.0941E+01 4.5070E+03 3.0890E+01!
FUNCTION GHSERCR 2.98150E+02 -8856.94+157.48*T-26.908*T*LN(T)
+.00189435*T**2-1.47721E-06*T**3+139250*T**(-1); 2.18000E+03 Y
-34869.344+344.18*T-50*T*LN(T)-2.88526E+32*T**(-9); 6.00000E+03 N !
FUNCTION GPCRBCC 298.14 0.0; 6000 N !
FUNCTION GPCRLIQ 298.14 0.0; 6000 N !
FUNCTION GCRFCC 298.14 +7284+.163*T+GHSERCR#; 6000 N !
FUNCTION GHSERVV 2.98150E+02 -7930.43+133.346053*T-24.134*T*LN(T)
-.003098*T**2+1.2175E-07*T**3+69460*T**(-1); 7.90000E+02 Y
-7967.842+143.291093*T-25.9*T*LN(T)+6.25E-05*T**2-6.8E-07*T**3;
2.18300E+03 Y
-41689.864+321.140783*T-47.43*T*LN(T)+6.44389E+31*T**(-9);
4.00000E+03 N !
FUNCTION GLIQV 2.98150E+02 +12833.687+123.890501*T-24.134*T*LN(T)
-.003098*T**2+1.2175E-07*T**3+69460*T**(-1)-5.19136E-22*T**7;
7.90000E+02 Y
+12796.275+133.835541*T-25.9*T*LN(T)+6.25E-05*T**2-6.8E-07*T**3
-5.19136E-22*T**7; 2.18300E+03 Y
-19617.51+311.055983*T-47.43*T*LN(T); 4.00000E+03 N !
FUNCTION UN_ASS 298.15 0; 300 N !
FUNCTION V001 298.15 0; 300 N !
FUNCTION V002 298.15 0; 300 N !
FUNCTION V003 298.15 0; 300 N !
FUNCTION V004 298.15 0; 300 N !
TYPE_DEFINITION % SEQ *!
DEFINE_SYSTEM_DEFAULT ELEMENT 2 !
DEFAULT_COMMAND DEF_SYS_ELEMENT VA /- !
PHASE LIQUID % 1 1.0 !
CONSTITUENT LIQUID :CR,V : !
PARAMETER L(LIQUID,CR,V;0) 2.98150E+02 -20224.86+V001#;
6.00000E+03 N REF3 !
PARAMETER L(LIQUID,CR,V;1) 2.98150E+02 -13469.9+V002#;
6.00000E+03 N REF3 !
PHASE BCC_A2 %& 2 1 3 !
TYPE_DEFINITION & GES A_P_D BCC_A2 MAGNETIC -1.0 4.00000E-01 !
CONSTITUENT BCC_A2 :CR,V : VA : !
PARAMETER G(BCC_A2,CR:VA;0) 2.98150E+02 +GHSERCR#; 6.00000E+03 N
REF1 !
PARAMETER TC(BCC_A2,CR:VA;0) 2.98150E+02 -311.5;
6.00000E+03 N REF1 !
PARAMETER BMAGN(BCC_A2,CR:VA;0) 2.98150E+02 -.008; 6.00000E+03 N
REF1 !
PARAMETER L(BCC_A2,CR,V:VA;0) 298.14 -8253.85-3.61592*T+V003#;
6000 N REF3 !
PARAMETER L(BCC_A2,CR,V:VA;1) 298.14 +7494.82-8.69424*T+V004#;
6000 N REF3 !
PARAMETER L(BCC_A2,CR,V:VA;2) 298.14 -17599.07+10.13142*T;
6000 N REF3 !
PARAMETER L(BCC_A2,CR,V:VA;3) 298.14 +1042.58;
6000 N REF3 !
PARAMETER G(BCC_A2,V:VA;0) 2.98150E+02 +GHSERVV#;
4.00000E+03 N REF1 !
PHASE HCP_A3 %) 2 1 .5 !
TYPE_DEFINITION ) GES A_P_D HCP_A3 MAGNETIC -3.0 2.80000E-01 !
CONSTITUENT HCP_A3 :CR,V : VA : !
PARAMETER G(HCP_A3,CR:VA;0) 2.98150E+02 -4418.94+157.48*T
-26.908*T*LN(T)+.00189435*T**2-1.47721E-06*T**3+139250*T**(-1);
2.18000E+03 Y
-30431.344+344.18*T-50*T*LN(T)-2.88526E+32*T**(-9);
6.00000E+03 N REF1 !
PARAMETER TC(HCP_A3,CR:VA;0) 2.98150E+02 -1109;
6.00000E+03 N REF1 !
PARAMETER BMAGN(HCP_A3,CR:VA;0) 2.98150E+02 -2.46;
6.00000E+03 N REF1 !
PARAMETER G(HCP_A3,V:VA;0) 2.98150E+02 -3930.43+135.746053*T
-24.134*T*LN(T)-.003098*T**2+1.2175E-07*T**3+69460*T**(-1);
7.90000E+02 Y
-3967.842+145.691093*T-25.9*T*LN(T)+6.25E-05*T**2-6.8E-07*T**3;
2.18300E+03 Y
-37689.864+323.540783*T-47.43*T*LN(T)+6.44389E+31*T**(-9);
4.00000E+03 N REF1 !
PARAMETER L(HCP_A3,CR,V:VA;0) 2.98150E+02 +5000;
6.00000E+03 N REF3 !
"""
In [2]:
from pycalphad import Database, Model, equilibrium, calculate, variables as v
from pycalphad.codegen.callables import build_phase_records
from pycalphad.core.utils import point_sample
In [3]:
import numpy as np
rng = np.random.RandomState(1769)
dbf = Database(TDB)
comps = [v.Species('CR'), v.Species('V'), v.Species('VA')]
phases = list(dbf.phases.keys())
conds = {v.N: 1, v.T: 2500, v.P: 101325, v.X('V'): 0.25}
param_count = 2000
param_names = ['V001', 'V002', 'V003', 'V004']
models = {name: Model(dbf, comps, name, parameters=param_names) for name in phases}
parameters = {'V001': rng.normal(loc=1000, scale=1000,size=param_count), 'V002': rng.normal(scale=1000,size=param_count),
'V003': rng.normal(scale=100,size=param_count), 'V004': rng.normal(scale=100,size=param_count)}
# Specific value of parameter doesn't matter here for build_phase_records, since we'll override it
phase_records = build_phase_records(dbf, comps, phases, conds, models,
parameters={key: value[0] for key, value in parameters.items()})
# This step is not optimized. There are better ways to get a dof sample array
dof_Y = {name: calculate(dbf, comps, [name], N=1, P=1e5, T=2500).Y.values.squeeze() for name in phases}
dof_X = {name: calculate(dbf, comps, [name], N=1, P=1e5, T=2500).X.values.squeeze() for name in phases}
In [4]:
results = {}
param_array = np.column_stack(tuple(parameters[n] for n in param_names))
for name in phases:
my_dof = np.concatenate((np.broadcast_to(np.array([1, 1e5, 2500])[np.newaxis], (dof_Y[name].shape[0], 3)),
np.array(dof_Y[name])), axis=1)
result = np.zeros((dof_Y[name].shape[0], parameters['V001'].size))
phase_records[name].obj_parameters_2d(result, my_dof, param_array)
results[name] = result
In [5]:
%%timeit
from pycalphad.core.hyperplane import hyperplane
target_composition = np.atleast_1d([0.7,0.3])
target_comp_idx = np.array([1], dtype=np.uint64)
target_chempot_idx = np.array([],dtype=np.uint64)
chempot_results = {}
for name in phases:
tmp = np.zeros((param_count, 2))
d = np.ascontiguousarray(dof_X[name])
for sample_idx in range(param_count):
chemical_potentials = np.zeros(2)
result_fractions = np.zeros(3)
result_simplex = np.zeros(3, dtype=np.int32)
energies = np.ascontiguousarray(results[name][:, sample_idx])
result = hyperplane(d, energies,target_composition, chemical_potentials,1.0, target_chempot_idx, target_comp_idx,result_fractions, result_simplex)
tmp[sample_idx, :] = chemical_potentials
chempot_results[name] = tmp
In [ ]: