How to use the pycalphad.variables.P function in pycalphad

def run_test():
    dbf = Database()
    dbf.elements = frozenset(['A'])
    dbf.add_phase('TEST', {}, [1])
    dbf.add_phase_constituents('TEST', [['A']])
    # add THETA parameters here
    dbf.add_parameter('THETA', 'TEST', [['A']], 0, 334.)
    conds = {v.T: np.arange(1.,800.,1), v.P: 101325}
    res = calculate(dbf, ['A'], 'TEST', T=conds[v.T], P=conds[v.P],
                    model=EinsteinModel, output='testprop')
    #res_TE = calculate(dbf, ['A'], 'TEST', T=conds[v.T], P=conds[v.P],
    #                model=EinsteinModel, output='einstein_temperature')
    import matplotlib.pyplot as plt
    plt.scatter(res['T'], res['testprop'])
    plt.xlabel('Temperature (K)')
    plt.ylabel('Molar Heat Capacity (J/mol-K)')
    plt.savefig('einstein.png')
    print(dbf.to_string(fmt='tdb'))

def get_multiphase_constraints(self, conds):
        fixed_chempots = [cond for cond in conds.keys() if isinstance(cond, v.ChemicalPotential)]
        multiphase_constraints = []
        for statevar in sorted(conds.keys(), key=str):
            if not is_multiphase_constraint(statevar):
                continue
            if isinstance(statevar, v.Composition):
                multiphase_constraints.append(Symbol('NP') * self.moles(statevar.species))
            elif statevar == v.N:
                multiphase_constraints.append(Symbol('NP') * (sum(self.moles(spec) for spec in self.nonvacant_elements)))
            elif statevar in [v.T, v.P]:
                return multiphase_constraints.append(S.Zero)
            else:
                raise NotImplementedError
        return multiphase_constraints

continue
            build_piecewise_matrix(expr, cur_exponents, float(lowlim), float(highlim), output_matrix, symbol_matrix)
    elif isinstance(sympy_obj, Symbol):
        symbol_matrix.append([low_temp, high_temp] + list(cur_exponents) + [str(sympy_obj)])
    elif isinstance(sympy_obj, Add):
        sympy_obj = sympy_obj.expand()
        for arg in sympy_obj.args:
            build_piecewise_matrix(arg, cur_exponents, low_temp, high_temp, output_matrix, symbol_matrix)
    elif isinstance(sympy_obj, Mul):
        new_exponents = np.array(cur_exponents)
        remaining_argument = S.One
        for arg in sympy_obj.args:
            if isinstance(arg, Pow):
                if arg.base == v.T:
                    new_exponents[1] += int(arg.exp)
                elif arg.base == v.P:
                    new_exponents[0] += int(arg.exp)
                else:
                    raise NotImplementedError
            elif arg == v.T:
                new_exponents[1] += 1
            elif arg == v.P:
                new_exponents[0] += 1
            elif arg == log(v.T):
                new_exponents[3] += 1
            elif arg == log(v.P):
                new_exponents[2] += 1
            else:
                remaining_argument *= arg
        if not isinstance(remaining_argument, Mul):
            build_piecewise_matrix(remaining_argument, new_exponents, low_temp, high_temp,
                                   output_matrix, symbol_matrix)

collected_argument = S.One
            piecewise_elem = None
            for arg in sympy_obj.args:
                if isinstance(arg, Piecewise):
                    piecewise_elem = arg
                elif isinstance(arg, (Float, Integer, Rational)):
                    collected_argument *= float(arg)
                else:
                    collected_argument *= arg
            remaining_argument = Piecewise(*[(collected_argument * expr, cond) for expr, cond in piecewise_elem.args])
        else:
            for arg in sympy_obj.args:
                if isinstance(arg, Pow):
                    if arg.base == v.T:
                        new_exponents[1] += int(arg.exp)
                    elif arg.base == v.P:
                        new_exponents[0] += int(arg.exp)
                    else:
                        raise NotImplementedError
                elif arg == v.T:
                    new_exponents[1] += 1
                elif arg == v.P:
                    new_exponents[0] += 1
                elif arg == sympy_log(v.T):
                    new_exponents[3] += 1
                elif arg == sympy_log(v.P):
                    new_exponents[2] += 1
                else:
                    remaining_argument *= arg
        if not isinstance(remaining_argument, Mul):
            build_piecewise_matrix(remaining_argument, new_exponents, low_temp, high_temp,
                                   output_matrix, symbol_matrix, param_symbols)

