Kimball Model of EthylAcetoacetate

We will write a simple Mathematica version 10 code to study EthylAcetoacetate. Our object is to compute the energy components of its groundstate, given a 3D structure from the PDB library, and to plot the different electronic constituents. All computations are transparent and operations explained in detail. This is a HTML copy of the interactive notebook.(ES 27 March 2013/05 October. 2014)

Input and Definitions

The coordinates are read from the Wolfram Mathematica ChemicalData repository (in pm = picometer). We use atomic units here, the universally applied system of theoretical chemistry and (micro)physics, see NIST. Length data are in Bohr: 1 a0 = 0.52917721 Å = 52.917721 pm (whose reciprocal is 0.01889726 bohr/pm, used below), electric charges in ± electron charges, and energies in Hartrees: 1 Eh = 2 Rydberg = 627.5095 kcal/mol = 2625.50 kJ/mol.

Sort array of atoms, show nuclear charges, number of “heavy” atoms, number of hydrogens

Compute Kimball radii from distance matrix

Select σ bonded pairs by a distance criterion

Build vector of core radii and proton excentricities from CH4, NH3, H2O gauge molecules (cnofhydb.pas)

From lower triangle of symmetrical matrix of distances subtract vector of core radii

Transpose to upper triangular matrix, show, then subtract vector of core radii from the other side, which leaves diameter of bonding clouds between “heavy” atoms and radius of H-clouds:

Determine total number of bonds and radii of bonding clouds between “heavy” atoms

Form an array, to be used later

Compute kinetic energy terms of σ bonding and core clouds. Kimballs ansatz 1.125/R^{2} per electron.

Total kinetic energy except for π-clouds and lone pairs

Determine connectivity matrix:

Localize double bonds and positions of π-clouds:

Transform the triangle of every target atom with two of its neighbors into the xy-plane and attach π-clouds above and below the plane to the target. Then backtransform the triangle plus π-clouds into the molecular coordinate array. This is done by the subroutine PItrans.m, that can be downloaded from the URL given below for reading it in. This subroutine is very fast!

Localize lone pairs, compute orientation

Transform the triangle of every target atom with two of its neighbors into the xy-plane and attach lone pair(s). Then backtransform the triangle plus lone pair(s) into molecular coordinate array. One side of the triangle is parallel to x-axis, the target atom at origin. Subroutines: XOtrans.m, and XOYtrans.m are very fast!

Bonding clouds: Atom pair, radius of cloud; statistics of sizes

Plot of Ethylacetoacetate and its partial constituents

Preparation of interaction matrices (w-w, w-n, check size arrays):

Computation of energy components (check charges):

Politzerratio

Results (energies in [Eh] Hartree)