<commit_before>#!/usr/bin/python
# -*- coding: utf-8 -*-

################################################################################
#
#	RMG - Reaction Mechanism Generator
#
#	Copyright (c) 2002-2009 Prof. <NAME> (<EMAIL>) and the
#	RMG Team (<EMAIL>)
#
#	Permission is hereby granted, free of charge, to any person obtaining a
#	copy of this software and associated documentation files (the 'Software'),
#	to deal in the Software without restriction, including without limitation
#	the rights to use, copy, modify, merge, publish, distribute, sublicense,
#	and/or sell copies of the Software, and to permit persons to whom the
#	Software is furnished to do so, subject to the following conditions:
#
#	The above copyright notice and this permission notice shall be included in
#	all copies or substantial portions of the Software.
#
#	THE SOFTWARE IS PROVIDED 'AS IS', WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
#	IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
#	FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
#	AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
#	LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
#	FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
#	DEALINGS IN THE SOFTWARE.
#
################################################################################

"""
A module for working with the thermodynamics of chemical species. This module
seeks to provide functionality for answering the question, "Given a species,
what are its thermodynamics?"
"""

import quantities as pq
import constants
import data
import math
import scipy
from scipy import linalg
from scipy import optimize

################################################################################

class ThermoData:
	"""
	A base class for all forms of thermodynamic data used by RMG. The common
	attributes are:
	
	========= ============================================================
	Attribute Meaning
	========= ============================================================
	`Trange`  a 2-tuple containin the minimum and maximum temperature in K
	`comment` a string describing the source of the data
	========= ============================================================
	"""
	
	def __init__(self, Trange=None, comment=''):
		self.Trange = Trange or (0.0, 0.0)
		self.comment = comment
	
	def isTemperatureValid(self, T):
		"""
		Return :data:`True` if the temperature `T` in K is within the valid
		temperature range of the thermodynamic data, or :data:`False` if not.
		"""
		if self.Trange is None:
			return True
		else:
			return self.Trange[0] <= T and T <= self.Trange[1]

################################################################################

