Applicability of DFT model in reactive distillation
-
Maciej Staszak
Abstract
The density functional theory (DFT) applicability to reactive distillation is discussed. Brief modeling techniques description of distillation and rectification with chemical reaction is provided as a background for quantum method usage description. The equilibrium and nonequilibrium distillation models are described for that purpose. The DFT quantum theory is concisely described. The usage of DFT in the modeling of reactive distillation is described in two parts. One of the fundamental and very important component of distillation modeling is vapor-liquid equilibrium description for which the DFT quantum approach can be used. The representative DFT models, namely COSMO-RS (Conductor like Screening Model for Real Solvents), COSMOSPACE (COSMO Surface Pair Activity Coefficient) and COSMO-SAC (SAC – segment activity coefficient) approaches are described. The second part treats the way in which the chemical reaction is described by means of quantum DFT method. The intrinsic reaction coordinate (IRC) method is described which is used to find minimum energy path of substrates to products transition. The DFT is one of the methods which can be used for that purpose. The literature data examples are provided which proves that IRC method is applicable for chemical reaction kinetics description.
1 Introduction
This article is divided into three main parts, the first introducing to the traditional concepts of reactive distillation modeling and clarifying the prospective need for quantum calculations. The second provides more details about density functional theory usage in vapor–liquid equilibria. This element of distillation modeling description is very important for accurate simulation due to the fundamental impact of vapor–liquid equilibria on the components distribution among the phases during the process. The third part presents quantum approach to chemical reaction kinetics description. For reactive distillation, it is another very important aspect which also has large effect on the accuracy and the realism of the results obtained.
There exist many highly specialized software tools available for engineers which enable fast and reliable design calculation. ChemCad, Hysys, UniSim are a few examples of them. Up to now no quantum method is adapted or included in commercial packages for distillation or reactive distillation design. But it is only a matter of time when such methods become wider available. Such a progression is visible in the field of vapor–liquid equilibria, where the COSMO methods are used to simulate distillation. Up to now only in the research field such quantum approach can be met, for cases with complex equilibria like three-phase distillation [1], or for the cases when the solvent screening is very difficult by traditional methods [2, 3].
One of the key advantages that a designer might find useful when performing calculations based on quantum approaches is that there should not exist a requirement to perform experiments prior to calculations. Performing ab initio (from the beginning) quantum calculations that prepares model for design phase, enables vapor–liquid equilibria and chemical reaction kinetics to be accounted for without additional experimental measurements. The disadvantage might be the high quantum calculation workload which for ab initio methods (like, e. g. MP2) is unavoidable. The solution proposed to that problem is the idea of databank creation where all fundamental molecular quantities (charge density surfaces, activation energies, etc.) are already calculated. Developers of one of the leading quantum software tools existing on the market, COSMOlogic, made the efforts to create such database available. Constantly the work of researchers is devoted to develop such databanks for wider use.
2 Distillation process
Distillation is well-known industrial process for separating different kind of liquid mixtures. This introductory article will give a concise description of rectification modeling which is the base for further text concerning molecular description of thermodynamic equilibrium.
The different tendency of the mixture components to gather between liquid and vapor phase is the key basis of this physical process. Industrial distillation is often referred as rectification, which is typically conducted in column apparatus where the inner, working elements providing phase mixing are trays or packing. Illustrative sketch of rectifying column is presented in Figure 1.
Distillation column schematics with designation of flows used in the text.
Example of cubic model lattice with coordination number z = 6.
Intrinsic reaction coordinate presenting energy states of substrate, transition state and reaction product.
This is also typical that distillation is used to intensify chemical reaction especially in the case of reversible, equilibrium one. In this case the equilibrium shift toward higher conversion is realized by removing the reaction products from the reaction mixture. The most typical design and simulation approach used in many modern software packages is to divide the column height into sections. Depending on its construction the section may refer to trays or package segment. Every section is then balanced based on the approach used. Such a description is referred to as “tray-by-tray” model and consists of several equations in one set. Besides the steady-state or dynamics approach that can be used in such models are divided by the way the interphase mass transfer is treated. The most popular is the so-called equilibrium approach in which the phases are assumed to be perfectly mixed at every column segment and being in equilibrium condition. The second so-called nonequilibrium and more complex approach is to refer for mass transfer at every stage assuming that stage is not at equilibrium conditions.
The model is built for N segments, which can be numbered from top to bottom, but this is not compulsory. At the top, the condenser constitutes unique equation and at the bottom the reboiler constitutes also specific equation, both of them describing the mass and energy balances around them. It depends on the level of complexity and details demanded of how to treat these auxiliary units. For this text, they will be treated as simple as possible due to the subject of the article.
In the most extended case when taking into consideration also the flow of chemical reaction, the model must be capable to account for chemical reaction kinetics. This is also typical that the mass of vapor phase is much smaller than liquid and is neglected during calculations.
2.1 Equilibrium approach
The equilibrium approach is sometimes referred to as efficiency modeling due to the fact that the actual mixing process in real system does not produce exactly equilibrium conditions. Consequently, the equilibrium approach, being more optimistic in the way that the smaller number of segments would be enough for assumed design requirements, must be corrected by the use of efficiencies calculated for every segment.
The transient total mass balance for given ith segment reads (the segments are assumed to be numbered from top to bottom of the column):
In the above formulation, the L is the liquid phase flowrate, V is the vapor phase flowrate and F is the feed stream flowrate. The flowrate may use molar or mass base of units’ description. The m is the amount of mass (using moles or kilograms) contained at ith stage. When considering chemical reaction the molar units system is most suitable due to the natural way of chemical reaction kinetics definition:
The ΔRi,j is the source of ith specie due to chemical reaction. The variable ri,j is the rate of component jth creation/consumption due to the chemical reaction at ith stage, while ri is the density of the mixture at ith stage. In the case of nonreactive flow the ΔRi,j is equal zero. The Δν is the molar change due to the reaction stoichiometry.
The total mass balance is used in the formulation for jth component mass balance at ith stage given below:
The variable xi,j refers to molar (mass) fractions of ith component at jth stage. The subscript F refers to feed composition. In the description above the feed is assumed to be liquid but it is not a process limitation but only a simplification applied here. This is typical that column contains one or a few only feed streams or in the case of batch distillation it may contain no feeds at all.
The energy balance is the important part of the balance equations which allows to estimate the temperatures along the column. The energy balance for ith segment can be presented by the relation:
The variable H refers to stream enthalpy for given phase at given stage. The ΔHr is the reaction enthalpy and the Δℜi reads:
where ri is the overall rate of chemical reaction at ith stage.
The above equations form the differential equations set describing whole column at its every stage. The key element of the distillation process description is the vapor–liquid equilibrium (VLE). The usual way of estimating the equilibrium is to use some of the well-known and established methods for equilibrium conditions estimations. In general the equilibrium relation is presented in the form:
This formulation relates the concentration of jth component in the vapor to its concentration in the liquid on every stage. The equilibrium constant, being in fact a parameter, Kj is the subject of calculation by several methods, which will be discussed in detail in Section 2.3.
The tray efficiency or overall column efficiency must be applied to account for realism of the distillation process. The simplest definition is the overall efficiency Eo which is defined according to the equation below:
In the above Neq is the number of trays used by equilibrium approach, Nact is the actual tray number. Such overall efficiency can be applied to some specified column section or to column as a whole. The exact calculation of such efficiency is a complex task and for general cases some estimations are only be proposed [4, 5] but for selected cases the procedures are designed, e. g. for alkane mixtures fractionation [6].
The most typical formulation for tray efficiency is the Murphree efficiency (EM) which is given by:
In the above formulation, the
The model presented is sufficient for calculating basic mass and energy balance which gives as a result the components concentrations and temperature profiles along the column. To be able to describe the hydraulics of the process, additional equations and procedures must be added. By applying the column geometry by defining the tray sizing, the total volume for selected stage can be calculated. Based on mass balance solution from the above equations the amount of liquid at every stage is calculated mi. In the case when the plates are constructed with weirs the Francis formula for liquid flow can be used. The excess amount of liquid at given plate can be calculated from the geometric plate volume and liquid mass mi. This gives the possibility to estimate the liquid flow over the weir to the tray below. For the vapor phase flowrates the different approach can be used. The mechanistic calculation for pressure evolution along the column can be calculated by applying the pressure drop calculation for given trays or packing. The reboiler pressure increase (due to heating and liquid vaporization) is then the driving parameter to estimate the vapor flows along the column. Calculation of the actual pressure drop between reboiler and a tray above and then tray-to-tray pressure drops along the whole column enables to perform the estimation of the actual vapor flowrates. The additional data which must be given prior to this algorithm of calculation is top pressure. The detailed correlations for pressure drop for different flow regimes and different types of trays and weirs or packings are given by Kister [11].
2.2 Nonequilibrium approach
The nonequilibrium approach is a modification to the equilibrium description in the way that additional mass fluxes are considered, namely component interphase fluxes. The fundamental elements of nonequilbrium distillation approach are:
material balances,
energy balances,
equilibrium relations,
mass and energy transfer models
The phases in this approach are balanced by independent equations. The mass balance model formulation for jth component for liquid phase then reads:
and for the vapor phase:
The mass transfer streams for ith component for both phases
The total mass balance is defined for both phases and for liquid is given by:
Consequently, for vapor phase it reads:
The
The energy balance for the liquid phase can be written as follows:
and for the vapor phase consequently:
In the above equations, the energy streams
The equilibrium is assumed only to exist at the interface between phases. The value of
2.3 Vapor–liquid equilibrium
The equilibrium and mass transfer approaches both require a method of VLE estimation. The choice of specific method is based mainly on the type of distilled mixture components and the secondary on the actual pressure-temperature range.
The traditional methodology is to use one of the following thermodynamic approaches:
activity methods
equation of state methods
special methods
quantum and molecular methods
A very short description of them is presented below (as they are not based on any quantum or molecular mechanic method) and the latter will be discussed in details later in the text.
2.3.1 Activity methods
This approach is developed from the principia relating excess functions with activity coefficients. For the so-called regular liquid solutions (characterized by nonzero heat of mixing while entropy obeys ideal mixing rule) the activity coefficient of ith component is related to the excess Gibbs free energy
The typical models that are based on that definition are: Margules one and two parameter models [19, 20], van Laar [21], Wilson [22], NRTL (NonRandom TwoLiquids) [23], UNIFAC [24], UNIQUAC [25]. The models mentioned are the subject of extensions which provide additional useful properties like three-phase predictions (modified Wilson [26], UNIFAC LLE [27]) or polymer solutions property estimation (UNIFAC for polymers [28]).
2.3.2 Equation of state methods
Describing thermodynamic state functions by single relation is the focus of research which began from the introduction Clapeyron equation of state. Relating the amount of matter along with space occupied for given pressure and temperature is the main objective of this methodology. The inaccuracy of Clapeyron equation led to further study which resulted in Van der Waals introduction of attraction and volumetric term into the equation of state. Further modifications allowed more accurate prediction of phase behavior by cubic equations by Redlich and Kwong [29], which was then later corrected by Soave [30]. Modifications, which aim into more accurate description, were further introduced (e. g. Peng-Robinson [31] with its modifications). Further development leads to quartic equations of state which are reported to be of increased accuracy [32].
At the same time, the statistical associating fluid theory gave rise to so-called SAFT equations of state family. The SAFT equations of state are based on the equilibrium description by residual Helmholtz energy which originates from:
model of hard spheres effect
dispersion effects
molecular chains effects
association effects
The resulting family of equations are typically presented as the summation of all contributing Helmholtz energies:
Each of the residual part is calculated by models discussed elsewhere [33, 34, 35, 36, 37].
3 Short outline of quantum modeling by DFT
3.1 Density functional theory
In the chemical calculation field, the highest accuracy for the molecular system can be achieved by solving the Schrödinger equation by specified approximation scheme without the need to rely on direct experimental measurements. The wave function of a molecular system can be found by solving the quantum equations directly using so-called ab initio methods. On the other hand, the density functional theory (DFT) can be used in cases where knowledge about wave function is not of highest importance. The DFT concept is utilized when the electron charge density is sufficient for describing the required molecular properties. The dominant role of DFT methods in the field of quantum calculations is due to their computational efficiency and high accuracy. The DFT theory was subject of computational research since the seventies of last century, but the fastest development occurred since nineties. From the applicative point of view the requirement to solve electron charge density as a function of position only (three spatial coordinates) is much an advantage over the ab initio, or Born–Oppenheimer [38, 39, 40] approximation. In the latter, the wave function is not only a function of spatial coordinates position but also of spin coordinate of each electron, which for system containing N electrons results in the 4N dimension space to be considered. The specific advantage is also sometimes recognized by the fact that wave function is not measureable while electron density is a property which can be directly measured by, e. g., X ray diffraction experiments. In fact, the discussion still holds if the wave function is a real physical entity or if it is only a mathematical concept which is appropriate for particle and molecular system description. The DFT is developed without using any variable parameters and thus in its origin it is an ab initio type method. In the DFT approach any ground state property can be described by the electron charge density [41]. The total energy of the N-particle system can be represented by energy components like kinetic energy of the moving electrons, nuclear-electron attraction potential energy, the repulsion energy of electron–electron system, and exchange correlation which describe other interactions between electrons. This is represented by the equation which reads [42]:
The most problematic and computationally challenging term in the above is the last component which represents the exchange correlation Exc[ρ]. The mathematical approach for computation of this term is discussed in detail elsewhere [43, 44, 45, 46, 47, 48, 49, 50]
One of the widest used numerical approximation to this term is the B3LYP [51] energy functional (Becke, three parameter Lee-Yang-Parr). Such combinative approach to evaluating hybrid functional approximations led to improved description of many properties of molecular systems. The increased quality in atomization energies, vibration frequencies and bond lengths over the description by simple ab initio methods are its main advantage. The B3LYP correlation was demonstrated to be of comparable accuracy to preceding ab initio quantum methods [52] like e. g. Møller–Plesset perturbation theory [53].
3.2 Basis sets
The exact solution to Schrödinger equation is not possible except hydrogen and helium atoms. Numerical approaches rely on specified and defined approximation. The wave function is formulated as a vector span over infinite dimension space, which is one of the reasons why the numerical solution is not able to reproduce it in exact way. The general method of approximating the wave function is to use functional basis sets. They are finite sets of orthonormal functions which form a solution to the quantum problem which is solved. The resultant representation is then an approximation to the actual orbitals by linear combination of such functions. The physically most appropriate basis sets are Slater-type orbitals [54] (STO), which represents the solution to any atoms with one electron (hydrogen-like atoms). The most widespread functional basis sets are formed by linear combination of Gaussian-type functions [55] (GTO). The Pople [56] basis sets are represented by symbols X-YZg in which X is the number of Gaussian functions covering each core atomic orbital basis function, the variables Y and Z represent the valence orbitals by corresponding number of Gaussian primitive functions. The basis sets are subject to intensive studies and are discussed in literature along with their improving and developing [56, 57, 58, 59, 60].
This very short outline does not exhaust the topic of quantum modeling and is provided to give background for next chapters.
4 Quantum approach to VLE calculations
4.1 Continuum solvation models
The dielectric continuum solvation models are the methods of quantum calculations which provide the way to accurately describe fluid properties. The advantage of this family of methods lies in the fact that the solvent is represented by the mean continuous field rather than by system of explicitly positioned molecules. Consideration of solvent-solute molecules interaction by the averaged dielectric field significantly reduces the computational workload. In the case of explicit solvent molecules interactions, the time of calculations is often unaffordable long. Such approach enables the liquid phase calculations with the effectiveness and robustness comparable to gas phase quantum models. The quantum calculations must be conducted prior to the use of the results obtained from continuum solvation models in VLE estimations. Therefore, there is significant preparation stage at which the solvent molecules must be characterized by means of quantum calculated charge density surfaces estimations, which may take quite a long time. However, this preparation must be done one time only and in fact the user may have choose to use some of the existing databases of molecules [61].
Regarding the type of approach that is represented by such models they are classified as activity methods. The advantage of using the continuum solvation models lies in the fact that there is no need to use any group contribution parameters (like in the case of UNIFAC method), nor they require any parameters adjustment (like binary interaction parameters in Wilson or NRTL methods).
4.1.1 COSMO-RS
The early works devoted to VLE calculations by the use of quantum approach were works of Klamt [62, 63] and Taylor [64]. The continuum solvation model (COSMO) and its extension to real solvents by COSMO-RS model implementation was included in simulation software ChemSep [65] that was used to predict and simulate the process of distillation unit. COSMO-RS stands for Conductor like Screening Model for Real Solvents and is an equilibrium approach of quantum chemistry which is capable of predicting the chemical potentials in mixtures. This method is aimed at screening charge density σ on the surface of molecules to estimate the chemical potential µ of each mixture component. The implicit solvation approach represents the fluid space by continuous medium.
The COSMO-RS is a two-stage model [66]. At the first stage quantum calculations are done for every component of the mixture analyzed. The COSMO calculations are used to obtain the molecule environment imitated by virtual conductor. In such an environment, the molecules of solute generate a polarization charge density σ at the surface of the molecule which produces a back-charge acting on the molecule itself and generating higher polarized electron density in comparison with vacuum. During this quantum calculation step the energetic optimum is searched for the molecule structure contained in the virtual solvent. For the COSMO-RS model the standard quantum method used is DFT with defined basis set. Such calculation can be performed by the use of many software tools for example Gaussian [67], Turbomole [68] and Jaguar [69].
The second stage involves statistical thermodynamics calculations based on the previously structure and charge density calculated results. The liquid system is considered to be molecules ensemble, span over a lattice. The fluid particles positioned on the lattice can be factorized into three components, which reads:
The Z0 factor is the entropic contribution to the system approximated by:
which enumerates permutations of identical particles. The factor ZC is called the combinatorial factor. It represents the partition sum of molecules ensemble which interact only due to steric constraints. The Staverman–Guggenheim [70, 71] equation is used to estimate this factor:
The above formulation is also popular in many equilibrium activity or group contribution models like UNIQUAC or UNIFAC. The variables xi, Φi, Θi are mole, volume and surface fractions of ith component, respectively. The variables ri and qi are relative volume and relative area. The parameter z is the coordination number of the lattice depending on the fluid spatial model applied, which is commonly assumed to be 10. For example cubic model of molecular packing results in z = 6 (Figure 2), while for hexagonal model z = 12.
The ZR factor represents the residual contribution and is caused by non-steric molecular, electrostatic and hydrogen bond interactions. It is the most important factor for liquids. The formulation for this part is derived from the assumption of surface–pair interaction models that such residual interactions can be outlined as local pairwise interactions of surface segments. The partition factor ZR is given by:
where P counts all possible total pairs of segments. The μ and ν represent different types of surface segments and εμν is the energy of interaction of μν pair. The function pμν(P) is the total number of pairs of kind μν. Several models exist to approximate its value since it is quite complex task to evaluate all different configurations considering the number of pairs of type μν.
4.1.2 COSMOSPACE
The COSMO Surface Pair Activity Coefficient (COSMOSPACE) model is utilized to estimate the partition sum Z. By assuming that all of the molecules in the ensemble have different surfaces for example due to vibrational effects, and considering M number of segments created by contact pairs it follows that M/2 pairs can be formed. Every pair can be placed in different lattice segments, which consequently results in M! different placements. In the light of the above the formulation for Z reads:
It is significant to realize that the expression above is formulated for only one segment type ν, but for very large number of this type segments. The ν(i;P) indicates the segment of type ν residing on site i in placement P.
Finally, two parts contribute to COSMOSPACE model activity coefficients of ith component in the mixture:
where the combinatorial part is estimated by modified [72] Stavermann–Guggenheim expression:
in which variables xi, Φi,
in the above ccomb is adjustable combinatorial parameter.
and
In the above ri and qi are relative volume and surface area respectively. The residual contribution to the eq. (24) is given by:
The residual activity coefficient is a function of
where the interaction parameter τμν for physical consistency is given by symmetric matrix which elements read:
in which u is the interaction energy between groups μν, μμ and νν. Equation (30) is the general equation of COSMOSPACE model and is solved using iterative procedure. The model depends only on a few adjusted parameters for each element to be modeled and dos not rely on any functional group. Consequently, any parameterized variable is totally general and can be used to calculate properties of any compound.
The first step of the calculation procedure is the preparation the molecular model by the use of some specific software tools for modeling three-dimensional chemical structures [73, 74]. The second calculation step is to perform DFT computation with dielectric continuum solvation model, which can be implemented in some software tools [75]. These computations return the screening charge density on surface of the molecule and its total energy. In the next stage, molecules of solvent and solute are taken into account as an ensemble of pairwise interacting surfaces. The particular type of the intermolecular interactions (i. e., electrostatic interactions and hydrogen bonds) is expressed by the screening charge densities of the contacting surface segment types.
Comparison [76] of the COSMO-RS model with more traditional models like UNIFAC, modified UNIFAC and ASOG proved it to be best for systems with alkyl halides or aromatics as solutes in water. The model was less successful for nonaqueous mixtures but good results were obtained for the mixtures that contained alkyl halides, ethanol solutions of alkanes, and ketones in alkanes.
4.1.3 COSMO-SAC
Essentially similar in origin to the COSMO-RS, the COSMO-SAC model [77], is theoretically rather different. COSMO-SAC (SAC – segment activity coefficient) resolves some problems that exist in the COSMO-RS model. The model is based on group contribution solvation (GCS) model [78] where the activity coefficients are estimated from the solvation free energy of molecules in a solution. The COSMOS-SAC equation for activity coefficients reads:
In the above formulation the σm is the segment with charge density at a fixed position in the solution, Γs(σ) is the segment activity coefficient,
and pi(σ) is the probability of finding a segment with a surface charge density σ, ni is the number of surface segment in which the molecule contributes.
The segment activity coefficient is defined by the relation:
in which the
The max and min functions in the formulation (35) indicate the larger or smaller values from the arguments in brackets, respectively. The nonelectrostatic contribution is assumed to be constant and cancel out in the formulation because the nonelectrostatic energy is assumed to be constant.
The computational procedure which is performed during COSMO-SAC calculations can be divided into four steps. The first step is to obtain the charge density σ profile for each component in a mixture using a quantum chemistry package which has COSMO model implemented [79]. For example, using the DFT with Becke-Perdew (BP) version of correlation VWN−BP (Vosko-Wilk-Nusair [80]) for many electron system of the spin-polarized homogeneous electron gas functional at the DNP (double numeric with polarization functions) basis set level, the equilibrium geometry of the molecules in the ideal gas phase is established. After that, the COSMO calculation is performed to obtain the ideal screening charges on the molecular surface for each molecule. The screening charge densities σ* as the result of the COSMO calculations are averaged to give the “apparent” charge density σ on a standard surface segment using the following expression [64]:
In the above formulation, the radius of the standard surface segment reff is given by:
while the radius of nth segment reads:
The variable dmn is the distance between m and n segments. The constant α can be derived based on electrostatics [81, 82]. The charge density profile is calculated from:
where Ai is the surface area of ith molecule and the probability of finding a segment with a surface charge density σ in pure liquid ith specie is:
The second step of the calculation process is the calculation of restoring the free energy of the solute in the mixture
The third step the molecular volume and surface area acquired from the COSMO computation of each compound are normalized to a standard volume and surface area of a functional group to give the r and q parameters, which are then substituted into the Staverman−tavermanrs, which are then eq. (33).
The final fourth step consists of calculating the activity coefficient. The calculated free energies and the Staverman−Guggenheim contribution from the previous two steps are used in the COSMO-SAC model (32) to obtain the value of activity coefficient.
The quantum mechanics is required to estimate only the σ profiles for the molecules in the first step. This calculation must be done only once for each molecule regardless of the mixture in which the component is to be used.
The DFT quantum calculations become constantly more popular in the engineering area despite such methods are not widely incorporated in modern software flowsheeting tools. Comparable, and in some areas better than for traditional VLE algorithms, results show that this tool can become a wide used approach in the design area and in a simulation software as well. The requirement to perform lengthy quantum calculations for new compounds need not be an obstacle because in the case of traditional approaches the new components binary interactions parameters also need to be estimated, typically by the means of experimental measurements.
5 Quantum approach to chemical reaction kinetics description
The rectification is often used as a multifunctional process together with chemical reaction [83, 84]. Such process intensification is advantageous especially for equilibrium reactions where removing the products from reacting mixture shift the reaction toward higher conversion. During the design phase, the key part of the modeling is to achieve correct and accurate kinetics description of the reaction. This enables the designer to create project which matches the design criteria set.
5.1 Intrinsic reaction coordinate
The intrinsic reaction coordinate is a minimum energy path solving method that connects substrates, the transition state and products of the reaction. Any parametrization of reaction path s is called reaction coordinate that can be given by:
The minimum energy path or reaction path is a line in coordinate space, which connects two minima by passing the saddle point, the transition structure of a potential energy surface (Figure 3). The energy of the saddle point is assumed to be the highest value which is placed along the reaction path.
There exist several models for tracking the reaction path, e. g. Jasien and Shepard [85], Elber and Karplus [86], Gonzalez and Schlegel [87] models. The IRC can be solved by starting at the transition state and following the steepest descent pathway down to the substrate and product minima according to formulation [88]:
where s is the arc length along the path, x is the coordinate vector, and g is the gradient of the potential energy surface (PES). Due to usually very high stiffness (large value of maximal to minimal eigenvalues ratio of the matrix representing the equations) of the above equation some special numerical techniques must be used to solve it.
The kinetic reaction constants which are of most importance can be estimated using Eyirng equation that resembles Arrhenius equation. The Eyring equation reads:
where kB is the Boltzmann’s constant, h is the Planck’s constant and ΔGAct is the Gibbs activation energy.
The activation energy is the difference of energies found for different states of the molecular system. In the case of chemical reaction this is the difference between substrate and product energetic states at equilibrium and energy state of transition structure (high energy state). The starting structure for vibrational analysis is the high energy state, when performing frequency analysis the zero point energies of the molecular system are found. The zero point energy corresponds to a state of minimum possible energy. For example the zero point energy for hydrogen atom after interpolation to absolute zero temperature is given by hν/2. The gradient search is performed both ways to find substrate and product assuming correct transition structure is provided. Several methods are capable to calculate the optimized structures and among them is DFT approach.
The typical approach to IRC is firstly to find the equilibrium and transition states of the reaction system. The transition state in the simplest case can be found by manipulating the equilibrium resulting structures by changing for example bond angles, to obtain a structure of about two to four times higher energy. This is the point where actual IRC calculations are to be performed. The computation of IRC paths is done by two-way direction method. Both ways consist on relaxation of the barrier (transition) structure energy to the equilibrium states. Many software tools report the zero-state energy which by the use of formulation (44) is the base for the kinetic reaction rate description.
Another method of finding the minimum energy path or reaction path is the application of Newton trajectory approach rather than steepest descent method [89]. The starting point is a pathway described by means of geometric definition which is considered as an reaction path. Only properties of the potential energy surfaces are taken into account, and no transient behavior of the molecule is taken into account. This idea can be generalized that any gradient direction g(x) selected over the potential surface is fixed:
where r is selected unit vector of the direction of the search and the corresponding curve is Newton trajectory. A curve fits the search direction r when the gradient of potential energy surface stays parallel to the gradient r at every point along the curve x(t).
Generally, the IRC can be defined by a variational integral, which in general form reads:
which depends on n continuously differentiable functions x(t) in an n-dimensional configuration space. The IRC is frequently named as the minimum energy path. Moreover there are other reaction path models. The term minimum energy path is used for the whole category of these pathways.
where min! stands for minimization of the functional [90], a and b are the parameters of substrate and product [91]. The function E(x) is n-dimensional potential energy surface. Different approach to reaction path is proposed by introduction the path length L in the denominator as [92]:
In the formulations (47) and (48) for the non-local variational integral the function l(x’(t)) is given by formula:
Several methods, that are proposed to calculate minimum of certain integral, as a result give energy path which can be different. Such that the results satisfy the relations provided, but the question still exist which path is the true reaction path between substrates and products.
There are literature reports illustrating the usefulness of the IRC method to describe the reaction mechanism and kinetics. The calculations of chemical reaction rate of F- + CH3OOH reported by literature [93] are promising and the accuracy of calculated value with comparison with the experimental value is reported to be excellent. Value determined from the simulation is k = (1.70 ± 0.07) × 10−9 cm3/molecules and the experimental k = 1.23 × 10−9 cm3/molecules.). As an example of reaction which could be conducted in distillation column and is analyzed by quantum DFT calculations is the esterification reaction [94]. The methanol and acetic acid (and also its halides) are taken into account and corresponding kinetic properties are calculated. The values are of reasonable accuracy when comparing with available experimental data. The reaction of gas phase decarboxylations of the β-keto carboxylic acids XCOCH2COOH (X = H, OH, and CH3) [95] is another example. The optimization of structures was applied at the MP2/6-31G* level of theory. The predicted values energy barrier are reported to be with good agreement with experimental data obtained for various solvents. Another example were authors present good agreement of calculated data with experimentally measured data is a OH radical reaction with hydrofluorocarbon [96]. They used typical quantum approach MP2 level theory together with B2LYP functional using Pople basis set 6-31G∗ and 6-311++G∗∗.
6 Summary
The presented quantum approaches which may utilize the DFT as a quantum tool are intended to predict vapor–liquid equilibria and chemical reactions rates. Both areas are of fundamental importance for reactive distillation simulations. Although the presented approaches are popular in computational chemistry rather than in industrial design in chemical engineering field, they can be useful for a mixture systems which are poorly examined experimentally. The ability of COSMO methods family to predict highly complex, non-ideal behavior of mixtures, including three-phase systems is large advantage. The traditional group contribution method like UNIFAC allows to create user component from predefined functional group. But for the case when no functional group is presented in the UNIFAC database, one cannot consider the component to be non-ideal, which may generate undesirable loss in calculation accuracy. On the other side the COSMO methods family do not rely on any functional groups and needs to use some quantum approach (MP2, DFT etc.) to estimate the charge density surfaces of the molecules of the mixture components.
The design of multifunctional reactors like in the presented case the reactive distillation columns requires not only vapor–liquid equilibria knowledge but also chemical reactions rates. The typical approach relies on experimental results. Such results, depending on the complexity of the reactive system, contain typically the values of activation energy, the frequency factors, etc. To obtain pilot design of reactive distillation system which can be required to estimate the approximate costs of the equipment and the process, there is no need to perform experiments. Such experiments besides the costs may also require quite complex equipment and analytic hardware. For the initial estimation of the design it can be quite useful to rely first on the result from quantum chemical reaction rate estimation, which in fact proved to be quite accurate. Incorporating such a method is a question of future designers’ tool development although.
Funding statement: The work was supported by grant 03/32/DSPB/0707.
References
[1] Gutiérrez JP, Meindersma GW, De Haan AB. COSMO-RS-based ionic-liquid selection for extractive distillation processes. Ind Eng Chem Res 2012;51(35):11518–11529.10.1021/ie301506nSearch in Google Scholar
[2] Scheffczyk J, Redepenning C, Jens CM, Winter B, Leonhard K, Marquardt W, et al. Massive, automated solvent screening for minimum energy demand in hybrid extraction–distillation using COSMO-RS. Chem Eng Res Des Part B 2016;115:433–442.10.1016/j.cherd.2016.09.029Search in Google Scholar
[3] Jongmans MT, Schuur B, De Haan AB. Ionic liquid screening for ethylbenzene/styrene separation by extractive distillation. Ind Eng Chem Res 2011;50(18):10800–10810.10.1021/ie2011627Search in Google Scholar
[4] Biddulph MW, Kalbassi MA. New column for measurement of multicomponent distillation design efficiencies. Chem Eng Res Des 1990;68(5):453–456.Search in Google Scholar
[5] Lockett MJ. Distillation tray fundamentals. Cambridge, England and New York: Cambridge University Press, 1986.Search in Google Scholar
[6] Klemola KT, Ilme JK. Distillation efficiencies of an industrial-scale i-butane/n-butane fractionator. Ind Eng Chem Res 1996;35(12):4579–4586.10.1021/ie960390rSearch in Google Scholar
[7] AIChE. Bubble tray design manual: prediction of fractionation efficiency. New York: AIChE, 1958.Search in Google Scholar
[8] Chan H, Fair JR. Prediction of point efficiencies on sieve trays. 1. Binary systems. Ind Eng Chem Proc Des Deu 1984;23:814–819.10.1021/i200027a032Search in Google Scholar
[9] Chan H, Fair JR. Prediction of point efficiencies on sieve trays. 2. Multicomponent systems. Ind Eng Chem Proc Des Deu 1984;23:820–827.10.1021/i200027a033Search in Google Scholar
[10] Zuiderweg FJ. Sieve trays – A view of the state of the art. Chem Eng Sci 1982;37:1441–1464.10.1016/0009-2509(82)80001-8Search in Google Scholar
[11] Kister H. Distillation design. New York: McGraw-Hill, 1992Search in Google Scholar
[12] Sundmacher K, Hoffmann U. Development of a new catalytic distillation process for fuel ethers via a detailed nonequilibrium model. Chem Eng Sci 1996;51(10):2359–2368.10.1016/0009-2509(96)00092-9Search in Google Scholar
[13] Peng J, Edgar TF, Eldridge RB. Dynamic rate-based and equilibrium models for a packed reactive distillation column. Chem Eng Sci 2003;58(12):2671–2680.10.1016/S0009-2509(03)00103-9Search in Google Scholar
[14] Higler AP, Taylor R, Krishna R. Nonequilibrium modelling of reactive distillation: multiple steady states in MTBE synthesis. Chem Eng Sci 1999;54(10):1389–1395.10.1016/S0009-2509(99)00056-1Search in Google Scholar
[15] Baur R, Higler AP, Taylor R, Krishna R. Comparison of equilibrium stage and nonequilibrium stage models for reactive distillation. Chem Eng J 2000;76(1):33–47.10.1016/S1385-8947(99)00114-XSearch in Google Scholar
[16] Taylor R, Krishna R. Modelling reactive distillation. Chem Eng Sci 2000;55(22):5183–5229.10.1016/S0009-2509(00)00120-2Search in Google Scholar
[17] Higler A, Taylor R, Krishna R. Modeling of a reactive separation process using a nonequilibrium stage model. Comp Chem Eng 1998;22:S111–S118.10.1016/S0098-1354(98)00044-1Search in Google Scholar
[18] Katariya AM, Kamath RS, Moudgalya KM, Mahajani SM. Non-equilibrium stage modeling and non-linear dynamic effects in the synthesis of TAME by reactive distillation. Comp Chem Eng 2008;32(10):2243–2255.10.1016/j.compchemeng.2007.11.009Search in Google Scholar
[19] Prausnitz JM. Molecular thermodynamics of fluid-phase equilibria. Englewood Cliffs, New York: Prentice-Hall, 1969.Search in Google Scholar
[20] Pitzer KS, Simonson JM. Thermodynamics of multicomponent, miscible, ionic systems: theory and equations. J Phys Chem 1986;90(13):3005–3009.10.1021/j100404a042Search in Google Scholar
[21] Nickmand Z, Aghamiri SF. A New modified van Laar model for associating mixtures. J Dispersion Sci Technol 2010;31(12):1638–1647.10.1080/01932690903297413Search in Google Scholar
[22] Hanks RW, Tan RL, Christensen JJ. Limits on the simultaneous correlation of gE and hE data by the NRTL, LEMF and Wilson’s equations. Thermochimica Acta 1978;23(1):41–55.10.1016/0040-6031(78)85110-7Search in Google Scholar
[23] Arlt W, Onken U. Liquid-liquid equilibria of organic compounds: measurement, correlation and prediction. Chem Eng Commun 1982;15(1–4):207–213.10.1080/00986448208911070Search in Google Scholar
[24] Eng R, Sandler SI. Vapor-liquid equilibria for three aldehyde/hydrocarbon mixtures. J Chem Eng Data 1984;29(2):156–161.10.1021/je00036a017Search in Google Scholar
[25] Malanowzki S, Skjold-Jørgensen S, Rasmussen P, Fredenslund A. Simultaneous representation of binary VLE, LLE and HE data using the UNIQUAC model. Chem Eng Sci 1981;36(10):1727–1730.10.1016/0009-2509(81)80018-8Search in Google Scholar
[26] Nagata I, Tamura K, Yamada T. Correlation of liquid-liquid equilibria in aqueous and organic systems using a modified Wilson model. J Solution Chem 1996;25:567–587.10.1007/BF00973086Search in Google Scholar
[27] Yee D, Simonetty J, Tassios D. Prediction of vapor-liquid equilibrium from ternary liquid-liquid equilibrium data. Ind Eng Chem Process Des Dev 1983;22(1):123–129.10.1021/i200020a020Search in Google Scholar
[28] Surana RK, Danner RP, De Haan AB, Beckers N. New technique to measure high-pressure and high-temperature polymer-solvent vapor-liquid equilibrium. Fluid Phase Equilib 1997;139(1):361–370.10.1016/S0378-3812(97)00172-6Search in Google Scholar
[29] Redlich O, Kwong JN. On the thermodynamics of solutions. V. An equation of state fugacities of gaseous solutions. Chem Rev 1949;44(1):233–244.10.1021/cr60137a013Search in Google Scholar
[30] Soave G. Equilibrium constants from a modified Redlich-Kwong equation of state. Chem Eng Sci 1972;27(6):1197–1203. ISSN 0009-2509.10.1016/0009-2509(72)80096-4Search in Google Scholar
[31] Peng DY, Robinson DB. A new two-constant equation of state. Ind Eng Chem Fundam 1976;15:59–64.10.1021/i160057a011Search in Google Scholar
[32] Zhi Y, Meiren S, Jun S, Lee H. A new quartic equation of state. Fluid Phase Equilib 2001;187–188:275–298.10.1016/S0378-3812(01)00542-8Search in Google Scholar
[33] Huang SH, Radosz M. Equation of state for small, large, polydisperse, and associating molecules. Ind Eng Chem Res 1990;29:2284–2294.10.1021/ie00107a014Search in Google Scholar
[34] Chapman WG, Gubbins KE, Jackson G, Radosz M. New reference equation of state for associating liquids. Ind Eng Chem Res 1990;29:1709–1721.10.1021/ie00104a021Search in Google Scholar
[35] Fu YH, Sandler S. A simplified SAFT equation of state for associating compounds and mixtures. Ind Eng Chem Res 1995;34:1897–1909.10.1021/ie00044a042Search in Google Scholar
[36] Huang SH, Radosz M. Equation of state for small, large, polydisperse, and associating molecules: extension to fluid mixtures. Ind Eng Chem Res 1991;30:1994–2005.10.1021/ie00056a050Search in Google Scholar
[37] Mansoori GA, Carnahan NF, Starling KE, Leland TW. Equilibrium thermodynamic properties of the mixture of hard spheres. J Chem Phys 1971;4:1523–1525.10.1063/1.1675048Search in Google Scholar
[38] Born M, Oppenheimer R. Zur quantentheorie der molekeln. Annalen der Physik 1927;389(20):457–484.10.1002/andp.19273892002Search in Google Scholar
[39] Sutcliffe BT, Woolley RG. On the quantum theory of molecules. J Chem Phys 2012;137:22–31.10.1063/1.4755287Search in Google Scholar PubMed
[40] Sutcliffe BT, Woolley RG. Comment on “On the quantum theory of molecules” [J. Chem. Phys. 137, 22A544 (2012)]. J Chem Phys 2014;140:3–4.10.1063/1.4861897Search in Google Scholar PubMed
[41] Hohenberg P, Kohn W. Inhomogeneous electron gas. Phys Rev 1964;136 :B864–B870.10.1103/PhysRev.136.B864Search in Google Scholar
[42] Kohn W, Sham LJ. Self-consistent equations including exchange and correlation effects. Phys Rev 1965;140:A1133–A1140.10.1103/PhysRev.140.A1133Search in Google Scholar
[43] Perdew JP, Chevary JA, Vosko SH, Jackson KA, Pederson MR, Singh DJ, et al. Atoms, molecules, solids, and surfaces: applications of the generalized gradient approximation for exchange and correlation. Phys Rev B 1992;46:6671–6676. Erratum Phys Rev B 1993, 48, 4978.10.1103/PhysRevB.46.6671Search in Google Scholar
[44] Perdew JP, Burke K, Ernzerhof M. Generalized gradient approximation made simple [Phys. Rev. Lett. 77, 3865 (1996)]. Phys Rev Lett 1997;78:1396–1402.10.1103/PhysRevLett.78.1396Search in Google Scholar
[45] Becke AD. Density-functional exchange-energy approximation with correct asymptotic behavior. Phys Rev A 1988;38:3098–3100.10.1103/PhysRevA.38.3098Search in Google Scholar
[46] Śmiga S, Buksztel A, Grabowski I. Chapter 7 – density-dependent exchange–correlation potentials derived from highly accurate ab initio calculations. In: Hoggan P, editor(s). Advances in quantum chemistry. Cambridge, Massachusetts: Academic Press, 2014Search in Google Scholar
[47] Cohen AJ, Mori-Sánchez P, Yang W. Challenges for density functional theory. Chem Rev 2012;112(1):289–320.10.1021/cr200107zSearch in Google Scholar PubMed
[48] Bilc DI, Orlando R, Shaltaf R, Rignanese GM, Íñiguez J, Ghosez P. Hybrid exchange-correlation functional for accurate prediction of the electronic and structural properties of ferroelectric oxides. Phys Rev B 2008;77:165107–165129.10.1103/PhysRevB.77.165107Search in Google Scholar
[49] Baer R, Neuhauser D. Density functional theory with correct long-range asymptotic behavior. Phys Rev Lett 2005;94:043002–5.10.1103/PhysRevLett.94.043002Search in Google Scholar PubMed
[50] Sen KD, Luque FJ. First-order correlation-kinetic contribution to Kohn–Sham exchange charge density function in atoms, using quantal density functional theory approach. Int J Quantum Chem 2005;101:231–238.10.1002/qua.20262Search in Google Scholar
[51] Becke AD. Density‐functional thermochemistry. III. The role of exact exchange. J Chem Phys 1993;98:5648–5652.10.1063/1.464913Search in Google Scholar
[52] Kim K, Jordan KD. Comparison of density functional and MP2 calculations on the water monomer and dimer. J Phys Chem 1994;98(40):10089–10094.10.1021/j100091a024Search in Google Scholar
[53] Møller C, Plesset MS. Note on an approximation treatment for many-electron systems. Phys Rev 1934;46:618–622.10.1103/PhysRev.46.618Search in Google Scholar
[54] Davidson E, Feller D. Basis set selection for molecular calculations. Chem Rev 1986;86(4):681–696.10.1021/cr00074a002Search in Google Scholar
[55] Gill PM. Molecular integrals over Gaussian basis functions”. Adv Quantum Chem 1994;25:141–205.10.1016/S0065-3276(08)60019-2Search in Google Scholar
[56] Ditchfield R, Hehre WJ, Pople JA. Self-consistent molecular-orbital methods. IX. An extended Gaussian-type basis for molecular-orbital studies of organic molecules. J Chem Phys 1971;54(2):724–728.10.1063/1.1674902Search in Google Scholar
[57] Manninen P, Vaara J. Systematic Gaussian basis-set limit using completeness-optimized primitive sets. A case for magnetic properties. J Comput Chem 2006;27(4):434–445.10.1002/jcc.20358Search in Google Scholar PubMed
[58] Kresse G, Furthmüller J. Efficiency of ab-initio total energy calculations for metals and semiconductors using a plane-wave basis set. Comput Materials Sci 1996;6(1):15–50.10.1016/0927-0256(96)00008-0Search in Google Scholar
[59] Zakharov M. Performance of numerical atom-centered basis sets in the ground-state correlated calculations of noncovalent interactions: water and methane dimer cases. Int J Quantum Chem 2013;113(15):1899–1918.10.1002/qua.24407Search in Google Scholar
[60] Ema I, García De La Vega JM, Ramírez G, López R, Fernández Rico J, Meissner H, et al. Polarized basis sets of Slater-type orbitals: H to Ne atoms. J Comput Chem 2003;24:859–868.10.1002/jcc.10227Search in Google Scholar PubMed
[61] Klamt A, Eckert F. Prediction of vapor liquid equilibria using COSMOtherm. Fluid Phase Equilib 2004;217(1):53–57.10.1016/j.fluid.2003.08.018Search in Google Scholar
[62] Klamt A. Conductor-like screening model for real solvents: a new approach to the quantitative calculation of solvation phenomena. J Phys Chem 1995;99(7):2224–2235.10.1021/j100007a062Search in Google Scholar
[63] Klamt A, Jonas V, Bürger T, Lohrenz JC. Refinement and Parametrization of COSMO-RS. J Phys Chem A 1998;102(26):5074–5085.10.1021/jp980017sSearch in Google Scholar
[64] Taylor R, Kooijman HA, Klamt A, Eckert F. Distillation simulation with COSMO−RS. Presented at the Distillation & Absorption Conference, Baden-Baden, Germany, 2002.Search in Google Scholar
[65] Kooijman HA, Taylor R. ChemSep – another software system for the simulation of separation processes. CACHE News 1992;35:1–9.Search in Google Scholar
[66] Klamt A, Krooshof GJ, Taylor R. COSMOSPACE: alternative to conventional activity-coefficient models. AIChE J 2002;48:2332–2349.10.1002/aic.690481023Search in Google Scholar
[67] Hehre WJ, Lathan WA, Ditchfield R, Newton MD, Pople JA. Gaussian 70. Quantum chemistry program exchange, Program No. 237, 1970.Search in Google Scholar
[68] Ahlrichs R, Bär M, Häser M, Horn H, Kölmel C. Electronic structure calculations on workstation computers: the program system turbomole. Chem Phys Lett 1989;162(3):165–169.10.1016/0009-2614(89)85118-8Search in Google Scholar
[69] Bochevarov AD, Harder E, Hughes TF, Greenwood JR, Braden DA, Philipp DM, et al. Jaguar: a high-performance quantum chemistry software program with strengths in life and materials sciences. Int J Quantum Chem 2013;113(18):2110–2142.10.1002/qua.24481Search in Google Scholar
[70] Staverman A. The entropy of high polymer solutions. Generalization of formulae. Recueil des Travaux Chimiques des Pays-Bas 1950;69(2):163–174.10.1002/recl.19500690203Search in Google Scholar
[71] Guggenheim EA. Mixtures: the theory of the equilibrium properties of some simple classes of mixtures, solutions and alloys. Wotton-under-Edge: Clarendon Press, 1952Search in Google Scholar
[72] Bosse D, Bart HJ. Viscosity calculations on the basis of Eyring’s absolute reaction rate theory and COSMOSPACE. Ind Eng Chem Res 2005;44(22):8428–8435.10.1021/ie048797gSearch in Google Scholar
[73] ACD/Labs Toronto. ACD/Structure Elucidator, version 15.01. Toronto, ON, Canada: Advanced Chemistry Development, Inc., 2015Search in Google Scholar
[74] Evans David A. History of the Harvard ChemDraw Project. Angewandte Chemie International Edition. 2014 8 11;53(42):11140–11145. DOI: 10.1002/anie.20140582010.1002/anie.201405820Search in Google Scholar PubMed
[75] TURBOMOLE GmbH. TURBOMOLE: program package for ab initio electronic structure calculations. Leverkusen, Germany: COSMOlogic, 2017Search in Google Scholar
[76] Putnam R, Taylor R, Klamt A, Eckert F, Schiller M. Prediction of infinite dilution activity coefficients using COSMO-RS. Ind Eng Chem Res 2003;42(15):3635–3641.10.1021/ie020974vSearch in Google Scholar
[77] Lin ST, Sandler SI. A priori phase equilibrium prediction from a segment contribution solvation model. Ind Eng Chem Res 2002;41(5):899–913.10.1021/ie001047wSearch in Google Scholar
[78] Lin ST, Sandler SI. Infinite dilution activity coefficients from Ab Initio solvation calculations. AIChE J 1999;45:2606–2618.10.1002/aic.690451217Search in Google Scholar
[79] Accelrys Software Inc.. Cerius2, Dmol3, version 4.0. San Diego, CA: Molecular Simulations Inc., 1999Search in Google Scholar
[80] Vosko SH, Wilk L, Nusair M. Accurate spin-dependent electron liquid correlation energies for local spin density calculations: a critical analysis. Can J Phys 1980;58(8):1200–1212.10.1139/p80-159Search in Google Scholar
[81] Klamt A, Schüürmann GJ. COSMO: a new approach to dielectric screening in solvents with explicit expressions for the screening energy and its gradient. Chem Soc Perkin Trans 1993;2:799–805.10.1039/P29930000799Search in Google Scholar
[82] Klamt A. Conductor-like screening model for real solvents: a new approach to the quantitative calculation of solvation phenomena. J Phys Chem 1995;99:2224–2230.10.1021/j100007a062Search in Google Scholar
[83] Mansouri SS, Sales-Cruz M, Huusom JK, Woodley JM, Gani R. Integrated process design and control of reactive distillation processes. IFAC-PapersOnLine 2015;48(8):1120–1125.10.1016/j.ifacol.2015.09.118Search in Google Scholar
[84] Pérez-Cisneros ES, Mena-Espino X, Rodríguez-López V, Sales-Cruz M, Viveros-García T, Lobo-Oehmichen R. An integrated reactive distillation process for biodiesel production. Comp Chem Eng 2016;91:233–246.10.1016/j.compchemeng.2016.01.008Search in Google Scholar
[85] Jasien PG, Shepard R. A general polyatomic potential energy surface fitting method. Int J Quantum Chem 1988;34:183–198.10.1002/qua.560340822Search in Google Scholar
[86] Elber R, Karplus M. A method for determining reaction paths in large molecules: application to myoglobin. Chem Phys Lett 1987;139(5):375–380.10.1016/0009-2614(87)80576-6Search in Google Scholar
[87] Gonzalez C, Schlegel HB. Reaction path following in mass-weighted internal coordinates. J Phys Chem 1990;94(14):5523–5527.10.1021/j100377a021Search in Google Scholar
[88] Hratchian HP, Schlegel HB. Using Hessian updating to increase the efficiency of a Hessian based predictor-corrector reaction path following method. J Chem Theory Comput 2005;1(1):61–69.10.1021/ct0499783Search in Google Scholar PubMed
[89] Hirsch M, Quapp W. Reaction pathways and convexity of the potential energy surface: application of newton trajectories. J Math Chem 2004;36(4):307–340.10.1023/B:JOMC.0000044520.03226.5fSearch in Google Scholar
[90] Zeidler E. Nonlinear functional analysis and its applications, II/B: nonlinear monotone operators. New York: Springer-Verlag, 1990.Search in Google Scholar
[91] Sevick EM, Bell AT, Theodorou DN. A chain of states method for investigating infrequent event processes occurring in multistate, multidimensional systems. J Chem Phys 1993;98(4):3196–3212.10.1063/1.464093Search in Google Scholar
[92] Pratt LR. A statistical method for identifying transition states in high dimensional problems. J Chem Phys 1986;85(9):5045–5048.10.1063/1.451695Search in Google Scholar
[93] López JG, Vayner G, Lourderaj U, Addepalli SV, Kato S, DeJong WA, et al. A direct dynamics trajectory study of F- + CH3OOH reactive collisions reveals a major non-IRC reaction path. J Am Chem Soc 2007;129(32):9976–9985.10.1021/ja0717360Search in Google Scholar PubMed
[94] Lawal MM, Govender T, Maguire GE, Honarparvar B, Kruger HG. Mechanistic investigation of the uncatalyzed esterification reaction of acetic acid and acid halides with methanol: a DFT study. J Mol Model 2016;22:235–246.10.1007/s00894-016-3084-zSearch in Google Scholar PubMed
[95] Huang CL, Wu CC, Lien MH. Ab initio studies of decarboxylations of the β-Keto Carboxylic Acids XCOCH2COOH (X = H, OH, and CH3). J Phys Chem A 1997;101(42):7867–7873.10.1021/jp9712664Search in Google Scholar
[96] Ali MA, Rajakumar B. Kinetics of OH radical reaction with CF3CHFCH2F (HFC-245eb) between 200 and 400 K: G3MP2, G3B3 and transition state theory calculations. J Mol Struct 2010;949(1–3):73–81.10.1016/j.theochem.2010.03.006Search in Google Scholar
© 2017 Walter de Gruyter GmbH, Berlin/Boston
