# Analyzing ab initio infrared spectra and electronic properties of polyethylenimine water complexes in the solid state

G. Herlem, T. Gharbi, N. Ben Sedrine

► To cite this version:

G. Herlem, T. Gharbi, N. Ben Sedrine. Analyzing ab initio infrared spectra and electronic properties of polyethylenimine water complexes in the solid state. Journal of Molecular Structure: THEOCHEM, 2010, 945 (1-3), pp.27-32. 10.1016/j.theochem.2010.01.009. hal-00535895

HAL Id: hal-00535895

https://hal.science/hal-00535895v1

Submitted on 29 Apr 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Analyzing *ab initio* infrared spectra and electronic properties of polyethylenimine water complexes in the solid state

Guillaume Herlem<sup>a,*</sup>, Tijani Gharbi<sup>a</sup>, Nebiha Ben Sedrine<sup>b</sup>

We focused our study on the hydrated polyethylenimine (PEI) crystals which are characterized by the presence of 0.5, 1.5 or 2 water molecules per ethylenimine monomer unit. We propose a theoretical study in the framework of periodic linear combination of atomic orbitals (LCAO) using different Hamiltonians based on density functional theory (DFT) and hybrid Hartree-Fock (HF) - DFT level of theory. Once hydro-gen atoms added when missing and redundant atoms withdrawn in the PEI structures, *ab initio* calculations yielded relaxed geometries at the ground state. The closest structural parameters to the experimental values are observed for hybrid Hamiltonians and the use of the 6-31++G* basis set reduces this difference. We provide the IR vibration modes at the *Γ* point, with assignments of the bands. The electronic study reveals that these PEI hydrates have an insulating character with a quasi constant band gap whatever the number of water molecules. From the projected density of states (PDOS), we have found that the nitrogen orbital has the major contribution to the valence band, while the conduction band is mainly constructed from hydrogen bondings.

## 1. Introduction

PEI is widely used in different industrial, agricultural and biological applications [1–3]. Due to its high cationic state and its ability to complexe DNA, the linear form of PEI is very attractive for gene transfection. In fact, this polymer takes many forms: linear or branched, hydrated or anhydrous, solid or liquid depending on its molecular weight and its synthesis way: organic, electrochemical [4–8]. In addition, water molecules can unroll the doubled stranded PEI structure when anhydrous leading to linear zigzag polymer chains. Few theoretical works are devoted to PEI especially concerning vibration and electronic properties of this polymer. The reason probably comes from the crystallographic file structures containing atomic redundant positions or missing hydrogen. Consequently it is necessary to correct these structures in order to investigate some relevant properties of the hydrated PEI crystal structures at a sufficient level of theory.

In a previous paper, we studied the anhydrous linear PEI (also known as crystalline PEI or L-PEI) at the B3LYP/6-21G level of theory [9]. In this work, we propose to study the effect of the presence of water molecules and additional intermolecular interactions on the infrared spectra and electronic properties of the hydrate complexes. We present periodic *ab initio* calculations on the crystalline

PEI hydrates at different level of theory using two different basis sets: 6-21G and 6-31++G* [10].

## 2. Geometry and computational details

The starting PEI hydrate structures obtained from X-ray diffraction were taken from the literature and belong to the monoclinic C2/c space group [12–14]. The lattice parameters of the hemihydrate PEI are *a* = 10.89 Å, *b* = 9.52 Å, *c* = 7.31 Å, *β* = 127.6°, whereas those of the sesquihydrate PEI are *a* = 11.55 Å, *b* = 9.93 Å, *c* = 7.36 Å, *β* = 104.5°, and those of the dihydrate PEI are *a* = 13.26 Å, *b* = 4.61 Å, *c* = 7.36 Å, *β* = 101.0°. The hemihydrate PEI structure contains eight monomers –CH₂–CH₂–NH– and four water molecules, the sesquihydrate form has eight monomer units and 12 water molecules in its unit cell and the dihydrate four monomers and eight water molecules.

The electronic structures and properties of the three hydrates were performed with the CRYSTAL06 code [15,16] within the LCAO method using gaussian-type function 6-21G or 6-31++G* basis sets. Four different DFT Hamiltonians were used to construct Bloch sums as basis functions for the crystalline orbitals. Self-consistent field (SCF) calculations were carried out to solve the DFT Hamiltonian where exchange-correlation functionals are separated in an exchange component (keyword EXCHANGE) and a correlation component (keyword CORRELAT). LDA is a local functional of the electron density only, PBE functional takes into account the

