Quantitative electronic Lewis structure derived from nuclear coordinates of a molecule:

Crambin

We will write a simple Mathematica version 10 code to study the small protein “Crambin”, , consisting of 46 amino acids, 644 atoms.

Sequence: TTCCPSIVARSNFNVCRLPGTPEALCATYTGCIIIPGATCPGDYAN. The three carboxyl groups of glutamic, aspartic acids as well as the C-end are protonated. Hence we have a neutral molecule with an even number of electrons.

Our object is to compute energy components of crambin, given a 3D structure from 1CRN.pdb, and to plot the different electronic constituents.

Crambin’s structure was first established by:

M.M.TEETER,W.A.HENDRICKSON, “HIGHLY ORDERED CRYSTALS OF THE PLANT SEED PROTEIN CRAMBIN”, J.MOL.BIOL.V.127,219,1979

Several corrections have been applied by the same and other authors, here included to 16-APR-87.

In 2013, J.Chem.Phys. 139,134101, Ch.Riplinger, B. Sandhoefer, A. Hansen, and F. Neese, describe the first successful quantum chemical treatment of crambin with the high accuracy coupled cluster method with double and triple excitations: CCSD(T). This was only possible by their invention of a nearly linear scaling procedure called DLPNO-CCSD(T). The run of crambin lasted 30 days on a Linux system on 4 octo-core Intel Xeon CPU’s with 256 GB main memory, using one core only! Without their development the normal, “canonical”, CCSD(T) computation would have lasted several thousand years. We use the coordinates of these authors, published as Supplement to the cited paper, based on the structure by Teeter et al., but protonated.

Here we apply Kimball’s model to get a rough idea about energy components and electron distribution of crambin. A complete Hellmann-Feynman electrostatic analysis is included at the end of this notebook.

All computations are transparent and annotated. The run lasts about 31 sec on a i7-4690/3.6GHz CPU, and 55 sec with added Hellmann-Feynman analysis. (ES 24 November 2017/16 June 2013/20 June 2017).

Input and Definitions

The coordinates are read in Å. We are using atomic units, the universally applied system of theoretical chemistry and (micro) physics, see NIST. Length data are in Bohr : 1 a0 = 0.52917721 Å = 52.917721 pm; electric charges in ± electron charges, and energies in Hartrees : 1 Eh = 2 Rydberg = 627.5095 kcal/mol = 2625.50 kJ/mol.

Normal Input for a structure given as a table with rows: |Atom_symbol x y z|

Normal Input for a structure downloaded from Wolfram ChemData repository

Analyze the atomic constituents

Unique nuclei | C | H | N | O | S |

··· number | 202 | 317 | 55 | 64 | 6 |

··· nuclear charge | 6 | 1 | 7 | 8 | 16 |

··· core cloud charge | -2 | 0 | -2 | -2 | -10 |

··· valence electrons | 4 | 1 | 5 | 6 | 6 |

··· full Lewis shell | 8 | 2 | 8 | 8 | 8 |

··· Core radii | 0.252 | 0 | 0.212 | 0.178 | 0.579 |

Analyze Lewis structure

Compute Kimball radii from distance matrix, show core radii derived from CH4, NH3, H2O gauge molecules (cnofhydb.pas), (cnofhydb.ex_ to be renamed into runnable cnofhydb.exe after download), H eccentricities, and number of σ bonds.

Distance Matrix :

Nuclear repulsion

Determine bonded pairs by a distance criterion

Bonded atom pairs: distances

Subtract proton eccentricities

Subtract core radii

Show radii determined

Summary of Lewis properties

Compute kinetic energy terms, bonding clouds, core clouds:

Total kinetic energy except for π - clouds and lone pairs

Determine connectivity matrix:

Localize double bonds and positions of π-clouds (PItrans.m)

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 back transform the π-clouds into the molecular coordinate array.

Localize lone pairs, compute size and orientation:

Subroutines: XOtrans.m XOYtrans.m CNCtrans.m LpyrNtrans.m

Transform the triangle of every target atom with two of its neighbors into the xy-plane and attach lone pair(s). Then

back transform the lone pair(s) into the molecular coordinate array. See one of the subroutines. LpyrNtrans puts the

base atoms of a pyramid into the xy plane and attaches LP’s as needed, then moves these back into the molecule.

σ Bonding clouds: Connected atom pair, radius of cloud

Plot molecule and its electronic partial constituents

Add coordinates of π-clouds and lone pairs. Prepare interaction matrices:

Compute energy components

Interactions for i not j

Interactions for i equals j

Kinetic energy of π clouds and lone pairs

Add components of Ne[10] cores; Politzer ratio

Results (energies in [Eh] Hartree)

Hellman-Feynman electrostatic forces

Compute core polarizations to balance remaining HF-forces

Total polarization energy in Eh. It is to be added to Etot of the crambin molecule, but is usually neglected in regard to its small size compared to Etot ~ -18400 Eh.