Chemical Reactivity Theory (CRT) Study of the Melanoidin M8: Local Conceptual Density Functional Theory Descriptors

Juan Frau¹, Norma Flores-Holguín², Daniel Glossman-Mitnik^1,2*

¹Departament de Química, Universitat de les Illes Balears, Palma de Mallorca, Spain.
²Laboratorio Virtual NANOCOSMOS, Departamento de Medio Ambiente y Energía, Centro de Investigación en Materiales Avanzados, Chihuahua, México.

DOI: 10.4236/cmb.2018.82004                

Abstract

This study evaluated a fixed long-range corrected range-separated hybrid (RSH) density functional associated with the Def2TZVP basis set alongside the Solvation Model based on Density (SMD) for the computation of the structure, molecular properties and chemical reactivity of the M8 intermediate melanoidin pigment. The preference of the active sites pertinent to radical, nucleophilic and electrophilic attacks is made through linking them with the electrophilic and nucleophilic Parr functions, Fukui function indices, and condensed Dual Descriptor. This study showed that the MN12SX density functional is the most suitable one for predicting the chemical reactivity of this system.

Keywords

Melanoidin M8, Conceptual DFT, Chemical Reactivity, Dual Descriptor, Parr Function

1. Introduction

The Maillard reaction was first observed by the French chemist, Louis-Camille Maillard. It involves a complex network of non-enzymatic reactions resulting from the initial condensation between an available amino group and a carbonyl-containing moiety, usually a reducing sugar. When a protein reacts with a reducing sugar, a compound known as a Schiff base is formed, and the accumulation of these compounds causes degradation of proteins, leading to the production of advanced glycation endproducts (AGEs), which in turn cause the build-up of fibrils in the brain. When biological molecules are left to react with reducing sugars under physiological conditions, the same reaction takes place through a process called glycation. In recent times, our emphasis has been on understanding how glycation occurs as well as on the chemical reactivity that the reducing carbohydrates have with the amino acids and peptides that take part in the process, which are frequently associated with a number of diseases like Alzheimer, Parkinson, and diabetes.

The lightweight intermediate colored products are known as colored Maillard reaction products (CMRPs) [1] . The CMRPs are attracting interest from the food and other industries because it has been hypothesized that they can be used as prospective photosensitizers, whereby they can be used as antioxidants as well as colorants that can be applied in dye-sensitized solar cells to produce substitute energy. These prospective applications elucidate the rationale behind the focus on their chemical properties, amid their molecular reactivity, attracting increasing interest. Conceptual density functional theory (DFT), also known as chemical reactivity theory is a highly influential tool used to analyze, interpret, and predict the effects generated by chemical reactions. This topic was first studied by Parr et al. [2] . A number of useful concepts were developed after their work using density of molecular systems by employing DFT. These concepts allow researchers to qualitatively envisage how the chemical reactivity will occur in a certain system. In addition, there is a possibility that these concepts can be quantified. Collectively, these concepts are known as Conceptual DFT descriptors. Thus, Melanodin M8 is amenable to be studied by analyzing their molecular reactivity properties.

One of these isolated molecules is called Melanodin M8 which has interesting properties as a colored molecule in water and we believe that it could be of interest to study their molecular reactivity by using the ideas of Conceptual DFT [2] [3] [4] [5] [6] , in the same way of our previous works [7] [8] [9] [10] . Thus, in this computational study we will assess a powerful density functional in calculating the molecular properties and structure of the Melanoidin M8 pigment in water. Following the same ideas of previous works, we will consider a fixed RSH functional instead of the optimally-tuned RSH density functionals that have attained great success [11] - [19] .

2. Theoretical Background

The theoretical background of this study is similar to the previous conducted research presented [7] [8] [9] [10] , and will be shown here for complete purposes, because this research is a component of a major project that it is in progress.

If we consider the KID (for Koopmans in DFT) procedure presented in our previous works [7] [8] [9] [10] together with a finite difference approximation, then the global reactivity descriptors can be written as:

Electronegativity $\chi =-\frac{1}{2}\left(I+A\right)\approx \frac{1}{2}\left({ϵ}_L+{ϵ}_H\right)$ [20] [21]

Global Hardness $\eta =\left(I-A\right)\approx \left({ϵ}_L-{ϵ}_H\right)$ [20] [21]

