**Numerical Solution of Linear Ordinary Differential Equations in Quantum Chemistry by Spectral Method**

Masoud Saravi1 and Seyedeh-Razieh Mirrajei2 *1Islamic Azad University, Nour Branch, Nour, 2Education Office of Amol, Amol, Iran* 

#### **1. Introduction**

The problem of the structure of hydrogen atom is the most important problem in the field of atomic and molecular structure. Bahr's treatment of the hydrogen atom marked the beginning of the old quantum theory of atomic structure, and wave mechanics had its inception in Schrodinger 's first paper, in which he gave the solution of the wave equation for the hydrogen atom. Since the most differential equations concerning physical phenomenon could not be solved by analytical method hence, the solutions of the wave equation are based on polynomial (series) methods. Even if we use series method, some times we need an appropriate change of variable, and even when we can, their closed form solution may be so complicated that using it to obtain an image or to examine the structure of the system is impossible. For example, if we consider Schrodinger equation, i.e.,

$$
\varphi'' + (2mEh^{-2} - a^2x^2)\varphi = 0,
$$

we come to a three-term recursion relation, which work with it takes, at least, a little bit time to get a series solution. For this reason we use a change of variable such as

$$
\varphi = e^{-ax^4/2} f(x),
$$

or when we consider the orbital angular momentum, it will be necessary to solve

$$\frac{d^2\mathbf{s}}{d\theta^2} + \cot\theta\frac{d\mathbf{s}}{d\theta} + \left(\frac{c}{h^2} - \frac{m^2}{\sin^2\theta}\right)\mathbf{s} = \mathbf{0}.$$

As we can observe, working with this equation is tedious. Another two equations which occur in the hydrogen atom wave equations, are Legendre and Laguerre equations, which can be solved only by power series methods.

In next section, after a historical review of spectral methods we introduce Clenshaw method, which is a kind of spectral method, and then solve such equations in last section. But, first of all, we put in mind that this method can not be applied to atoms with more electrons. With

Numerical Solution of Linear Ordinary Differential

consider the following differential equation:

( )

*i*

very smooth basis functions, as given below

(2), we define the residual term by ( ) *Nr x* as follows

are bounded variation.

**3. Clenshaw method** 

where 0

*M*

unknowns 0 1 , ,..., . *<sup>N</sup> aa a*

of a Chebyshev series expansion.

Consider the following differential equation:

Equations in Quantum Chemistry by Spectral Method 5

is analytic on [,] *a b* if is infinitely differentiable and with all its derivatives on this interval

In this section, we are going to introduce Clenshaw method. For this reason, first we

*i*

( ) ( ), [ 1,1],

*M i L f xD* , and *if* , *i Mf* 0,1,..., , , are known real functions of , *<sup>i</sup> <sup>x</sup> <sup>D</sup>* denotes

*th <sup>i</sup>* order of differentiation with respect to *<sup>x</sup>* , is a linear functional of rank *M* and *<sup>M</sup> <sup>C</sup>* .

Here (3) can be initial, boundary or mixed conditions. The basis of spectral methods to solve this class of equations is to expand the solution function, *y* , in (2) and (3) as a finite series of

> 0 () () *N N n n n*

where, <sup>0</sup> ( ) *<sup>N</sup> T x <sup>n</sup>* is sequence of Chebyshev polynomials of the first kind. By replacing *Ny* in

In spectral methods, the main target is to minimize ( ) *Nr x* , throughout the domain as much as possible with regard to (3), and in the sense of pointwise convergence. Implementation of these methods leads to a system of linear equations with *N* 1 equations and *N* 1

are assumed to be elements of a complete set of orthonormal functions. The approximate

the number of independent boundary constraints *BuN* 0 that must be applied. Here we are going to use a Tau method developed by Clenshaw for the solution of linear ODE in terms

The Tau method was invented by Lanczos in 1938. The expansion functions

solution is assumed to be expanded in terms of those functions as

*M i Ly f x D y f x x* (2)

*y C* , (3)

*y x aT x* , (4)

( ) *N N r x Ly f* . (5)

,

*u a* where *m* is

*N m N nn n*

 1 *<sup>n</sup>* ( 1,2,3,...) *n*

0

*M*

the increasing complexity of the atom, the labour of making calculations increases tremendously. In these cases, one can use variation or perturbation methods for overcoming such problems.

#### **2. Historical review**

Spectral methods arise from the fundamental problem of approximation of a function by interpolation on an interval, and are very much successful for the numerical solution of ordinary or partial differential equations. Since the time of Fourier (1882), spectral representations in the analytic study of differential equations have been used and their applications for numerical solution of ordinary differential equations refer, at least, to the time of Lanczos.

Spectral methods have become increasingly popular, especially, since the development of Fast transform methods, with applications in problems where high accuracy is desired.

Spectral methods may be viewed as an extreme development of the class of discretization schemes for differential equations known generally as the *method of weighted residuals* (MWR) (Finlayson and Scriven (1966)). The key elements of the MWR are the trial functions (also called expansion approximating functions) which are used as basis functions for a truncated series expansion of the solution, and the test functions (also known as weight functions) which are used to ensure that the differential equation is satisfied as closely as possible by the truncated series expansion. The choice of such functions distinguishes between the three most commonly used spectral schemes, namely, Galerkin, Collocation(also called Pseudospectral) and Tau version. The Tau approach is a modification of Galerkin method that is applicable to problems with non-periodic boundary conditions. In broad terms, Galerkin and Tau methods are implemented in terms of the expansion coefficients, where as Collocation methods are implemented in terms of physical space values of the unknown function.

The basis of spectral methods to solve differential equations is to expand the solution function as a finite series of very smooth basis functions, as follows

$$\text{l.} \mathcal{y}\_N\{\mathbf{x}\} = \sum\_{n=0}^N a\_n \phi\_n(\mathbf{x}) \,. \tag{1}$$

in which, one of our choice of *<sup>n</sup>* , is the eigenfunctions of a singular Sturm-Liouville problem. If the solution is infinitely smooth, the convergence of spectral method is more rapid than any finite power of *1/N.* That is the produced error of approximation (1), when *N* , approaches zero with exponential rate. This phenomenon is usually referred to as "spectral accuracy". The accuracy of derivatives obtained by direct, term by term differentiation of such truncated expansion naturally deteriorates. Although there will be problem but for high order derivatives truncation and round off errors may deteriorate, but for low order derivatives and sufficiently high-order truncations this deterioration is negligible. So, if the solution function and coefficient functions of the differential equation are analytic on[,] *a b* , spectral methods will be very efficient and suitable. We call function *y*

is analytic on [,] *a b* if is infinitely differentiable and with all its derivatives on this interval are bounded variation.

### **3. Clenshaw method**

4 Quantum Chemistry – Molecules for Innovations

the increasing complexity of the atom, the labour of making calculations increases tremendously. In these cases, one can use variation or perturbation methods for overcoming

Spectral methods arise from the fundamental problem of approximation of a function by interpolation on an interval, and are very much successful for the numerical solution of ordinary or partial differential equations. Since the time of Fourier (1882), spectral representations in the analytic study of differential equations have been used and their applications for numerical solution of ordinary differential equations refer, at least, to the

Spectral methods have become increasingly popular, especially, since the development of Fast transform methods, with applications in problems where high accuracy is desired.

Spectral methods may be viewed as an extreme development of the class of discretization schemes for differential equations known generally as the *method of weighted residuals* (MWR) (Finlayson and Scriven (1966)). The key elements of the MWR are the trial functions (also called expansion approximating functions) which are used as basis functions for a truncated series expansion of the solution, and the test functions (also known as weight functions) which are used to ensure that the differential equation is satisfied as closely as possible by the truncated series expansion. The choice of such functions distinguishes between the three most commonly used spectral schemes, namely, Galerkin, Collocation(also called Pseudospectral) and Tau version. The Tau approach is a modification of Galerkin method that is applicable to problems with non-periodic boundary conditions. In broad terms, Galerkin and Tau methods are implemented in terms of the expansion coefficients, where as Collocation methods are implemented in terms of physical space values of the unknown

The basis of spectral methods to solve differential equations is to expand the solution

 0 () () *N N n n n*

problem. If the solution is infinitely smooth, the convergence of spectral method is more rapid than any finite power of *1/N.* That is the produced error of approximation (1), when *N* , approaches zero with exponential rate. This phenomenon is usually referred to as "spectral accuracy". The accuracy of derivatives obtained by direct, term by term differentiation of such truncated expansion naturally deteriorates. Although there will be problem but for high order derivatives truncation and round off errors may deteriorate, but for low order derivatives and sufficiently high-order truncations this deterioration is negligible. So, if the solution function and coefficient functions of the differential equation are analytic on[,] *a b* , spectral methods will be very efficient and suitable. We call function *y*

*y x ax* , (1)

*<sup>n</sup>* , is the eigenfunctions of a singular Sturm-Liouville

function as a finite series of very smooth basis functions, as follows

such problems.

time of Lanczos.

function.

in which, one of our choice of

**2. Historical review** 

In this section, we are going to introduce Clenshaw method. For this reason, first we consider the following differential equation:

$$Ly = \sum\_{0}^{M} f\_{M-i} \{ \mathbf{x} \} D^{i} y = f \{ \mathbf{x} \}, \mathbf{x} \in [-1, 1], \tag{2}$$

*y C* , (3)

where 0 ( ) *M i M i L f xD* , and *if* , *i Mf* 0,1,..., , , are known real functions of , *<sup>i</sup> <sup>x</sup> <sup>D</sup>* denotes

*th <sup>i</sup>* order of differentiation with respect to *<sup>x</sup>* , is a linear functional of rank *M* and *<sup>M</sup> <sup>C</sup>* .

Here (3) can be initial, boundary or mixed conditions. The basis of spectral methods to solve this class of equations is to expand the solution function, *y* , in (2) and (3) as a finite series of very smooth basis functions, as given below

$$\text{If } \mathbf{y}\_N \{ \mathbf{x} \} = \sum\_{n=0}^N a\_n T\_n \{ \mathbf{x} \} \tag{4}$$

where, <sup>0</sup> ( ) *<sup>N</sup> T x <sup>n</sup>* is sequence of Chebyshev polynomials of the first kind. By replacing *Ny* in (2), we define the residual term by ( ) *Nr x* as follows

$$\mathbf{r}\_N \{ \mathbf{x} \} = \mathbf{L} \mathbf{y}\_N - \mathbf{f} \ . \tag{5}$$

In spectral methods, the main target is to minimize ( ) *Nr x* , throughout the domain as much as possible with regard to (3), and in the sense of pointwise convergence. Implementation of these methods leads to a system of linear equations with *N* 1 equations and *N* 1 unknowns 0 1 , ,..., . *<sup>N</sup> aa a*

