Molten Salt Research with various ionic cholrides with alkali earth metals

molten_salt_research.py

from __future__ import division
from thermo.chemical import Chemical, Mixture
from matplotlib import pyplot as plt
import thermo
import math
import decimal
import sys

'''
LR
	- needs polynomial fitted dielectric constants to a range of temperature and dielectric constants of pure subgroups (Dortmund Data Bank)
	- needs charge of ions
	- molar mas of the solvent
	- density of solvent 
	- volume fractions 
	- molarity of ions
	- relative dieletric constat of mixture
'''

'''
MR
	[*]- ionic strength (needs charge of ions, molarity of ions) (mol/kg)
		- thermo.electrochem.ionic_strength([molarities of each ion], [charges of each ion])
		- mole_fractions = ms3.zs
			- CaCl2_molarity
				- Ca+2_molarity_amount = 1 / thermo.elements.molecular_weight({"Ca": CaCl2.atoms["Ca"]}) * 1000
				- Ca+2_molarity =  2_molarity_amount  * mole_fractions[0]
	[*] - molar mass of mix (mole fraction * molar mass of solvent)
	- stoichiometric coefficients
	- volume fractions of solvents

	- seperate out all the subgroups
		- ions
		- solvents
		- groups

	- use ion equation if only ions are present
		- use group equation if only groups are present?



	[*]-calc and cache
		- b func
		- db_func_di
'''
subgroup_objects = None
subgroup_keys = None
b_inc = None
db_dI = None

def drange(x, y, jump):
  while x < y:
	yield float(x)
	x += decimal.Decimal(jump)

def compute_ionic_stregth_of_ions_in_mix(PeriodicTable, list_of_chemicals):
	mols = []
	chrgs = []
	molarities = {}
	charges = {}
	for molecule in list_of_chemicals:
		for k,v in molecule.atoms.iteritems():
			temp_moles = 1 / float(thermo.elements.molecular_weight({k: v}) * 1000)
			if k in molarities:
				molarities[k] += temp_moles
			else:
				molarities[k] = temp_moles


	for k,v in molarities.iteritems():
		mols.append(molarities[k])
		charges[k] = Chemical(PeriodicTable[k].name).charge
		chrgs.append(charges[k])

	return thermo.electrochem.ionic_strength(mols, chrgs), charges, molarities

