Lattice Energy Calculation via Born-Lande Equation. Uses Sympy with Custom Units. Outputs LaTex

LatticeEnergy_Born-Land-Equation.py

#!/usr/bin/env python3
from pprint import pprint

import sympy
from sympy import physics as sympy_physics
# https://github.com/sympy/sympy/tree/master/sympy
# /usr/local/lib/python3.8/site-packages/sympy/
# https://docs.sympy.org/latest/tutorial/printing.html
# https://docs.sympy.org/latest/modules/physics/units/examples.html?highlight=set_quantity_dimension#dimensional-analysis
# https://github.com/sympy/sympy/wiki/Unit-systems

import chempy
# https://bjodah.github.io/chempy/latest/#working-with-units

# https://chem.libretexts.org/Bookshelves/Inorganic_Chemistry/Modules_and_Websites_(Inorganic_Chemistry)/Crystal_Lattices/Thermodynamics_of_Lattices/The_Born-Lande'_equation
# https://chem-is-you.blogspot.com/2013/01/solid-state-chemistry-and-standard.html


def build_avogadros_constant():
	avogadro_number = sympy_physics.units.quantities.Quantity( "avogadros_number" ,
		abbrev="avo's num" ,
		latex_repr="\\left(\\frac{6.02214076\\cdot 10^{23}}{mol}\\right)"  ,
		)
	sympy_physics.units.systems.SI.set_quantity_dimension( avogadro_number , sympy_physics.units.Dimension( ( float( 6.02214076e23 ) / sympy_physics.units.mol ) ) )
	sympy.physics.units.systems.SI.set_quantity_scale_factor( avogadro_number , float( 6.02214076e23 ) )
	return avogadro_number

# \left(1.602176634\cdot \:10^{-19}\right)
def build_elementary_charge_constant():
	elementary_charge = sympy_physics.units.quantities.Quantity( "elementary_charge" ,
		abbrev="eCharge" ,
		latex_repr="\\left(1.602176634\\cdot 10^{-19}\\cdot Coulomb\\right)"  ,
		)
	sympy_physics.units.systems.SI.set_quantity_dimension( elementary_charge , sympy_physics.units.Dimension( sympy_physics.units.charge ) )
	sympy.physics.units.systems.SI.set_quantity_scale_factor( elementary_charge , ( float( 1.602176634e-19 ) * sympy_physics.units.coulomb ) )
	return elementary_charge

def build_pi_constant():
	pi_constant = sympy_physics.units.quantities.Quantity( "pi_constant" ,
		abbrev="pi" ,
		latex_repr="\\pi"  ,
		)
	sympy_physics.units.systems.SI.set_quantity_dimension( pi_constant , sympy_physics.units.Dimension( sympy.pi ) )
	sympy.physics.units.systems.SI.set_quantity_scale_factor( pi_constant , sympy.pi )
	return pi_constant

def build_positive_charge_symbol( charge_value=1 ):
	positive_charge_symbol = sympy_physics.units.quantities.Quantity( f"|z^+{charge_value}|" ,
		abbrev=f"|z^+{charge_value}|" ,
		latex_repr="\\left|z^{+"+ str(charge_value) + "}\\right|"  ,
		)
	sympy_physics.units.systems.SI.set_quantity_dimension( positive_charge_symbol , sympy_physics.units.Dimension( "z^+1" ) )
	sympy.physics.units.systems.SI.set_quantity_scale_factor( positive_charge_symbol , sympy.Abs( charge_value ) )
	return positive_charge_symbol

def build_negative_charge_symbol( charge_value=1 ):
	negative_charge_symbol = sympy_physics.units.quantities.Quantity( f"|z^+{charge_value}|" ,
		abbrev=f"|z^-{charge_value}|" ,
		latex_repr="\\left|z^{-"+ str(charge_value) + "}\\right|"  ,
		)
	sympy_physics.units.systems.SI.set_quantity_dimension( negative_charge_symbol , sympy_physics.units.Dimension( "z^-1" ) )
	sympy.physics.units.systems.SI.set_quantity_scale_factor( negative_charge_symbol , sympy.Abs( charge_value ) )
	return negative_charge_symbol

