janaf.py - code to calculate janaf thermodynamics (v1.1)

#!/usr/bin/python

import sys, scipy, os, re, math, random
from scipy.optimize import minimize

if len(sys.argv) != 5:
  print "Usage: janaf.py <file> <minimum_temp> <common_temp> <maximum_temp>"
  sys.exit(1)

R = 8.314459848	# J/K/mol

elements = {
  "H" : 1.00794,
  "Li" : 6.941,
  "B" : 10.811,
  "Be" : 9.012182,
  "C" : 12.0107,
  "N" : 14.0067,
  "O" : 15.9994049,
  "F" : 18.998403,
  "Na" : 22.98976,
  "Mg" : 24.3050,
  "Al" : 26.98153851111,
  "Si" : 28.0855,
  "P" : 30.97696,
  "S" : 32.065,
  "Cl" : 35.453,
  "K" : 39.0983,
  "Ca" : 40.078,
  "Cr" : 51.9962,
  "Fe" : 55.845,
  "Co" : 58.93319,
  "Zn" : 65.38,
  "Ga" : 69.723,
  "Br" : 79.904,
  "I" : 126.09044,
}
low = float(sys.argv[2])
common = float(sys.argv[3])
high = float(sys.argv[4])
reference = 298.15
zero = 0.0

def read(category):
  lines = open(sys.argv[1]).readlines()
  state = "Searching"
  data = []
  data2 = []
  index = 0
  name = "UNKNOWN"
  for line in lines:
    if state == "Searching":
      if re.match('^\s*"name":\s*"%s",\s*\n$' % category, line):
        state = "values"
      else:
        match = re.match('^\s*"chemicalFormula":\s*"(.*)",\s*\n$', line)
        if match:
          name = match.group(1)
    elif state == "values":
      if re.match('^\s*"name":\s*"Temperature",\s*\n$', line):
        state = "temperatures"
      else:
        match = re.match('^\s*"value":\s*"(-?\d*\.?\d+)"\n$', line)
        if match:
          data.append(float(match.group(1)))
    elif state == "temperatures":
      if re.match('^\s*"name":\s*".*",\s*\n$', line):
        break
      else:
        match = re.match('^\s*"value":\s*"(-?\d*\.?\d+)"\n$', line)
        if match:
          data2.append( (float(match.group(1)), data[index]) )
          index += 1
    else:
      print "Unknown state!"
      sys.exit(1)
  return data2, name

tmp = read("Specific heat capacity")
specific_heat = tmp[0]
name = tmp[1]
enthalpy = read("Relative standard enthalpy")[0]
entropy = read("Standard entropy")[0]

molweight = 0.0;
keys = elements.keys()
keys.sort(key=len, reverse=True)
element_string = "|".join(keys)
matches = re.findall('(%s)(\d*)' % element_string, name)
for match in matches:
  count = 1
  if match[1] != "":
    count = float(match[1])
  molweight += elements[match[0]] * count

low_point = len(specific_heat)
low_specific_heat = None
for i in range(len(specific_heat)):
  if specific_heat[i][0] >= low - 0.000001:
    low_point = i
    low_specific_heat = specific_heat[i][1]
    break

common_point = len(specific_heat)
common_specific_heat = None
for i in range(len(specific_heat)):
  if specific_heat[i][0] >= common - 0.000001:
    common_point = i
    common_specific_heat = specific_heat[i][1]
    break

high_point = len(specific_heat)
high_specific_heat = None
for i in range(len(specific_heat)):
  if specific_heat[i][0] >= high - 0.000001:
    high_point = i
    high_specific_heat = specific_heat[i][1]
    break

dataset_low  = specific_heat[low_point:common_point] + ([specific_heat[common_point]] * 9)
dataset_high = ([specific_heat[common_point]] * 10) + specific_heat[common_point:high_point]

entropy_low = len(entropy)
entropy_common = len(entropy)
entropy_high = len(entropy)
for i in range(len(specific_heat)):
  if entropy_low == len(entropy):
    if entropy[i][0] >= low - 0.000001:
      entropy_low = entropy[i][1]
  elif entropy_common == len(entropy):
    if entropy[i][0] >= common - 0.000001:
      entropy_common = entropy[i][1]
  else:
    if entropy[i][0] >= high - 0.000001:
      entropy_high = entropy[i][1]
      break