The Tau method was invented by Lanczos in 1938. The expansion functions *<sup>n</sup>* ( 1,2,3,...) *n* are assumed to be elements of a complete set of orthonormal functions. The approximate solution is assumed to be expanded in terms of those functions as 1 , *N m N nn n u a* where *m* is the number of independent boundary constraints *BuN* 0 that must be applied. Here we are going to use a Tau method developed by Clenshaw for the solution of linear ODE in terms of a Chebyshev series expansion.

Consider the following differential equation:

Numerical Solution of Linear Ordinary Differential

 cos , 0,1,..., . *<sup>k</sup> k x k N*

and using <sup>1</sup> ( ) cos( cos ), *Tx i x <sup>i</sup>* we get

"

and nodes

where, notation

we get

Equations in Quantum Chemistry by Spectral Method 7

 

*k ik <sup>P</sup> N NN*

> 

 

 

(1) (1) (1)

m p odd

(2) (2) 2 2 (2) (2)

m p even

*y x a T x a pp m a m N a a*

<sup>1</sup> ( ) ( ), ( ) , 0,1,..., 2 , 0 ,

*mm m p N N*

() () () () ( ) ( ) ( ),

( 1) ,

(1) .

*a T xT x a T xT x a T xT x S x* (10)

*yx a T x a pa m N a*

<sup>2</sup> ( ) ( ), , 0,1,..., 1, 0,

*mm m p N*

*k ik <sup>R</sup> N NN*

*k ik <sup>P</sup> N NN*

(cos( ))cos( ),

(cos( ))cos( ) .

*k ik <sup>Q</sup> N NN* (9)

*c*

(11)

1

*k k P T N NN* ,

(cos( )) (cos( ))

(cos( ))cos( ),

means first and last terms become half .Therefore, we will have :

(cos( ))cos( ),

 

*<sup>N</sup>* That is, we put:

 "

*k*

*N*

0

 "

*k*

 "

*k*

*N*

 "

 "

*k*

0

*k*

*N*

0

*N*

0

*N*

0

*i i*

*i*

*i*

Now, substituting (4) and (9) in equations (6), and using the fact that

*N N*

*m m p m*

*c*

0 1

*i N i i*

*N i i*

*i*

0

 0

(2) (1)

*a T*

*a T*

*im i m im i m imi m*

*i m i m i m*

*N N N N N N*

0 0 0 0 0 0

*mp m*

*N*

2

*i*

*i*

