Update prismspf coupling

ammber/BinarySystems.py

@@ -3,20 +3,23 @@
 from scipy.interpolate import CubicSpline
 import numpy as np
 from pycalphad import calculate
-from ammber.utils import write_binary_isothermal_parabolic_parameters
+from ammber.utils import add_to_dict
 class BinaryIsothermalDiscretePhase:
     "Class representing a single phase at one temperature described by a set of points in composition-Gibbs free energy space"
-    def __init__(self, xdata, Gdata):
+    def __init__(self, name, xdata, Gdata):
         """
         Constructor. Initializes the phase object from points in composition-Gibbs free energy space.
 
         Parameters
         ----------
+        name : string
+            name of phase to be added
         xdata : list, array
             compositions
         Gdata : list, array
             Gibbs free energies corresponding to compositions
         """
+        self.name = name
         self.xdata, ordering = np.unique(xdata, return_index=True)
         self.Gdata = Gdata[ordering]
 
@@ -36,7 +39,7 @@ def resample(self, xpoints):
         newx = xpoints[np.logical_and(xpoints >= self.xdata[0], xpoints <= self.xdata[-1])]
         if len(newx)>=3:
             spl = CubicSpline(self.xdata, self.Gdata)
-            return BinaryIsothermalDiscretePhase(newx, spl(newx))
+            return BinaryIsothermalDiscretePhase(self.name, newx, spl(newx))
         # Assume a line compound, don't resample
         return self
 
@@ -82,14 +85,17 @@ def free_energy(self, x):
             print("Attempted to evaluate free energy of a compound outside its range. This will cause errors.")
             return None
 
+
 class BinaryIsothermal2ndOrderPhase:
     "Class representing a single phase at one temperature described by a 2nd order polynomial (parabola)."
-    def __init__(self, fmin=0.0, kwell=1.0, cmin=0.5, discrete=None, kwellmax=1e9):#todo
+    def __init__(self, name, fmin=0.0, kwell=1.0, cmin=0.5, discrete=None, kwellmax=1e9):#todo
         """
         Constructor.
 
         Parameters
         ----------
+        name : string
+            name of phase to be added
         cmin : float
             composition of the parabolic free energy minimum
         fmin : float
@@ -104,10 +110,12 @@ def __init__(self, fmin=0.0, kwell=1.0, cmin=0.5, discrete=None, kwellmax=1e9):#
         BinaryIsothermal2ndOrderPhase object
         """
         if discrete is None:
+            self.name = name
             self.fmin = fmin
             self.kwell = kwell
             self.cmin = cmin
         else:
+            self.name = discrete.name
             self.fit_phase(discrete.xdata, discrete.Gdata, kwellmax=kwellmax)
 
     def fit_phase(self, xdata, Gdata, kwellmax=1e9):
@@ -177,40 +185,52 @@ def discretize(self, xdata=None, xrange=(1e-14, 1.0 - 1e-14), npts=1001):
         """
         if xdata is None:
             xdata = np.linspace(*xrange, npts)
-        return BinaryIsothermalDiscretePhase(xdata, self.free_energy(xdata))
+        return BinaryIsothermalDiscretePhase(self.name, xdata, self.free_energy(xdata))
 
 
 class BinaryIsothermalDiscreteSystem:
     "Class representing a seet of phases at one temperature described by a set of points in composition-Gibbs free energy space"
-    def __init__(self):
-        "Constructor. Initializes with an empty set of phases."
+    def __init__(self, component="", solution_component=""):
+        """
+        Constructor. Initializes with an empty set of phases.
+
+        Parameters
+        ----------
+        component : string
+            name of the x-component
+        solution_component : string
+            name of the (1-x)-component
+        """
+        self.component = component
+        self.solution_component = solution_component
         self.phases = {}
         self.lc_hull = None
         self._is_equilibrium_system = False
 
-    def fromTDB(self, db, elements, component, temperature, phase_list=None):
+    def fromTDB(self, db, component, solution_component, temperature, phase_list=None):
         """
         Constructs discrete phase from all the phases specified in a pycalphad database
 
         Parameters
         ----------
         db : pycalphad database