<footnote><sup>a</sup>Institut Femto-ST/Department of Optic UMR CNRS 6174, Bat. E, University of Franche-Comté, UFR Sciences & Techniques, 16 Route de Gray, 25030 Besançon Cedex, France</footnote>

<footnote><sup>b</sup>Laboratoire de Photovoltaïque et de Semiconducteurs, Centre de Recherche et de Technologie de l'Energie, BP. 95, Hammam Lif 2050, Tunisia</footnote>

<footnote>* Corresponding author.
E-mail address: guillaume.herlem@univ-fcomte.fr (G. Herlem).</footnote>
<table><caption>**Table 1**
Global keywords defining hybrid exchange-correlation functional completely.</caption><thead><tr><th rowspan="2">Name</th><th colspan="4">Keywords</th></tr><tr><th>EXCHANGE</th><th>CORRELAT</th><th>HYBRID</th><th>NON-LOCAL</th></tr></thead><tbody><tr><td>LDA</td><td>LDA</td><td>PZ</td><td>-</td><td>-</td></tr><tr><td>PBE</td><td>PBE</td><td>PBE</td><td>-</td><td>-</td></tr><tr><td>B3LYP</td><td>BECKE</td><td>LYP</td><td>20</td><td>0.9-0.81</td></tr><tr><td>PBE0</td><td>PBE</td><td>PBE</td><td>25</td><td>-</td></tr></tbody></table>

gradient of the electron density (generalized gradient approximation, GGA) while B3LYP [17,18] and PBE0 are hybrid ones. Part of the exact HF exchange is incorporated into the exchange functional through the keyword HYBRID. The resulting exchange functional is then a linear combination of HF, local and gradient-corrected exchange term (Table 1). NON-LOCAL allows modifying the relative weight of the local and non-local part of both the exchange and the correlation potential.

The SCF controls used for numerical accuracy and convergence were TOLINTEG, the thresholds for bielectronic integral truncations (6, 6, 6, 6, 12) corresponding to the Coulomb and exchange overlap and TOLDEE, the convergence threshold on total energy (10⁻⁹). This level of numerical accuracy did not bring any basis set linear dependence. Higher threshold values for bielectronic integral truncations (7, 7, 7, 7, 14) led to the same results as for standard parameters. An improved convergence rate was obtained with FMIXING by mixing Fock/Kohn-Sham matrix elements of 30% and LEVSHIFT, shifting the eigenvalue levels of 1 Hartree before the diagonalization and removing it after. Sampling of *k* points in reciprocal space according to the Monkhorst-Pack method in the irreducible Brillouin zone was defined by SHRINK, the shrinking factor in reciprocal space (4) and for Gilat net (4) for all the systems.

# 3. Results and discussion

## 3.1. Geometry optimization

The optimization of atomic coordinates with fixed lattice parameters and minima localization on the potential energy surface was performed using a quasi-Newton optimization scheme combined with a quadratic linear search by means of the BFGS algorithm [16]. The direction of the step is computed at each cycle, while the length is determined by linear minimization along an extrapolated quadratic polynomial, and the test on total energy change between two steps was set to 10⁻⁹ Hartree.

Prior to the study of the different PEI hydrates, various functions were tested on the PEI hemihydrate lattice in order to find the most suited level of theory for these systems. In this purpose, some bond lengths and angles were compared after geometry optimization in Table 2. As the LDA functional approximates the energy of the true density by the energy of a local constant density, it fails in situations where the density undergoes rapid changes such as PEI hydrates. An improvement was made by considering the gradient of the electron density, the so-called generalized gradient approximation (GGA) [19]. In addition, the structure of the hemihydrate at the B3LYP level of theory (reference) was compared in

<table><caption>**Table 3**
Structural comparison between two PEI hemihydrates conformations: a reference and an alternative one optimized with a different level of theory than the reference. The results are the root mean square deviation (RMSD).</caption><thead><tr><th>Reference</th><th>Alternative</th><th>RMSD</th></tr></thead><tbody><tr><td rowspan="3">Hemihydrate B3LYP/6-21G</td><td>Hemihydrate LDA/6-21G</td><td>0.066</td></tr><tr><td>Hemihydrate PBE/6-21G</td><td>0.023</td></tr><tr><td>Hemihydrate PBE0/6-21G</td><td>0.015</td></tr><tr><td rowspan="3">Hemihydrate B3LYP/6-31++G*</td><td>Hemihydrate LDA/6-31++G*</td><td>0.087</td></tr><tr><td>Hemihydrate PBE/6-31++G*</td><td>0.027</td></tr><tr><td>Hemihydrate PBE0/6-31++G*</td><td>0.013</td></tr></tbody></table>