Electrophilicity $\omega =\frac{{\mu }^2}{2\eta }=\frac{{\left(I+A\right)}^2}{4\left(I-A\right)}\approx \frac{{\left({ϵ}_L+{ϵ}_H\right)}^2}{4\left({ϵ}_L-{ϵ}_H\right)}$ [22]

Electrodonating Power ${\omega }^{-}=\frac{{\left(3I+A\right)}^2}{16\left(I-A\right)}\approx \frac{{\left(3{ϵ}_H+{ϵ}_L\right)}^2}{16\eta }$ [23]

Electroaccepting Power ${\omega }^{+}=\frac{{\left(I+3A\right)}^2}{16\left(I-A\right)}\approx \frac{{\left({ϵ}_H+3{ϵ}_L\right)}^2}{16\eta }$ [23]

Net Electrophilicity $\Delta {\omega }^{\pm }={\omega }^{+}-\left(-{\omega }^{-}\right)={\omega }^{+}+{\omega }^{-}$ [24]

where ${ϵ}_H$ and ${ϵ}_L$ are the energies of the highest occupied and the lowest unoccupied molecular orbitals (HOMO and LUMO), respectively.

Applying the same ideas, the definitions for the local reactivity descriptors are:

Nucleophilic Fukui Function $f^{+}\left(r\right)={\rho }_{N+1}\left(r\right)-{\rho }_N\left(r\right)$ [20]

Electrophilic Fukui Function $f^{-}\left(r\right)={\rho }_N\left(r\right)-{\rho }_{N-1}\left(r\right)$ [20]

Dual Descriptor $\Delta f\left(r\right)={\left(\frac{\partial f\left(r\right)}{\partial N}\right)}_{\upsilon \left(r\right)}$ [25] - [31]

Nucleophilic Parr Function $P^{-}\left(r\right)={\rho }_s^{rc}\left(r\right)$ [32] [33]

Electrophilic Parr Function $P^{+}\left(r\right)={\rho }_s^{ra}\left(r\right)$ [32] [33]

where ${\rho }_{N+1}\left(r\right)$ , ${\rho }_N\left(r\right)$ , and ${\rho }_{N-1}\left(r\right)$ are the electronic densities at point $r$ for the system with N + 1, N, and N − 1 electrons, respectively, and ${\rho }_s^{rc}\left(r\right)$ and ${\rho }_s^{ra}\left(r\right)$ ) are related to the atomic spin density (ASD) at the r atom of the radical cation or anion of a given molecule, respectively [34] .

3. Results and Discussion

Following the lines of our previous work [7] [8] [9] , the computational studies were performed with the Gaussian 09 [35] series of programs with the MN12SX density functional as implemented in the computational package. The basis set used in this work was Def2SVP for geometry optimization and frequencies, while Def2TZVP was considered for the calculation of the electronic properties [36] [37] . All the calculations were performed using water as a solvent by integral equation formalism-polarized continuum model (IEF-PCM) computa- tions according to the SMD solvation model [38] .