-        elements : list (string)
-            The space of compositions of interest (Binary). Element abbreviations must be all-caps.
-        componenet : string
-            Element abbreviation for element corresponding to x=1.
+        component : string
+            Element abbreviation for element corresponding to x. Element abbreviations must be all-caps.
+        solution_component : string
+            Element abbreviation for element corresponding to (1-x). Element abbreviations must be all-caps.
         phase_list : list (string)
             If specified, only listed phases will be constructed, otherwise, all available phases will be constructed.
         """
-        if "VA" not in elements:
-            elements.append("VA")
+        elements = [component, solution_component, "VA"]
+        self.component = component
+        self.solution_component = solution_component
         if phase_list is None:
             phase_list = list(db.phases.keys())
             cal = calculate(db, elements, phase_list, P=101325, T=temperature, output='GM', pdens=1001)
             phase_list = set(cal.Phase.data[0][0][0])
         for phase_name in phase_list:
             result = calculate(db, elements, phase_name, P=101325, T=temperature, output='GM', pdens=1001)
-            self.phases[phase_name] = BinaryIsothermalDiscretePhase(
+            self.phases[phase_name] = BinaryIsothermalDiscretePhase(phase_name,
                 result.X.sel(component=component).data[0][0][0],
                 result.GM.data[0][0][0])
 
@@ -227,7 +247,7 @@ def add_phase(self, name, xdata, Gdata):
         Gdata : list, array
             Gibbs free energies corresponding to compositions
         """
-        self.phases[name] = BinaryIsothermalDiscretePhase(xdata, Gdata)
+        self.phases[name] = BinaryIsothermalDiscretePhase(name, xdata, Gdata)
 
     def get_lc_hull(self, recalculate=False):#todo
         """
@@ -268,13 +288,14 @@ def get_equilibrium_system(self, miscibility_gap_threshold=0.1):#todo
         """
         lc_hull = self.get_lc_hull()
   #      print(lc_hull)
-        equilibrium_system = BinaryIsothermalDiscreteSystem()
+        equilibrium_system = BinaryIsothermalDiscreteSystem(component=self.component,
+            solution_component=self.solution_component)
         equilibrium_system._is_equilibrium_system = True
         for phase in self.phases:
             mask = np.logical_and(np.isin(self.phases[phase].xdata, lc_hull[:, 0]),
                                   np.isin(self.phases[phase].Gdata, lc_hull[:, 1]))
             if np.count_nonzero(mask) > 0:
