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”, Crambin_HFsipi6b_1.png, 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|

Crambin_HFsipi6b_2.gif

Crambin_HFsipi6b_3.png

Normal Input for a structure downloaded from Wolfram ChemData repository

Crambin_HFsipi6b_4.png

Analyze the atomic constituents

Crambin_HFsipi6b_5.png

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

Crambin_HFsipi6b_6.png

Crambin_HFsipi6b_7.png

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 :

Crambin_HFsipi6b_8.gif

Nuclear repulsion

Crambin_HFsipi6b_9.gif

Crambin_HFsipi6b_10.png

Determine bonded pairs by a distance criterion

Crambin_HFsipi6b_11.gif

Crambin_HFsipi6b_12.png

Bonded atom pairs:  distances

Crambin_HFsipi6b_13.png

Crambin_HFsipi6b_14.png

Subtract proton eccentricities

Crambin_HFsipi6b_15.png

Crambin_HFsipi6b_16.png

Crambin_HFsipi6b_17.gif

Crambin_HFsipi6b_18.png

Subtract core radii

Crambin_HFsipi6b_19.gif

Crambin_HFsipi6b_20.png

Show radii determined

Crambin_HFsipi6b_21.gif

Crambin_HFsipi6b_22.png

Summary of Lewis properties

Crambin_HFsipi6b_23.gif

Crambin_HFsipi6b_24.png

Crambin_HFsipi6b_25.png

Crambin_HFsipi6b_26.png

Compute kinetic energy terms, bonding clouds, core clouds:

Crambin_HFsipi6b_27.gif

Crambin_HFsipi6b_28.png

Crambin_HFsipi6b_29.png

Crambin_HFsipi6b_30.png

Total kinetic energy except for π - clouds and lone pairs

Crambin_HFsipi6b_31.png

Crambin_HFsipi6b_32.png

Crambin_HFsipi6b_33.gif

Crambin_HFsipi6b_34.png

Crambin_HFsipi6b_35.png

Determine connectivity matrix:

Crambin_HFsipi6b_36.gif

Crambin_HFsipi6b_37.png

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.

Crambin_HFsipi6b_38.png

Crambin_HFsipi6b_39.png

Crambin_HFsipi6b_40.gif

Crambin_HFsipi6b_41.png

Crambin_HFsipi6b_42.png

Crambin_HFsipi6b_43.png

Crambin_HFsipi6b_44.png

Crambin_HFsipi6b_45.png

Crambin_HFsipi6b_46.png

Crambin_HFsipi6b_47.png

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.

Crambin_HFsipi6b_48.png

Crambin_HFsipi6b_49.png

Crambin_HFsipi6b_50.png

Crambin_HFsipi6b_51.png

Crambin_HFsipi6b_52.gif

Crambin_HFsipi6b_53.png

Crambin_HFsipi6b_54.png

Crambin_HFsipi6b_55.png

Crambin_HFsipi6b_56.png

Crambin_HFsipi6b_57.png

Crambin_HFsipi6b_58.png

Crambin_HFsipi6b_59.png

σ Bonding clouds: Connected atom pair, radius of cloud

Crambin_HFsipi6b_60.gif

Crambin_HFsipi6b_61.png

Crambin_HFsipi6b_62.png

Crambin_HFsipi6b_63.png

Crambin_HFsipi6b_64.png

Crambin_HFsipi6b_65.png

Crambin_HFsipi6b_66.png

Graphics:Cloud radii (bohr)

Plot molecule and its electronic partial constituents

Crambin_HFsipi6b_68.gif

Crambin_HFsipi6b_69.gif

Graphics:Core skeleton

Graphics:σ skeleton

Graphics:π-clouds on skeleton

Graphics:Lone Pairs

Graphics:H atoms

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

Crambin_HFsipi6b_75.png

Crambin_HFsipi6b_76.png

Crambin_HFsipi6b_77.png

Compute energy components

Interactions for i not j

Crambin_HFsipi6b_78.gif

Interactions for i equals j

Crambin_HFsipi6b_79.png

Kinetic energy of π clouds and lone pairs

Crambin_HFsipi6b_80.png

Crambin_HFsipi6b_81.png

Crambin_HFsipi6b_82.png

Add components of Ne[10] cores; Politzer ratio

Crambin_HFsipi6b_83.png

Crambin_HFsipi6b_84.png

Results (energies in [Eh] Hartree)

Crambin_HFsipi6b_85.png

Crambin_HFsipi6b_86.png

Crambin_HFsipi6b_87.png

Crambin_HFsipi6b_88.png

Crambin_HFsipi6b_89.png

Crambin_HFsipi6b_90.png

Crambin_HFsipi6b_91.png

Crambin_HFsipi6b_92.png

Crambin_HFsipi6b_93.png

Crambin_HFsipi6b_94.png

Crambin_HFsipi6b_95.png

Crambin_HFsipi6b_96.png

Crambin_HFsipi6b_97.png

Crambin_HFsipi6b_98.png

Hellman-Feynman electrostatic forces

Crambin_HFsipi6b_99.gif

Crambin_HFsipi6b_100.png

Graphics:Residual Hellmann-Feynman force ellipsoids, core polarization neglected

Compute core polarizations to balance remaining HF-forces

Crambin_HFsipi6b_102.png

Graphics:Residual Hellmann-Feynman force ellipsoids (magenta) core polarization included

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.

Crambin_HFsipi6b_104.png

Crambin_HFsipi6b_105.png

Created with the Wolfram Language