def MR_LIQUAC_IONS(chemgroups, PeriodicTable, list_of_chemicals, subgroup_data=None, interaction_data=None):
	global subgroup_objects, subgroup_keys, b_inc, db_dI
	
	subgroups = subgroup_data
	interactions = interaction_data

	group_counts = {}
	for groups in chemgroups:
		for group, count in groups.items():
			if group in group_counts:
				group_counts[group] += count
			else:
				group_counts[group] = count

	# compute ionic stregth of ions in mix 
	ionic_strength, molarities, charges = compute_ionic_stregth_of_ions_in_mix(PeriodicTable, list_of_chemicals)


	# compute b_func and db_dI

	if (subgroup_objects is None) or (subgroup_keys is None) or (b_inc is None) or (db_dI is None):
		subgroup_objects = {}
		subgroup_keys = {}
		b_inc = {i: {j: 0 for j in group_counts} for i in group_counts}
		db_dI = {i: {j: 0 for j in group_counts} for i in group_counts}
		for i in group_counts:
			for j in group_counts:
				main1 = subgroup_data[i].main_group_id
				main2 = subgroup_data[j].main_group_id
				main_c2 = subgroup_data[j].group

				if i not in subgroup_objects:
					main_c1 = subgroup_data[i].group
					c1_chem = Chemical(main_c1)
					subgroup_objects[i] = {"charge": c1_chem.charge, "key": c1_chem.atoms.keys()[0]}
					subgroup_keys[c1_chem.atoms.keys()[0]] = i

				if j not in subgroup_objects:
					main_c1 = subgroup_data[j].group
					c2_chem = Chemical(main_c1)
					subgroup_objects[j] = {"charge": c2_chem.charge, "key": c2_chem.atoms.keys()[0]}
					subgroup_keys[c2_chem.atoms.keys()[0]] = j
				
				a, b, c = 0, 0, 0
				
				try:
					a, b, c = interaction_data[main1][main2]
				except:
					pass

				a2 = 0
				try:
					a2 = 0.25 /  abs(subgroup_objects[i]["charge"] - subgroup_objects[j]["charge"])
				except:
					pass
				

				if ionic_strength != 0:
					sqrt_is = math.sqrt(ionic_strength)
					exp_term = exp((a * sqrt_is) + (a2 * ionic_strength))
					b_inc[i][j] =  b + (c * exp_term)
					db_dI[i][j] = ((a / (2 * sqrt_is)) + a2) * (c * exp_term)
				else:
					b_inc[i][j] = b
					db_dI[i][j] = 0


			


	term_sums = []
	for groups in chemgroups:
		term_3_sum = 0
		term_4_sum = 0
		for group, count in groups.items():
			temp_4_sum = 0
			charge_factor = ((charges[subgroup_objects[group]["key"]] ** 2) / 2.0)
			for j in group_counts:
				term_3_sum += molarities[subgroup_objects[j]["key"]] * b_inc[j][subgroup_keys[subgroup_objects[group]["key"]]]
				for z in group_counts:
					if (charges[subgroup_objects[j]["key"]] > 0) and (charges[subgroup_objects[z]["key"]] < 0):
						temp_4_sum += molarities[subgroup_objects[j]["key"]] * molarities[subgroup_objects[z]["key"]] * db_dI[j][z]
		
			temp_4_sum *= charge_factor
			term_4_sum += temp_4_sum
		term_sums.append(term_3_sum + term_4_sum)

	return term_sums


def calc_vap_pressure(K, delta_Hvap, A):
	r = 8.3144598
	exp_val = - delta_Hvap / (r * K)
	return A * math.exp(exp_val)


r = 8.3144598 #J / K / mol
KNO3 = Chemical('Potassium Nitrate')
NaNO2 = Chemical('Sodium Nitrite')
#NaNO2.Tb = 
NaNO3 = Chemical('Sodium Nitrate')
MgBr2 = Chemical('MgBr2')
MgBr2.Tb = 1520
K3BO3 = Chemical('Potassium borate')
K3BO3.Tb = 1848

CBr4 = Chemical("carbon tetrabromide")
AlB2 = Chemical('aluminium boride')

mmHg_to_Pa = 133.322365
kelvin_offset = 298.15
CaCl2 = Chemical('CaCl2')
#https://m.chemicalbook.com/ChemicalProductProperty_EN_CB3465202.htm
CaCl2.Psat_310 = 0.01 * mmHg_to_Pa
CaCl2.Psat_1602 = (8 * (10**4))
CaCl2.delta_Hvap = math.log(CaCl2.Psat_310 / CaCl2.Psat_1602) / ((310 - 1602) / (310 * 1602)) * r
CaCl2.A_const = math.exp(math.log(CaCl2.Psat_310) + (CaCl2.delta_Hvap / (r * 310)))

MgCl2 = Chemical('MgCl2')
MgCl2.Psat_310 = 0.01 * mmHg_to_Pa
MgCl2.Psat_1387 = (1 * (10**5))
MgCl2.delta_Hvap = math.log(MgCl2.Psat_310 / MgCl2.Psat_1387) / ((310 - 1387 ) / (310 * 1387)) * r
MgCl2.A_const = math.exp(math.log(MgCl2.Psat_310) + (MgCl2.delta_Hvap / (r * 310)))

#ln(P1/P2)=(Hvap/R)(T1-T2/T1xT2)