Table 3 to the one optimized with other functionals (alternative). The EnsembleMatch tool supplied in the Chimera's tool menu was used to calculate the superposition and root mean square deviation (RMSD) values for every equivalent atom between two molecules containing the same set of atoms [20]. The lower the RMSD value, the higher the structures match. As expected, the structure optimized with B3LYP/6-31++G* is the closest to the one obtained with PBE0/6-31++G*. The RMSD between hemihydrates LDA/6-31++G* and B3LYP/6-31++G* is greater than that between hemihydrates LDA/6-21G and B3LYP/6-21G denoting a noticeable improvement that comes from the choice of the 6-31++G* basis set over 6-21G. The 6-21 split valence basis set, which is usually adequate for covalent crystals and near to the variational limit in crystals, seems to be too small to render efficiently the intermolecular hydrogen bonds [11].

From these results, the B3LYP/6-31++G* level of theory was applied to the other PEI hydrates. Atomic coordinates from our calculations are given as Supplementary data. From the three starting structures, some corrections were essential because some atomic positions were redundant or wrong (sesqui and dihydrate) or missing (hydrogen atoms in hemihydrate). A key to reliable geometry optimization, comes from the polymer chains keeping their planar all-trans zigzag conformation in the solid state that has been achieved at the B3LYP/6-31++G* level of theory. In all the hydrates, the N···H···O intermolecular bonds are shorter than the N···H···O ones, as it was supposed by the authors of the experimental structure diffractions [12–14]. We find in the three hydrates that d(N···H) < d(O···H), it could be a consequence that d(O···H) < d(N···H). In return, a good agreement is found between the optimized atomic parameters and the experimental ones especially the bonds as shown in Table 4. The use of the medium split valence basis set 6-31++G* proposed for molecules gave geometrical features (Fig. 1) that match well the experimental values [12–14].

## 3.2. Phonon calculations at Γ center of the Brillouin zone

Phonon frequencies were obtained via the response function calculation of the total energy second derivative with respect to phonon perturbation. Details of the computation of the full dynamical matrix at Γ can be found elsewhere [21,22]. The geometries of the three hydrates are at the ground state, since only three imaginary frequencies between 0 and −10 cm⁻¹ corresponding to the rotation and translation modes were obtained in each lattice. The whole set of the calculated modes together with their related infra-

<table><caption>**Table 2**
Comparison of bond length and angle calculated using LDA, PBE, B3LYP and PBE0 Hamiltonians, with the X-ray diffraction data from [12–14] for PEI hemihydrate.</caption><thead><tr><th></th><th>LDA/6-21G</th><th>LDA/6-31++G*</th><th>PBE/6-21G</th><th>PBE/6-31++G*</th><th>B3LYP/6-21G</th><th>B3LYP/6-31++G*</th><th>PBE0/6-21G</th><th>PBE0/6-31++G*</th><th>Expt</th></tr></thead><tbody><tr><td>C–C</td><td>1.509</td><td>1.512</td><td>1.517</td><td>1.520</td><td>1.517</td><td>1.519</td><td>1.515</td><td>1.518</td><td>1.539</td></tr><tr><td>C–N</td><td>1.484</td><td>1.457</td><td>1.493</td><td>1.465</td><td>1.487</td><td>1.462</td><td>1.482</td><td>1.457</td><td>1.467</td></tr><tr><td>C–C–N</td><td>111.5</td><td>112.4</td><td>109.0</td><td>111.2</td><td>109.5</td><td>110.7</td><td>110.0</td><td>111.2</td><td>110.7</td></tr><tr><td>C–N–C</td><td>107.3</td><td>109.7</td><td>107.3</td><td>109.7</td><td>108.6</td><td>110.6</td><td>108.6</td><td>110.7</td><td>110.7</td></tr></tbody></table>
<table><caption>**Table 4**
Structural parameters of the three PEI hydrates optimized at the B3LYP/6-31++G* level of theory compared to the experimental ones.</caption><thead><tr><th rowspan="2"></th><th colspan="2">Sesquihydrate</th><th colspan="2">Hemihydrate</th><th colspan="2">Dihydrate</th></tr><tr><th>B3LYP/6-31++G*</th><th>Experiment<fnref n="1" /></th><th>B3LYP/6-31++G*</th><th>Experiment<fnref n="1" /></th><th>B3LYP/6-31++G*</th><th>Experiment<fnref n="1" /></th></tr></thead><tbody><tr><td>C-H</td><td>1.103</td><td>1.092</td><td>1.103</td><td>-</td><td>1.103</td><td>1.095</td></tr><tr><td>N-H</td><td>1.031</td><td>0.993</td><td>1.033</td><td>-</td><td>1.032</td><td>0.986</td></tr><tr><td>O-H</td><td>0.978</td><td>0.981</td><td>1.004</td><td>-</td><td>0.999</td><td>0.984</td></tr><tr><td>C-C</td><td>1.528</td><td>1.538</td><td>1.519</td><td>1.54</td><td>1.520</td><td>1.539</td></tr><tr><td>C-N</td><td>1.468</td><td>1.473</td><td>1.461</td><td>1.50</td><td>1.474</td><td>1.467</td></tr><tr><td>H-C-H</td><td>106.4</td><td>110.1</td><td>106.8</td><td>-</td><td>109.9</td><td>110.2</td></tr><tr><td>C-C-N</td><td>110.5</td><td>110.4</td><td>110.7</td><td>110.7</td><td>109.0</td><td>110.7</td></tr><tr><td>C-N-C</td><td>112.5</td><td>110.1</td><td>110.6</td><td>110.7</td><td>111.3</td><td>110.7</td></tr><tr><td>N···H</td><td>1.757</td><td>-</td><td>1.897</td><td>-</td><td>1.718</td><td>-</td></tr><tr><td>O···H</td><td>2.000</td><td>-</td><td>2.048</td><td>-</td><td>2.178</td><td>-</td></tr></tbody></table>

