Kimball Model of EthylCyanoAcetate

We will write a simple Mathematica version 10 code to study Ethylcyanoacetate.Our object is to compute its energy components 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 an HTML copy of the original Mathematica Notebook EtCyanoAc_n1.nb, which can be downloaded and interactively used.

(ES 27 March 2013/26 Sept. 2014)

Input and Definitions

The stored coordinates are read from the Wolfram Mathematica ChemicalData repository (coordinates in pm = picometer) into the t1 array. 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 = 2’625.50 kJ/mol. Other units are indicated explicitely.

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

Compute Kimball radii from distance matrix

Select sigma 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 some arrays, to be used later

Compute kinetic energy of sigma 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 two subroutines PItrans.m, and PICNtrans.m, that can be downloaded from the URL’s given below for reading them in.

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. Some such transformations shown below for perpendicular and in-plane oxygen lone pairs. One side of the triangle is parallel to x-axis, the target atom at origin. Subroutines: NCNtrans.m, ODtrans.m, and OBWtrans.m.

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

Plot of EthylCyanoacetate and its partial constituents

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

Compute energy components (check charges):

Politzerratio

Results (energies in [Eh] Hartree)