janaf.py - code to calculate janaf thermodynamics, broken

#!/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

low = float(sys.argv[2])
common = float(sys.argv[3])
high = float(sys.argv[4])

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]

elements = {
  "Al" : 26.98153851111,
  "O" : 15.9994049,
  "H" : 1.00794,
}
molweight = 0.0;
matches = re.findall('(Al|O|H)(\d*)', name)
for match in matches:
  count = 1
  if match[1] != "":
    count = float(match[1])
  molweight += elements[match[0]] * count

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

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

high_point = len(specific_heat)
for i in range(len(specific_heat)):
  if specific_heat[i][0] >= high - 0.000001:
    high_point = i
    high_value = 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 evaluate(a, dataset, verbose = False):
  total_error = 0.0
  for i in dataset:
    T = i[0]
    x = R * ((((a[4]*T + a[3])*T + a[2])*T + a[1])*T + a[0])
    if x < 0:
      return 999999999999999.9
    error = math.pow(math.fabs(x - i[1]), 2)
    total_error += error
    if verbose:
      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]))

#H2
#high_results = (0, [2.99142, 0.000700064, -5.63383e-08, -9.23158e-12, 1.58275e-15])
#low_results  = (0, [3.29812, 0.000824944, -8.14302e-07, -9.47543e-11, 4.13487e-13])
#print "Should be: -835.034 -1.35511 / -1012.52 -3.29409"

#HCl
#high_results = (0, [2.75534, 0.00147358, -4.97125e-07, 8.10866e-11, -5.07206e-15]);
#low_results  = (0, [3.33853, 0.00126821, -3.66692e-06, 4.70399e-09, -1.83601e-12]);
#print "Should be: -11918.1 6.51512 / -12131.5 3.19355"

#SiC
#high_results = (0, [5.02427, -0.000492089, 3.10931e-07, -6.90134e-11, 5.21573e-15]);
#low_results  = (0, [2.42781, 0.00955193, -2.79663e-06, -1.36001e-08, 9.19632e-12]);
#print "Should be: 85310.3 -2.4788 / 85465.1 9.17925"

low_enthalpy = None
for entry in enthalpy:
  if entry[0] >= common - 0.000001:
    if entry[0] <= common + 0.000001:
      low_enthalpy = entry[1]
      break
    else:
      print "ERROR - enthalpy-low 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)

low_enthalpy *= 1000	# kJ -> J
high_enthalpy *= 1000	#

T = common
a = low_results[1]
calculated_low_enthalpy  = R * ((((a[4]*T/5 + a[3]/4)*T + a[2]/3)*T + a[1]/2)*T + a[0])*T
T = high
a = high_results[1]
calculated_high_enthalpy = R * ((((a[4]*T/5 + a[3]/4)*T + a[2]/3)*T + a[1]/2)*T + a[0])*T

#print calculated_high_enthalpy, high_enthalpy, (calculated_high_enthalpy - high_enthalpy)

low_enthalpy_offset  = (low_enthalpy - calculated_low_enthalpy) / R
high_enthalpy_offset = (high_enthalpy - calculated_high_enthalpy) / R

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

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

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

T1 = common
T2 = low
a = low_results[1]
calculated_low_entropy  = R*(a[4]*(math.pow(T1, 4) - math.pow(T2, 4))/4 + a[3] * (math.pow(T1, 3) - math.pow(T2, 3))/3 + a[2] * (math.pow(T1, 2) - math.pow(T2, 2))/2 + a[1] * (T1 - T2)+a[0] * math.log(T1 / T2))
T1 = high
T2 = common
a = high_results[1]
calculated_high_entropy = R*(a[4]*(math.pow(T1, 4) - math.pow(T2, 4))/4 + a[3] * (math.pow(T1, 3) - math.pow(T2, 3))/3 + a[2] * (math.pow(T1, 2) - math.pow(T2, 2))/2 + a[1] * (T1 - T2)+a[0] * math.log(T1 / T2))

low_entropy = common_entropy_val - low_entropy_val
high_entropy = high_entropy_val - common_entropy_val

#print calculated_high_entropy, high_entropy, (calculated_high_entropy - high_entropy)

low_entropy_offset  = (low_entropy - calculated_low_entropy) / R / -(1.0/common - 1.0/low)
high_entropy_offset = (high_entropy - calculated_high_entropy) / R / -(1.0/high - 1.0/common)

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 "}"