-                equilibrium_system.phases[phase] = BinaryIsothermalDiscretePhase(
+                equilibrium_system.phases[phase] = BinaryIsothermalDiscretePhase(phase,
                     self.phases[phase].xdata[mask],
                     self.phases[phase].Gdata[mask])
         # break apart phases that phase separate
@@ -299,12 +320,12 @@ def get_equilibrium_system(self, miscibility_gap_threshold=0.1):#todo
                     #print(start,end)
                     if end - start >= 1:
                         # phase is represented by multiple points
-                        equilibrium_system.phases[phase +'_'+ str(i)] = BinaryIsothermalDiscretePhase(
+                        equilibrium_system.phases[phase +'_'+ str(i)] = BinaryIsothermalDiscretePhase(phase +'_'+ str(i),
                             phase_temp.xdata[start:end], phase_temp.Gdata[start:end])
                     elif(start == end):
                         # phase is represented by single point
                         print("phase "+phase+'_'+str(i)+" is represented by a single point, a finer discretization may be needed")
-                        equilibrium_system.phases[phase +'_'+ str(i)] = BinaryIsothermalDiscretePhase(
+                        equilibrium_system.phases[phase +'_'+ str(i)] = BinaryIsothermalDiscretePhase(phase +'_'+ str(i),
                             phase_temp.xdata[start], phase_temp.Gdata[start])
         return equilibrium_system
 
@@ -363,7 +384,8 @@ def resample_near_equilibrium(self, x, xdist=0.001, npts=101):
         BinaryIsothermalDiscreteSystem object
         """
         equilibrium = self.get_equilibrium(x)
-        new_system = BinaryIsothermalDiscreteSystem()
+        new_system = BinaryIsothermalDiscreteSystem(component=self.component,
+            solution_component=self.solution_component)
         for phase in equilibrium:
             if phase in self.phases:
                 new_system.phases[phase] = self.phases[phase].resample_near_xpoint(
@@ -378,12 +400,16 @@ def resample_near_equilibrium(self, x, xdist=0.001, npts=101):
 
 class BinaryIsothermal2ndOrderSystem:
     "Class representing a seet of phases at one temperature described by a second order polynomial (parabola)"
-    def __init__(self, phases=None):
+    def __init__(self, component="", solution_component="", phases=None):
         """
         Constructor.
 
         Parameters
         ----------
+        component : string
+            name of the x-component
+        solution_component : string
+            name of the (1-x)-component
         phases : dict {string phase_name : BinaryIsothermal2ndOrderPhase phase}
             (optional) composition to be sampled
         """
@@ -402,8 +428,10 @@ def from_discrete(self, discrete_system, kwellmax=1e6):
         kwellmax : float
             (optional) kwell to be used when fitting line compounds
         """
+        self.component = discrete_system.component
+        self.solution_component = discrete_system.solution_component
         for phase_name in discrete_system.phases.keys():
-            self.phases[phase_name] = BinaryIsothermal2ndOrderPhase()
+            self.phases[phase_name] = BinaryIsothermal2ndOrderPhase(phase_name)
             self.phases[phase_name].fit_phase(discrete_system.phases[phase_name].xdata,
                 discrete_system.phases[phase_name].Gdata, kwellmax=kwellmax)
 
@@ -441,7 +469,8 @@ def to_discrete(self, xdata=None, xrange=(1e-14, 1.0 - 1e-14), npts=1001):
         -------
         BinaryIsothermalDiscretePhase object
         """
-        discrete_system = BinaryIsothermalDiscreteSystem()
+        discrete_system = BinaryIsothermalDiscreteSystem(component=self.component,
+            solution_component=self.solution_component)
         for phase_name in self.phases.keys():
             discrete_system.phases[phase_name] = self.phases[phase_name].discretize(
                 xdata=xdata, xrange=xrange, npts=npts)

ammber/utils.py

@@ -1,78 +1,50 @@
-def write_binary_isothermal_parabolic_parameters(binaryIsothermalSys, output_file, template_file, phases=None, component="comp_1"):
+def add_to_dict(binaryIsothermalSys, dict, add_templates=False, c0={}, Vm=None, convert_to_volumetric_energy=True):
     """
-    Creates a parameters file from a BinaryIsothermal2ndOrderSystem
-
+    Adds the parameters of a BinaryIsothermal2ndOrderSystem to a dictionary
     Parameters
     ----------
     binaryIsothermalSys : BinaryIsothermal2ndOrderSystem
         system to draw parameters from
-    output_file : string
-        path to file to create
-    template_file : string
-        path to template parameter file if needed
-    phases : list (string)
-        (optional) list of specific phases to copy to parameters
-    component : string
-        (optional) name of the x-component
+    dict : dict
+        dictionary to add parameters to
+    add_templates : bool
+        add templates or additional values needed for kinetics simulation
+    c0 : dict
+        phase-name keys, initial composition values. only used if add_templates==True
+    Vm : float
+        nominal number-volume (eg. molar volume) of system
+    convert_to_volumetric_energy : bool
+        a setting for AMMBER kinetics. true if binaryIsothermalSys composition axes are number fraction
     """
-    if phases is None:
-        phase_list = binaryIsothermalSys.phases.keys()
-    else:
-        phase_list = phases
-    cmin = {}
-    fmin = {}
-    kwell= {}
-    for phase_name in phase_list:
-        cmin[phase_name] = {component : binaryIsothermalSys.phases[phase_name].cmin}
-        fmin[phase_name] = binaryIsothermalSys.phases[phase_name].fmin
-        kwell[phase_name] = {component : binaryIsothermalSys.phases[phase_name].kwell}
-    _PRISMS_parabolic_parameters(cmin=cmin, fmin=fmin, kwell=kwell, 
-                            phases = phase_list ,comps=[component],
-                            template_file=template_file, output_file=output_file)
-
-
-def _PRISMS_parabolic_parameters(cmin, fmin, kwell, phases, comps, output_file, template_file):
-    ## Step 2: Make necessary strings
-    num_phases = len(phases)
-    num_comps = len(comps)+1
-    sf = 3
-    fmin_str = ""
-    cmin_str = ""
-    kwell_str = ""
-    phase_names_str = ""
-    comp_names_str = ""
-    for phase in phases:
-        fmin_str += f'{fmin[phase]:.{sf}e}' + ", "
-        phase_names_str += phase + ", " 
-        for comp in comps:
-            cmin_str += f'{cmin[phase][comp]:.{sf}f}' + ","
-            kwell_str += f'{kwell[phase][comp]:.{sf}e}' + ","
-        cmin_str +=  " "
-        kwell_str += " "
-
-    for comp in comps:
-        comp_names_str += comp + ", "
-
-    ## Step 3: read template and write prm
-    import re
-    with open(template_file, 'r') as f_in, open(output_file, 'w') as f_out:
-            for line in f_in:
-                if re.match(r'^set Model constant num_phases\b', line):
-                    f_out.write("set Model constant num_phases = " + str(num_phases) + ", INT\n")
-                elif re.match(r'^set Model constant num_comps\b', line):
-                    f_out.write("set Model constant num_comps = " + str(num_comps) + ", INT\n")
-                elif re.match(r'^set Model constant fWell\b', line):
-                    f_out.write("set Model constant fWell = " + fmin_str + "DOUBLE ARRAY\n")
-                elif re.match(r'^set Model constant kWell\b', line):
-                    f_out.write("set Model constant kWell = " + kwell_str + "DOUBLE ARRAY\n")
-                elif re.match(r'^set Model constant cmin\b', line):
-                    f_out.write("set Model constant cmin = " + cmin_str + "DOUBLE ARRAY\n")
-                elif re.match(r'^set Model constant phase_names\b', line):
-                    f_out.write("set Model constant phase_names = " + phase_names_str + "STRING ARRAY\n")
-                elif re.match(r'^set Model constant comp_names\b', line):
-                    f_out.write("set Model constant comp_names = " + comp_names_str + "STRING ARRAY\n")
-                else:
-                    f_out.write(line)
-
-    print("Model Parameters written to " + output_file + "\n")
+    dict["solution_component"] = binaryIsothermalSys.solution_component
+    dict["components"] = [binaryIsothermalSys.component]
+    dict["convert_fractional_to_volumetric_energy"] = convert_to_volumetric_energy
+    if "phases" not in dict:
+        dict["phases"] = {}
+    comp = binaryIsothermalSys.component
+    for phase_name in binaryIsothermalSys.phases.keys():
+        phase = binaryIsothermalSys.phases[phase_name]
+        if phase_name not in dict:
+            dict["phases"][phase_name] = {}
+
+        if comp not in dict["phases"][phase_name]:
+            dict["phases"][phase_name][comp] = {}
+        dict["phases"][phase_name][comp]["k_well"] = phase.kwell
+        dict["phases"][phase_name][comp]["c_min"] = phase.cmin
+        dict["phases"][phase_name][comp]["f_min"] = phase.fmin
 
+        if add_templates:
+            c0_phase_keys = list(c0.keys())
+            if "c0" not in dict["phases"][phase_name][comp]:
+                dict["phases"][phase_name][comp]["c0"] = c0[phase_name] if phase_name in c0_phase_keys else -1.0
+            for phase_prop in ["mu_int", "D", "sigma"]:
+                if phase_prop not in dict["phases"][phase_name]:
+                    dict["phases"][phase_name][phase_prop] = -1.0
+    if add_templates:
+        if "l_int" not in dict:
+            dict["l_int"] = -1.0
+        if "Vm" not in dict:
+            dict["Vm"] = Vm if Vm is not None else -1.0
+        if "order_parameters" not in dict:
+            dict["order_parameters"] = list(binaryIsothermalSys.phases.keys())
+