enthalpy_low = len(enthalpy)
enthalpy_common = len(enthalpy)
enthalpy_high = len(enthalpy)
for i in range(len(specific_heat)):
  if enthalpy_low == len(enthalpy):
    if enthalpy[i][0] >= low - 0.000001:
      enthalpy_low = enthalpy[i][1]
  elif enthalpy_common == len(enthalpy):
    if enthalpy[i][0] >= common - 0.000001:
      enthalpy_common = enthalpy[i][1]
  else:
    if enthalpy[i][0] >= high - 0.000001:
      enthalpy_high = enthalpy[i][1]
      break

def calculate_specific_heat(T, a):
  return R * ((((a[4]*T + a[3])*T + a[2])*T + a[1])*T + a[0])

def evaluate(a, dataset):
  total_error = 0.0
#  print a
  for i in dataset:
    T = i[0]
    x = calculate_specific_heat(T, a)
    if x < 0:
      return 999999999999999.9
    error = math.pow(math.fabs(x - i[1]), 2) / len(dataset)
    total_error += error
#    print "  %f: %6f/%6f (%f)" % (T, x, i[1], math.fabs(x - i[1]))
  return total_error

def optimize(dataset):
  params = [1e1, 8e-3, 6e-6, 4e-9, 2e-12]
  best = evaluate(params, dataset)
  i = 0
  while i < 60 * math.sqrt(best):
#    print "Cycle %d" % i
    trial = list(params)
    while True:
      if random.random() < math.pow(1.0 / (len(trial) + 1), 0.4):
        break
      choice = random.randrange(len(trial))
      trial[choice] += trial[choice] * math.pow(random.random()/(random.random() + 0.0000001), 0.3) - trial[choice] * math.pow(random.random()/(random.random() + 0.0000001), 0.3)
      trial[choice] *= math.pow(random.random()/(random.random() + 0.2), 0.3)
    res = minimize(evaluate, list(trial), method='Nelder-Mead', args=(dataset,), options={"maxiter": 1000})
    val = evaluate(res.x, dataset)
#    print "%f, %f" % (val, best)
    if val < best:
      best = val
      params = list(res.x)
#      print "", params
    i += 1
  return best, params

low_results = optimize(dataset_low)
high_results = optimize(dataset_high)

sys.stderr.write("Errors: low=%f, high=%f\n" % (low_results[0], high_results[0]))
if low_results[0] >= 0.03:
  sys.stderr.write("  Warning: low-Cp error is very high - consider raising Tlow.\n")
elif low_results[0] >= 0.01:
  sys.stderr.write("  Warning: low-Cp error is somewhat high - consider raising Tlow.\n")
if high_results[0] >= 0.03:
  sys.stderr.write("  Warning: high-Cp error is very high - consider raising Tcommon or lowering Thigh.\n")
elif high_results[0] >= 0.01:
  sys.stderr.write("  Warning: high-Cp error is somewhat high - consider raising Tcommon or lowering Thigh.\n")

zero_enthalpy = None
for entry in enthalpy:
  if entry[0] >= zero - 0.000001:
    if entry[0] <= zero + 0.000001:
      zero_enthalpy = entry[1]
      break
    else:
      print "ERROR - enthalpy-zero not found: %f, %f" % (zero, entry[0])
      sys.exit(1)

low_enthalpy = None
for entry in enthalpy:
  if entry[0] >= low - 0.000001:
    if entry[0] <= low + 0.000001:
      low_enthalpy = entry[1]
      break
    else:
      print "ERROR - enthalpy-low not found: %f, %f" % (low, entry[0])
      sys.exit(1)

reference_enthalpy = None
for entry in enthalpy:
  if entry[0] >= reference - 0.000001:
    if entry[0] <= reference + 0.000001:
      reference_enthalpy = entry[1]
      break
    else:
      print "ERROR - enthalpy-reference not found: %f, %f" % (reference, entry[0])
      sys.exit(1)