$$\begin{aligned} P\{\mathbf{x}\} \mathbf{y''} + Q\{\mathbf{x}\} \mathbf{y'} + R\{\mathbf{x}\} \mathbf{y} &= S\{\mathbf{x}\}, \mathbf{x} \in \{-1, 1\}, \\ \mathbf{y}(-1) &= \alpha \end{aligned} \tag{6}$$
  $\mathbf{y}(-1) = \alpha$   $\mathbf{y}(1) = \boldsymbol{\beta}$   $\mathbf{z}$ 

First, for an arbitrary natural number *N* , we suppose that the approximate solution of equations (6) is given by (4). Our target is to find 0 1 ( , ,..., )*<sup>t</sup> <sup>N</sup> a aa a* . For this reason, we put

$$\begin{aligned} P\{\mathbf{x}\} & \equiv \sum\_{i=0}^{N} \tilde{\boldsymbol{\varphi}}\_{i} T\_{i}\{\mathbf{x}\}, \\ Q\{\mathbf{x}\} & \equiv \sum\_{i=0}^{N} \boldsymbol{\gamma}\_{i} T\_{i}\{\mathbf{x}\}, \\ R\{\mathbf{x}\} & \equiv \sum\_{i=0}^{N} \boldsymbol{\lambda}\_{i} T\_{i}\{\mathbf{x}\}. \end{aligned} \tag{7}$$

Using this fact that the Chebyshev expansion of a function <sup>2</sup> ( 1,1) *<sup>w</sup> u L* is

 1 1 0 <sup>2</sup> ( ) ( ); ˆ ˆ () () () *kk k <sup>k</sup> k k u x u T x u u x T x w x dx c* , we can find coefficients , *i i* and *<sup>i</sup>* as

follows:

$$\begin{aligned} \xi\_i &= \frac{2}{\pi c\_i} \int\_{-1}^{1} \frac{P\{\mathbf{x}\} T\_i\{\mathbf{x}\}}{\sqrt{1 - \mathbf{x}^2}} d\mathbf{x} \\ \mathcal{I}\_i &= \frac{2}{\pi c\_i} \int\_{-1}^{1} \frac{Q\{\mathbf{x}\} T\_i\{\mathbf{x}\}}{\sqrt{1 - \mathbf{x}^2}} d\mathbf{x} \\ \mathcal{A}\_i &= \frac{2}{\pi c\_i} \int\_{-1}^{1} \frac{R\{\mathbf{x}\} T\_i\{\mathbf{x}\}}{\sqrt{1 - \mathbf{x}^2}} d\mathbf{x} \end{aligned} \tag{8}$$

where, *<sup>c</sup>*<sup>0</sup> <sup>2</sup> and 1 *<sup>i</sup> c* for *<sup>i</sup>* 1.

To compute the right-hand side of (8) it is sufficient to use an appropriate numerical integration method. Here, we use ( 1) *N* - point Gauss-Chebyshev-Lobatto quadrature

$$\mathbf{x}\_{j} = \cos \frac{\pi \ j}{N}, \mathbf{w}\_{j} = \frac{\pi}{\tilde{c}\_{j} N}, \ 0 \le j \le N \ \mathbf{y}\_{j}$$

where *c c* <sup>0</sup> *<sup>N</sup>* <sup>2</sup> and <sup>1</sup> *<sup>j</sup> <sup>c</sup>* for *j N* 1,2,..., 1 .

Note that, for simplicity of the notation, these points are arranged in descending order, namely, 1 10 ... *N N x x xx* , with weights

$$\mathbf{w}\_k = \frac{\pi}{N} \quad , \quad 1 \le k \le N - 1 \quad ,$$

$$= \frac{\pi}{2N} \quad , \quad k = \mathbf{0} \text{ }, k = N \quad ,$$

First, for an arbitrary natural number *N* , we suppose that the approximate solution of

*i N i i*

0

( ) ( ),

*Px T x*

*N i i*

Using this fact that the Chebyshev expansion of a function <sup>2</sup> ( 1,1) *<sup>w</sup> u L* is

*i*

*i*

*i*

1

*i*

*i*

*i*

*i*

*P xT x dx*

*i*

*QxT x dx*

*i*

*RxT x dx*

1 2

2 ()() 1 2 ()() 1 <sup>2</sup> ()() , 1

*c x*

*c x*

*c x*

To compute the right-hand side of (8) it is sufficient to use an appropriate numerical integration method. Here, we use ( 1) *N* - point Gauss-Chebyshev-Lobatto quadrature

Note that, for simplicity of the notation, these points are arranged in descending order,

*w kN <sup>k</sup> <sup>N</sup>*

, 1 1,

*k kN*

, 0, , <sup>2</sup>

*N*

1 2

1 2

1

1

*u x u T x u u x T x w x dx c*

<sup>2</sup> ( ) ( ); ˆ ˆ () () () *kk k <sup>k</sup>*

1 1

 

*k k*

where, *<sup>c</sup>*<sup>0</sup> <sup>2</sup> and 1 *<sup>i</sup> c* for *<sup>i</sup>* 1.

 cos , , 0 *j j j <sup>j</sup> x w jN N cN* ,

where *c c* <sup>0</sup> *<sup>N</sup>* <sup>2</sup> and <sup>1</sup> *<sup>j</sup> <sup>c</sup>* for *j N* 1,2,..., 1 .

namely, 1 10 ... *N N x x xx* , with weights

0

follows:

*i N i i*

0

( ) ( ),

*Qx T x*

*i*

0

( ) ( ).

*Rx T x*

( ) ( ) ( ) ( ) , ( 1,1) ,

*y y* (6)

, we can find coefficients

*<sup>N</sup> a aa a* . For this reason, we put

(7)

 ,*i i* and

(8)

*<sup>i</sup>* as

 

*Px y Qx y Rx y Sx x*

equations (6) is given by (4). Our target is to find 0 1 ( , ,..., )*<sup>t</sup>*

( 1) , (1) .

and nodes cos , 0,1,..., . *<sup>k</sup> k x k N <sup>N</sup>* That is, we put:

$$\xi\_i \equiv \frac{\pi}{N} \sum\_{k=0}^{N} P\{\cos(\frac{k\pi}{N})\} T\_i(\cos(\frac{k\pi}{N})) $$

and using <sup>1</sup> ( ) cos( cos ), *Tx i x <sup>i</sup>* we get

$$\zeta\_i \equiv \frac{\pi}{N} \sum\_{k=0}^{N} P\{\cos(\frac{k\pi}{N})\} \cos(\frac{\pi ik}{N}).$$

where, notation " means first and last terms become half .Therefore, we will have :

$$\begin{aligned} \dot{\xi}\_{l} & \equiv \frac{\pi}{N} \sum\_{k=0}^{N} \text{P} (\cos(\frac{k\pi}{N})) \cos(\frac{\pi lk}{N}), \\\\ \gamma\_{i} & \equiv \frac{\pi}{N} \sum\_{k=0}^{N} \text{Q} (\cos(\frac{k\pi}{N})) \cos(\frac{\pi ik}{N}), \\\\ \lambda\_{i} & \equiv \frac{\pi}{N} \sum\_{k=0}^{N} \text{R} (\cos(\frac{k\pi}{N})) \cos(\frac{\pi ik}{N}). \end{aligned} \tag{9}$$

Now, substituting (4) and (9) in equations (6), and using the fact that

$$\mathbf{y}'(\mathbf{x}) \equiv \sum\_{m=0}^{N} a\_m^{(1)} T\_m \{ \mathbf{x} \} \ , \ a\_m^{(1)} = \frac{2}{c\_m} \sum\_{p=m+1}^{N} p a\_p \ , \ m = 0, 1, \ldots, N-1, \ a\_N^{(1)} = 0,$$

$$\mathbf{m} + \mathbf{p} = \text{odd}$$

$$\mathbf{y}''(\mathbf{x}) \approx \sum a\_m^{(2)} T\_m \{ \mathbf{x} \} \ , \ a\_m^{(2)} = \frac{1}{c\_m} \sum\_{p=m+2}^{N} p (p^2 - m^2) a\_p \ , \ m = 0, 1, \ldots, N-2 \ , \ a\_{N-1}^{(2)} = a\_N^{(2)} = \mathbf{0} \ ,$$

$$\mathbf{m} + \mathbf{p} = \text{even}$$

we get

$$\sum\_{l=0}^{N} \sum\_{m=0}^{N} \xi\_l a\_m^{(2)} T\_l \{ \mathbf{x} \} T\_m \{ \mathbf{x} \} + \sum\_{l=0}^{N} \sum\_{m=0}^{N} \gamma\_l a\_m^{(1)} T\_l \{ \mathbf{x} \} T\_m \{ \mathbf{x} \} + \sum\_{l=0}^{N} \sum\_{m=0}^{N} \lambda\_l a\_m T\_l \{ \mathbf{x} \} T\_m \{ \mathbf{x} \} = \mathbf{S} \{ \mathbf{x} \}, \tag{10}$$

$$\begin{aligned} \sum\_{i=0}^{N} a\_i T\_i (-1) &= \alpha \ , \\ \sum\_{i=0}^{N} a\_i T\_i (1) &= \beta \ . \end{aligned} \tag{11}$$

Numerical Solution of Linear Ordinary Differential

As we expected when *N* increases, errors decrease. *Example 2*. Consider Legendre's equation given by

�(�) = 1 − 3��. The results for *N=*4, 6, 10 were:

We end this section by solving Laguerre's equation.

Suppose �=2 and boundary conditions are given by �(−1) <sup>=</sup> �

Numer. Anal, Vol.1, pp. 193-213, 1981.

Springer- Verlag,NewYork,1988.

University Press, 1985.

SIAM,Philadelphia,1982.

exact solution is a polynomial.

The exact solution is (�) = 1 − 2� � ��

*Example 3*. Consider

respectively.

numerical result.

**5. References** 

function was singular.

Equations in Quantum Chemistry by Spectral Method 9

(1−��)��� − 2��� � �(��1)� = 0. As we know, this equation for �=2, and boundary conditions �(±1) = −2 has solution

5.5511 × 10���, 2.2204 × 10���, 2.7756 × 10���. Since our solution is a polynomial then for ��3, we come to a solution with error very closed to zero. If such cases you find the error is not zero but closed to it, is because of rounding error. We must put in our mind that the results by this method will be good if the

���� � (1 − �)�� � �� = 0.

Here we have again a polynomial solution, so we expect a solution with very small error. We examined for different values of *N* such as *N*=2, 3 and get the results 0 and 3 × 10���,

Results in these examples show the efficiency of Clenshaw method for obtaining a good

In case of singularity, one can use pseudo-spectral method. Some papers also modified pseudo-spectral method and overcome the problem of singularity even if the solution

Babolian. E, Bromilow. T. M, England. R, Saravi. M, '*A modification of pseudo-spectral method* 

Babolian. E, Delves. L .M, *A fast Galerkin scheme for linear integro-differential equations*, IMAJ.

Canuto. C, Hussaini. M. Y, Quarteroni. A, Zang. T. A, *Spectral Methods in Fluid Dynamics,*

Delves. L. M, Mohamed. J. L, *Computational methods for integral equations*, Cambridge

Gottlieb. D, Orszag. S. A, *Numerical Analysis of Spectral Methods, Theory and Applications*,

*for olving linear ODEs with singularity',* AMC 188 (2007) 1260-1266.

2� .

��,�(1) = − �

� .

Now, we multiply both sides of (10) by <sup>2</sup> 2 ( ) 1 *j j T x c x* , and integrate from -1 to 1, we obtain

$$\frac{2}{\pi\pi}\sum\_{j}^{N}\sum\_{i=0}^{N}\left[\xi\_{i}a\_{m}^{(2)} + \gamma\_{i}a\_{m}^{(1)} + \lambda\_{i}a\_{m}\right]\int\_{-1}^{1}\frac{T\_{i}\{\mathbf{x}\}T\_{m}\{\mathbf{x}\}T\_{j}\{\mathbf{x}\}}{\sqrt{1-\mathbf{x}^{2}}}d\mathbf{x}$$

$$=\frac{2}{\pi c\_{j}}\int\_{-1}^{1}\frac{S(\mathbf{x})T\_{j}\{\mathbf{x}\}}{\sqrt{1-\mathbf{x}^{2}}}d\mathbf{x}, j=0,1,...,N-2,\tag{12}$$

where,

$$\int \int\_{-1}^{1} \frac{T\_i(\mathbf{x}) T\_m(\mathbf{x}) T\_j(\mathbf{x})}{\sqrt{1 - \mathbf{x}^2}} d\mathbf{x} = \begin{cases} \pi & \text{, } i = m = j = 0 \\\\ \frac{\pi}{2} \delta\_{i, m} & \text{, } i + m > 0 \text{, } j = 0 \\\\ \frac{\pi}{4} (\delta\_{j, i + m} + \delta\_{j, \left| i - m \right|}) & \text{, } j > 0 \text{ , } \end{cases} \tag{13}$$

with, , <sup>1</sup> *i j* ,when *i j* , and zero when *i j* .

We can also compute the integrals in the right-hand side of (12) by the method of numerical integration using *N* 1 -point Gauss-Chebyshev-Lobatto quadrature. Therefore, substituting (13) in (12) and using the fact that ( 1) ( 1) , *<sup>i</sup> Ti* equations (12) and (11) make a system of *N* 1 equations for *N* 1 unknowns, 0 1 , ,..., *<sup>N</sup> aa a* , hence we can find 0 1 ( , ,..., )*<sup>t</sup> <sup>N</sup> aa a* from this system.

#### **4. Numerical examples**

As we mentioned the important problem in the field of atomic and molecular structure, is solution of wave equation for hydrogen atom. In this section we will solve Schrodinger, Legendre and Laguerre equations, which occur in the hydrogen atom wave equations, by Clenshaw method and observe the power of this method comparing with usual numerical methods such as Euler's or Runge-Kutta's methods. We start with Schrodinger's equation.

*Example 1*. Let us consider

$$
\varphi'' + (2mEh^{-2} - a^2 \mathfrak{x}^2)\varphi = 0.
$$

Assume � = 2, ����� = −1, with �(0) = 1, �(1) = �. The exact solution is �(�) = ���� .

Here interval is chosen as [0,1], but using change of variable such as � = ��� � we can transfer interval [0,1] to [-1,1].

We solve this equation by Clenshaw method and compare the results for different values of *N*. The results for *N*=4, 7, 10, 13, respectively, were:

$$1.660 \times 10^{-2}, 4.469 \times 10^{-5}, 5.901 \times 10^{-8}, 7.730 \times 10^{-11}.$$

As we expected when *N* increases, errors decrease.

*Example 2*. Consider Legendre's equation given by

$$(1 - x^2)\mathbf{y}'' - 2x\mathbf{y}' + \lambda(\lambda + 1)\mathbf{y} = 0.1$$

As we know, this equation for �=2, and boundary conditions �(±1) = −2 has solution �(�) = 1 − 3��. The results for *N=*4, 6, 10 were:

$$15.5511 \times 10^{-17}, 2.2204 \times 10^{-16}, 2.7756 \times 10^{-17}.$$

Since our solution is a polynomial then for ��3, we come to a solution with error very closed to zero. If such cases you find the error is not zero but closed to it, is because of rounding error. We must put in our mind that the results by this method will be good if the exact solution is a polynomial.

We end this section by solving Laguerre's equation.

*Example 3*. Consider

8 Quantum Chemistry – Molecules for Innovations

*T x c x*

, and integrate from -1 to 1, we obtain

(12)

(13)

.

we can transfer

�

1

*imj*

*j*

, 0,

*im j*

*imj*

*T xT xT x*

 <sup>2</sup> 2 ( ) 1 *j*

> 

 

,

*i m*

 

We can also compute the integrals in the right-hand side of (12) by the method of numerical integration using *N* 1 -point Gauss-Chebyshev-Lobatto quadrature. Therefore, substituting (13) in (12) and using the fact that ( 1) ( 1) , *<sup>i</sup> Ti* equations (12) and (11) make a system of *N* 1 equations for *N* 1 unknowns, 0 1 , ,..., *<sup>N</sup> aa a* , hence we can find

As we mentioned the important problem in the field of atomic and molecular structure, is solution of wave equation for hydrogen atom. In this section we will solve Schrodinger, Legendre and Laguerre equations, which occur in the hydrogen atom wave equations, by Clenshaw method and observe the power of this method comparing with usual numerical methods such as Euler's or Runge-Kutta's methods. We start with Schrodinger's equation.

��� + (2����� − ����)� = 0.

We solve this equation by Clenshaw method and compare the results for different values of

1.660 × 10��, 4.469 × 10��, 5.901 × 10��, 7.730 × 10���.

Assume � = 2, ����� = −1, with �(0) = 1, �(1) = �. The exact solution is �(�) = ����

Here interval is chosen as [0,1], but using change of variable such as � = ���

*N*. The results for *N*=4, 7, 10, 13, respectively, were:

<sup>2</sup> () () , 0,1,..., 2,

*aaa dx*

*dx j N*

 

*ji m jim*

( ), 0 , <sup>4</sup>

, ,

, 0, 0, <sup>2</sup>

*j*

*c x*

*j*

*S xT x*

*im im im*

 1 (2) (1) <sup>1</sup> <sup>2</sup> 0 0 <sup>2</sup> () () () [ ]

> 1

> > =

*c x*

1 2

1

 

Now, we multiply both sides of (10) by

 1

1 2

1

, <sup>1</sup> *i j* ,when *i j* , and zero when *i j* .

() () ()

*T xT xT x imj dx x*

where,

with, 

0 1 ( , ,..., )*<sup>t</sup>*

*<sup>N</sup> aa a* from this system.

**4. Numerical examples** 

*Example 1*. Let us consider

interval [0,1] to [-1,1].

*N N*

*j*

*ji m*

$$x\mathcal{y}'' + (1-x)y' + \lambda y = 0.$$

Suppose �=2 and boundary conditions are given by �(−1) <sup>=</sup> � ��,�(1) = − � � .

The exact solution is (�) = 1 − 2� � �� 2� .

Here we have again a polynomial solution, so we expect a solution with very small error. We examined for different values of *N* such as *N*=2, 3 and get the results 0 and 3 × 10���, respectively.

Results in these examples show the efficiency of Clenshaw method for obtaining a good numerical result.

In case of singularity, one can use pseudo-spectral method. Some papers also modified pseudo-spectral method and overcome the problem of singularity even if the solution function was singular.

#### **5. References**


**0**

**2**

*Brazil*

**at CCSD(T) Level**

Nelson Henrique Morgon *Universidade Estadual de Campinas*

**Composite Method Employing Pseudopotential**

Thermochemical data are among the most fundamental and useful information of chemical species which can be used to predict chemical reactivity and relative stability. Thus, it is not surprising that an important goal of computational chemistry is to predict thermochemical parameters with reasonable accuracy (Morgon, 1995a). Reliability is a critical feature of any theoretical model, and for practical purposes the model should be efficient in order to be widely applicable in estimating the structure, energy and other properties of systems, as isolated ions, atoms, molecules(Ochterski et al., 1995), or gas phase reactions(Morgon, 2008a).

For instance gas phase reactions between molecules and ions, and molecules and electrons are known to be important in many scientifically and technologically environments. On the cosmic scale, the chemistry that produces molecules in interstellar clouds is dominated by ion-molecule reactions. Shrinking down to our own planet, the upper atmosphere is a plasma, and contains electrons and various positive ions. Certain anthropogenic chemical compounds (including SF6 and perfluorocarbons) can probably not be destroyed within the troposhere or stratosphere, but may be removed by reactions with ions or electrons in the ionosphere. Recent years have seen a massive growth in the industrial use of plasmas, particularly in the fabrication of microelectronic devices and components. The chemistry within the plasma, much of which involves ion-molecule and electron-molecule reactions, determines the species that etch the surface, and hence the outcome and rate of an etching process. Much of the chemistry that is often labelled as 'organic' or 'inorganic' involves ion-molecule reactions, usually carried out in the presence of a solvent. For instance, S*N*2 reactions, such as OH<sup>−</sup> + CH3Cl, fall into this category. To gain a clearer picture of how these (gas phase) reactions occur, it is advantageous to study them removed from the (very great) perturbations due to

So, the need for thorough studies of ion-molecule and electron-molecule reactions are thus well established, ranging from the astrophysical origins of molecules, through the survival of the earth's atmosphere, to modelling the plasmas that underpin many advanced processing technologies. There is intrinsic interest too in the studies, as they help to explore the nature and progress of binary encounters between molecules and ions, and molecules and electrons. At the most basic level answers are needed to the following questions - how fast does a reaction proceed? and what are the products of the reaction? What determines which

**1. Introduction**

the solvent.

What is the importance of these studies?


## **Composite Method Employing Pseudopotential at CCSD(T) Level**

Nelson Henrique Morgon *Universidade Estadual de Campinas Brazil*

#### **1. Introduction**

10 Quantum Chemistry – Molecules for Innovations

Lanczos. C, *Trigonometric interpolation of empirical and analytical functions,* J. Math. Phys. 17

Levine. Ira N, *Quantum Chemistry*, 5th ed, City University of NewYork, Prentice-Hall

Pauling. L, Wilson. E.B, *Quantum Mechanic*, McGraw-Hill Book Company, 1981.

(1938) 123-129.

Publication, 2000.

Thermochemical data are among the most fundamental and useful information of chemical species which can be used to predict chemical reactivity and relative stability. Thus, it is not surprising that an important goal of computational chemistry is to predict thermochemical parameters with reasonable accuracy (Morgon, 1995a). Reliability is a critical feature of any theoretical model, and for practical purposes the model should be efficient in order to be widely applicable in estimating the structure, energy and other properties of systems, as isolated ions, atoms, molecules(Ochterski et al., 1995), or gas phase reactions(Morgon, 2008a).

What is the importance of these studies?

For instance gas phase reactions between molecules and ions, and molecules and electrons are known to be important in many scientifically and technologically environments. On the cosmic scale, the chemistry that produces molecules in interstellar clouds is dominated by ion-molecule reactions. Shrinking down to our own planet, the upper atmosphere is a plasma, and contains electrons and various positive ions. Certain anthropogenic chemical compounds (including SF6 and perfluorocarbons) can probably not be destroyed within the troposhere or stratosphere, but may be removed by reactions with ions or electrons in the ionosphere. Recent years have seen a massive growth in the industrial use of plasmas, particularly in the fabrication of microelectronic devices and components. The chemistry within the plasma, much of which involves ion-molecule and electron-molecule reactions, determines the species that etch the surface, and hence the outcome and rate of an etching process. Much of the chemistry that is often labelled as 'organic' or 'inorganic' involves ion-molecule reactions, usually carried out in the presence of a solvent. For instance, S*N*2 reactions, such as OH<sup>−</sup> + CH3Cl, fall into this category. To gain a clearer picture of how these (gas phase) reactions occur, it is advantageous to study them removed from the (very great) perturbations due to the solvent.

So, the need for thorough studies of ion-molecule and electron-molecule reactions are thus well established, ranging from the astrophysical origins of molecules, through the survival of the earth's atmosphere, to modelling the plasmas that underpin many advanced processing technologies. There is intrinsic interest too in the studies, as they help to explore the nature and progress of binary encounters between molecules and ions, and molecules and electrons. At the most basic level answers are needed to the following questions - how fast does a reaction proceed? and what are the products of the reaction? What determines which

of anions. Gaussian-n theories (G1, G2, G3, and G4) (Curtiss et al., 1997; 1998; 2000; 2007) have given good results for properties like proton and electron affinities, enthalpies of formation, atomization energies, and ionization potentials. These theories are a composite technique in which a sequence of well-defined ab initio molecular orbital calculations is performed to arrive at a total energy of a given molecular species. There are other techniques that have been demonstrated to predict accurate thermochemical properties of chemical species, and are alternative to the Gaussian-n methods: the Correlation Consistent Composite Approach (ccCA)(DeYonker et al., 2006), the Multireference Correlation Consistent Composite Approach (MR-ccCA)(Oyedepo & Wilson, 2010), the Complete Basis Set Methods (CBS) and its versions: CBS-4M, CBS-Lq, CBS-Q, CBS-QB3, CBS-APNO(Montgomery Jr. et al., 2000; Nyden & Petersson, 1981; Ochterski et al., 1996; Peterson et al., 1991), and Weizmann Theories (W1 to W4)(Boese et al., 2004; Karton et al., 2006; Martin & De Oliveira, 1999; Parthiban & Martin,

Composite Method Employing Pseudopotential at CCSD(T) Level 13

Recently, we have implemented and tested a pseudopotential to be used with the G3 theory for molecules containing first-, second-, and non-transition third-row atoms (G3CEP) (Pereira et al., 2011). The final average total absolute deviation using this methodology and the all-electron G3 were 5.39 kJ mol−<sup>1</sup> and 4.85 kJ mol−1, respectively. Depending on the size of the molecules and the type of atoms considered, the CPU time was drastically decreased.

In this chapter we have developed a computational model similar to version of the G2(MP2, SVP) theory (Curtiss et al., 1996). Both theories are based on the additivity approximations to estimate the high level energy for the extended function basis set. While G2(MP2,SVP) is based on the additivity approximation to estimate the QCISD(T) energy for the extended 6-311+G (3df,2p) basis set: E[QCISD(T)/ 6-311+G(3df,2p)] E[QCISD(T)/6-31G(d)] + E[MP2/6-311+G(3df,2p)] - E[MP2/6-31G(d)]. Our methodology employs CCSD(T) energies in addition to the the valence basis sets adapted for pseudopotential (ECP) (Stevens et al., 1984) using the Generator Coordinate Method (GCM) procedure (Mohallem & Dreizler, 1986;

The present methodology which relies on small basis sets (representation of the core electrons by ECP) and an easier and simpler way for correcting the valence region (mainly of anionic systems) appears as an interesting alternative for the calculation of thermochemical data such

The GCM has been very useful in the study of basis sets(Morgon, 1995a;b; 2006; 2008b; 2011; Morgon et al., 1997). It considers the monoelectronic functions *ψ*(1) as an integral transform,

where *f*(*α*) and *φ*(*α*, 1) are the weight and generator functions respectively (gaussian functions are used in this work), and *α* is the generator. The existence of the weight functions (graphical display of the linear combination of basis functions) is an essential condition for the use of GCM. Analysis of the behavior of the weight functions by the GCM permits the

*f*(*α*) *φ*(*α*, 1) *dα* (1)

as electron and proton affinities or heat of formation for larger systems.

*<sup>ψ</sup>*(1) = <sup>∞</sup>

0

2001).

**2. Computational methods**

Mohallem & Trsic, 1985).

**2.1 Development of basis sets**

reactions occur? and what products are formed? Beyond these may come questions about the detailed dynamics of the reaction, such as how changing the energy of the reactants may influence the progress and outcome of the reaction.

Many powerful experimental techniques have been developed to give the basic data of reaction rate coefficients and products (usually just the identification of the ion product). These results are part of the raw data needed to understand and model the complex chemistry occurring in the diverse environments identified above. There is much information that is not directly available from the experimental data. This includes identification of the neutral products of a reaction, knowledge of the thermochemistry of the reaction, and characterization of the pathway that connects reactants to products. By invoking some general rules, the experimental observations can be used to provide partial answers. Thus the fast flow techniques that are used to provide much of the experimental data on ion-molecule reactions can only detect the occurrence of very rapid reactions (k >= 10−<sup>12</sup> cm3 s−1), which places an upper bound on the exothermicity of the reaction of +20 kJ mol−1. In cases where there are existing reliable enthalpies of formation of each of the species in a proposed pathway to an observed ion product, this rule can test whether the suggested neutral products may be correct. In other cases, where the enthalpy of formation of just one of the species involved in a reaction (usually the product ion) is unknown, the observation of a specific reaction pathway can be used to place a bound on the previously unkown enthalpy of formation. Finally for reactions which are known to be exothermic, if the experimenal rate coefficient is observed to be less than the capture theory rate coefficient, then it is usual to conclude that there must be some bottleneck or barrier to the reaction.

What can theoretical calculations add to the experimental data?

Three important and fundamental gas-phase thermochemical properties from a theoretical and experimental point of view are the standard heat of formation (Δ*<sup>f</sup> H<sup>o</sup> gas*), the electron (EA) and proton (PA) affinities. Thus, it is not surprising that an important goal of computational chemistry is to predict such thermochemical parameters with reasonable accuracy, which can be useful in the gas phase reaction studies. Proton transfer reactions are also of great importance in chemistry and in biomolecular processes of living organisms(Ervin, 2001). Absolute values of proton affinities are not always easy to obtain and are often derived from relative measurements with respect to reference molecules. Relative proton affinities are usually measured by means of high pressure mass spectrometry, with triple quadrupole and ion trap mass spectrometers (Mezzache et al., 2005) or using ion mobility spectrometry (Tabrizchi & Shooshtari, 2003). The importance and utility of the EA extend well beyond the regime of gas-phase ion chemistry. A survey of examples illustrates the diversity of areas in which electron affinities play a role: silicon, germanium clusters, interstellar chemistry, microelectronics, and so on.

The standard heat of formation, which measures the thermodynamic stability, is useful in the interpretation of the mechanisms of chemical reactions (Badenes et al., 2000).

On the other hand, theoretical calculations represent one attempt to study absolute values of electron or proton affinity and other thermochemical properties (Smith & Radom, 1991). However, accurate calculations of these properties require sophisticated and high level methods, and great amount of computational resources. This is particularly true for atoms of the 2nd, 3rd, ..., periods and for calculating properties like the proton and electron affinity of anions. Gaussian-n theories (G1, G2, G3, and G4) (Curtiss et al., 1997; 1998; 2000; 2007) have given good results for properties like proton and electron affinities, enthalpies of formation, atomization energies, and ionization potentials. These theories are a composite technique in which a sequence of well-defined ab initio molecular orbital calculations is performed to arrive at a total energy of a given molecular species. There are other techniques that have been demonstrated to predict accurate thermochemical properties of chemical species, and are alternative to the Gaussian-n methods: the Correlation Consistent Composite Approach (ccCA)(DeYonker et al., 2006), the Multireference Correlation Consistent Composite Approach (MR-ccCA)(Oyedepo & Wilson, 2010), the Complete Basis Set Methods (CBS) and its versions: CBS-4M, CBS-Lq, CBS-Q, CBS-QB3, CBS-APNO(Montgomery Jr. et al., 2000; Nyden & Petersson, 1981; Ochterski et al., 1996; Peterson et al., 1991), and Weizmann Theories (W1 to W4)(Boese et al., 2004; Karton et al., 2006; Martin & De Oliveira, 1999; Parthiban & Martin, 2001).

Recently, we have implemented and tested a pseudopotential to be used with the G3 theory for molecules containing first-, second-, and non-transition third-row atoms (G3CEP) (Pereira et al., 2011). The final average total absolute deviation using this methodology and the all-electron G3 were 5.39 kJ mol−<sup>1</sup> and 4.85 kJ mol−1, respectively. Depending on the size of the molecules and the type of atoms considered, the CPU time was drastically decreased.

#### **2. Computational methods**

2 Will-be-set-by-IN-TECH

reactions occur? and what products are formed? Beyond these may come questions about the detailed dynamics of the reaction, such as how changing the energy of the reactants may

Many powerful experimental techniques have been developed to give the basic data of reaction rate coefficients and products (usually just the identification of the ion product). These results are part of the raw data needed to understand and model the complex chemistry occurring in the diverse environments identified above. There is much information that is not directly available from the experimental data. This includes identification of the neutral products of a reaction, knowledge of the thermochemistry of the reaction, and characterization of the pathway that connects reactants to products. By invoking some general rules, the experimental observations can be used to provide partial answers. Thus the fast flow techniques that are used to provide much of the experimental data on ion-molecule reactions can only detect the occurrence of very rapid reactions (k >= 10−<sup>12</sup> cm3 s−1), which places an upper bound on the exothermicity of the reaction of +20 kJ mol−1. In cases where there are existing reliable enthalpies of formation of each of the species in a proposed pathway to an observed ion product, this rule can test whether the suggested neutral products may be correct. In other cases, where the enthalpy of formation of just one of the species involved in a reaction (usually the product ion) is unknown, the observation of a specific reaction pathway can be used to place a bound on the previously unkown enthalpy of formation. Finally for reactions which are known to be exothermic, if the experimenal rate coefficient is observed to be less than the capture theory rate coefficient, then it is usual to conclude that there must be

Three important and fundamental gas-phase thermochemical properties from a theoretical

and proton (PA) affinities. Thus, it is not surprising that an important goal of computational chemistry is to predict such thermochemical parameters with reasonable accuracy, which can be useful in the gas phase reaction studies. Proton transfer reactions are also of great importance in chemistry and in biomolecular processes of living organisms(Ervin, 2001). Absolute values of proton affinities are not always easy to obtain and are often derived from relative measurements with respect to reference molecules. Relative proton affinities are usually measured by means of high pressure mass spectrometry, with triple quadrupole and ion trap mass spectrometers (Mezzache et al., 2005) or using ion mobility spectrometry (Tabrizchi & Shooshtari, 2003). The importance and utility of the EA extend well beyond the regime of gas-phase ion chemistry. A survey of examples illustrates the diversity of areas in which electron affinities play a role: silicon, germanium clusters, interstellar chemistry,

The standard heat of formation, which measures the thermodynamic stability, is useful in the

On the other hand, theoretical calculations represent one attempt to study absolute values of electron or proton affinity and other thermochemical properties (Smith & Radom, 1991). However, accurate calculations of these properties require sophisticated and high level methods, and great amount of computational resources. This is particularly true for atoms of the 2nd, 3rd, ..., periods and for calculating properties like the proton and electron affinity

interpretation of the mechanisms of chemical reactions (Badenes et al., 2000).

*gas*), the electron (EA)