def _map_coord_to_variable(coord):
    """
    Map a coordinate to a StateVariable object.

    Parameters
    ----------
    coord : str
        Name of coordinate in equilibrium object.

    Returns
    -------
    pycalphad StateVariable
    """
    vals = {'T': v.T, 'P': v.P}
    if coord.startswith('X_'):
        return v.X(coord[2:])
    elif coord in vals:
        return vals[coord]
    else:
        return coord

print('Phases:', end=' ')
    for name in progressbar(active_phases, desc='Initialize (1/3)', unit='phase', disable=not pbar):
        mod = models[name]
        if isinstance(mod, type):
            models[name] = mod = mod(dbf, comps, name)
        variables = sorted(mod.energy.atoms(v.StateVariable).union({key for key in conditions.keys() if key in [v.T, v.P]}), key=str)
        site_fracs = sorted(mod.energy.atoms(v.SiteFraction), key=str)
        maximum_internal_dof = max(maximum_internal_dof, len(site_fracs))
        # Extra factor '1e-100...' is to work around an annoying broadcasting bug for zero gradient entries
        #models[name].models['_broadcaster'] = 1e-100 * Mul(*variables) ** 3
        out = models[name].energy
        undefs = list(out.atoms(Symbol) - out.atoms(v.StateVariable))
        for undef in undefs:
            out = out.xreplace({undef: float(0)})
        callable_dict[name], grad_callable_dict[name], hess_callable_dict[name] = \
            build_functions(out, [v.P, v.T] + site_fracs)

        # Adjust gradient by the approximate chemical potentials
        hyperplane = Add(*[v.MU(i)*mole_fraction(dbf.phases[name], comps, i)
                           for i in comps if i != 'VA'])
        plane_obj, plane_grad, plane_hess = build_functions(hyperplane, [v.MU(i) for i in comps if i != 'VA']+site_fracs)
        mass_obj, mass_grad, mass_hess = build_functions(Add(*[mole_fraction(dbf.phases[name], comps, i)
                                                               for i in comps if i != 'VA']), site_fracs)
        phase_records[name.upper()] = PhaseRecord(variables=variables,
                                                  grad=grad_callable_dict[name],
                                                  hess=hess_callable_dict[name],
                                                  plane_grad=plane_grad,
                                                  plane_hess=plane_hess,
                                                  mass_obj=mass_obj,
                                                  mass_grad=mass_grad,
                                                  mass_hess=mass_hess)
        if verbose:

for x in symbols_to_fit:
        if isinstance(dbf.symbols[x], sympy.Piecewise):
            logging.debug('Replacing {} in database'.format(x))
            dbf.symbols[x] = dbf.symbols[x].args[0].expr

    # construct the models for each phase, substituting in the SymPy symbol to fit.
    logging.log(TRACE, 'Building phase models (this may take some time)')
    import time
    t1 = time.time()
    phases = sorted(dbf.phases.keys())
    parameters = dict(zip(symbols_to_fit, [0]*len(symbols_to_fit)))
    models = instantiate_models(dbf, comps, phases, parameters=parameters)
    if make_callables:
        eq_callables = build_callables(dbf, comps, phases, models, parameter_symbols=symbols_to_fit,
                            output='GM', build_gradients=True, build_hessians=True,
                            additional_statevars={v.N, v.P, v.T})
    else:
        eq_callables = None
    t2 = time.time()
    logging.log(TRACE, 'Finished building phase models ({:0.2f}s)'.format(t2-t1))
    logging.log(TRACE, 'Getting non-equilibrium thermochemical data (this may take some time)')
    t1 = time.time()
    thermochemical_data = get_thermochemical_data(dbf, comps, phases, datasets, weight_dict=data_weights, symbols_to_fit=symbols_to_fit)
    t2 = time.time()
    logging.log(TRACE, 'Finished getting non-equilibrium thermochemical data ({:0.2f}s)'.format(t2-t1))
    logging.log(TRACE, 'Getting equilibrium thermochemical data (this may take some time)')
    t1 = time.time()
    eq_thermochemical_data = get_equilibrium_thermochemical_data(dbf, comps, phases, datasets, parameters, data_weight_dict=data_weights)
    t2 = time.time()
    logging.log(TRACE, 'Finished getting equilibrium thermochemical data ({:0.2f}s)'.format(t2-t1))
    logging.log(TRACE, 'Getting ZPF data (this may take some time)')
    t1 = time.time()