common_enthalpy = None
for entry in enthalpy:
  if entry[0] >= common - 0.000001:
    if entry[0] <= common + 0.000001:
      common_enthalpy = entry[1]
      break
    else:
      print "ERROR - enthalpy-common not found: %f, %f" % (common, entry[0])
      sys.exit(1)

high_enthalpy = None
for entry in enthalpy:
  if entry[0] >= high - 0.000001:
    if entry[0] <= high + 0.000001:
      high_enthalpy = entry[1]
      break
    else:
      print "ERROR - enthalpy-high not found: %f, %f" % (high, entry[0])
      sys.exit(1)

zero_enthalpy *= 1000
low_enthalpy *= 1000
reference_enthalpy *= 1000
common_enthalpy *= 1000
high_enthalpy *= 1000

def calculate_enthalpy(T, a):
  return R * (a[4]*math.pow(T, 5)/5 + a[3]*math.pow(T, 4)/4 + a[2]*math.pow(T, 3)/3 + a[1]*math.pow(T, 2)/2 + a[0]*T)

low_enthalpy_offset = 0.0
high_enthalpy_offset = 0.0
last_temperature = None
last_calculated_enthalpy = None
last_enthalpy_difference = None
last_enthalpy = None

for entry in enthalpy:
  temperature = entry[0]
  cur_enthalpy = entry[1] * 1000
  calculated_enthalpy = None
  enthalpy_difference = None
  if temperature >= low - 0.000001:
    if temperature <= common + 0.000001:
      calculated_enthalpy = calculate_enthalpy(temperature, low_results[1])
      enthalpy_difference = cur_enthalpy - calculated_enthalpy
      if last_enthalpy_difference != None:
        average_enthalpy_difference = (last_enthalpy_difference + enthalpy_difference) / 2
        low_enthalpy_offset += average_enthalpy_difference * (temperature - last_temperature) / R
    if temperature >= common - 0.000001:
      calculated_enthalpy = calculate_enthalpy(temperature, high_results[1])
      enthalpy_difference = cur_enthalpy - calculated_enthalpy
      average_enthalpy_difference = (last_enthalpy_difference + enthalpy_difference) / 2
      if temperature >= common + 0.000001:
        high_enthalpy_offset += average_enthalpy_difference * (temperature - last_temperature) / R
        if temperature >= high - 0.000001:
          break
  last_temperature = temperature
  last_calculated_enthalpy = calculated_enthalpy
  last_enthalpy_difference = enthalpy_difference
  last_enthalpy = cur_enthalpy

low_enthalpy_offset /= (common - low)
high_enthalpy_offset /= (high - common)

zero_entropy = None
for entry in entropy:
  if entry[0] >= zero - 0.000001:
    if entry[0] <= zero + 0.000001:
      zero_entropy = entry[1]
      break
    else:
      print "ERROR - entropy-zero not found: %f, %f" % (zero, entry[0])
      sys.exit(1)

low_entropy = None
for entry in entropy:
  if entry[0] >= low - 0.000001:
    if entry[0] <= low + 0.000001:
      low_entropy = entry[1]
      break
    else:
      print "ERROR - entropy-low not found: %f, %f" % (low, entry[0])
      sys.exit(1)

reference_entropy = None
for entry in entropy:
  if entry[0] >= reference - 0.000001:
    if entry[0] <= reference + 0.000001:
      reference_entropy = entry[1]
      break
    else:
      print "ERROR - entropy-high not found: %f, %f" % (reference, entry[0])
      sys.exit(1)

common_entropy = None
for entry in entropy:
  if entry[0] >= common - 0.000001:
    if entry[0] <= common + 0.000001:
      common_entropy = entry[1]
      break
    else:
      print "ERROR - entropy-high not found: %f, %f" % (common, entry[0])
      sys.exit(1)

high_entropy = None
for entry in entropy:
  if entry[0] >= high - 0.000001:
    if entry[0] <= high + 0.000001:
      high_entropy = entry[1]
      break
    else:
      print "ERROR - entropy-high not found: %f, %f" % (high, entry[0])
      sys.exit(1)

