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

Lorazepam

All computations are transparent and annotated. The run lasts about 1 sec on a i7-4690 CPU. (ES 16 June 2013/ 17 December 2016).

Structure from Wolfram ChemData.

Input and Definitions

The coordinates are read in pm. 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

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 frame.

σ 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) (RHF/6-31G(d) opt for comparison)

Hellmann-Feynman force analysis

The force vectors are shown below as “force ellipsoids”. If these forces really existed, the molecule would be in an excited vibrational state and vibrate causing “vibrational ellipsoids” proportional to those pictured, similar to “thermal ellipsoids” shown in Tutorial6a. It is easy to see that these ellipsoids do not exist but are artifacts of the computation to be remedied below. If any ground state molecule had permanent, non vanishing forces to move nuclei, it would spontaneously vibrate much more vigorously than what happens in thermal equilibrium. Hence, this is a perfect demonstration of a perpetuum mobile, cooling the environment by driving the molecular “machine” which could produce work and radiate infrared light. The artifact of the computation is a simplification: we have placed all heavy nuclei in the center of their core clouds. This is only correct if the core and its nucleus are in a centrosymmetric location like C[1s2] in methane. In all other environments location of nucleus and center of core cloud are not exactly coincident! Hence, the Hellmann-Feynman electrostatic theorem cannot be demonstrated without core polarization, neither in Kimballs model nor in any QC computation. That is, of course, well known, but barely mentioned in any chemistry text (except QC). Often a second correction is necessary: Peaks of internuclear electron density are not on the straight line connecting two nuclei (‘bent bonds’). See Didactics issues where this is looked at in the context of ‘ring strain’ and its release within Kimballs model.

For “heavy” nuclei a tiny shift (< 0.0001 to 0.001 a0) of the core cloud is computed. This separates its center from that of the nucleus and thus introduces a polarization force. If applied in the reverse direction that balances the residual electrostatic force. In this notebook all nuclei except protons sit in the center of their core clouds. Introducing polarization does not change the energy terms by more than 0.001 to 0.01 Eh (or 0.1 Eh in large molecules like crambin) and, hence, is usually neglected. Next we apply this correction to our computation if the coincidence is not required by symmetry. Compare the plots and error vectors below!

Residual HF forces after polarization of the core clouds:

Polarization vectors for every heavy nucleus, i.e. radius vector of every shifted nucleus inside its core cloud, x,y,z components shown. This balances exactly the artificial residual HF-force for C1 to O21, see table above

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

We have not corrected the small HF-forces remaining for H nuclei. They could arise from small errors in the experimental structure. Proton locations are the least exact in X-ray crystallographic structure determinations. Furthermore, a second order correction for proton locations in an X-H bond is again required in our computation, if protons do not sit on the straight line X - center bonding cloud - proton, which is assumed. This is only correct if the three players [X, center X-H cloud, H] are in a cylinder symmetric surrounding around the X-H axis. Because of the small nuclear charge the effect of a lower symmetry does not count at all in the energy computation. For the HF analysis we let the errors stand uncorrected!