red (IR) intensities and the normal modes identified thanks to the ANALYSIS keyword are listed in Supplementary data. Hemihydrate, sesquihydrate and dihydrate PEI present 114, 150 and 84 infrared and/or Raman active vibration modes respectively. The resulting irreducible representations of the vibration modes at the Γ point are given as follows:

$$ \Gamma_{\text{hemihydrate}} = 28A_g \oplus 29B_g \oplus 28A_u \oplus 29B_u $$

$$ \Gamma_{\text{sesquihydrate}} = 37A_g \oplus 38B_g \oplus 37A_u \oplus 28B_u $$

$$ \Gamma_{\text{dihydrate}} = 42A_g \oplus 42A_u $$

The A<sub>g</sub> and B<sub>g</sub> modes are RAMAN active, whereas the A<sub>u</sub> and B<sub>u</sub> modes are IR active. As it can be seen, the three lower frequency modes correspond to the three translation modes of the C2/c space group with the irreducible representation A<sub>u</sub> ⊕ 2B<sub>u</sub> (2B<sub>u</sub> correspond to the two degenerated translation modes given by the calculation procedure). In fact, we obtained three imaginary frequencies between 0 and -10 cm<sup>-1</sup> (or), small enough for our systems to be considered as stable. These values result from the exchange and correlation functional that cannot reach the translation invariance property because it would require a grid with a very high resolution.

Computed IR spectra with normalized intensities (FWHM = 30 cm<sup>-1</sup>) in the frequency regions 400-2000 and 2750-3750 cm<sup>-1</sup> are gathered Fig. 2a and b, respectively. We can notice that the computed infrared spectra are quite different each others as this is the case for the experimental spectra. The explanation has to be found on the level of the increasing hydrogen bond strength with increasing water content, since the hydrates have identical polymeric chemistry. The three different hydrogen bonds N-H···O, O-H···O and O-H···N, formed between the polymer chain have a prominent role in the cohesion of the PEI hydrate lattices. Linear polyamine can adopt much conformational behavior. However, the all-trans linear zigzag chain, which is strongly determined by the occurrence of intermolecular hydrogen bonds in the solid state, has the lower energy conformation [23,24]. This kind of close contact gives rise to infinite polymeric chains here.

The definitive assignments of the vibration modes from our calculations complete the characteristic bands given elsewhere and allow raising the doubt on some bands [25-28]. In case of overlapping bands during the hydration process of the crystalline PEI with the mixing of different phases, the bands at 876 and 918 cm<sup>-1</sup> are clearly a superposition of hemi-, sesqui- and dihydrate. In spite of our assignment, only high resolution spectra could decide between one and the other.