def calculate_entropy(T, a):
  if T == 0.0:
    return 1.0
  return R * (a[4]*math.pow(T, 4)/4 + a[3]*math.pow(T, 3)/3 + a[2]*math.pow(T, 2)/2 + a[1]*T + a[0]*math.log(T))

low_entropy_offset = 0.0
high_entropy_offset = 0.0
last_temperature = None
last_calculated_entropy = None
last_entropy_difference = None

for entry in entropy:
  temperature = entry[0]
  cur_entropy = entry[1]
  calculated_entropy = None
  entropy_difference = None
  if temperature >= low - 0.000001:
    if temperature <= common + 0.000001:
      calculated_entropy = calculate_entropy(temperature, low_results[1])
      entropy_difference = cur_entropy - calculated_entropy
      if last_entropy_difference != None:
        average_entropy_difference = (last_entropy_difference + entropy_difference) / 2
        low_entropy_offset += average_entropy_difference * (temperature - last_temperature) / R
    if temperature >= common - 0.000001:
      calculated_entropy = calculate_entropy(temperature, high_results[1])
      entropy_difference = cur_entropy - calculated_entropy
      average_entropy_difference = (last_entropy_difference + entropy_difference) / 2
      if temperature >= common + 0.000001:
        high_entropy_offset += average_entropy_difference * (temperature - last_temperature) / R
        if temperature >= high - 0.000001:
          break
  last_temperature = temperature
  last_calculated_entropy = calculated_entropy
  last_entropy_difference = entropy_difference

low_entropy_offset /= (common - low)
high_entropy_offset /= (high - common)

sys.stderr.write("Specific heat checks, low:  %9.2f %9.2f / %9.2f %9.2f\n" % (calculate_specific_heat(low, low_results[1]), low_specific_heat, calculate_specific_heat(common, low_results[1]), common_specific_heat))
sys.stderr.write("Specific heat checks, high: %9.2f %9.2f / %9.2f %9.2f\n" % (calculate_specific_heat(common, high_results[1]), common_specific_heat, calculate_specific_heat(high, high_results[1]), high_specific_heat))
sys.stderr.write("Enthalpy checks, low:       %9.2f %9.2f / %9.2f %9.2f\n" % (calculate_enthalpy(low, low_results[1]) + low_enthalpy_offset * R, low_enthalpy, calculate_enthalpy(common, low_results[1]) + low_enthalpy_offset * R, common_enthalpy))
sys.stderr.write("Enthalpy checks, low:       %9.2f %9.2f / %9.2f %9.2f\n" % (calculate_enthalpy(common, high_results[1]) + high_enthalpy_offset * R, common_enthalpy, calculate_enthalpy(high, high_results[1]) + high_enthalpy_offset * R, high_enthalpy))
sys.stderr.write("Entropy checks, low:        %9.2f %9.2f / %9.2f %9.2f\n" % (calculate_entropy(low, low_results[1]) + low_entropy_offset * R, low_entropy, calculate_entropy(common, low_results[1]) + low_entropy_offset * R, common_entropy))
sys.stderr.write("Entropy checks, low:        %9.2f %9.2f / %9.2f %9.2f\n" % (calculate_entropy(common, high_results[1]) + high_entropy_offset * R, common_entropy, calculate_entropy(high, high_results[1]) + high_entropy_offset * R, high_entropy))

print "%s" % name
print "{"
print "    specie"
print "    {"
print "        nMoles          1;"
print "        molWeight       %s;" % molweight
print "    }"
print "    thermodynamics"
print "    {"
print "        Tlow            %f;" % low
print "        Thigh           %f;" % high
print "        Tcommon         %f;" % common
print "        highCpCoeffs    ( %g %g %g %g %g %g %g );" % (high_results[1][0], high_results[1][1], high_results[1][2], high_results[1][3], high_results[1][4], high_enthalpy_offset, high_entropy_offset)
print "        lowCpCoeffs     ( %g %g %g %g %g %g %g );" % (low_results[1][0], low_results[1][1], low_results[1][2], low_results[1][3], low_results[1][4], low_enthalpy_offset, low_entropy_offset)
print "    }"
#print "    transport"
#print "    {"
#print "        As              1.67212e-06;"
#print "        Ts              170.672;"
#print "    }"
print "}"