remaining_argument = S.One
        for arg in sympy_obj.args:
            if isinstance(arg, Pow):
                if arg.base == v.T:
                    new_exponents[1] += int(arg.exp)
                elif arg.base == v.P:
                    new_exponents[0] += int(arg.exp)
                else:
                    raise NotImplementedError
            elif arg == v.T:
                new_exponents[1] += 1
            elif arg == v.P:
                new_exponents[0] += 1
            elif arg == log(v.T):
                new_exponents[3] += 1
            elif arg == log(v.P):
                new_exponents[2] += 1
            else:
                remaining_argument *= arg
        if not isinstance(remaining_argument, Mul):
            build_piecewise_matrix(remaining_argument, new_exponents, low_temp, high_temp,
                                   output_matrix, symbol_matrix)
        else:
            raise NotImplementedError(sympy_obj)
    else:
        raise NotImplementedError

import hashlib
from copy import deepcopy

# ast.Num is deprecated in Python 3.8 in favor of as ast.Constant
# Both are whitelisted for compatability across versions
_AST_WHITELIST = [ast.Add, ast.BinOp, ast.Call, ast.Constant, ast.Div,
                  ast.Expression, ast.Load, ast.Mult, ast.Name, ast.Num,
                  ast.Pow, ast.Sub, ast.UAdd, ast.UnaryOp, ast.USub]

# Avoid symbol names clashing with objects in sympy (gh-233)
clashing_namespace = {}
clashing_namespace.update(_clash)
clashing_namespace['CC'] = Symbol('CC')
clashing_namespace['FF'] = Symbol('FF')
clashing_namespace['T'] = v.T
clashing_namespace['P'] = v.P
clashing_namespace['R'] = v.R


def _sympify_string(math_string):
    "Convert math string into SymPy object."
    # drop pound symbols ('#') since they denote function names
    # we detect those automatically
    expr_string = math_string.replace('#', '')
    # sympify doesn't recognize LN as ln()
    expr_string = \
        re.sub(r'(?

tuple
        Tuple of (Gibbs energies, phases, phase fractions, compositions, site fractions, chemical potentials)

    Notes
    -----
    Assumes that potentials are fixed and there is just a 1d composition grid.
    Minimizes the use of Dataset objects.

    """
    calc_opts = calc_opts or {}
    conditions = _adjust_conditions(conditions)

    # 'calculate' accepts conditions through its keyword arguments
    if 'pdens' not in calc_opts:
        calc_opts['pdens'] = 500
    grid = calculate(dbf, comps, phases, T=conditions[v.T], P=conditions[v.P],
                     parameters=parameters, fake_points=True, output='GM',
                     callables=callables, model=model, **calc_opts)

    # assume only one independent component
    indep_comp_conds = [c for c in conditions if isinstance(c, v.X)]
    num_indep_comp = len(indep_comp_conds)
    if num_indep_comp != 1:
        raise ConditionError(
            "Convex hull independent components different than one.")
    max_num_phases = (len(indep_comp_conds) + 1,)  # Gibbs phase rule
    comp_grid = build_composition_grid(comps, conditions)
    calc_grid_shape = comp_grid.shape[:-1]
    num_comps = comp_grid.shape[-1:]

    grid_energy_values = grid.GM.values.squeeze()
    grid_composition_values = grid.X.values.squeeze()