In this study, the molecular structure of the M8 melanoidin pigment was obtained from PubChem (https://pubchem.ncbi.nlm.nih.gov), a website that serves as the public storage area for information related to chemical substances and the biological activities associated with them. The selection of the most stable conformer constitutes the pre-optimization of the final system. The selection was performed through random sampling that involved molecular mechanics methods and the addition of a variety of torsional angles using the general MMFF94 force field [39] [40] [41] [42] [43] . This was done using the Marvin View 17.15 program that comprises an advanced chemical viewer favorable for single and multiple reactions, structures and chemical queries (https://www.chemaxon.com). Subsequently, the structure assumed by the resulting lower-energy conformer was reoptimized n the same way as in our recent works on melanoidins [7] [8] [9] [10] , where it has been found that the model chemistry formed by the connection between the MN12SX density functional and the Def2TZVP basis set is the best for justifying the fulfilling of the KID procedure for the Melanoidin M8 molecule, both in the presence of water and dioxane. Thus, instead of presenting the comparison of the values of the orbital energies with the ionization potential I and the electron affinity A for the different functionals, we are showing the results for the global descriptors in Table 1.

Having verified that the MN12SX/Def2TZVP model chemistry is a good choice for the calculation of the global reactivity descriptors, we now present in Figure 1 the molecular structure of Melanoidin M8 through a schematic representation generated through the Chemcraft molecular analysis program [44] , and the calculated bond lengths and bond angles are shown in Table 2 and Table 3.

The calculations of the condensed Fukui functions and Dual descriptor are done using the Chemcraft molecular analysis program to extract the Mulliken and NPA atomic charges [44] beginning with single-point energy calculations involving the MN12SX density functional that uses the Def2TZVP basis set in line with the SMD solvation model with water utilized as a solvent.

Considering the potential application of the Melanoidin M8 molecule as an antioxidant, it is of interest to gain insight into the active sites for radical attack. A graphical representation of the condensed radical Fukui function $f^0$ is presented in Figure 2.

A graphical representation of the Dual Descriptor $\Delta f\left(r\right)$ of the Melanoidin M8 molecule calculated with the MN12SX/Def2TZVP model chemistry in water is presented in Figure 3.

From the results for the Dual Descriptor $\Delta f\left(r\right)$ in Figure 3, it can be concluded that for the M8 melanoidin molecule, C11 will be the preferred site for a nucleophilic attack and that this atom will act as an electrophilic species in

Table 1. Global reactivity descriptors for the M8 intermediate melanoidin pigment calculated with the MN12SX density functional.

Table 2. Calculated bond lengths (in Å) of the M8 intermediate melanoidin pigment calculated with the MN12SX density functional.

Table 3. Calculated bond angles (in ˚) of the M8 intermediate melanoidin pigment calculated with the MN12SX density functional.

Figure 1. A schematic representation of the optimized structure of the M8 intermediate melanoidin pigment calculated with the MN12SX density functional in connection with the Def2TZVP basis set showing the numbering of the atoms.

Figure 2. A graphical representation of the radical Fukui function $f^0$ of the M8 intermediate melanoidin pigment in water where the red regions represent positive values of the magnitude and the blue zones represent the negative ones.

a chemical reaction. In turn, it can be appreciated that N1 and C3 will be prone to electrophilic attacks and that these atomic sites will act as nucleophilic species in chemical reactions where the M8 molecule is involved.

Finally, the condensed electrophilic and nucleophilic Parr functions $P_k^{+}$ and $P_k^{-}$ over the atoms of the Melanoidin M8 molecule in water have been calculated by extracting the Mulliken (or MPA) and Hirshfeld (or CM5) atomic charges using the Chemcraft molecular analysis program [44] starting from single-point energy calculations of the ionic species with the MN12SX density functional using the Def2TZVP basis set in the presence of the solvents according to the SMD solvation model. The maximum value of $P_k^{+}$ is located

Figure 3. A graphical representation of the Dual Descriptor $\Delta f\left(r\right)$ of the Melanoidin M8 molecule calculated with the MN12SX/Def2TZVP model chemistry in water where the red regions represent positive values of the magnitude and the blue zones represent the negative ones.

over C11 (0.4771 for MPA and 0.2479 for CM5) and the maximum value of $P_k^{-}$ is located over N1 (0.5011 for MPA and 0.3597 for CM5) and C3 (0.5158 for MPA and 0.3475 for HPA) for the calculation in the presence of water. Thus, there is a nice agreement between the results describing the local reactivity of the Melanoidin M8 molecule provided by Fukui functions, Dual Descriptor $\Delta f\left(r\right)$ (or its condensed counterpart) and the Parr functions.

4. Conclusions

It is not the objective of Computational Chemistry to reproduce the known results of experiments but to predict in advance the behavior of the molecular systems for those cases where the experimental results are unknown or very difficult to obtain. As a matter of fact, this is the case for the subject of this computational study. A powerful fixed RSH density functional, namely MN12SX, was examined to establish whether they fulfill the empirical KID procedure. The assessment was done by comparing the values from HOMO and LUMO calculations to those that the ΔSCF technique for the Melanoidin M8 molecule generates. It is an intermediate melanoidin pigment that is of both academic as well as industrial interest. This study observed that the range-separated hybrid meta-NGA MN12SX density functional tended to be the most suited to meeting this goal in agreement with our previous studies on the chemical reactivity of melanoidins. In this case, the MN12SX density functional emerged as a suitable alternative to those density functionals that are tuned using a gap-fitting procedure. It also exhibited desirable prospects of how it would benefit future studies in understanding the chemical reactivity of colored melanoidins with larger molecular weights interacting with reducing sugars react with proteins and peptides.

From the results of this work, it has become evident that it is easy to predict the sites of interaction of the melanodin M8 pigment under investigation. This would involve having DFT-based reactivity descriptors including Parr functions and Dual Descriptor calculations. Evidently, the descriptors were useful in characterizing and describing the preferred reactive sites. They were also useful in comprehensively explaining the reactivity of the molecule.

Acknowledgements

This work has been partially supported by CIMAV, SC and Consejo Nacional de Ciencia y Tecnologa (CONACYT, Mexico) through Grant 219566-2014 for Basic Science Research. Daniel Glossman-Mitnik conducted this work while a Visiting Lecturer at the University of the Balearic Islands from which support is gratefully acknowledged. This work was cofunded by the Ministerio de Economía y Competitividad (MINECO) and the European Fund for Regional Development (FEDER) (CTQ2014-55835-R).

Conflicts of Interest

The authors declare no conflicts of interest.

References

[1]	Nursten, H., Ed. (2005) The Maillard Reaction—Chemistry, Biochemistry and Implications. The Royal Society of Chemistry, Cambridge, UK.
[2]	Parr, R. and Yang, W. (1989) Density-Functional Theory of Atoms and Molecules. Oxford University Press, New York.
[3]	Mineva, T., Sicilia, E. and Russo, N. (1998) Density-Functional Approach to Hardness Evaluation and Its Use in the Study of the Maximum Hardness Principle, Journal of the American Chemical Society, 120, 9053-9058.
[4]	Mineva, T., Russo, N. and Sicilia, E. (1998) Solvation Effects on Reaction Profiles by the Polarizable Continuum Model Coupled with the Gaussian Density Functional Method. Journal of Computational Chemistry, 19, 290-299. https://doi.org/10.1002/(SICI)1096-987X(199802)19:3<290::AID-JCC3>3.0.CO;2-O
[5]	De Luca, G., Sicilia, E., Russo, N. and Mineva, T. (2002) On the Hardness Evaluation in Solvent for Neutral and Charged Systems. Journal of the American Chemical Society, 124, 1494-1499.
[6]	Sicilia, E., Russo, N. and Mineva, T. (2001) Correlation between Energy, Polarizability, and Hardness Profiles in the Isomerization Reaction of HNO and ClNO. The Journal of Physical Chemistry A, 105, 442-450. https://doi.org/10.1021/jp002350d
[7]	Frau, J. and Glossman-Mitnik, D. (2018) Molecular Reactivity and Absorption Properties of Melanoidin Blue-G1 through Conceptual DFT. Molecules, 23, 559-515. https://doi.org/10.3390/molecules23030559
[8]	Frau, J. and Glossman-Mitnik, D. (2018) Conceptual DFT Study of the Local Chemical Reactivity of the Dilysyldipyrrolones A and B Intermediate Melanoidins. Theoretical Chemistry Accounts, 137, 1210. https://doi.org/10.1007/s00214-018-2244-x
[9]	Frau, J. and Glossman-Mitnik, D. (2018) Conceptual DFT Study of the Local Chemical Reactivity of the Colored BISARG Melanoidin and Its Protonated Derivative. Frontiers in Chemistry, 6, 1-9.
[10]	Frau, J. and Glossman-Mitnik, D. (2018) Molecular Reactivity of some Maillard Reaction Products Studied through Conceptual DFT. Contemporary Chemistry, 1, 1-14.
[11]	Karolewski, A., Stein, T., Baer, R. and Kümmel, S. (2011) Communication: Tailoring the Optical Gap in Light-Harvesting Molecules. The Journal of Chemical Physics, 134, Article ID: 151101.
[12]	Karolewski, A., Kronik, L. and Kümmel, S. (2013) Using Optimally Tuned Range Separated Hybrid Functionals in Ground-State Calculations: Consequences and Caveats. The Journal of Chemical Physics, 138, Article ID: 204115. https://doi.org/10.1063/1.4807325
[13]	Koppen, J.V., Hapka, M., Szczeniak, M.M. and Chalasinski, G. (2012) Optical Absorption Spectra of Gold Clusters Au(n) (n = 4, 6, 8,12, 20) from Long-Range Corrected Functionals with Optimal Tuning. The Journal of Chemical Physics, 137, Article ID: 114302.
[14]	Kronik, L., Stein, T., Refaely-Abramson, S. and Baer, R. (2012) Excitation Gaps of Finite-Sized Systems from Optimally Tuned Range-Separated Hybrid Functionals. Journal of Chemical Theory and Computation, 8, 1515-1531.
[15]	Kuritz, N., Stein, T., Baer, R. and Kronik, L. (2011) Charge-Transfer-Like π → π* Excitations in Time-Dependent Density Functional Theory: A Conundrum and Its Solution. Journal of Chemical Theory and Computation, 7, 2408-2415. https://doi.org/10.1021/ct2002804
[16]	Refaely-Abramson, S., Baer, R. and Kronik, L. (2011) Fundamental and Excitation Gaps in Molecules of Relevance for Organic Photovoltaics From an Optimally Tuned Range-Separated Hybrid Functional. Physical Review B, 84, Article ID: 075144.
[17]	Stein, T., Kronik, L. and Baer, R. (2009) Prediction of Charge-Transfer Excitations in Coumarin-Based Dyes Using a Range-Separated Functional Tuned from First Principles. The Journal of Chemical Physics, 131, Article ID: 244119.
[18]	Stein, T., Kronik, L. and Baer, R. (2009) Reliable Prediction of Charge Transfer Excitations in Molecular Complexes Using Time-Dependent Density Functional Theory. Journal of the American Chemical Society, 131, 2818-2820.
[19]	Sun, H. and Autschbach, J. (2014) Electronic Energy Gaps for π-Conjugated Oligomers and Polymers Calculated with Density Functional Theory. Journal of Chemical Theory and Computation, 10, 1035-1047.
[20]	Parr, R. and Yang, W. (1984) Density Functional Approach to the Frontier-Electron Theory of Chemical Reactivity. Journal of the American Chemical Society, 106, 4049-4050.
[21]	Geerlings, P., De Proft, F. and Langenaeker, W. (2003) Conceptual Density Functional Theory. Chemical Reviews, 103, 1793-1873.
[22]	Parr, R., Szentpaly, L. and Liu, S. (1999) Electrophilicity Index. Journal of the American Chemical Society, 121, 1922-1924. https://doi.org/10.1021/ja983494x
[23]	Gázquez, J., Cedillo, A. and Vela, A. (2007) Electrodonating and Electroaccepting Powers. Journal of Physical Chemistry A, 111, 1966-1970.
[24]	Chattaraj, P., Chakraborty, A. and Giri, S. (2009) Net Electrophilicity. Journal of Physical Chemistry A, 113, 10068-10074.
[25]	Morell, C., Grand, A. and Toro-Labbé, A. (2005) New Dual Descriptor for Chemical Reactivity. Journal of Physical Chemistry A, 109, 205-212.
[26]	Morell, C., Grand, A. and Toro-Labbé, A. (2006) Theoretical Support for Using the Descriptor. Chemical Physics Letters, 425, 342-346. https://doi.org/10.1016/j.cplett.2006.05.003
[27]	Cárdenas, C., Rabi, N., Ayers, P., Morell, C., Jaramillo, P. and Fuentealba, P. (2009) Chemical Reactivity Descriptors for Ambiphilic Reagents: Dual Descriptor, Local Hypersoftness, and Electrostatic Potential. Journal of Physical Chemistry A, 113, 8660-8667.
[28]	Toro-Labbé, A. (2007) Theoretical Aspects of Chemical Reactivity. Elsevier Science, Amsterdam.
[29]	Ayers, P., Morell, C., De Proft, F. and Geerlings, P. (2007) Understanding the Woodward-Hoffmann Rules by Using Changes in Electron Density. Chemistry—A European Journal, 13, 8240-8247.
[30]	Morell, C., Ayers, P., Grand, A., Gutiérrez-Oliva, S. and Toro-Labbé, A. (2008) Rationalization of the Diels-Alder Reactions through the Use of the Dual Reactivity Descriptor . Physical Chemistry Chemical Physics, 10, 7239-7246.
[31]	Morell, C., Hocquet, A., Grand, A. and Jamart-Grégoire, B. (2008) A Conceptual DFT Study of Hydrazino Peptides: Assessment of the Nucleophilicity of the Nitrogen Atoms by Means of the Dual Descriptor . Journal of Molecular Structure: THEOCHEM, 849, 46-51.
[32]	Domingo, L.R., Pérez, P. and Sáez, J. (2013) Understanding the Local Reactivity in Polar Organic Reactions through Electrophilic and Nucleophilic Parr Functions. RSC Advances, 3, 1486-1494.
[33]	Chamorro, E., Pérez, P. and Domingo, L.R. (2013) On the Nature of Parr Functions to Predict the Most Reactive Sites along Organic Polar Reactions. Chemical Physics Letters, 582, 141-143.
[34]	Domingo, L.R., Ríos-Gutiérrez, M. and Pérez, P. (2016) Applications of the Conceptual Density Functional Theory Indices to Organic Chemistry Reactivity. Molecules, 21, 748.
[35]	Frisch, M.J., Trucks, G.W., Schlegel, H.B., Scuseria, G.E., Robb, M.A., Cheeseman, J.R., Scalmani, G., Barone, V., Mennucci, B., Petersson, G.A., Nakatsuji, H., Caricato, M., Li, X., Hratchian, H.P., Izmaylov, A.F., Bloino, J., Zheng, G., Sonnenberg, J.L., Hada, M., Ehara, M., Toyota, K., Fukuda, R., Hasegawa, J., Ishida, M., Nakajima, T., Honda, Y., Kitao, O., Nakai, H., Vreven, T., Montgomery, J.A., Peralta, J.E., Ogliaro, F., Bearpark, M., Heyd, J.J., Brothers, E., Kudin, K.N., Staroverov, V.N., Kobayashi, R., Normand, J., Raghavachari, K., Rendell, A., Burant, J.C., Iyengar, S.S., Tomasi, J., Cossi, M., Rega, N., Millam, J.M., Klene, M., Knox, J.E., Cross, J.B., Bakken, V., Adamo, C., Jaramillo, J., Gomperts, R., Stratmann, R.E., Yazyev, O., Austin, A.J., Cammi, R., Pomelli, C., Ochterski, J.W., Martin, R.L., Morokuma, K., Zakrzewski, V.G., Voth, G.A., Salvador, P., Dannenberg, J.J., Dapprich, S., Daniels, A.D., Farkas, O., Foresman, J.B., Ortiz, J.V., Cioslowski, J. and Fox, D.J. (2018) Gaussian 09 Revision D.01. Gaussian Inc., Wallingford.
[36]	Weigend, F. and Ahlrichs, R. (2005) Balanced Basis Sets of Split Valence, Triple Zeta Valence and Quadruple Zeta Valence Quality for H to Rn: Design and Assessment of Accuracy. Physical Chemistry Chemical Physics, 7, 3297-3305. https://doi.org/10.1039/b508541a
[37]	Weigend, F. (2006) Accurate Coulomb-Fitting Basis Sets for H to R. Physical Chemistry Chemical Physics, 8, 1057-1065.
[38]	Marenich, A., Cramer, C. and Truhlar, D. (2009) Universal Solvation Model Based on Solute Electron Density and a Continuum Model of the Solvent Defined by the Bulk Dielectric Constant and Atomic Surface Tensions. Journal of Physical Chemistry B, 113, 6378-6396.
[39]	Halgren, T.A. (1996) Merck Molecular Force Field. I. Basis, Form, Scope, Parameterization, and Performance of MMFF94. Journal of Computational Chemistry, 17, 490-519.
[40]	Halgren, T.A. (1996) Merck Molecular Force Field. II. MMFF94 van der Waals and Electrostatic Parameters for Intermolecular Interactions. Journal of Computational Chemistry, 17, 520-552.
[41]	Halgren, T.A. (1996) MMFF VI. MMFF94s Option for Energy Minimization Studies. Journal of Computational Chemistry, 20, 720-729.
[42]	Halgren, T.A. and Nachbar, R.B. (1996) Merck Molecular Force Field. IV. Conformational Energies and Geometries for MMFF94. Journal of Computational Chemistry, 17, 587-615.
[43]	Halgren, T.A. (1996) Merck Molecular Force field. V. Extension of MMFF94 Using Experimental Data, Additional Computational Data, and Empirical Rules. Journal of Computational Chemistry, 17, 616-641.
[44]	Zhurko, G. and Zhurko, D. (2012) Chemcraft Program Revision 1.6. http://www.chemcraft.com/