'''
molten_salt_init_mult = max([0, min([1.0, 1.0])])
ms1 = Mixture(
	[NaNO3.CAS, KNO3.CAS, NaNO2.CAS],
	ws=[
		16.1/100.00 * molten_salt_init_mult,
		54.7/100.00 * molten_salt_init_mult,
		29.2/100.00 * molten_salt_init_mult],
	T=124+kelvin_offset
)
'''
NaNO3_KNO2_NaNO2_chem_group = {
	204: 1,
	203: 2,
	208: 2,
	211: 1
}
MgBr2_chem_group = {
	210: 1,
	207: 2
}
CaCl2_chem_group = {
	205: 1,
	206: 2
}
MgCl2_chem_group = {
	206: 2,
	210: 1
}

chemgroups = [
	CaCl2_chem_group,
	MgCl2_chem_group
]
'''
Basis functions to determine Energy level at deg K:
	- https://www.youtube.com/watch?v=43Kd3yUG5io
	-https://www.youtube.com/watch?v=pu4uL7deCNw&list=PLkNVwyLvX_TFBLHCvApmvafqqQUHb6JwF
		-CaF2
			- determine basis functions for CaF2 (describe varying amplitude of electron density)
				- sp… orbital basis functions
			- determine halitonian for for electron
				- H = Kintic Energy + Nuclear Energy + Electron - Electron Reuplsion
					- hartree - fock
						- slater determinenant (many electronss)
		-MgF2

		- orbitals 
			- Mg2+ = 1s^2 + 2s^2 + 2p^6 + 3s^2 (valence)
			- Ca2+ = 1s^2 + 2s^2 + 2p^6 + 3s^2 + 3p^6 + 4s^2 (valence)
			- F- = 1s^2 + 2s^2 + 2p^5 (valence)
			
	- Kohm- Shan theory
		- spin appoximations * (local spin dnsity appoximation) + gradient density correction (Becke-Lee-Yang-Par)
			- exchnage denisty Energy formula
				- (9 * alpha / 8) * ((3 / pi) ^ (1 / 3)) * rho_func ^(1/3)(r) (use dirac alpha)

	Ideally Use the MN12‐L dft b

	B3LYP
		E_xc = (1 -a) * E_x(lsda) + a * E_x(HF) + b * delta_E_x(B) + (1-c) * E_c(lsda) + c * E_c(lyp)
		a = 0.20
		b = 0.72
		c = 0.81

	density f
		-function of ocupied orbitals
			- orbital basis functions

	
		
'''