But, from our calculations (Fig. 2a) unambiguous crystal phase identification can be carried out from characteristic regions with strong bands such as 1700-1800 cm<sup>-1</sup> for the dihydrate, 1300-1400 cm<sup>-1</sup> for the sesquihydrate and 1200-1300 cm<sup>-1</sup> for the hemihydrate.

<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /><figcaption>**Fig. 1.** Optimized crystal lattices of (a) hemihydrates, (b) sesquihydrate and (c) dihydrate.</figcaption></figure>

In Fig. 2c, we compare our computed IR spectrum of the PEI dihydrate (FWHM = 30 cm<sup>-1</sup>) with the experimental data by Hashida et al. [27]. Although all peak positions of the theoretical

<footnote n="1">¹ Experimental values are taken from Refs. [12,13].</footnote>
<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /><figcaption>**Fig. 2.** Calculated IR spectra of the PEI hydrates at the ground state between (1) 400 and 2000 cm⁻¹, (2) between 2750 and 3750 cm⁻¹ and (3) computed IR spectrum of dehydrate compared to the experimental data by Hashida et al. [26].</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Fig. 3.** Calculated electronic band structure and total density of states of PEI (a) hemihydrates from band 40 to 60, (b) sesquihydrate from band 60 to 90 and (c) dehydrate from band 30 to 50 [29].</figcaption></figure>

spectrum are slightly shifted from the experimental ones, our computed spectrum without scaling factor correction matches well all the bands. In addition, this good agreement between our calculations and experience confirm the dihydrate geometry we propose.

The anharmonicity for X–H stretching modes (where X = C, N and O atoms) was also investigated via the CRYSTAL06 code. As expected, we have found a decrease of all the frequency values of the stretching modes in the case of the anharmonic oscillator: C–H, N–H and O–H frequencies are slightly decreased of about 240, 180 and 100 cm⁻¹, respectively, due to the intermolecular hydrogen bond strength.

### 3.3. Electronic properties

Both the electronic band structures and the corresponding total density of states (TDOS, dot line) are displayed in Fig. 3 [29]. The computed band structure clearly reveals a band gap of 6.5 eV with the Fermi level at the top of the valence band which corresponds to an expected insulating behavior. The difference it exists with L-PEI may come from the crystal lattices, the way the polymer chains are stacked (helix for L-PEI and linear zigzag chains for hydrates) but above all from the water molecule presence.

However, it should be recalled that theoretically derived band gaps in insulators and semiconductors may be underestimated

using DFT calculations. This does not affect the accuracy of the description of the total energy and related properties of crystals.

The PDOS on the left of the band structures corresponding to hydrogen, nitrogen, carbon, and oxygen atom contributions show that the nitrogen orbital has the major contribution to the valence band with a few oxygen and hydrogen orbital contributions, while the conduction band is mainly constructed from hydrogen and oxygen orbitals. We can see also that water molecule influence is preponderant for conduction bands whereas it is lower for valence bands. The valence bandwidths are narrow with sharp TDOS, indicating small overlap between neighbors. On the other hand, the lowest unoccupied bands in the conduction band and corresponding to hydrogen atom projection have prominent contribution to the PDOS. The resulting broad bands in the PDOS for all the hydrates come from the hydrogen bond overlapping with neighbors.

## 4. Conclusions

In summary, we have presented a detailed structural, vibrational and electronic analysis of the PEI hydrate lattices based on periodic *ab initio* calculations. B3LYP/6-31++G* yields a relaxed system in the ground state for the three systems. The original hydrate structures were corrected from missing or redundant atomic positions. It was then possible to calculate the IR spectra with band intensities, and to make some corrections for the anharmonic stretching modes where a better assignment of bands can be made due to various hydrogen bonding. We have found that the water molecule presence induces important differences between anhydrous L-PEI and hydrates, especially in the crystal structure and electronic point of views.

## Acknowledgment

One of the authors (G.H.) is grateful to the CINES, FRANCE for supercomputing facilities and the kind help from the staff concerning the cib2438 project.

## Appendix A. Supplementary data

Supplementary data associated with this article can be found, in the online version, at doi:10.1016/j.theochem.2010.01.009.

## References

[1] D. Mandler, A. Kaminsky, I. Willner, Electrochim. Acta 37 (1992) 2765-2767.

[2] A. Kaminsky, I. Willner, D. Mandler, J. Electrochem. Soc. 40 (3) (1993) L215-L227.