influence the progress and outcome of the reaction.

some bottleneck or barrier to the reaction.

microelectronics, and so on.

What can theoretical calculations add to the experimental data?

and experimental point of view are the standard heat of formation (Δ*<sup>f</sup> H<sup>o</sup>*

In this chapter we have developed a computational model similar to version of the G2(MP2, SVP) theory (Curtiss et al., 1996). Both theories are based on the additivity approximations to estimate the high level energy for the extended function basis set. While G2(MP2,SVP) is based on the additivity approximation to estimate the QCISD(T) energy for the extended 6-311+G (3df,2p) basis set: E[QCISD(T)/ 6-311+G(3df,2p)] E[QCISD(T)/6-31G(d)] + E[MP2/6-311+G(3df,2p)] - E[MP2/6-31G(d)]. Our methodology employs CCSD(T) energies in addition to the the valence basis sets adapted for pseudopotential (ECP) (Stevens et al., 1984) using the Generator Coordinate Method (GCM) procedure (Mohallem & Dreizler, 1986; Mohallem & Trsic, 1985).

The present methodology which relies on small basis sets (representation of the core electrons by ECP) and an easier and simpler way for correcting the valence region (mainly of anionic systems) appears as an interesting alternative for the calculation of thermochemical data such as electron and proton affinities or heat of formation for larger systems.

#### **2.1 Development of basis sets**

The GCM has been very useful in the study of basis sets(Morgon, 1995a;b; 2006; 2008b; 2011; Morgon et al., 1997). It considers the monoelectronic functions *ψ*(1) as an integral transform,

$$
\psi(1) = \int\_0^\infty f(\alpha) \,\phi(\alpha, 1) \,d\alpha \tag{1}
$$

where *f*(*α*) and *φ*(*α*, 1) are the weight and generator functions respectively (gaussian functions are used in this work), and *α* is the generator. The existence of the weight functions (graphical display of the linear combination of basis functions) is an essential condition for the use of GCM. Analysis of the behavior of the weight functions by the GCM permits the

**2.2 Molecular calculations**

represents a significant challenge.

electrons.

In many problems to be addressed by electronic structure methodology, high accuracy is of crucial importance. In order to obtain energies that may approach chemical accuracy (≈10 kJ.mol−1) calculations must take account of electron correlation. The 'ideal' methodology would have been a multi-reference configuration interaction all electron calculation with several large, flexible basis sets to enable extrapolation to the complete basis set limit. Depend on the size of the systems, performing accurate calculations (methodology and basis sets)

Composite Method Employing Pseudopotential at CCSD(T) Level 15

Morgon *et al.* (Morgon, 1998; Morgon et al., 1997; Morgon & Riveros, 1998) have been developing techniques to tackle such problems. These are centered around the use of effective core potentials, in which the inner electrons are represented by an effective potential derived from calculations on atoms. The electronic wavefunction itself then only contains the outer

**(a)** optimization of the molecular geometries and vibrational analysis are carried out at HF/B0 level. The harmonic frequencies confirm that the stationary points correspond to minima

**(c)** at the MP2 equilibrium geometry corrections to the total energies are performed at higher level of theory. First, this is carried out at CR-CCSD[T]/B0 level (Completely Renormalized Coupled-Cluster with Single and Double and Perturbative Triple excitation) (Kowalski & Piecuch, 2000) (or at CCSD(T)/B0 level for EA calculations), and later by addition of extra

Thus, these results coupled to additive approximations for the energy yield an effective

≈ *E*[*CR* − *CCSD*[*T*]/(*B*1)] = *E*[*CR* − *CCSD*[*T*]/(*B*0)] +

The CR-CCSD[T] (Completely Renormalized Coupled-Cluster with Single and Double and Perturbative Triple excitation) methodology refers to size-extensive left eigenstate completely renormalized (CR) coupled-cluster (CC) singles (S), doubles (D), and noniterative triples (T). This approach is abbreviated as CR-CCL and is appropriately described by Piecuch(Piecuch

An alternative model was developed for the study of heat of formation. This model employs valence basis sets aug-CCpVnZ (n = 2, 3, and 4) (Dunning, 1989). These basis sets were adapted to ECP using the GCM and are identified by ECP+ACCpVnZm (m = modified). The energies are obtained through the extrapolations to the complete basis set limit (CBS) using Peterson mixed exponential/Gaussian function extrapolation scheme (Feller & Peterson,

where x = 2, 3, and 4 come from ECP+ACCpV2Zm, ECP+ACCpV3Zm and ECP+ACCpV4Zm

where *scal* (0.89) is the scaling factor on the vibrational frequencies.

+ *E*[*MP*2/(*B*1)] − *E*[*MP*2/(*B*0)] + *ZPE*[*HF*/(*B*0)] ∗ *scal* (3)

*<sup>E</sup>*(*MP*2) = *ECBS* <sup>+</sup> *B exp*[−(*<sup>x</sup>* <sup>−</sup> <sup>1</sup>)] + *C exp*[−(*<sup>x</sup>* <sup>−</sup> <sup>1</sup>)2] (4)

The procedure to the molecular calculations employing this methodology is:

and are used to compute the zero-point energies; **(b)** further optimization is carried out at MP2/B0 level;

functions (*s, p, d,* and *f*) at MP2/B1 level.

calculation at a high level of theory,

et al., 2002) and Ge(Ge et al., 2007).

1999).

energies, respectively.

atomic basis set to be adapted in such a way as to yield a better description of the core electrons (represented by ECP) and the valence orbitals (corrected by addition of the extra diffuse functions), in the molecular environment. With the exception of some simple systems the analytical expression of the weight functions is unknown. Thus, an analytical solution of the integral transform (Eq. 1) is not viable in most cases, and suggests the need of numerical techniques to solve Eq. 1 (Custodio, Giordan, Morgon & Goddard, 1992; Custodio, Goddard, Giordan & Morgon, 1992). The solution can be carried out by an appropriate choice of discrete points on the generate coordinate, represented by:

$$\mu\_{i,(k)} = \exp[\Omega\_{o,(k)} + (i - 1) \cdot \Delta \Omega\_{(k)}]\_\prime \quad i = 1, 2, 3, \dots \, N\_{(k)} \tag{2}$$

The discretization of the set is defined by the following parameters: an initial value (Ω*o*), an increment (ΔΩ), and by the number of primitives used (N) for a given orbital *k* (*s, p, d,* ... ). The search for the best representation is obtained using the total energy of the electronic ground state as the minimization criterion.

The SIMPLEX search method (Nelder & Mead, 1965) can be adapted to the any electronic structure program to provide the minimum energy of the ground state of the atom corresponding to the optimized discretization parameters.

The basic procedure consists of the following steps:


To the heavy atoms *f* type polarization functions are not available in these valence basis sets for this kind of ECP(Stevens et al., 1984). So, it was need to define the value of these *f* functions for Br and I atoms. The determination of the best value was carried out considering the smaller difference between the PAs (experimental and theoretical) values considering the Br− and I− anions. The *f* exponent values found are 0.7 and 0.3 for Br and I atoms, respectively.

In fact two sets of basis functions are used, a small basis and a larger basis, with extra diffuse and polarization functions, B0 and B1, respectivally. Calculations with basis B1 are naturally much more expensive than those employing basis B0, so it is important to have computational schemes that perform the minimum number of calculations using basis B1.

For instance, the B0 basis set is defined as: (31) for H; (311/311) for C, O, F, S, and Cl; and (411/411) Br and I. For more refined energy calculations (B1 basis set are used), this set was augmented with additional diffuse and polarization functions (*p* for H and *s*, *p*, *d*, and *f* for heavy atoms) to yield a (311/11) set for H; (311/311/11/1) for C, O, F, S, and Cl; and (411/411/11/1) for Br and I atoms.

#### **2.2 Molecular calculations**

4 Will-be-set-by-IN-TECH

atomic basis set to be adapted in such a way as to yield a better description of the core electrons (represented by ECP) and the valence orbitals (corrected by addition of the extra diffuse functions), in the molecular environment. With the exception of some simple systems the analytical expression of the weight functions is unknown. Thus, an analytical solution of the integral transform (Eq. 1) is not viable in most cases, and suggests the need of numerical techniques to solve Eq. 1 (Custodio, Giordan, Morgon & Goddard, 1992; Custodio, Goddard, Giordan & Morgon, 1992). The solution can be carried out by an appropriate choice of discrete

The discretization of the set is defined by the following parameters: an initial value (Ω*o*), an increment (ΔΩ), and by the number of primitives used (N) for a given orbital *k* (*s, p, d,* ... ). The search for the best representation is obtained using the total energy of the electronic

The SIMPLEX search method (Nelder & Mead, 1965) can be adapted to the any electronic structure program to provide the minimum energy of the ground state of the atom

**(a)** search of the optimum discretized parameter set for the atoms using the GCM for variation on the generator coordinate space. The core electrons are represented by a pseudopotential. The discretization parameters (Ω*<sup>o</sup>* and ΔΩ) are defined with conjunction

**(b)** the minimum energy criterion is observed and the characteristics of the atomic orbital

**(c)** extra functions (polarizaton or diffuse funcions - *s, p, d* and *f*) are obtained by observing the convergent behavior of the weight functions of the outer atomic orbitals (*s* and *p*). These extra functions are needed for the correct description of the electronic distribution in an

To the heavy atoms *f* type polarization functions are not available in these valence basis sets for this kind of ECP(Stevens et al., 1984). So, it was need to define the value of these *f* functions for Br and I atoms. The determination of the best value was carried out considering the smaller difference between the PAs (experimental and theoretical) values considering the Br− and I−

In fact two sets of basis functions are used, a small basis and a larger basis, with extra diffuse and polarization functions, B0 and B1, respectivally. Calculations with basis B1 are naturally much more expensive than those employing basis B0, so it is important to have computational

For instance, the B0 basis set is defined as: (31) for H; (311/311) for C, O, F, S, and Cl; and (411/411) Br and I. For more refined energy calculations (B1 basis set are used), this set was augmented with additional diffuse and polarization functions (*p* for H and *s*, *p*, *d*, and *f* for heavy atoms) to yield a (311/11) set for H; (311/311/11/1) for C, O, F, S, and Cl; and

anions. The *f* exponent values found are 0.7 and 0.3 for Br and I atoms, respectively.

schemes that perform the minimum number of calculations using basis B1.

*αi*,(*k*) = *exp*[Ω*o*,(*k*) + (*i* − 1) · ΔΩ(*k*)], *i* = 1, 2, 3, . . . *N*(*k*) (2)

points on the generate coordinate, represented by:

ground state as the minimization criterion.

with this ECP;

weight functions are analyzed;

(411/411/11/1) for Br and I atoms.

corresponding to the optimized discretization parameters.

The basic procedure consists of the following steps:

anion (diffuse character of electronic cloud).

In many problems to be addressed by electronic structure methodology, high accuracy is of crucial importance. In order to obtain energies that may approach chemical accuracy (≈10 kJ.mol−1) calculations must take account of electron correlation. The 'ideal' methodology would have been a multi-reference configuration interaction all electron calculation with several large, flexible basis sets to enable extrapolation to the complete basis set limit. Depend on the size of the systems, performing accurate calculations (methodology and basis sets) represents a significant challenge.

Morgon *et al.* (Morgon, 1998; Morgon et al., 1997; Morgon & Riveros, 1998) have been developing techniques to tackle such problems. These are centered around the use of effective core potentials, in which the inner electrons are represented by an effective potential derived from calculations on atoms. The electronic wavefunction itself then only contains the outer electrons.

The procedure to the molecular calculations employing this methodology is:


Thus, these results coupled to additive approximations for the energy yield an effective calculation at a high level of theory,

$$\approx E[CR-CCSD[T]/(B1)] = E[CR-CCSD[T]/(B0)] + $$

$$+ E[MP2/(B1)] - E[MP2/(B0)] + 2PE[HF/(B0)] \approx \text{scal} \tag{3}$$

where *scal* (0.89) is the scaling factor on the vibrational frequencies.

The CR-CCSD[T] (Completely Renormalized Coupled-Cluster with Single and Double and Perturbative Triple excitation) methodology refers to size-extensive left eigenstate completely renormalized (CR) coupled-cluster (CC) singles (S), doubles (D), and noniterative triples (T). This approach is abbreviated as CR-CCL and is appropriately described by Piecuch(Piecuch et al., 2002) and Ge(Ge et al., 2007).

An alternative model was developed for the study of heat of formation. This model employs valence basis sets aug-CCpVnZ (n = 2, 3, and 4) (Dunning, 1989). These basis sets were adapted to ECP using the GCM and are identified by ECP+ACCpVnZm (m = modified). The energies are obtained through the extrapolations to the complete basis set limit (CBS) using Peterson mixed exponential/Gaussian function extrapolation scheme (Feller & Peterson, 1999).

$$E(MP2) = E\_{CBS} + B \exp[-(\mathbf{x} - \mathbf{1})] + C \exp[-(\mathbf{x} - \mathbf{1})^2] \tag{4}$$

where x = 2, 3, and 4 come from ECP+ACCpV2Zm, ECP+ACCpV3Zm and ECP+ACCpV4Zm energies, respectively.

better description of the electrons in the molecular environment. The analysis of the weight functions of the outermost atomic orbitals suggested the need for improvements of the basis sets for the heavy atoms. Observing the plots of the weight functions of the outermost orbitals it is possible to establish the best fit of the basis sets. Using the atom of Cl as an example, the representation of the weight function of the atomic orbital 3s is shown in Fig. 1. The continuous line represents the plot of the weight function of the original primitive basis set for the all electron system. The dashed line represents the same weight function obtained using the pseudopotential. This figure shows that the 10 inner functions (large *α*) have a contribution close to zero towards the description of the weight function. The vertical solide line cutoff of the basis set indicates precisely where the pseudopotential starts to represent the

Composite Method Employing Pseudopotential at CCSD(T) Level 17


Fig. 1. Weight functions for the 3s atomic orbital of Cl in systems with all electrons (ae, continuous line), with the pseudopotential (pp, dashed line), with the addition of one diffuse function plus pseudopotential (dif+pp, dot-dashed line) and the cutoff point line represented

Proton affinity is a very sensitive property of the electronic structure and it is appropriate to test our methodology. Table 1 shows the results of the acidity (in kJ.mol−1) for a set of anions systems. These results were obtained using the B0 and B1 basis sets and Eq. 3. A comparison is also presented with experimental results. One can observe that our theoretical results are


by the vertical solid line.

**3.2 Proton affinity**

0.00

0.50

f(α)

1.00

core atomic region.

For this electronic property, molecular calculations consist of:


$$\begin{aligned} \approx E[\text{CR} - \text{CCL} / \text{ECP} + \text{ACC5Z}m] &= E[\text{CR} - \text{CCL} / \text{ECP} + \text{ACC2Z}m] + \\ &+ E[\text{MP2} / \text{ECP} + \text{ACC5Z}m] - E[\text{MP2} / \text{ECP} + \text{ACC2Z}m] - \\ &+ \text{ZPE}[\text{HF} / \text{ECP} + \text{ACC2Z}m] \* \text{scal} + E(\text{HLC}) \end{aligned} \tag{5}$$

where *scal* is the scaling factor on ZPE.

The method also includes an empirical higher-level correction (HLC) term. This term is given by either Eq. 6 or Eq. 7 depending on whether the species is a molecule or an atom:

$$HLC\_{mlec.} = -C \cdot n\_{\beta} - D \cdot (n\_{\alpha} - n\_{\beta}) \tag{6}$$

$$HLC\_{atom} = -A \cdot n\_{\beta} - B \cdot (n\_{\alpha} - n\_{\beta}) \tag{7}$$

In these equations, n*<sup>α</sup>* and n*<sup>β</sup>* are the numbers of *α* and *β* electrons, respectively. The parameters A (4.567mH), B (2.363mH), C (4.544mH), and D (2.337mH) were obtained by fitting to the experimental data of heats of formation.

It is important to note that since it only depends on the number of *α* and *β* valence electrons, the HLC cancels entirely from most reaction energies, except when the reactions involve a mixture of atoms and molecules (as in heats of formation and bond dissociation energies) and/or when spin is not conserved (Lin et al., 2009).

It should be also noted that the absolute values of the calculated energies have no real significance, as no common energy zero for different atoms has been used. This arises from how the effective core potentials are constructed. When differences are formed, the differences between the zeros cancel, to leave for instance the difference in energy between the products and reactants of a reaction.

Additionally, an alternative approach for molecular geometry optimization and harmonic frequencies calculation can be considered through the use of the DFT (B3LYP, M06, ...).

#### **3. Results and discussion**

#### **3.1 Basis sets**

The existence of the weight functions (graphic representations of the linear combination of the atomic orbitals) is the fundamental condition to use the GCM. The analysis of the behavior of the weight functions by the GCM allows the fitting of the atomic basis sets in order to get a 6 Will-be-set-by-IN-TECH

**(a)** optimization of the molecular geometries and vibrational analysis are carried out at HF/ECP+ACC2Zm level. The harmonic frequencies are employed to characterize the local

**(c)** at the MP2 equilibrium geometry corrections to the total energies are performed at higher level of theory. First, at CR-CCL/ECP+ACC2Zm level, and calculations by addition of extra functions (*s, p, d,* and *f*) at MP2/ECP+ACC3Zm and MP2/ECP+ACC4Zm levels. The

**(d)** Finally, the results are coupled through additive approximations, and the energy

≈ *E*[*CR* − *CCL*/*ECP* + *ACC*5*Zm*] = *E*[*CR* − *CCL*/*ECP* + *ACC*2*Zm*] +

+ *ZPE*[*HF*/*ECP* + *ACC*2*Zm*] ∗ *scal* + *E*(*HLC*)

The method also includes an empirical higher-level correction (HLC) term. This term is given

In these equations, n*<sup>α</sup>* and n*<sup>β</sup>* are the numbers of *α* and *β* electrons, respectively. The parameters A (4.567mH), B (2.363mH), C (4.544mH), and D (2.337mH) were obtained by

It is important to note that since it only depends on the number of *α* and *β* valence electrons, the HLC cancels entirely from most reaction energies, except when the reactions involve a mixture of atoms and molecules (as in heats of formation and bond dissociation energies)

It should be also noted that the absolute values of the calculated energies have no real significance, as no common energy zero for different atoms has been used. This arises from how the effective core potentials are constructed. When differences are formed, the differences between the zeros cancel, to leave for instance the difference in energy between the products

Additionally, an alternative approach for molecular geometry optimization and harmonic frequencies calculation can be considered through the use of the DFT (B3LYP, M06, ...).

The existence of the weight functions (graphic representations of the linear combination of the atomic orbitals) is the fundamental condition to use the GCM. The analysis of the behavior of the weight functions by the GCM allows the fitting of the atomic basis sets in order to get a

by either Eq. 6 or Eq. 7 depending on whether the species is a molecule or an atom:

+ *E*[*MP*2/*ECP* + *ACC*5*Zm*] − *E*[*MP*2/*ECP* + *ACC*2*Zm*] − (5)

*HLCmolec*. = −*C* · *n<sup>β</sup>* − *D* · (*n<sup>α</sup>* − *nβ*) (6)

*HLCatom* = −*A* · *n<sup>β</sup>* − *B* · (*n<sup>α</sup>* − *nβ*) (7)

For this electronic property, molecular calculations consist of:

**(b)** further optimization is carried out at MP2/ECP+ACC2Zm level;

corresponds to an effective calculation at a high level of theory,

E[MP2/ECP+ACC5Zm] is estimated throught the Eq. 4.

minima and to compute the zero-point energies;

where *scal* is the scaling factor on ZPE.

fitting to the experimental data of heats of formation.

and/or when spin is not conserved (Lin et al., 2009).

and reactants of a reaction.

**3. Results and discussion**

**3.1 Basis sets**

better description of the electrons in the molecular environment. The analysis of the weight functions of the outermost atomic orbitals suggested the need for improvements of the basis sets for the heavy atoms. Observing the plots of the weight functions of the outermost orbitals it is possible to establish the best fit of the basis sets. Using the atom of Cl as an example, the representation of the weight function of the atomic orbital 3s is shown in Fig. 1. The continuous line represents the plot of the weight function of the original primitive basis set for the all electron system. The dashed line represents the same weight function obtained using the pseudopotential. This figure shows that the 10 inner functions (large *α*) have a contribution close to zero towards the description of the weight function. The vertical solide line cutoff of the basis set indicates precisely where the pseudopotential starts to represent the core atomic region.

Fig. 1. Weight functions for the 3s atomic orbital of Cl in systems with all electrons (ae, continuous line), with the pseudopotential (pp, dashed line), with the addition of one diffuse function plus pseudopotential (dif+pp, dot-dashed line) and the cutoff point line represented by the vertical solid line.

#### **3.2 Proton affinity**

Proton affinity is a very sensitive property of the electronic structure and it is appropriate to test our methodology. Table 1 shows the results of the acidity (in kJ.mol−1) for a set of anions systems. These results were obtained using the B0 and B1 basis sets and Eq. 3. A comparison is also presented with experimental results. One can observe that our theoretical results are

System EA*<sup>a</sup>*

with experimental values.

Br 3.36 (-)*<sup>c</sup>* - Br2 2.53 (-0.11) 2.42 CBr3 2.49 (0.08) 2.57 ± 0.12 CCl3 2.22 (-0.05) 2.17 ± 0.10 CF3 1.74 (0.08) 1.82 ± 0.05 CF3COO 4.56 (-0.10) 4.46 ± 0.18 CH2=CH 0.71 (-0.04) 0.67 ± 0.02 CH2=CHCHBr 0.96 (-) - CH2=CHCHF 0.50 (-) - CH2=CCH3 0.69 (-0.13) 0.56 ± 0.19 CH2=NO2 2.40 (0.07) 2.48 ± 0.01 CH2Br 0.91 (-0.12) 0.79 ± 0.14 CH2BrCOO 3.99 (-0.01) 3.98 ± 0.16 CH2Cl 0.70 (0.04) 0.74 ± 0.16 CH2ClCOO 4.02 (-0.11) 3.91 ± 0.16 CH2F 0.17 (0.08) 0.25 ± 0.18 CH2FCOO 4.55 (-0.76) 3.79 ± 0.16 CH3 0.08 (0.00) 0.08 ± 0.03 CH3CH2O 1.75 (-0.03) 1.71 CH3CH2S 2.02 (-0.07) 1.95 CH3O 1.53 (0.05) 1.57 CH3S 1.75 (0.12) 1.87 CHBr2 1.74 (-0.03) 1.71 ± 0.08 CHBr2COO 4.37 (-0.11) 4.26 ± 0.16 CHClBr 1.64 (-0.17) 1.47 ± 0.04 CHF2COO 4.25 (-0.10) 4.15 ± 0.16 Cl 3.60 (0.01) 3.61 Cl2 2.55 (-0.15) 2.40 ± 0.20 F 3.34 (0.06) 3.40 F2 2.92 (0.08) 3.00 ± 0.07 Formamide-N,N-dimethyl 0.65 (-) - Formamide 3.34 (-) - HS 2.27 (0.05) 2.32 NH2 0.70 (0.07) 0.77 ± 0.01 OH 1.77 (0.06) 1.83 iPrO 1.87 (0.01) 1.87 nPrO 1.74 (0.05) 1.79 ± 0.03 *<sup>δ</sup>rmsd* 0.15 *<sup>a</sup>*The higher level correction (HLC) from G3 theory was added to the final energy. *<sup>b</sup>*Experimental values from NIST Webbase - http://webbook.nist.gov/chemistry/. *<sup>c</sup>*(EA*Exp*. - EA*Calc*.). Table 2. Electron Affinity (eV) calculated with the method given by Eq. 3, and comparison

Composite Method Employing Pseudopotential at CCSD(T) Level 19

*Calc*. EA*<sup>b</sup>*

*Exp*.


very close to the experimental errors (well within 5 kJ.mol−1) with root mean square deviation of 4.14 kJ mol−1.

*<sup>δ</sup>rmsd* 4.14 *<sup>a</sup>*Experimental values from NIST Webbase - http://webbook.nist.gov/chemistry/. *<sup>b</sup>*(PA*Exp*. - PA*Calc*.).

Table 1. Proton Affinity (kJ.mol−1) calculated with the method given by Eq. 3 and comparison with experimental values.

#### **3.3 Electron affinity**

In the Table 2 are the electron affinities calculated with the CCSD(T)/B1 energy from Eq. 3. It also shows a comparison between our results and the experimental values, where experimental data are available. The root mean square deviation (*δrmsd*) calculated is 0.15 eV.

The use of pseudopotential is competitive, mainly in systems containing S, Cl, and Br atoms. The computational time is almost constant for analogous systems with Cl and Br atoms, because in these cases we have an equal number of outer electrons. In calculations involving all electrons the computational performance is totally different and increases with the number of electrons. The CR-CCSD[T]/ECP computational demand is decreased by 10% when compared with all-electron calculations. For molecules containing Cl, Br or I atom the time is drastically decreased (Morgon, 2006).

#### **3.4 Heat of formation**

In the Table 3 are the heats of formation calculated with the CCR-CL/ ECP+ACC5Zm from Eq. 6. It also shows a comparison between our results and the experimental values, where experimental data are available. The average error using this methodology with respect to experimental results is closer to 10 kJ.mol−1.

8 Will-be-set-by-IN-TECH

very close to the experimental errors (well within 5 kJ.mol−1) with root mean square deviation

System PA*Calc*. PA*<sup>a</sup>*

F− 1556.17 (-2.17) 1554.0 CH2F<sup>−</sup> 1710.21 ( 0.79) 1711 ± 17 Cl− 1395.56 (-0.56) 1395.0 CH2Cl<sup>−</sup> 1658.18 (-1.18) 1657 ± 13 Br<sup>−</sup> 1354.80 (-1.80)*<sup>b</sup>* 1353.0 <sup>±</sup> 8.80 CHClBr<sup>−</sup> 1563.4 (0.82) 1560.0 ± 13.00 CH2Br<sup>−</sup> 1640.54 (2.46) 1643 ± 13

− 1315.10 1315.0

Table 1. Proton Affinity (kJ.mol−1) calculated with the method given by Eq. 3 and

In the Table 2 are the electron affinities calculated with the CCSD(T)/B1 energy from Eq. 3. It also shows a comparison between our results and the experimental values, where experimental data are available. The root mean square deviation (*δrmsd*) calculated is 0.15

The use of pseudopotential is competitive, mainly in systems containing S, Cl, and Br atoms. The computational time is almost constant for analogous systems with Cl and Br atoms, because in these cases we have an equal number of outer electrons. In calculations involving all electrons the computational performance is totally different and increases with the number of electrons. The CR-CCSD[T]/ECP computational demand is decreased by 10% when compared with all-electron calculations. For molecules containing Cl, Br or I atom the time is

In the Table 3 are the heats of formation calculated with the CCR-CL/ ECP+ACC5Zm from Eq. 6. It also shows a comparison between our results and the experimental values, where experimental data are available. The average error using this methodology with respect to

<sup>−</sup> 1617.47 (-1.47) 1616 ± 21 Acetamide,2,2,2-trichloride 1493,6 1436 ± 8.8 Formamide-N,N-dimethyl 1671.34 (-1.34) 1670.0 ± 17.00 Formamide 1510.13 (-5.13) 1505.0 ± 8.80 methyl cyclopentanol ( 2.88) 1559.00 ± 8.40 2,methyl cyclopentanol ( 2.99) 1558.00 ± 8.40 CH2CHCHI<sup>−</sup> 1555.24 (-4.24) 1551.0 ± 8.8 *<sup>δ</sup>rmsd* 4.14 *<sup>a</sup>*Experimental values from NIST Webbase - http://webbook.nist.gov/chemistry/. *<sup>b</sup>*(PA*Exp*. - PA*Calc*.).

*Exp*.

of 4.14 kJ mol−1.

I

CH2I

comparison with experimental values.

drastically decreased (Morgon, 2006).

experimental results is closer to 10 kJ.mol−1.

**3.3 Electron affinity**

**3.4 Heat of formation**

eV.


*<sup>δ</sup>rmsd* 0.15 *<sup>a</sup>*The higher level correction (HLC) from G3 theory was added to the final energy. *<sup>b</sup>*Experimental values from NIST Webbase - http://webbook.nist.gov/chemistry/. *<sup>c</sup>*(EA*Exp*. - EA*Calc*.).

Table 2. Electron Affinity (eV) calculated with the method given by Eq. 3, and comparison with experimental values.

**5. Acknowledgments**

*Phys.* 120: 4129–4141.

Ditter, G. & Niemann, U. (1982). *Phillips J. Res.* 37: 1. Dunning, T. H. (1989). *J. Chem. Phys.* 90: 1007. Ervin, M. K. (2001). *Chem. Rev.* 101: 391.

Feller, D. & Peterson, K. A. (1999). *J. Chem. Phys.* 110: 8384.

Kowalski, K. & Piecuch, P. (2000). *J. Chem. Phys.* 18: 113.

(2005). *J. Mass Spectrom.* 40: 1300.

Morgon, N. H. (1995a). *J. Phys. Chem. A* 99: 17832. Morgon, N. H. (1995b). *J. Phys. Chem.* 99: 17832. Morgon, N. H. (1998). *J. Phys. Chem. A* 102: 2050. Morgon, N. H. (2006). *Int. J. Quantum Chem.* 106: 2658. Morgon, N. H. (2008a). *J. Braz. Chem. Soc.* 19: 74. Morgon, N. H. (2008b). *Int. J. Quantum Chem.* 108: 2454. Morgon, N. H. (2011). *Int. J. Quantum Chem.* 111: 1555–1561.

Nelder, J. A. & Mead, R. (1965). *Computer J.* 7: 308.

Ge, Y., Gordon, M. S. & Piecuch, P. (2007). *J. Chem. Phys.* 174: 174106.

Martin, J. M. L. & De Oliveira, G. (1999). *J. Chem. Phys.* 111: 1843–1856.

Mohallem, J. R. & Trsic, M. (1985). *Z. Phys. A - Atoms and Nuclei* 322: 538.

Morgon, N. H. & Riveros, J. M. (1998). *J. Phys. Chem. A* 102: 10399.

Nyden, M. R. & Petersson, G. A. (1981). *J. Chem. Phys.* 75: 1843.

109: 7764–7776.

42: 411.

112: 6532.

119: 1708.

**6. References**

I would like to thank the computational facilities of Chemistry Institute at UNICAMP and the financial support from Conselho Nacional de Desenvolvimento Científico e Tecnológico

Composite Method Employing Pseudopotential at CCSD(T) Level 21

Badenes, M. P., Tucceri, M. E. & Cobos, C. J. (2000). *Zeitschrift für Physikalische Chemie* 214: 1193. Boese, A. D., Oren, M., Atasoylu, O., Martin, J. M. L., Kállay, M. & Gauss, J. (2004). *J. Chem.*

Curtiss, L. A., Raghavachari, K., Redfern, P. C. & Pople, J. A. (1997). *J. Chem. Phys.* 106: 1063. Curtiss, L. A., Raghavachari, K., Redfern, P. C. & Pople, J. A. (1998). *J. Chem. Phys.*

Curtiss, L. A., Raghavachari, K., Redfern, P. C. & Pople, J. A. (2000). *J. Chem. Phys.* 112: 7374.

Custodio, R., Giordan, M., Morgon, N. H. & Goddard, J. D. (1992). *Int. J. Quantum Chem.*

Karton, A., Rabinovich, E., Martin, J. M. L. & Ruscic, B. (2006). *J. Chem. Phys.* 125(14): 144108.

Mezzache, S., Bruneleau, N., Vekey, K., Afonso, C., Karoyan, P. Fournier, F. & Tabet, J.-C.

Montgomery Jr., J. A., Frisch, M. J., Ochterski, J. W. & Petersson, G. A. (2000). *J. Chem. Phys.*

Morgon, N. H., Argenton, A. B., Silva, M. L. P. & Riveros, J. M. (1997). *J. Am. Chem. Soc.*

Lin, C. Y., Hodgson, J. L., Namazian, M. & Coote, M. L. (2009). *J. Phys. Chem. A* 113: 3690.

Mohallem, J. R. & Dreizler, R. M. Trsic, M. (1986). *Int. J. Quantum Chem. - Symp.* 20: 45.

Custodio, R., Goddard, J. D., Giordan, M. & Morgon, N. H. (1992). *Can. J. Chem.* 70: 580. DeYonker, N. J., Cundari, T. R. & Wilson, A. K. (2006). *J. Chem. Phys.* 124: 114104.

Curtiss, L. A., Redfern, P. C. & Raghavachari, K. (2007). *J. Chem. Phys.* 126: 084108. Curtiss, L. A., Redfern, P. C., Smith, B. J. & Radom, L. (1996). *J. Chem. Phys.* 104: 5148.

(CNPq) and Fundação de Amparo à Pesquisa de São Paulo (FAPESP).


*<sup>a</sup>*Experimental values from NIST Webbase - http://webbook.nist.gov/chemistry/. *<sup>b</sup>* Ref. (Ditter & Niemann, 1982). *<sup>c</sup>*(Δ*<sup>f</sup> Hgas Calc*. - Δ*<sup>f</sup> Hgas Exp*.).

Table 3. Heats of formation (in kJ mol−1) with the method given by Eq. 6, and comparison with experimental values.

#### **4. Conclusions**

The proton and electron affinities and the heats of formation of some simple systems obtained by the procedure outlined in this paper are in very good agreement with experimental values. These results can be compared with those obtained by sophisticated and computationally more expensive calculations. ECP-based methods have been shown to be powerful, and of affordable computational cost for the systems addressed in this work. This is due to three features:


The use of adapted basis functions for atoms by the Generator Coordinate Method along with the use of the pseudopotential allows a high quality calculation at a lower computational cost.

The present methodology - Eqs. 3 and 6, which relies on small basis sets (representation of the core electrons by ECP) and an easier and simpler way for correcting the valence region (mainly of anionic systems) appears as an interesting alternative for the calculation of thermochemical data such as electron and proton affinities and enthalpies of formation for larger systems.

The CCSD(T)/B1 method have been shown to be powerful, and of affordable computational cost for the systems containing atoms of the 2nd and 3rd periods.

#### **5. Acknowledgments**

I would like to thank the computational facilities of Chemistry Institute at UNICAMP and the financial support from Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq) and Fundação de Amparo à Pesquisa de São Paulo (FAPESP).

#### **6. References**

10 Will-be-set-by-IN-TECH

OF C∞*<sup>v</sup>* <sup>2</sup>Π 110.999 (2.21)*<sup>c</sup>* 108.78 OF2 C2*<sup>v</sup>* 1A1 28.582 (4.06) 24.52 OCl C∞*<sup>v</sup>* <sup>2</sup>Π 103.482 (2.26) 101.22 OCl2 C2*<sup>v</sup>* 1A1 92.341 (4.48) 87.86 SF C∞*<sup>v</sup>* <sup>2</sup>Π 8.661 (-4.31) 12.97 SF2 C2*<sup>v</sup>* 1A1 -289.573 (7.08) -296.65 SF3 C*<sup>s</sup>* 2A' -512.605 (-9.57) -503.03 SF4 C2*<sup>v</sup>* 1A1 -775.362 (-12.20) -763.16 SF5 C4*<sup>v</sup>* 2A1 -887.770 (20.68) -908.45

SCl C∞*<sup>v</sup>* <sup>2</sup>Π 153.843 (-2.62) 156.47 SCl2 C2*<sup>v</sup>* 1A1 -19.262 (-1.69) -17.57 SCl3 C*<sup>s</sup>* 2A' 22.346 (-14.47) 36.82*<sup>b</sup>* SCl4 C2*<sup>v</sup>* 1A1 -15.301 (-12.38) -2.92*<sup>b</sup>* SCl5 C4*<sup>v</sup>* 2A1 -38.76 (-2.78) -35.98*<sup>b</sup>*

*<sup>a</sup>*Experimental values from NIST Webbase - http://webbook.nist.gov/chemistry/. *<sup>b</sup>* Ref. (Ditter & Niemann, 1982). *<sup>c</sup>*(Δ*<sup>f</sup> Hgas Calc*. - Δ*<sup>f</sup> Hgas Exp*.). Table 3. Heats of formation (in kJ mol−1) with the method given by Eq. 6, and comparison

The proton and electron affinities and the heats of formation of some simple systems obtained by the procedure outlined in this paper are in very good agreement with experimental values. These results can be compared with those obtained by sophisticated and computationally more expensive calculations. ECP-based methods have been shown to be powerful, and of affordable computational cost for the systems addressed in this work. This is due to three

The use of adapted basis functions for atoms by the Generator Coordinate Method along with the use of the pseudopotential allows a high quality calculation at a lower computational cost. The present methodology - Eqs. 3 and 6, which relies on small basis sets (representation of the core electrons by ECP) and an easier and simpler way for correcting the valence region (mainly of anionic systems) appears as an interesting alternative for the calculation of thermochemical data such as electron and proton affinities and enthalpies of formation for larger systems.

The CCSD(T)/B1 method have been shown to be powerful, and of affordable computational

*gas* calc Δ*<sup>f</sup> H<sup>o</sup>*

1A1*<sup>g</sup>* -1218.986 (1.48) -1220.47

1A1*<sup>g</sup>* -87.54 (-4.7) -82.84*<sup>b</sup>*

*gas* exp*<sup>a</sup>*

Molecule Point Group Ground State Δ*<sup>f</sup> H<sup>o</sup>*

SF6 O*<sup>h</sup>*

SCl6 O*<sup>h</sup>*

1) the number of steps employed during the calculations, 2) the smaller basis sets used in our methodology, and

cost for the systems containing atoms of the 2nd and 3rd periods.

with experimental values.

**4. Conclusions**

3) the use of ECP.

features:

Badenes, M. P., Tucceri, M. E. & Cobos, C. J. (2000). *Zeitschrift für Physikalische Chemie* 214: 1193. Boese, A. D., Oren, M., Atasoylu, O., Martin, J. M. L., Kállay, M. & Gauss, J. (2004). *J. Chem. Phys.* 120: 4129–4141.

Curtiss, L. A., Raghavachari, K., Redfern, P. C. & Pople, J. A. (1997). *J. Chem. Phys.* 106: 1063.

Curtiss, L. A., Raghavachari, K., Redfern, P. C. & Pople, J. A. (1998). *J. Chem. Phys.* 109: 7764–7776.


**Part 2** 

**Electronic Structures and Molecular Properties** 

Ochterski, J. W., Petersson, G. A. & Montgomery Jr., J. A. (1996). *J. Chem. Phys.* 104: 2598.

Ochterski, J. W., Petersson, G. A. & Wiberg, K. B. (1995). *J. Amer. Chem. Soc.* 117: 11299.

Oyedepo, G. A. & Wilson, A. K. (2010). *J. Chem. Phys.* 114: 8806.

Parthiban, S. & Martin, J. M. L. (2001). *J. Chem. Phys.* 114: 6014–6029.

Pereira, D. H., Ramos, A. F., Morgon, N. H. & Custodio, R. (2011). *J. Chem. Phys.* 135: 034106.

Peterson, K. A., Tensfeldt, T. G. & Montgomery Jr., J. A. (1991). *J. Chem. Phys.* 94: 6091.

Piecuch, P., Kucharski, S. A., Kowalski, K. & Musial, M. (2002). *Comp. Phys. Commun.* 149: 71.

Smith, B. J. & Radom, L. (1991). *J. Phys. Chem.* 95: 10549.

Stevens, W. J., Basch, H. & Krauss, M. (1984). *J. Chem. Phys.* 81: 6026.

Tabrizchi, M. & Shooshtari, S. (2003). *J. Chem. Thermodyn.* 35: 863.