class ThermoGAData(ThermoData):
	"""
	A set of thermodynamic parameters as determined from Benson's group
	additivity data. The attributes are:
	
	========= ========================================================
	Attribute Meaning
	========= ========================================================
	`H298`    the standard enthalpy of formation at 298 K in J/mol
	`S298`    the standard entropy of formation at 298 K in J/mol*K
	`Cp`      the standard heat capacity in J/mol*K at 300, 400, 500,
		      600, 800, 1000, and 1500 K
	========= ========================================================
	"""
	# I think putting a comment with a colon just before the thing is defined
	# makes it show up in the documentation with autodoc. (or is that just in constants.py?)
	#: The list of temperatures at which heat capacity is defined = [300,400,500,600,800,1000,1500]
	CpTlist = list(pq.Quantity([300.0, 400.0, 500.0, 600.0, 800.0, 1000.0, 1500.0], 'K').simplified)
	CpTlist = [float(T) for T in CpTlist]
	# refer to it as ThermoGAData.CpTlist, even when calling it from methods within this Class.
	
	def __init__(self, H298=0.0, S298=0.0, Cp=None, comment='', index=''):
		"""Initialize a set of group additivity thermodynamic data."""
		ThermoData.__init__(self, Trange=(298.0, 2500.0), comment=comment)
		self.H298 = H298
		self.S298 = S298
		self.Cp = Cp or [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
		self.index = index
		self.__cache_T=None
		self.__cache = dict()
	
	def __add__(self, other):
		"""
		Add two sets of thermodynamic data together. All parameters are
		considered additive. Returns a new :class:`ThermoGAData` object that is
		the sum of the two sets of thermodynamic data.
		"""
		new = ThermoGAData()
		new.H298 = self.H298 + other.H298
		new.S298 = self.S298 + other.S298
		new.Cp = [self.Cp[i] + other.Cp[i] for i in range(len(self.Cp))]
		if self.comment == '': new.comment = other.comment
		elif other.comment == '': new.comment = self.comment
		else: new.comment = self.comment + '+ ' + other.comment
		new.index = self.index + '+' + other.index
		return new
	
	def __repr__(self):
		string = 'ThermoGAData(H298=%s, S298=%s, Cp=%s, index="%s")'%(self.H298,self.S298,self.Cp,self.index)
		return string
	
	def __str__(self):
		"""
		Return a string summarizing the thermodynamic data.
		"""
		string = ''
		string += 'Enthalpy of formation: %s J/mol\n' % (self.H298)
		string += 'Entropy of formation: %s J/mol*K\n' % (self.S298)
		string += 'Heat capacity (J/mol*K): '
		for T, Cp in zip(ThermoGAData.CpTlist, self.Cp):
			string += '%s(%sK) ' % (Cp,T)
		string += '\n'
		string += 'Index: %s\tComment: %s' % (self.index, self.comment)
		return string
	
	def getHeatCapacity(self, T):
		"""
		Return the heat capacity in J/mol*K at temperature `T` in K.
		"""
		if not self.isTemperatureValid(T):
			raise data.TemperatureOutOfRangeException('Invalid temperature for heat capacity estimation from group additivity.')
		if T < 300.0:
			return self.Cp[0]
		elif T > ThermoGAData.CpTlist[-1]:
			return self.Cp[-1]
		else:
			for Tmin, Tmax, Cpmin, Cpmax in zip(ThermoGAData.CpTlist[:-1], \
					ThermoGAData.CpTlist[1:], self.Cp[:-1], self.Cp[1:]):
				if Tmin <= T and T <= Tmax:
					return (Cpmax - Cpmin) * ((T - Tmin) / (Tmax - Tmin)) + Cpmin
	
	def getEnthalpy(self, T):
		"""
		Return the enthalpy in J/mol at temperature `T` in K.
		"""	
		if self.__cache_T==T:
			try: return self.__cache['H']
			except KeyError: pass # 'H' isn't cached yet
		else:  # temperature has changed
			self.__cache = dict()
			self.__cache_T = T
			
		H = self.H298
		if not self.isTemperatureValid(T):
			raise data.TemperatureOutOfRangeException('Invalid temperature for enthalpy estimation from group additivity.')
		for Tmin, Tmax, Cpmin, Cpmax in zip(ThermoGAData.CpTlist[:-1], \
				ThermoGAData.CpTlist[1:], self.Cp[:-1], self.Cp[1:]):
			if T > Tmin:
				slope = (Cpmax - Cpmin) / (Tmax - Tmin)
				intercept = (Cpmin * Tmax - Cpmax * Tmin) / (Tmax - Tmin)
				if T < Tmax:	H += 0.5 * slope * (T*T - Tmin*Tmin) + intercept * (T - Tmin)
				else:			H += 0.5 * slope * (Tmax*Tmax - Tmin*Tmin) + intercept * (Tmax - Tmin)
		if T > ThermoGAData.CpTlist[-1]:
			H += self.Cp[-1] * (T - ThermoGAData.CpTlist[-1])
		self.__cache['H'] = H
		return H
	
	def getEntropy(self, T):
		"""
		Return the entropy in J/mol*K at temperature `T` in K.
		"""
		import math
		
		if self.__cache_T==T:
			try: return self.__cache['S']
			except KeyError: pass # 'S' isn't cached yet
		else:  # temperature has changed
			self.__cache = dict()
			self.__cache_T = T
			
		S = self.S298
		if not self.isTemperatureValid(T):
			raise data.TemperatureOutOfRangeException('Invalid temperature for entropy estimation from group additivity.')
		for Tmin, Tmax, Cpmin, Cpmax in zip(ThermoGAData.CpTlist[:-1], \
				ThermoGAData.CpTlist[1:], self.Cp[:-1], self.Cp[1:]):
			if T > Tmin:
				slope = (Cpmax - Cpmin) / (Tmax - Tmin)
				intercept = (Cpmin * Tmax - Cpmax * Tmin) / (Tmax - Tmin)
				if T < Tmax:	S += slope * (T - Tmin) + intercept * math.log(T/Tmin)
				else:			S += slope * (Tmax - Tmin) + intercept * math.log(Tmax/Tmin)
		if T > ThermoGAData.CpTlist[-1]:
			S += self.Cp[-1] * math.log(T / ThermoGAData.CpTlist[-1])
		self.__cache['S'] = S
		return S
	
	def getFreeEnergy(self, T):
		"""
		Return the Gibbs free energy in J/mol at temperature `T` in K.
		"""
		if self.__cache_T==T:
			try: return self.__cache['G']
			except KeyError: pass # 'G' isn't cached yet
		else:  # temperature has changed
			self.__cache = dict()
			self.__cache_T = T
			
		if not self.isTemperatureValid(T):
			raise data.TemperatureOutOfRangeException('Invalid temperature for free energy estimation from group additivity.')
		G = self.getEnthalpy(T) - T * self.getEntropy(T)
		self.__cache['G'] = G
		return G
	
	def fromDatabase(self, data, comment):
		"""
		Process a list of numbers `data` and associated description `comment`
		generated while reading from a thermodynamic database.
		"""
		
		if len(data) != 12:
			raise Exception('Invalid list of thermo data; should be a list of numbers of length 12.')
		
		H298, S298, Cp300, Cp400, Cp500, Cp600, Cp800, Cp1000, Cp1500, \
			dH, dS, dCp = data
		
		self.H298 = float(pq.Quantity(H298, 'kcal/mol').simplified)
		self.S298 = float(pq.Quantity(S298, 'cal/(mol*K)').simplified)
		self.Cp = list(pq.Quantity([Cp300, Cp400, Cp500, Cp600, Cp800, Cp1000, Cp1500], 'cal/(mol*K)').simplified)
		for i in range(len(self.Cp)): self.Cp[i] = float(self.Cp[i])
		self.comment = comment
	
	def toXML(self, dom, root):
		"""
		Generate an XML representation of the thermodynamic data using the
		:data:`xml.dom.minidom` package. The `dom` and `root` parameters
		represent the DOM and the element in the DOM used as the parent of
		the generated XML.
		"""
		
		thermo = dom.createElement('thermodynamics')
		thermo.setAttribute('index', self.index)
		thermo.setAttribute('comment', self.comment)
		root.appendChild(thermo)
		
		enthalpy = dom.createElement('enthalpyOfFormation')
		enthalpy.setAttribute('temperature', '298.0 K')
		thermo.appendChild(enthalpy)
		data.createXMLQuantity(dom, enthalpy, self.H298, 'J/mol')
		
		entropy = dom.createElement('entropyOfFormation')
		entropy.setAttribute('temperature', '298.0 K')
		thermo.appendChild(entropy)
		data.createXMLQuantity(dom, entropy, self.S298, 'J/(mol*K)')
		
		for i, Cp in enumerate(self.Cp):
			heatCapacity = dom.createElement('heatCapacity')
			heatCapacity.setAttribute('temperature', '%s K' % (ThermoGAData.CpTlist[i]) )
			thermo.appendChild(heatCapacity)
			data.createXMLQuantity(dom, heatCapacity, Cp, 'J/(mol*K)')
	

################################################################################

class ThermoNASAPolynomial(ThermoData):
	"""
	A single NASA polynomial for thermodynamic data. The `coeffs` attribute
	stores the seven polynomial coefficients
	:math:`\\mathbf{a} = \\left[a_1\\ a_2\\ a_3\\ a_4\\ a_5\\ a_6\\ a_7 \\right]`
	from which the relevant thermodynamic parameters are evaluated via the
	expressions
	
	.. math:: \\frac{C_\\mathrm{p}(T)}{R} = a_1 + a_2 T + a_3 T^2 + a_4 T^3 + a_5 T^4
	
	.. math:: \\frac{H(T)}{RT} = a_1 + \\frac{1}{2} a_2 T + \\frac{1}{3} a_3 T^2 + \\frac{1}{4} a_4 T^3 + \\frac{1}{5} a_5 T^4 + \\frac{a_6}{T}
	
	.. math:: \\frac{S(T)}{R} = a_1 \\ln T + a_2 T + \\frac{1}{2} a_3 T^2 + \\frac{1}{3} a_4 T^3 + \\frac{1}{4} a_5 T^4 + a_7
	
	The above was adapted from `this page <http://www.me.berkeley.edu/gri-mech/data/nasa_plnm.html>`_.
	"""
	
	def __init__(self, T_range=None, coeffs=None, comment=''):
		ThermoData.__init__(self, Trange=T_range, comment=comment)
		self.coeffs = coeffs or (0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0)
	def __repr__(self):
                return "ThermoNASAPolynomial(%r,%r,'%s')"%(self.Trange,self.coeffs,self.comment)
	def getHeatCapacity(self, T):
		"""
		Return the heat capacity in J/mol*K at temperature `T` in K.
		"""
		import constants
		if not self.isTemperatureValid(T):
			raise data.TemperatureOutOfRangeException('Invalid temperature for heat capacity estimation from NASA polynomial.')
		c = self.coeffs
		T2 = T*T
		# Cp/R = a1 + a2 T + a3 T^2 + a4 T^3 + a5 T^4
#		HeatCapacityOverR = c[0] + c[1]*T + c[2]*T*T + c[3]*T*T*T + c[4]*T*T*T*T
		HeatCapacityOverR = c[0] + T*(c[1] + T*(c[2] + T*(c[3] + c[4]*T)))
		Cp = HeatCapacityOverR * constants.R
		return Cp
	
	def getEnthalpy(self, T):
		"""
		Return the enthalpy in J/mol at temperature `T` in K.
		"""
		import constants
		if not self.isTemperatureValid(T):
			raise data.TemperatureOutOfRangeException('Invalid temperature for enthalpy estimation from NASA polynomial.')
		c = self.coeffs
		T2 = T*T
		T4 = T2*T2
		# H/RT = a1 + a2 T /2 + a3 T^2 /3 + a4 T^3 /4 + a5 T^4 /5 + a6/T
		EnthalpyOverR = ( c[0]*T + c[1]*T2/2 + c[2]*T2*T/3 + c[3]*T4/4 +
						  c[4]*T4*T/5 + c[5] )
		H = EnthalpyOverR * constants.R
		return H
	
	def getEntropy(self, T):
		"""
		Return the entropy in J/mol*K at temperature `T` in K.
		"""
		import math
		import constants
		if not self.isTemperatureValid(T):
			raise data.TemperatureOutOfRangeException('Invalid temperature for entropy estimation from NASA polynomial.')
		# S/R  = a1 lnT + a2 T + a3 T^2 /2 + a4 T^3 /3 + a5 T^4 /4 + a7
		c = self.coeffs
		T2 = T*T
		T4 = T2*T2
		EntropyOverR = ( c[0]*math.log(T) + c[1]*T + c[2]*T2/2 +
					c[3]*T2*T/3 + c[4]*T4/4 + c[6] )
		S = EntropyOverR * constants.R
		return S
	
	def getFreeEnergy(self, T):
		"""
		Return the Gibbs free energy in J/mol at temperature `T` in K.
		"""
		if not self.isTemperatureValid(T):
			raise data.TemperatureOutOfRangeException('Invalid temperature for free energy estimation from NASA polynomial.')
		return self.getEnthalpy(T) - T * self.getEntropy(T)
	
	def toXML(self, dom, root):
### prime-like:
#   <polynomial>
#      <validRange>
#         <bound kind="lower" property="temperature" units="K">300.000</bound>
#         <bound kind="upper" property="temperature" units="K">1000</bound>
#      </validRange>
#      <coefficient id="1" label="a1">4.0733</coefficient>
#      <coefficient id="2" label="a2">0.011308</coefficient>
#      <coefficient id="3" label="a3">-1.6565e-005</coefficient>
#      <coefficient id="4" label="a4">1.1784e-008</coefficient>
#      <coefficient id="5" label="a5">-3.3006e-012</coefficient>
#      <coefficient id="6" label="a6">-19054.23</coefficient>
#      <coefficient id="7" label="a7">4.4083</coefficient>
#   </polynomial>
		pass
	

################################################################################

class ThermoNASAData(ThermoData):
	"""
	A set of thermodynamic parameters given by NASA polynomials. This class
	stores a list of :class:`ThermoNASAPolynomial` objects in the `polynomials`
	attribute. When evaluating a thermodynamic quantity, a polynomial that
	contains the desired temperature within its valid range will be used.
	"""
	
	def __init__(self, polynomials=None, comment='', Trange=None):
		ThermoData.__init__(self, Trange=Trange, comment=comment)
		self.polynomials = polynomials or []
	
	def addPolynomial(self, polynomial):
		if not isinstance(polynomial,ThermoNASAPolynomial):
			raise TypeError("Polynomial attribute should be instance of ThermoNASAPolynomial")
		self.polynomials.append(polynomial)
	
	def selectPolynomialForTemperature(self, T):
		for poly in self.polynomials:
			if poly.isTemperatureValid(T): break
		else:
			raise data.TemperatureOutOfRangeException("No polynomial found for T=%s" % T)
		return poly
	
	def __repr__(self):
		return "ThermoNASAData(%s, '%s')"%(repr(self.polynomials),self.comment)
	
	def toXML(self, dom, root):
### prime-like: 
# <thermodynamicPolynomials type="nasa7">
#   <referenceState>
#      <Tref units="K">298.15</Tref>
#      <Pref units="Pa">100000</Pref>
#   </referenceState>
#   <dfH units="J/mol">-145177.5731</dfH>
#   <polynomial> FROM NASA_polynomials </polynomial>
# </thermodynamicPolynomials>
		pass
	
	def getHeatCapacity(self, T):
		"""
		Return the heat capacity in J/mol*K at temperature `T` in K.
		"""
		poly = self.selectPolynomialForTemperature(T)
		return poly.getHeatCapacity(T)
	
	def getEnthalpy(self, T):
		"""
		Return the enthalpy in J/mol at temperature `T` in K.
		"""
		poly = self.selectPolynomialForTemperature(T)
		return poly.getEnthalpy(T)
	
	def getEntropy(self, T):
		"""
		Return the entropy in J/mol*K at temperature `T` in K.
		"""
		poly = self.selectPolynomialForTemperature(T)
		return poly.getEntropy(T)
	
	def getFreeEnergy(self, T):
		"""
		Return the Gibbs free energy in J/mol at temperature `T` in K.
		"""
		poly = self.selectPolynomialForTemperature(T)
		return poly.getFreeEnergy(T)
	

################################################################################

class ThermoWilhoitData(ThermoData):
	"""
	A set of thermodynamic parameters given by Wilhoit polynomials. 
	
	"""
	B = 500.0 # Kelvin. Default temperature.
	
	def __init__(self, cp0, cpInf, a0, a1, a2, a3, I, J, comment='', ):
		"""Initialise the Wilhoit polynomial. Trange is set to (0,9999.9)"""
		Trange = (0,9999.9) # Wilhoit valid over all temperatures
		ThermoData.__init__(self, Trange=Trange, comment=comment)
		self.cp0 = cp0
		self.cpInf = cpInf
		self.B = ThermoWilhoitData.B
		self.a0 = a0
		self.a1 = a1
		self.a2 = a2
		self.a3 = a3
		self.I = I
		self.J = J
	
	def __repr__(self):
		return "ThermoWilhoitData(%.2g,%.2g,%.2g,%.2g,%.2g,%.2g,%.2g,%.2g,%.2g'%s')"%(self.cp0, self.cpInf, self.B, self.a0, self.a1, self.a2, self.a3, self.I, self.J, self.comment)
	
	def toXML(self, dom, root):
		pass
	
	def getHeatCapacity(self, T):
		"""
		Return the heat capacity in J/mol*K at temperature `T` in K.
		"""
		y = T/(T+self.B)
		Cp = self.cp0+(self.cpInf-self.cp0)*y*y*( 1 + 
			(y-1)*(self.a0 + y*(self.a1 + y*(self.a2 + y*self.a3))) )
		return Cp
	
	def getEnthalpy(self, T):
		"""
		Return the enthalpy in J/mol at temperature `T` in K.
		"""
		#cp0 = self.cp0 * R
		#cpInf = self.cpInf * R
		cp0 = self.cp0
		cpInf = self.cpInf
		B = self.B
		a0 = self.a0
		a1 = self.a1
		a2 = self.a2
		a3 = self.a3
		I = self.I
		
		enthalpy = I + WilhoitInt0(cp0, cpInf, B, a0, a1, a2, a3, T)
		
		return enthalpy
	
	def getEntropy(self, T):
		"""
		Return the entropy in J/mol*K at temperature `T` in K.
		"""
		#cp0 = self.cp0 * R
		#cpInf = self.cpInf * R
		cp0 = self.cp0
		cpInf = self.cpInf
		B = self.B
		a0 = self.a0
		a1 = self.a1
		a2 = self.a2
		a3 = self.a3
		J = self.J
		
		entropy = J + WilhoitIntM1(cp0, cpInf, B, a0, a1, a2, a3, T)
		
		return entropy
	
	def getFreeEnergy(self, T):
		"""
		Return the Gibbs free energy in J/mol at temperature `T` in K.
		"""
		return self.getEnthalpy(T) - T*self.getEntropy(T)
	

###########
def convertGAtoWilhoit(GAthermo, atoms, rotors, linear):
	"""Convert a Group Additivity thermo instance into a Wilhoit thermo instance.
	
	Takes a `ThermoGAData` instance of themochemical data, and some extra information 
	about the molecule used to calculate high- and low-temperature limits of Cp.
	These are the number of atoms (integer `atoms`), the number of rotors (integer `rotors`)
	and whether the molecule is linear (boolean `linear`)
	Returns a `ThermoWilhoitData` instance.
	
	cf. <NAME>'s thesis, p. 185-186
	"""
	
	# get info from incoming group additivity thermo
	H298 = GAthermo.H298
	S298 = GAthermo.S298
	Cp_list = GAthermo.Cp
	T_list = ThermoGAData.CpTlist  # usually [300, 400, 500, 600, 800, 1000, 1500] but why assume?
	R = constants.R
	B = ThermoWilhoitData.B # Constant, set once in the class def.
	
	# convert from K to kK
	T_list = [t/1000. for t in T_list] 
	B = B/1000.  
	
	Cp_list = [x/R for x in Cp_list] # convert to Cp/R
	
	(cp0, cpInf) = CpLimits(atoms, rotors, linear) # determine the heat capacity limits (non-dimensional)
	
	#create matrices for linear least squares problem
	m = len(T_list)  # probably m=7
	# A = mx4
	# b = mx1
	# x = 4x1
	A = scipy.zeros([m,4])
	b = scipy.zeros([m,1])
	for i in range(m):
		y = T_list[i]/(T_list[i]+B)
		A[i,0] = (cpInf-cp0) * y*y*(y-1)
		A[i,1] = A[i,0] * y
		A[i,2] = A[i,1] * y
		A[i,3] = A[i,2] * y
		b[i] = Cp_list[i]-cp0 - y*y*(cpInf-cp0)
		
	#solve least squares problem A*x = b; http://docs.scipy.org/doc/scipy/reference/tutorial/linalg.html#solving-linear-least-squares-problems-and-pseudo-inverses
	x,resid,rank,sigma = linalg.lstsq(A,b)
	a0 = x[0]
	a1 = x[1]
	a2 = x[2]
	a3 = x[3]
	
	# this doesn't work yet:
	# err = rmsErrWilhoit(t, cp, cp0, cpInf, a0, a1, a2, a3) #(optional); display rmsError (dimensionless units)
	
	# scale everything back
	# T_list = [t*1000. for t in T_list] # not needed because not stored here
	# B = B*1000. # not needed because stored elsewhere
	# Cp_list = [x*R for x in Cp_list] # not needed becasue not stored
	
	# cp0 and cpInf should be in units of J/mol-K
	cp0 = cp0*R
	cpInf = cpInf*R
	
	# output comment; ****we could also include fitting accuracy ("err") in the output below
	comment = 'Fitted to GA data with Cp0=%2g and Cp_inf=%2g. '%(cp0,cpInf) + GAthermo.comment
	
	# first set I=J=0, then calculate what they should be 
	# by referring to H298, S298
	I = 0
	J = 0
	# create Wilhoit instance
	WilhoitThermo = ThermoWilhoitData( cp0, cpInf, a0, a1, a2, a3, I, J, comment=comment)
	# calculate correct I, J (integration constants for H, S, respectively)
	I = H298 - WilhoitThermo.getEnthalpy(298.15)
	J = S298 - WilhoitThermo.getEntropy(298.15)
	# update Wilhoit instance with correct I,J
	WilhoitThermo.I = I
	WilhoitThermo.J = J
	
	return WilhoitThermo

def CpLimits(atoms, rotors, linear):
	"""Calculate the zero and infinity limits for heat capacity.
	
	Input: number of atoms, number of rotors, linearity (`True` if molecule is linear)
	Output: Cp(0 K)/R, Cp(infinity)/R
	"""
	#(based off of lsfp_wilh1.f in GATPFit)

	if (atoms == 1): # monoatomic
		cp0 = 2.5
		cpInf = 2.5
	elif linear: # linear
		cp0	 = 3.5
		cpInf = 3*atoms - 1.5
	else: # nonlinear
		cp0	 = 4.0
		cpInf = 3*atoms - (2 + 0.5*rotors)
	return cp0, cpInf

def rmsErrWilhoit(self,t):
	#calculate the RMS error between the Wilhoit form and training data points; result will have same units as cp inputs; cp, cp0, and cpInf should agree in units (e.g. Cp/R); units of B and t should be consistent, based, for example on kK or K 
	m = len(t)
	rms = 0.0
	for i in range(m):
		err = cp[i]-getHeatCapacity(self,t[i])
		rms += err*err
	rms = rms/m
	rms = math.sqrt(rms)
	
	return rms

#this function already exists in getHeatCapacity
#def Wilhoit_Cp(t, cp0, cpInf, B, a0, a1, a2, a3):
#	#calculate Cp/R based on Wilhoit form;  result will have same units as cp0 and cpInf (e.g. Cp/R); cp is Cp/R; units of B and t should be consistent, based, for example on kK or K
#	y = t/(t+B)
#	cp = cp0+(cpInf-cp0)*y*y*(1+(y-1)*(a0+a1*y+a2*y*y+a3*y*y*y))
#	
#	return cp

################################################################################
def convertWilhoitToNASA(Wilhoit):
	"""Convert a Wilhoit thermo instance into a NASA polynomial thermo instance.
	
	Takes a `ThermoWilhoitData` instance of themochemical data.
	Returns a `ThermoNASAData` instance containing two `ThermoNASAPolynomial` 
	polynomials
	"""
	
	# Temperature ranges for resulting polynomials
	Tmin = 298.0
	Tintg = 1000.0
	Tmax = 6000.0
	
	#for now, do not allow tint to float so tint = Tint(guess)
	fixed = 1
	
	#for now, use weighting of 1/T
	weighting = 1
	
	# get info from incoming Wilhoit thermo
	cp0 = Wilhoit.cp0
	cpInf = Wilhoit.cpInf
	B = Wilhoit.B
	a0 = Wilhoit.a0
	a1 = Wilhoit.a1
	a2 = Wilhoit.a2
	a3 = Wilhoit.a3
	R = constants.R
	
	#scale the values to kK and dimensionless (for heat capacity)
	cp0 = cp0/R
	cpInf = cpInf/R
	B = B/1000
	Tmin = Tmin/1000
	Tintg = Tintg/1000
	Tmax = Tmax/1000
	
	# do your clever maths here
	#if we are using fixed tint, set tint equal to Tintg and do not allow tint to float
	if(fixed == 1):
		(b1, b2, b3, b4, b5, b6, b7, b8, b9, b10) = Wilhoit2NASA(cp0, cpInf, B, a0, a1, a2, a3, Tmin, Tmax, Tintg, weighting)
		err = TintOpt_objFun(Tintg, cp0, cpInf, B, a0, a1, a2, a3, Tmin, Tmax, weighting) #to print the objective function value
		tint = Tintg
	else:
		(b1, b2, b3, b4, b5, b6, b7, b8, b9, b10, tint) = Wilhoit2NASA_TintOpt(cp0, cpInf, B, a0, a1, a2, a3, Tmin, Tmax, Tintg, weighting)
		# rmsErr = rmsErrNASA(t, cp, b1, b2, b3, b4, b5, b6, b7, b8, b9, b10, tint) #this needs group data
		
	#restore to conventional units of K for Tint and units based on K rather than kK in NASA polynomial coefficients
	tint=tint*1000.
	Tmin = Tmin*1000
	Tmax = Tmax*1000
	b2 = b2/1000.
	b7 = b7/1000.
	b3 = b3/1000000.
	b8 = b8/1000000.
	b4 = b4/1000000000.
	b9 = b9/1000000000.
	b5 = b5/1000000000000.
	b10= b10/1000000000000.
	
	#for now, set H and S parameters equal to zero; this will need to be fixed****
	Hlow = 0.0
	Slow = 0.0
	Hhigh = 0.0
	Shigh = 0.0
	
	coeffs_low = (b1,b2,b3,b4,b5,Hlow,Slow)
	coeffs_high = (b6,b7,b8,b9,b10,Hhigh,Shigh)
	
	# could we include fitting accuracy in the expression below?
	# output comment
	comment = 'Fitted to Wilhoit data. '+Wilhoit.comment
		
	# create ThermoNASAPolynomial instances
	polynomial_low = ThermoNASAPolynomial( T_range=(Tmin,tint), comment=comment, coeffs=coeffs_low)
	polynomial_high = ThermoNASAPolynomial( T_range=(tint,Tmax), comment=comment, coeffs=coeffs_high)
	
	NASAthermo = ThermoNASAData( Trange=(Tmin,Tmax), polynomials=[polynomial_low,polynomial_high], comment=comment)
	return NASAthermo

################################################################################
def rmsErrNASA(t, cp, b1, b2, b3, b4, b5, b6, b7, b8, b9, b10, tint):
	#calculate the RMS error between the NASA polynomial and training data points in non-dimensional units (/R); cp is Cp/R; units of tint, t, and bi should be consistent, based, for example on kK or K 
	m = len(t)
	rms = 0.0
	for i in range(m):
		err = cp[i]-NASA_CpR(t[i],b1, b2, b3, b4, b5, b6, b7, b8, b9, b10, tint)
		rms += err*err
	rms = rms/m
	rms = math.sqrt(rms)
	
	return rms


def NASA_CpR(t, b1, b2, b3, b4, b5, b6, b7, b8, b9, b10, tint):
	#calculate Cp/R based on NASA polynomial; cp is Cp/R; units of tint, t, and bi should be consistent, based, for example on kK or K; does not take into account lower and upper temperature limits, and will extrapolate
	if(t < tint):
		cp = b1 + b2*t + b3*t*t + b4*t*t*t + b5*t*t*t*t
	else:
		cp = b6 + b7*t + b8*t*t + b9*t*t*t + b10*t*t*t*t
	<commit_msg>adding H,S consideration to Wilhoit to NASA conversion; in principle, it should work, but I still need to test it<commit_after>#!/usr/bin/python
# -*- coding: utf-8 -*-

################################################################################
#
#	RMG - Reaction Mechanism Generator
#
#	Copyright (c) 2002-2009 Prof. <NAME> (<EMAIL>) and the
#	RMG Team (<EMAIL>)
#
#	Permission is hereby granted, free of charge, to any person obtaining a
#	copy of this software and associated documentation files (the 'Software'),
#	to deal in the Software without restriction, including without limitation
#	the rights to use, copy, modify, merge, publish, distribute, sublicense,
#	and/or sell copies of the Software, and to permit persons to whom the
#	Software is furnished to do so, subject to the following conditions:
#
#	The above copyright notice and this permission notice shall be included in
#	all copies or substantial portions of the Software.
#
#	THE SOFTWARE IS PROVIDED 'AS IS', WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
#	IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
#	FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
#	AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
#	LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
#	FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
#	DEALINGS IN THE SOFTWARE.
#
################################################################################

"""
A module for working with the thermodynamics of chemical species. This module
seeks to provide functionality for answering the question, "Given a species,
what are its thermodynamics?"
"""

import quantities as pq
import constants
import data
import math
import scipy
from scipy import linalg
from scipy import optimize

################################################################################

class ThermoData:
	"""
	A base class for all forms of thermodynamic data used by RMG. The common
	attributes are:
	
	========= ============================================================
	Attribute Meaning
	========= ============================================================
	`Trange`  a 2-tuple containin the minimum and maximum temperature in K
	`comment` a string describing the source of the data
	========= ============================================================
	"""
	
	def __init__(self, Trange=None, comment=''):
		self.Trange = Trange or (0.0, 0.0)
		self.comment = comment
	
	def isTemperatureValid(self, T):
		"""
		Return :data:`True` if the temperature `T` in K is within the valid
		temperature range of the thermodynamic data, or :data:`False` if not.
		"""
		if self.Trange is None:
			return True
		else:
			return self.Trange[0] <= T and T <= self.Trange[1]

################################################################################

class ThermoGAData(ThermoData):
	"""
	A set of thermodynamic parameters as determined from Benson's group
	additivity data. The attributes are:
	
	========= ========================================================
	Attribute Meaning
	========= ========================================================
	`H298`    the standard enthalpy of formation at 298 K in J/mol
	`S298`    the standard entropy of formation at 298 K in J/mol*K
	`Cp`      the standard heat capacity in J/mol*K at 300, 400, 500,
		      600, 800, 1000, and 1500 K
	========= ========================================================
	"""
	# I think putting a comment with a colon just before the thing is defined
	# makes it show up in the documentation with autodoc. (or is that just in constants.py?)
	#: The list of temperatures at which heat capacity is defined = [300,400,500,600,800,1000,1500]
	CpTlist = list(pq.Quantity([300.0, 400.0, 500.0, 600.0, 800.0, 1000.0, 1500.0], 'K').simplified)
	CpTlist = [float(T) for T in CpTlist]
	# refer to it as ThermoGAData.CpTlist, even when calling it from methods within this Class.
	
	def __init__(self, H298=0.0, S298=0.0, Cp=None, comment='', index=''):
		"""Initialize a set of group additivity thermodynamic data."""
		ThermoData.__init__(self, Trange=(298.0, 2500.0), comment=comment)
		self.H298 = H298
		self.S298 = S298
		self.Cp = Cp or [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
		self.index = index
		self.__cache_T=None
		self.__cache = dict()
	
	def __add__(self, other):
		"""
		Add two sets of thermodynamic data together. All parameters are
		considered additive. Returns a new :class:`ThermoGAData` object that is
		the sum of the two sets of thermodynamic data.
		"""
		new = ThermoGAData()
		new.H298 = self.H298 + other.H298
		new.S298 = self.S298 + other.S298
		new.Cp = [self.Cp[i] + other.Cp[i] for i in range(len(self.Cp))]
		if self.comment == '': new.comment = other.comment
		elif other.comment == '': new.comment = self.comment
		else: new.comment = self.comment + '+ ' + other.comment
		new.index = self.index + '+' + other.index
		return new
	
	def __repr__(self):
		string = 'ThermoGAData(H298=%s, S298=%s, Cp=%s, index="%s")'%(self.H298,self.S298,self.Cp,self.index)
		return string
	
	def __str__(self):
		"""
		Return a string summarizing the thermodynamic data.
		"""
		string = ''
		string += 'Enthalpy of formation: %s J/mol\n' % (self.H298)
		string += 'Entropy of formation: %s J/mol*K\n' % (self.S298)
		string += 'Heat capacity (J/mol*K): '
		for T, Cp in zip(ThermoGAData.CpTlist, self.Cp):
			string += '%s(%sK) ' % (Cp,T)
		string += '\n'
		string += 'Index: %s\tComment: %s' % (self.index, self.comment)
		return string
	
	def getHeatCapacity(self, T):
		"""
		Return the heat capacity in J/mol*K at temperature `T` in K.
		"""
		if not self.isTemperatureValid(T):
			raise data.TemperatureOutOfRangeException('Invalid temperature for heat capacity estimation from group additivity.')
		if T < 300.0:
			return self.Cp[0]
		elif T > ThermoGAData.CpTlist[-1]:
			return self.Cp[-1]
		else:
			for Tmin, Tmax, Cpmin, Cpmax in zip(ThermoGAData.CpTlist[:-1], \
					ThermoGAData.CpTlist[1:], self.Cp[:-1], self.Cp[1:]):
				if Tmin <= T and T <= Tmax:
					return (Cpmax - Cpmin) * ((T - Tmin) / (Tmax - Tmin)) + Cpmin
	
	def getEnthalpy(self, T):
		"""
		Return the enthalpy in J/mol at temperature `T` in K.
		"""	
		if self.__cache_T==T:
			try: return self.__cache['H']
			except KeyError: pass # 'H' isn't cached yet
		else:  # temperature has changed
			self.__cache = dict()
			self.__cache_T = T
			
		H = self.H298
		if not self.isTemperatureValid(T):
			raise data.TemperatureOutOfRangeException('Invalid temperature for enthalpy estimation from group additivity.')
		for Tmin, Tmax, Cpmin, Cpmax in zip(ThermoGAData.CpTlist[:-1], \
				ThermoGAData.CpTlist[1:], self.Cp[:-1], self.Cp[1:]):
			if T > Tmin:
				slope = (Cpmax - Cpmin) / (Tmax - Tmin)
				intercept = (Cpmin * Tmax - Cpmax * Tmin) / (Tmax - Tmin)
				if T < Tmax:	H += 0.5 * slope * (T*T - Tmin*Tmin) + intercept * (T - Tmin)
				else:			H += 0.5 * slope * (Tmax*Tmax - Tmin*Tmin) + intercept * (Tmax - Tmin)
		if T > ThermoGAData.CpTlist[-1]:
			H += self.Cp[-1] * (T - ThermoGAData.CpTlist[-1])
		self.__cache['H'] = H
		return H
	
	def getEntropy(self, T):
		"""
		Return the entropy in J/mol*K at temperature `T` in K.
		"""
		import math
		
		if self.__cache_T==T:
			try: return self.__cache['S']
			except KeyError: pass # 'S' isn't cached yet
		else:  # temperature has changed
			self.__cache = dict()
			self.__cache_T = T
			
		S = self.S298
		if not self.isTemperatureValid(T):
			raise data.TemperatureOutOfRangeException('Invalid temperature for entropy estimation from group additivity.')
		for Tmin, Tmax, Cpmin, Cpmax in zip(ThermoGAData.CpTlist[:-1], \
				ThermoGAData.CpTlist[1:], self.Cp[:-1], self.Cp[1:]):
			if T > Tmin:
				slope = (Cpmax - Cpmin) / (Tmax - Tmin)
				intercept = (Cpmin * Tmax - Cpmax * Tmin) / (Tmax - Tmin)
				if T < Tmax:	S += slope * (T - Tmin) + intercept * math.log(T/Tmin)
				else:			S += slope * (Tmax - Tmin) + intercept * math.log(Tmax/Tmin)
		if T > ThermoGAData.CpTlist[-1]:
			S += self.Cp[-1] * math.log(T / ThermoGAData.CpTlist[-1])
		self.__cache['S'] = S
		return S
	
	def getFreeEnergy(self, T):
		"""
		Return the Gibbs free energy in J/mol at temperature `T` in K.
		"""
		if self.__cache_T==T:
			try: return self.__cache['G']
			except KeyError: pass # 'G' isn't cached yet
		else:  # temperature has changed
			self.__cache = dict()
			self.__cache_T = T
			
		if not self.isTemperatureValid(T):
			raise data.TemperatureOutOfRangeException('Invalid temperature for free energy estimation from group additivity.')
		G = self.getEnthalpy(T) - T * self.getEntropy(T)
		self.__cache['G'] = G
		return G
	
	def fromDatabase(self, data, comment):
		"""
		Process a list of numbers `data` and associated description `comment`
		generated while reading from a thermodynamic database.
		"""
		
		if len(data) != 12:
			raise Exception('Invalid list of thermo data; should be a list of numbers of length 12.')
		
		H298, S298, Cp300, Cp400, Cp500, Cp600, Cp800, Cp1000, Cp1500, \
			dH, dS, dCp = data
		
		self.H298 = float(pq.Quantity(H298, 'kcal/mol').simplified)
		self.S298 = float(pq.Quantity(S298, 'cal/(mol*K)').simplified)
		self.Cp = list(pq.Quantity([Cp300, Cp400, Cp500, Cp600, Cp800, Cp1000, Cp1500], 'cal/(mol*K)').simplified)
		for i in range(len(self.Cp)): self.Cp[i] = float(self.Cp[i])
		self.comment = comment
	
	def toXML(self, dom, root):
		"""
		Generate an XML representation of the thermodynamic data using the
		:data:`xml.dom.minidom` package. The `dom` and `root` parameters
		represent the DOM and the element in the DOM used as the parent of
		the generated XML.
		"""
		
		thermo = dom.createElement('thermodynamics')
		thermo.setAttribute('index', self.index)
		thermo.setAttribute('comment', self.comment)
		root.appendChild(thermo)
		
		enthalpy = dom.createElement('enthalpyOfFormation')
		enthalpy.setAttribute('temperature', '298.0 K')
		thermo.appendChild(enthalpy)
		data.createXMLQuantity(dom, enthalpy, self.H298, 'J/mol')
		
		entropy = dom.createElement('entropyOfFormation')
		entropy.setAttribute('temperature', '298.0 K')
		thermo.appendChild(entropy)
		data.createXMLQuantity(dom, entropy, self.S298, 'J/(mol*K)')
		
		for i, Cp in enumerate(self.Cp):
			heatCapacity = dom.createElement('heatCapacity')
			heatCapacity.setAttribute('temperature', '%s K' % (ThermoGAData.CpTlist[i]) )
			thermo.appendChild(heatCapacity)
			data.createXMLQuantity(dom, heatCapacity, Cp, 'J/(mol*K)')
	

################################################################################

class ThermoNASAPolynomial(ThermoData):
	"""
	A single NASA polynomial for thermodynamic data. The `coeffs` attribute
	stores the seven polynomial coefficients
	:math:`\\mathbf{a} = \\left[a_1\\ a_2\\ a_3\\ a_4\\ a_5\\ a_6\\ a_7 \\right]`
	from which the relevant thermodynamic parameters are evaluated via the
	expressions
	
	.. math:: \\frac{C_\\mathrm{p}(T)}{R} = a_1 + a_2 T + a_3 T^2 + a_4 T^3 + a_5 T^4
	
	.. math:: \\frac{H(T)}{RT} = a_1 + \\frac{1}{2} a_2 T + \\frac{1}{3} a_3 T^2 + \\frac{1}{4} a_4 T^3 + \\frac{1}{5} a_5 T^4 + \\frac{a_6}{T}
	
	.. math:: \\frac{S(T)}{R} = a_1 \\ln T + a_2 T + \\frac{1}{2} a_3 T^2 + \\frac{1}{3} a_4 T^3 + \\frac{1}{4} a_5 T^4 + a_7
	
	The above was adapted from `this page <http://www.me.berkeley.edu/gri-mech/data/nasa_plnm.html>`_.
	"""
	
	def __init__(self, T_range=None, coeffs=None, comment=''):
		ThermoData.__init__(self, Trange=T_range, comment=comment)
		self.coeffs = coeffs or (0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0)
	def __repr__(self):
                return "ThermoNASAPolynomial(%r,%r,'%s')"%(self.Trange,self.coeffs,self.comment)
	def getHeatCapacity(self, T):
		"""
		Return the heat capacity in J/mol*K at temperature `T` in K.
		"""
		import constants
		if not self.isTemperatureValid(T):
			raise data.TemperatureOutOfRangeException('Invalid temperature for heat capacity estimation from NASA polynomial.')
		c = self.coeffs
		T2 = T*T
		# Cp/R = a1 + a2 T + a3 T^2 + a4 T^3 + a5 T^4
#		HeatCapacityOverR = c[0] + c[1]*T + c[2]*T*T + c[3]*T*T*T + c[4]*T*T*T*T
		HeatCapacityOverR = c[0] + T*(c[1] + T*(c[2] + T*(c[3] + c[4]*T)))
		Cp = HeatCapacityOverR * constants.R
		return Cp
	
	def getEnthalpy(self, T):
		"""
		Return the enthalpy in J/mol at temperature `T` in K.
		"""
		import constants
		if not self.isTemperatureValid(T):
			raise data.TemperatureOutOfRangeException('Invalid temperature for enthalpy estimation from NASA polynomial.')
		c = self.coeffs
		T2 = T*T
		T4 = T2*T2
		# H/RT = a1 + a2 T /2 + a3 T^2 /3 + a4 T^3 /4 + a5 T^4 /5 + a6/T
		EnthalpyOverR = ( c[0]*T + c[1]*T2/2 + c[2]*T2*T/3 + c[3]*T4/4 +
						  c[4]*T4*T/5 + c[5] )
		H = EnthalpyOverR * constants.R
		return H
	
	def getEntropy(self, T):
		"""
		Return the entropy in J/mol*K at temperature `T` in K.
		"""
		import math
		import constants
		if not self.isTemperatureValid(T):
			raise data.TemperatureOutOfRangeException('Invalid temperature for entropy estimation from NASA polynomial.')
		# S/R  = a1 lnT + a2 T + a3 T^2 /2 + a4 T^3 /3 + a5 T^4 /4 + a7
		c = self.coeffs
		T2 = T*T
		T4 = T2*T2
		EntropyOverR = ( c[0]*math.log(T) + c[1]*T + c[2]*T2/2 +
					c[3]*T2*T/3 + c[4]*T4/4 + c[6] )
		S = EntropyOverR * constants.R
		return S
	
	def getFreeEnergy(self, T):
		"""
		Return the Gibbs free energy in J/mol at temperature `T` in K.
		"""
		if not self.isTemperatureValid(T):
			raise data.TemperatureOutOfRangeException('Invalid temperature for free energy estimation from NASA polynomial.')
		return self.getEnthalpy(T) - T * self.getEntropy(T)
	
	def toXML(self, dom, root):
### prime-like:
#   <polynomial>
#      <validRange>
#         <bound kind="lower" property="temperature" units="K">300.000</bound>
#         <bound kind="upper" property="temperature" units="K">1000</bound>
#      </validRange>
#      <coefficient id="1" label="a1">4.0733</coefficient>
#      <coefficient id="2" label="a2">0.011308</coefficient>
#      <coefficient id="3" label="a3">-1.6565e-005</coefficient>
#      <coefficient id="4" label="a4">1.1784e-008</coefficient>
#      <coefficient id="5" label="a5">-3.3006e-012</coefficient>
#      <coefficient id="6" label="a6">-19054.23</coefficient>
#      <coefficient id="7" label="a7">4.4083</coefficient>
#   </polynomial>
		pass
	

################################################################################

class ThermoNASAData(ThermoData):
	"""
	A set of thermodynamic parameters given by NASA polynomials. This class
	stores a list of :class:`ThermoNASAPolynomial` objects in the `polynomials`
	attribute. When evaluating a thermodynamic quantity, a polynomial that
	contains the desired temperature within its valid range will be used.
	"""
	
	def __init__(self, polynomials=None, comment='', Trange=None):
		ThermoData.__init__(self, Trange=Trange, comment=comment)
		self.polynomials = polynomials or []
	
	def addPolynomial(self, polynomial):
		if not isinstance(polynomial,ThermoNASAPolynomial):
			raise TypeError("Polynomial attribute should be instance of ThermoNASAPolynomial")
		self.polynomials.append(polynomial)
	
	def selectPolynomialForTemperature(self, T):
		for poly in self.polynomials:
			if poly.isTemperatureValid(T): break
		else:
			raise data.TemperatureOutOfRangeException("No polynomial found for T=%s" % T)
		return poly
	
	def __repr__(self):
		return "ThermoNASAData(%s, '%s')"%(repr(self.polynomials),self.comment)
	
	def toXML(self, dom, root):
### prime-like: 
# <thermodynamicPolynomials type="nasa7">
#   <referenceState>
#      <Tref units="K">298.15</Tref>
#      <Pref units="Pa">100000</Pref>
#   </referenceState>
#   <dfH units="J/mol">-145177.5731</dfH>
#   <polynomial> FROM NASA_polynomials </polynomial>
# </thermodynamicPolynomials>
		pass
	
	def getHeatCapacity(self, T):
		"""
		Return the heat capacity in J/mol*K at temperature `T` in K.
		"""
		poly = self.selectPolynomialForTemperature(T)
		return poly.getHeatCapacity(T)
	
	def getEnthalpy(self, T):
		"""
		Return the enthalpy in J/mol at temperature `T` in K.
		"""
		poly = self.selectPolynomialForTemperature(T)
		return poly.getEnthalpy(T)
	
	def getEntropy(self, T):
		"""
		Return the entropy in J/mol*K at temperature `T` in K.
		"""
		poly = self.selectPolynomialForTemperature(T)
		return poly.getEntropy(T)
	
	def getFreeEnergy(self, T):
		"""
		Return the Gibbs free energy in J/mol at temperature `T` in K.
		"""
		poly = self.selectPolynomialForTemperature(T)
		return poly.getFreeEnergy(T)
	

################################################################################

class ThermoWilhoitData(ThermoData):
	"""
	A set of thermodynamic parameters given by Wilhoit polynomials. 
	
	"""
	B = 500.0 # Kelvin. Default temperature.
	
	def __init__(self, cp0, cpInf, a0, a1, a2, a3, I, J, comment='', ):
		"""Initialise the Wilhoit polynomial. Trange is set to (0,9999.9)"""
		Trange = (0,9999.9) # Wilhoit valid over all temperatures
		ThermoData.__init__(self, Trange=Trange, comment=comment)
		self.cp0 = cp0
		self.cpInf = cpInf
		self.B = ThermoWilhoitData.B
		self.a0 = a0
		self.a1 = a1
		self.a2 = a2
		self.a3 = a3
		self.I = I
		self.J = J
	
	def __repr__(self):
		return "ThermoWilhoitData(%.2g,%.2g,%.2g,%.2g,%.2g,%.2g,%.2g,%.2g,%.2g'%s')"%(self.cp0, self.cpInf, self.B, self.a0, self.a1, self.a2, self.a3, self.I, self.J, self.comment)
	
	def toXML(self, dom, root):
		pass
	
	def getHeatCapacity(self, T):
		"""
		Return the heat capacity in J/mol*K at temperature `T` in K.
		"""
		y = T/(T+self.B)
		Cp = self.cp0+(self.cpInf-self.cp0)*y*y*( 1 + 
			(y-1)*(self.a0 + y*(self.a1 + y*(self.a2 + y*self.a3))) )
		return Cp
	
	def getEnthalpy(self, T):
		"""
		Return the enthalpy in J/mol at temperature `T` in K.
		"""
		#cp0 = self.cp0 * R
		#cpInf = self.cpInf * R
		cp0 = self.cp0
		cpInf = self.cpInf
		B = self.B
		a0 = self.a0
		a1 = self.a1
		a2 = self.a2
		a3 = self.a3
		I = self.I
		
		enthalpy = I + WilhoitInt0(cp0, cpInf, B, a0, a1, a2, a3, T)
		
		return enthalpy
	
	def getEntropy(self, T):
		"""
		Return the entropy in J/mol*K at temperature `T` in K.
		"""
		#cp0 = self.cp0 * R
		#cpInf = self.cpInf * R
		cp0 = self.cp0
		cpInf = self.cpInf
		B = self.B
		a0 = self.a0
		a1 = self.a1
		a2 = self.a2
		a3 = self.a3
		J = self.J
		
		entropy = J + WilhoitIntM1(cp0, cpInf, B, a0, a1, a2, a3, T)
		
		return entropy
	
	def getFreeEnergy(self, T):
		"""
		Return the Gibbs free energy in J/mol at temperature `T` in K.
		"""
		return self.getEnthalpy(T) - T*self.getEntropy(T)
	

###########
def convertGAtoWilhoit(GAthermo, atoms, rotors, linear):
	"""Convert a Group Additivity thermo instance into a Wilhoit thermo instance.
	
	Takes a `ThermoGAData` instance of themochemical data, and some extra information 
	about the molecule used to calculate high- and low-temperature limits of Cp.
	These are the number of atoms (integer `atoms`), the number of rotors (integer `rotors`)
	and whether the molecule is linear (boolean `linear`)
	Returns a `ThermoWilhoitData` instance.
	
	cf. <NAME>'s thesis, p. 185-186
	"""
	
	# get info from incoming group additivity thermo
	H298 = GAthermo.H298
	S298 = GAthermo.S298
	Cp_list = GAthermo.Cp
	T_list = ThermoGAData.CpTlist  # usually [300, 400, 500, 600, 800, 1000, 1500] but why assume?
	R = constants.R
	B = ThermoWilhoitData.B # Constant, set once in the class def.
	
	# convert from K to kK
	T_list = [t/1000. for t in T_list] 
	B = B/1000.  
	
	Cp_list = [x/R for x in Cp_list] # convert to Cp/R
	
	(cp0, cpInf) = CpLimits(atoms, rotors, linear) # determine the heat capacity limits (non-dimensional)
	
	#create matrices for linear least squares problem
	m = len(T_list)  # probably m=7
	# A = mx4
	# b = mx1
	# x = 4x1
	A = scipy.zeros([m,4])
	b = scipy.zeros([m,1])
	for i in range(m):
		y = T_list[i]/(T_list[i]+B)
		A[i,0] = (cpInf-cp0) * y*y*(y-1)
		A[i,1] = A[i,0] * y
		A[i,2] = A[i,1] * y
		A[i,3] = A[i,2] * y
		b[i] = Cp_list[i]-cp0 - y*y*(cpInf-cp0)
		
	#solve least squares problem A*x = b; http://docs.scipy.org/doc/scipy/reference/tutorial/linalg.html#solving-linear-least-squares-problems-and-pseudo-inverses
	x,resid,rank,sigma = linalg.lstsq(A,b)
	a0 = x[0]
	a1 = x[1]
	a2 = x[2]
	a3 = x[3]
	
	# this doesn't work yet:
	# err = rmsErrWilhoit(t, cp, cp0, cpInf, a0, a1, a2, a3) #(optional); display rmsError (dimensionless units)
	
	# scale everything back
	# T_list = [t*1000. for t in T_list] # not needed because not stored here
	# B = B*1000. # not needed because stored elsewhere
	# Cp_list = [x*R for x in Cp_list] # not needed because not stored
	
	# cp0 and cpInf should be in units of J/mol-K
	cp0 = cp0*R
	cpInf = cpInf*R
	
	# output comment; ****we could also include fitting accuracy ("err") in the output below
	comment = 'Fitted to GA data with Cp0=%2g and Cp_inf=%2g. '%(cp0,cpInf) + GAthermo.comment
	
	# first set I=J=0, then calculate what they should be 
	# by referring to H298, S298
	I = 0
	J = 0
	# create Wilhoit instance
	WilhoitThermo = ThermoWilhoitData( cp0, cpInf, a0, a1, a2, a3, I, J, comment=comment)
	# calculate correct I, J (integration constants for H, S, respectively)
	I = H298 - WilhoitThermo.getEnthalpy(298.15)
	J = S298 - WilhoitThermo.getEntropy(298.15)
	# update Wilhoit instance with correct I,J
	WilhoitThermo.I = I
	WilhoitThermo.J = J
	
	return WilhoitThermo

def CpLimits(atoms, rotors, linear):
	"""Calculate the zero and infinity limits for heat capacity.
	
	Input: number of atoms, number of rotors, linearity (`True` if molecule is linear)
	Output: Cp(0 K)/R, Cp(infinity)/R
	"""
	#(based off of lsfp_wilh1.f in GATPFit)

	if (atoms == 1): # monoatomic
		cp0 = 2.5
		cpInf = 2.5
	elif linear: # linear
		cp0	 = 3.5
		cpInf = 3*atoms - 1.5
	else: # nonlinear
		cp0	 = 4.0
		cpInf = 3*atoms - (2 + 0.5*rotors)
	return cp0, cpInf

def rmsErrWilhoit(self,t):
	#calculate the RMS error between the Wilhoit form and training data points; result will have same units as cp inputs; cp, cp0, and cpInf should agree in units (e.g. Cp/R); units of B and t should be consistent, based, for example on kK or K 
	m = len(t)
	rms = 0.0
	for i in range(m):
		err = cp[i]-getHeatCapacity(self,t[i])
		rms += err*err
	rms = rms/m
	rms = math.sqrt(rms)
	
	return rms

#this function already exists in getHeatCapacity
#def Wilhoit_Cp(t, cp0, cpInf, B, a0, a1, a2, a3):
#	#calculate Cp/R based on Wilhoit form;  result will have same units as cp0 and cpInf (e.g. Cp/R); cp is Cp/R; units of B and t should be consistent, based, for example on kK or K
#	y = t/(t+B)
#	cp = cp0+(cpInf-cp0)*y*y*(1+(y-1)*(a0+a1*y+a2*y*y+a3*y*y*y))
#	
#	return cp

################################################################################
def convertWilhoitToNASA(Wilhoit):
	"""Convert a Wilhoit thermo instance into a NASA polynomial thermo instance.
	
	Takes a `ThermoWilhoitData` instance of themochemical data.
	Returns a `ThermoNASAData` instance containing two `ThermoNASAPolynomial` 
	polynomials
	"""
	
	# Temperature ranges for resulting polynomials
	Tmin = 298.0
	Tintg = 1000.0
	Tmax = 6000.0
	
	#for now, do not allow tint to float so tint = Tint(guess)
	fixed = 1
	
	#for now, use weighting of 1/T
	weighting = 1
	
	# get info from incoming Wilhoit thermo
	cp0 = Wilhoit.cp0
	cpInf = Wilhoit.cpInf
	B = Wilhoit.B
	a0 = Wilhoit.a0
	a1 = Wilhoit.a1
	a2 = Wilhoit.a2
	a3 = Wilhoit.a3
	R = constants.R
	
	#scale the values to kK and dimensionless (for heat capacity)
	cp0 = cp0/R
	cpInf = cpInf/R
	B = B/1000
	Tmin = Tmin/1000
	Tintg = Tintg/1000
	Tmax = Tmax/1000
	
	# do your clever maths here
	#if we are using fixed tint, set tint equal to Tintg and do not allow tint to float
	if(fixed == 1):
		(b1, b2, b3, b4, b5, b6, b7, b8, b9, b10) = Wilhoit2NASA(cp0, cpInf, B, a0, a1, a2, a3, Tmin, Tmax, Tintg, weighting)
		err = TintOpt_objFun(Tintg, cp0, cpInf, B, a0, a1, a2, a3, Tmin, Tmax, weighting) #to print the objective function value
		tint = Tintg
	else:
		(b1, b2, b3, b4, b5, b6, b7, b8, b9, b10, tint) = Wilhoit2NASA_TintOpt(cp0, cpInf, B, a0, a1, a2, a3, Tmin, Tmax, Tintg, weighting)
		# rmsErr = rmsErrNASA(t, cp, b1, b2, b3, b4, b5, b6, b7, b8, b9, b10, tint) #this needs group data
		
	#restore to conventional units of K for Tint and units based on K rather than kK in NASA polynomial coefficients
	tint=tint*1000.
	Tmin = Tmin*1000
	Tmax = Tmax*1000
	b2 = b2/1000.
	b7 = b7/1000.
	b3 = b3/1000000.
	b8 = b8/1000000.
	b4 = b4/1000000000.
	b9 = b9/1000000000.
	b5 = b5/1000000000000.
	b10= b10/1000000000000.
	
	#set H and S parameters equal to zero; we will solve for the correct values below, based on the enthalpy/entropy with the parameters equal to zero 
	Hlow = 0.0
	Slow = 0.0
	Hhigh = 0.0
	Shigh = 0.0
	
	coeffs_low = (b1,b2,b3,b4,b5,Hlow,Slow)
	coeffs_high = (b6,b7,b8,b9,b10,Hhigh,Shigh)
	
	# could we include fitting accuracy in the expression below?
	# output comment
	comment = 'Fitted to Wilhoit data. '+Wilhoit.comment
		
	# create ThermoNASAPolynomial instances
	polynomial_low = ThermoNASAPolynomial( T_range=(Tmin,tint), comment=comment, coeffs=coeffs_low)
	polynomial_high = ThermoNASAPolynomial( T_range=(tint,Tmax), comment=comment, coeffs=coeffs_high)

        #for the low polynomial, we want the results to match the Wilhoit value at 298.15K
        #low polynomial enthalpy:
	polynomial_low.coeffs[5] = (Wilhoit.getEnthalpy(298.15) - polynomial_low.getEnthalpy(298.15))/constants.R
	#low polynomial entropy:
	polynomial_low.coeffs[6] = (Wilhoit.getEntropy(298.15) - polynomial_low.getEntropy(298.15))/constants.R
        #for the high polynomial, we want the results to match the low polynomial value at tint
	#high polynomial enthalpy:
	polynomial_high.coeffs[5] = (polynomial_low.getEnthalpy(tint) - polynomial_high.getEnthalpy(tint))/constants.R
	#high polynomial entropy:
	polynomial_high.coeffs[6] = (polynomial_low.getEntropy(tint) - polynomial_high.getEntropy(tint))/constants.R
	
	NASAthermo = ThermoNASAData( Trange=(Tmin,Tmax), polynomials=[polynomial_low,polynomial_high], comment=comment)
	return NASAthermo

################################################################################
def rmsErrNASA(t, cp, b1, b2, b3, b4, b5, b6, b7, b8, b9, b10, tint):
	#calculate the RMS error between the NASA polynomial and training data points in non-dimensional units (/R); cp is Cp/R; units of tint, t, and bi should be consistent, based, for example on kK or K 
	m = len(t)
	rms = 0.0
	for i in range(m):
		err = cp[i]-NASA_CpR(t[i],b1, b2, b3, b4, b5, b6, b7, b8, b9, b10, tint)
		rms += err*err
	rms = rms/m
	rms = math.sqrt(rms)
	
	return rms


def NASA_CpR(t, b1, b2, b3, b4, b5, b6, b7, b8, b9, b10, tint):
	#calculate Cp/R based on NASA polynomial; cp is Cp/R; units of tint, t, and bi should be consistent, based, for example on kK or K; does not take into account lower and upper temperature limits, and will extrapolate
	if(t < tint):
		cp = b1 + b2*t + b3*t*t + b4*t*t*t + b5*t*t*t*t
	else:
		cp = b6 + b7*t + b8*t*t + b9*t*t*t + b10*t*t*t*t
	
{'fraction_non_alphanumeric': 0.17072065266656594, 'fraction_numerical': 0.0401306843934129, 'mean_word_length': 3.8209213401311, 'pattern_counts': {'":': 0, '<': 97, '<?xml version=': 0, '>': 91, 'https://': 0, 'lorem ipsum': 0, 'www.': 2, 'xml': 2}, 'pii_count': 0, 'substrings_counts': 8, 'word_list_counts': {'cursed_substrings.json': 44, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '1960944', 'n_tokens_mistral': 22214, 'n_tokens_neox': 19376, 'n_words': 5758}