[3] B. Lakard, G. Herlem, M. de Labachelerie, W. Daniau, G. Martin, J.-C. Jeannot, L. Robert, B. Fahys, Biosens. Bioelect. 19 (2004) 595-606.

[4] P.A. Gembitskii, N.A. Kleshcheva, A.I. Chhmarin, D.S. Zuk, Polym. Sci. URSS (Engl. Trans.) 20 (1978) 2982-2984.

[5] C.R. Dick, G.E. Ham, J. Macromol. Sci. Chem. A4 (6) (1970) 1301-1311.

[6] G. Herlem, K. Reybier, A. Trokourey, B. Fahys, J. Electrochem. Soc. 147 (2) (2000) 597.

[7] B. Lakard, G. Herlem, A. Etcheberry, J. Morvan, B. Fahys, Surf. Sci. 502 (2002) 296-303.

[8] B. Lakard, J.-C. Jeannot, M. Spajer, G. Herlem, M. de Labachelerie, P. Blind, B. Fahys, Electrochim. Acta 50 (9) (2005) 1863-1869.

[9] G. Herlem, J. Lakard, J. Chem. Phys. 120 (19) (2004) 9376-9382.

[10] P.C. Hariharan, J.A. Pople, Theor. Chim. Acta 28 (1973) 213-222 (Basis sets were obtained from the Extensible Computational Chemistry Environment Basis Set Database, Version 02/25/04, as developed and distributed by the Molecular Science Computing Facility, Environmental and Molecular Sciences Laboratory which is part of the Pacific Northwest Laboratory, P.O. Box 999, Richland, Washington 99352, USA. http://www.emsl.pnl.gov/forms/basisform.html).

[11] M.A. Spackman, A.S. Mitchell, Phys. Chem. Chem. Phys. 3 (2001) 1518-1523.

[12] Y. Chatani, H.T. Tadokoro, H. Ikeda Saegusa, Macromolecules 14 (1981) 315-321.

[13] Y. Chatani, T. Kobatake, H. Tadokoro, Macromolecules 15 (1982) 170-176.

[14] Y. Chatani, T. Kobatake, H. Tadokoro, Macromolecules 16 (1983) 199-204.
[15] B. Civalleri, P. D'arco, R. Orlando, V.R. Saunders, R. Dovesi, Chem. Phys. Lett. 348 (2001) 131-138.

[16] R. Dovesi, V.R. Saunders, C. Roetti, R. Orlando, C.M. Zicovich-Wilson, F. Pascale, B. Civalleri, K. Doll, N.M. Harrison, J. Bush, Ph. D'Arco, M. Llunell, CRYSTAL06 User's Manual, University of Torino, Italy, 2007.

[17] A.D. Becke, Phys. Rev. A 38 (1988) 3098-3100.

[18] C. Lee, W. Yang, R.G. Parr, Phys. Rev. B 37 (1988) 785-789.

[19] A. St-Amant, W.D. Cornell, P.A. Kollman, T.A. Halgren, J. Comput. Chem. 16 (1995) 1483-1506.

[20] E.C. Meng, E.F. Pettersen, G.S. Couch, C.C. Huang, T.E. Ferrin, BMC Bioinform. 12 (7) (2006) 349.

[21] F. Pascale, C.M. Zicovich-Wilson, F. Lopez, B. Civalleri, R. Orlando, R. Dovesi, J. Comput. Chem. 25 (2004) 888-897.

[22] C.M. Zicovich-Wilson, F. Pascale, C. Roetti, V.R. Saunders, R. Orlando, R. Dovesi, J. Comput. Chem. 25 (2004) 1873-1881.

[23] M.P.P. Marques, L.A.E. Batista de Carvalho, J. Tomkinson, J. Phys. Chem. A 106 (2002) 2473-2482.

[24] S.S. York, S.E. Boesch, R.A. Wheeler, R. Frech, Macromolecules 36 (2003) 7348-7351.

[25] T. Hashida, K. Tashiro, Y. Inaki, Polymer 44 (2003) 1721-1724.

[26] T. Hashida, K. Tashiro, Polymer 48 (2007) 7614-7622.

[27] T. Hashida, K. Tashiro, S. Aoshima, Y. Inaki, Macromolecules 35 (2002) 4330-4336.

[28] T. Hashida, K. Tashiro, Y. Inaki, J. Polym, Science B 41 (2003) 2937-2948.

[29] B.G. Searle, Comput. Phys. Commun. 137 (2001) 25-32.