'''
	POSSIBLE EUTECTIC POINT
	CaCl2.CAS, MgCl2.CAS
	Chemical Potential 1: 0.0128150678975
	Chemical Potential 2: 0.0161191905054
	Temp: 407.461, mole frac 1: 0.417990272081, mole frac 2: 0.582009727919, weight frac 1: 0.455682 weight frac 2: 0.544318
	 ms3 = Mixture(
		[CaCl2.CAS, MgCl2.CAS],
		ws=[
			0.455682,
			0.544318
		],
		T=407.461
	 )

	range C: 1568.0
		partial_cacl2: 32017.9026159, partial_mgcl2: 64313.5936705, max_temp: 1976.0

ms3.ws=[0.455682,0.544318]
1 * ms3.Cps * (407.461 - 30 - kelvin_offset)
'''
# https://en.wikipedia.org/wiki/Activity_coefficient
# https://en.wikipedia.org/wiki/Chemical_equilibrium
# todo, need to loop over ws, T, the calcluate use coffs to determine where chemical potentials are equal to each other
#100, 1200, 5
#100, 120, 0.01
# drange(decimal.Decimal(407.46), decimal.Decimal(407.48), 0.001):
is_run_boil = True
list_of_chemicals = [CaCl2, MgCl2]
list_of_cas = [x.CAS for x in list_of_chemicals]
PeriodicTable = thermo.elements.PeriodicTable(thermo.elements.element_list).symbol_to_elements
boiling_point = 101325
range_max = 0.95
init_temp = 408
total_vapor_pressure = []
total_vapor_pressure1 = []
total_vapor_pressure2 = []
temperature = []
prev_data = [None,None]
for d_temp in drange(decimal.Decimal(init_temp), decimal.Decimal(4000), 1):
	if is_run_boil:

		ms3 = Mixture(
			[CaCl2.CAS, MgCl2.CAS],
			ws=[
				0.455682,
				0.544318
			],
			T=d_temp
		)
		mole_fractions = ms3.zs
		coeff_offset_CaCl2, coeff_offset_MgCl2 = thermo.unifac.UNIFAC(init_temp,
			xs=mole_fractions,
			chemgroups=chemgroups,
			subgroup_data=thermo.unifac.DOUFSG,
			interaction_data=thermo.unifac.DOUFIP2016,
			modified=True
		)

		coeff_offset_CaCl2_MR, coeff_offset_MgCl2_MR = MR_LIQUAC_IONS(
			chemgroups,
			PeriodicTable,
			list_of_chemicals,
			subgroup_data=thermo.unifac.DOUFSG, 
			interaction_data=thermo.unifac.DOUFIP2016
		)

		SR_MR_CaCl2 = math.exp(math.log(coeff_offset_CaCl2) + coeff_offset_CaCl2_MR)
		SR_MR_MgCl2 = math.exp(math.log(coeff_offset_MgCl2) + coeff_offset_MgCl2_MR)

		partial_cacl2 = calc_vap_pressure(d_temp, CaCl2.delta_Hvap, CaCl2.A_const)
		partial_mgcl2 = calc_vap_pressure(d_temp, MgCl2.delta_Hvap, MgCl2.A_const)

		yi_cacl2 = (mole_fractions[0] * partial_cacl2) / float(boiling_point)
		yi_mgcl2 = (mole_fractions[1] * partial_mgcl2) / float(boiling_point)

		fgacity_cacl2 = (mole_fractions[0] * SR_MR_CaCl2 * partial_cacl2) / (yi_cacl2 * boiling_point)
		fgacity_mgcl2 = (mole_fractions[1] * SR_MR_MgCl2 * partial_mgcl2) / (yi_mgcl2 * boiling_point)

		adj_cacl2 = (mole_fractions[0] * SR_MR_CaCl2 * partial_cacl2)
		adj_mgcl2 = (mole_fractions[1] * SR_MR_MgCl2 * partial_mgcl2)

		#print yi_cacl2, yi_mgcl2, adj_cacl2, adj_mgcl2

		prev_data[0] = "range C: %s" %  (d_temp - init_temp)
		prev_data[1] = "\tpartial_cacl2: %s, partial_mgcl2: %s, max_temp: %s" % (adj_cacl2, adj_mgcl2, d_temp)
		total_vapor_pressure.append(adj_cacl2+adj_mgcl2)
		total_vapor_pressure1.append(adj_cacl2)
		total_vapor_pressure2.append(adj_mgcl2)
		temperature.append(d_temp)

		if (adj_cacl2+adj_mgcl2) > (boiling_point * range_max):
			print prev_data[0]
			print prev_data[1]
			break
		#else:
		

		'''
		partial_cacl2 = calc_vap_pressure(d_temp, CaCl2.delta_Hvap, CaCl2.A_const) * mole_fractions[0]
		partial_mgcl2 = calc_vap_pressure(d_temp, MgCl2.delta_Hvap, MgCl2.A_const) * mole_fractions[1]

		ideal_vap_pressure = (partial_cacl2 + partial_mgcl2)
		print SR_MR_CaCl2, SR_MR_MgCl2
		
		if ideal_vap_pressure > (boiling_point * range_max):
			print "range C:", d_temp - init_temp
			print "\tpp1: %s, pp2: %s, max_temp: %s" % (partial_cacl2, partial_mgcl2, d_temp)
			break
		'''

		#print ideal_vap_pressure, d_temp, partial_cacl2, partial_mgcl2
	else:
		for cacl2_ws in drange(decimal.Decimal(0.45), decimal.Decimal(0.46), 0.000001):
			ms3 = Mixture(
				list_of_cas,
				ws=[
					cacl2_ws,
					1.0 - cacl2_ws
				],
				T=d_temp
			)
			mole_fractions = ms3.zs
			coeff_offset_CaCl2, coeff_offset_MgCl2 = thermo.unifac.UNIFAC(ms3.T,
				xs=mole_fractions,
				chemgroups=chemgroups,
				subgroup_data=thermo.unifac.DOUFSG,
				interaction_data=thermo.unifac.DOUFIP2016,
				modified=True
			)

			coeff_offset_CaCl2_MR, coeff_offset_MgCl2_MR = MR_LIQUAC_IONS(
				chemgroups,
				PeriodicTable,
				list_of_chemicals,
				subgroup_data=thermo.unifac.DOUFSG, 
				interaction_data=thermo.unifac.DOUFIP2016
			)

			SR_MR_CaCl2 = exp(log(coeff_offset_CaCl2) + coeff_offset_CaCl2_MR)
			SR_MR_MgCl2 = exp(log(coeff_offset_MgCl2) + coeff_offset_MgCl2_MR)

			total_cp1 = 0
			total_cp2 = 0
			if (CaCl2.phase == 'l'):
				CaCl2.calculate(T=CaCl2.Tm-0.01)
				total_cp1 += CaCl2.Cpsm * (CaCl2.Tm-0.01)
				try:
					CaCl2.calculate(T=d_temp)
					total_cp1 += CaCl2.Cplm * (d_temp)
				except:
					print "CaCl2", d_temp
					sys.exit()
			elif (CaCl2.phase == 's'):
				CaCl2.calculate(T=d_temp)
				total_cp1 += CaCl2.Cpsm * (d_temp-kelvin_offset)

			if (MgCl2.phase == 'l'):
				MgCl2.calculate(T=MgCl2.Tm-0.01)
				total_cp2 += MgCl2.Cpsm * (MgCl2.Tm-0.01)
				MgCl2.calculate(T=d_temp)
				total_cp2 += MgCl2.Cplm * (d_temp)
			elif (MgCl2.phase == 's'):
				MgCl2.calculate(T=d_temp)
				total_cp2 += MgCl2.Cpsm * (d_temp-kelvin_offset)

			cp1_real = (total_cp1 + (r*(d_temp)*math.log(mole_fractions[0] * SR_MR_CaCl2))) if cacl2_ws > 0 else 0
			cp2_real = (total_cp2 + (r*(d_temp)*math.log(mole_fractions[1] * SR_MR_MgCl2))) if (1.0 - cacl2_ws) > 0 else 0

			if (abs(cp1_real) < 0.1) and (abs(cp2_real) < 0.1):
				print "Chemical Potential 1: %s\nChemical Potential 2: %s\n\tTemp: %s, mole frac 1: %s, mole frac 2: %s, weight frac 1: %s weight frac 2: %s" % (cp1_real, cp2_real, (d_temp), mole_fractions[0], mole_fractions[1], cacl2_ws, 1.0 - cacl2_ws)


if is_run_boil:
	fig = plt.figure(figsize=(14.0, 8.0))
	try:
		close_plot = 'vapor_pressure.png'
		plt.xlabel('Temperature (K)')
		plt.ylabel('Combined Vapor Pressure (Pa)')
		plt.title('Net PNL Occurance Distrubiton')
		plt.grid(True)
		plt.plot(temperature, total_vapor_pressure, c='green', hold=True)
		plt.plot(temperature, total_vapor_pressure1, c='blue', hold=True)
		plt.plot(temperature, total_vapor_pressure2, c='red')
		axes = plt.gca()
		axes.set_ylim([math.floor(0), math.ceil(max(total_vapor_pressure))])
		plt.draw()
		plt.savefig(close_plot, format='png')
		plt.clf()
	except Exception,e:
		print e
		pass