def build_ro_distance_symbol( ro_distance=1 ):
	if ro_distance == 1:
		latex = "r_o"
	else:
		latex = f"\\left({ro_distance}\\cdot r_o\\right)"
	ro_distance_symbol = sympy_physics.units.quantities.Quantity( "ro" ,
		abbrev=f"{ro_distance}*ro" ,
		latex_repr=latex  ,
		)
	sympy_physics.units.systems.SI.set_quantity_dimension( ro_distance_symbol , sympy_physics.units.Dimension( "r_o" ) )
	sympy.physics.units.systems.SI.set_quantity_scale_factor( ro_distance_symbol , ro_distance )
	return ro_distance_symbol

def build_eo_distance_symbol( eo_distance=1 ):
	if eo_distance == 1:
		latex = "e_o"
	else:
		latex = f"\\left({eo_distance}\\cdot e_o\\right)"
	eo_distance_symbol = sympy_physics.units.quantities.Quantity( "eo" ,
		abbrev=f"{eo_distance}*ro" ,
		latex_repr=latex  ,
		)
	sympy_physics.units.systems.SI.set_quantity_dimension( eo_distance_symbol , sympy_physics.units.Dimension( "e_o" ) )
	sympy.physics.units.systems.SI.set_quantity_scale_factor( eo_distance_symbol , eo_distance )
	return eo_distance_symbol

def build_born_exponent_symbol( value=1 ):
	if value == 1:
		latex = "n"
	else:
		latex = f"{value}"
	borne_exponent_symbol = sympy_physics.units.quantities.Quantity( "born_exponent" ,
		abbrev="n" ,
		latex_repr=latex ,
		)
	sympy_physics.units.systems.SI.set_quantity_dimension( borne_exponent_symbol , sympy_physics.units.Dimension( "n" ) )
	sympy.physics.units.systems.SI.set_quantity_scale_factor( borne_exponent_symbol , value )
	return borne_exponent_symbol

def build_madelung_constant( value=1 ):
	if value == 1:
		latex = "M"
	else:
		latex = f"{value}"
	madelung_constant = sympy_physics.units.quantities.Quantity( "madelung_constant" ,
		abbrev="M" ,
		latex_repr=latex ,
		)
	sympy_physics.units.systems.SI.set_quantity_dimension( madelung_constant , sympy_physics.units.Dimension( "M" ) )
	sympy.physics.units.systems.SI.set_quantity_scale_factor( madelung_constant , value )
	return madelung_constant

def lattice_energy_via_born_land_equation( molecule="Na" ):
	elementary_charge_constant = build_elementary_charge_constant()
	avogadros_constant = build_avogadros_constant()
	pi_constant = build_pi_constant()
	madelung_constant = build_madelung_constant( 1 )
	z_positive = build_positive_charge_symbol( 4 )
	z_negative = build_negative_charge_symbol( 3 )
	ro_distance = build_ro_distance_symbol( 1 )
	eo_distance = build_eo_distance_symbol( 23 )
	borne_exponent = build_born_exponent_symbol( 1 )

	part_one = ( ( -1 * sympy.Pow( elementary_charge_constant , 2 ) ) * sympy.Abs( z_positive ) * sympy.Abs( z_negative ) ) / ( 4 * pi_constant * eo_distance * ro_distance )
	part_two = ( ( madelung_constant ) * ( avogadros_constant ) * ( 1 - ( 1 / borne_exponent ) ) )
	equation = sympy.Mul(
		( ( part_one ) * ( part_two ) ) ,
		evaluate=True
		#evaluate=False
		)
	return equation


if __name__ == "__main__":
	sympy.init_printing()
	equation = lattice_energy_via_born_land_equation()
	equation_latex = f"\\Delta H_{{lattice}} = {sympy.latex( equation )}"
	print( equation_latex )
	pprint( equation )


#Tomas Barfod Pulsing (feat. Nina K)