A data driven, Quantum Mechanical understanding of Chemistry

AKA Trying to make sense of 134k quantum calculations

by Benjamin Sanchez Lengeling (bsanchezlengeling@g.harvard.edu)

Load custom CSS styling


In [10]:
from IPython.core.display import HTML
import os
def css_styling():
    """Load default custom.css file from ipython profile"""
    base = os.getcwd()
    styles = "<style>\n%s\n</style>" % (open(os.path.join(base,'custom.css'),'r').read())
    return HTML(styles)
css_styling()


Out[10]:

Molecular party!

0 Intro

0.1 Motivation

This CS 209 project is an exploration of Chemistry using data science techniques. I want to use the methods learned in the class to apply them in my own context, which is Quantum Chemistry.

Quantum chemistry is a branch of theoretical chemistry which applies quantum mechanics to address problems in chemistry.

In Quantum Chemistry you can calculate a molecule on a computer, simulating and solving the Schrödinger equation to obtain molecular properties.

Some pretty awesome examples of Quantum Chemistry applications are:

  • Designing and finding better molecules and dyes, to construct more efficient solar cells
  • Finding out why things have the color they do with no previous knowledge of the world, for example why is a tomato red? You can simulate a photon (light) hitting the molecules on the frontier of a tomato and seeing what type of wavelength you get back.
  • Studying possible new materials like graphene, or carbon nanotubes.
  • Designing better of newer drugs by modeling the molecules when they come in "contact" with other molecules in our bodies.
  • Finding alternative reactions for common industrial products, typically you want these reactions to be more ecological (better byproducts), cheaper or more efficient.
  • Studying the composition of molecules and reactions that can occur easily in space, to better understand all the material world surrounding us right now, including life.

Due to advances in computing now we can calculate a large number of molecules in relatively small time frame, creating huge masses of data.

0.1 Motivation

The pharmaceutical industry has established for a long time now the practice of using molecular calculations to drive drug discovery. This has olny been done with classical mechanics and with very limited quantum treatment of molecules. There is now a trend of trying to adapt the tools and strategies of this industry to material science.

So a current open question is: How can we incorporate quantum effects/information into chemoinformatics?

Very little work has been done on this aspect, in part due to the computational constraint of the treatment of quantum effects and also due to the complexity of this information.

I want to put my newly earned skills to the test, on data I will most probably be working with during my graduate studies.

One aspect of this work is the research, there are two main projects/works that have tried to utilize data science in Quantum Chemistry:

  • The Harvard Clean Energy project, the largest database of quantum calculation in the world, all seeking to discover the new generation of plastic solar cell materials. http://cleanenergy.molecularspace.org/ , http://www.molecularspace.org/.
  • A recent 2014 Scientific Data article "Quantum chemistry structures and properties of 134 kilo molecules" (http://www.nature.com/articles/sdata201422). One of the first projects to grant access to a large database of quantum calculations. The article sets-up the rational behind the creation of the dataset, they plan to use machine learning to predict molecular properties.

The other aspect relates to accesibility in science, in using great tools as ipython notebooks, data repositories, open source libraries to accelerate and make more open the process of doing science. In Chemistry two side projects which are of great interest are:

  • OpenBabel (http://openbabel.org): library designed to support molecular modeling, chemistry, and many related areas, including interconversion of file formats and data
  • RDkit (http://www.rdkit.org/): library for Cheminformatics and Machine Learning.

0.3 Goals

  1. Cleanup the dataset and present it in a more user and sharable format (pandas dataframe, pickle, excel, json database).
  2. Visualize the dataset, explore all parameters of the data and present global trends.
  3. Find patterns within the data, how do these patterns reflect with chemistry concepts?

1 Data Wrangling

1.0 The Dataset

Dataset is hosted on http://dx.doi.org/10.6084/m9.figshare.978904, it is a few hundred MB's.

I also have pandas dataframe pickle in my dropbox:

https://www.dropbox.com/sh/6iysi4w0xmmevlt/AABrTLUFZJvrJPDeDCzuhxJYa?dl=0.

They reported computed geometric, energetic, electronic, and thermodynamic properties for 134k stable small organic molecules made up of CHONF atoms. These molecules correspond to the subset of all 133,885 species with up to nine heavy atoms (CONF) out of the GDB-17 chemical universe of 166 billion organic molecules.

All calculations have been carried out at the B3LYP/6-31G(2df,p) level of quantum chemistry. (Decent, but not the best)

Results of each calculation has been embeded within a xyz file in a very specific way. The data provided by the article is not very user friendly ( see: http://www.ch.imperial.ac.uk/rzepa/blog/?p=12803), a better format would have been a json or excel database.

We will parse out the information acoording to the readme.txt file and create for each molecule a row within a pandas dataframe. This dataframe can thn easily be exported to a variety of popular data formats.

1.1 Libraries


In [1]:
#Libraries
import numpy as np
import scipy as sp
import pandas as pd
import sklearn
import seaborn as sns
from matplotlib import pyplot as plt
import random
#chemo-informatics
import imolecule
import openbabel
import pybel
import rdkit
# visual stuff
%matplotlib inline
from IPython.core import display
from PIL import Image
#rdkit
from rdkit import Chem
from rdkit.Chem import Draw
from rdkit.Chem import AllChem
import rdkit.rdBase
from rdkit.Chem.MACCSkeys import GenMACCSKeys
from rdkit import DataStructs
from rdkit.DataStructs import BitVectToText
from rdkit.Chem.Draw import IPythonConsole
from rdkit.Chem import Draw
from rdkit import DataStructs
from rdkit.Chem import MCS as MCS
from rdkit.Chem import Descriptors as Descriptors
from rdkit.ML.Descriptors import MoleculeDescriptors
from rdkit.Chem import PandasTools as PandasTools
from rdkit.Chem import Descriptors as Descriptors
#awesome plot options
plt.style.use('fivethirtyeight')
plt.rcParams['figure.figsize'] = (12.0, 6.0)
# utility
from io import BytesIO


#functions for ipython

def display_pil_image(im):
    """Displayhook function for PIL Images, rendered as PNG."""

    b = BytesIO()
    im.save(b, format='png')
    data = b.getvalue()

    ip_img = display.Image(data=data, format='png', embed=True)
    return ip_img._repr_png_()

# register display func with PNG formatter:
png_formatter = get_ipython().display_formatter.formatters['image/png']
dpi = png_formatter.for_type(Image.Image, display_pil_image)

#functions for plotting
def sns_cycle_palette():
    """Cycle seaborn's color palete to have a variety of colors"""

    pal = sns.color_palette()
    pal=pal[1:]+[pal[0]]
    sns.set_palette(pal)
    return

1.2 Read XYZ files and load to DataFrame

We create a simple function to read an xyz file and return a row of data to be appended to a dataframe. The we iterate over all the xyz files, this took a few hours. We save the dataframe periodically unto a python pickle.


In [7]:
def readFormatedXYZ(xyzfile):
    fileData=[]
    with(open(xyzFile,'r')) as afile:
        numAtoms=int(afile.next())
        #1 num atoms
        fileData.append(numAtoms)
        datline=afile.next().split()
        # 2 molecular index within the database
        fileData.append(int(datline[1]))
        # 3-5 rotational constants
        fileData.append(float(datline[2]))
        fileData.append(float(datline[3]))
        fileData.append(float(datline[4]))
        # 6 dipole moment
        fileData.append(float(datline[5]))
        # 7 polarizability
        fileData.append(float(datline[6]))
        # 8-10 orbital energies (lumo,homo)
        fileData.append(float(datline[7]))
        fileData.append(float(datline[8]))
        fileData.append(float(datline[9]))
        # 11 spatial extent
        fileData.append(float(datline[10]))
        # 12 zpe
        fileData.append(float(datline[11]))
        # 13-14 Internal Energy
        fileData.append(float(datline[12]))
        fileData.append(float(datline[13]))
        # 15 Enthalpy Energy
        fileData.append(float(datline[14]))
        # 16 Free Energy
        fileData.append(float(datline[15]))
        # 17 Heat Capacity
        fileData.append(float(datline[16]))
        
        #read geometry
        atomCoords=np.zeros((numAtoms,3))
        atomList=[]
        for i in range(numAtoms):
            datline=afile.next().split()
            atomList.append(datline[0])
            atomCoords[i]=np.array([float(d) for d in datline[1:4]])
        fileData.append(atomList)
        fileData.append(atomCoords)
        # next frequencies
        fileData.append(np.array([float(f) for f in afile.next().split()]))
        #relaxed geometry smile
        fileData.append(afile.next().split()[1])
        #relaxed geometry InChI
        fileData.append( afile.next().split()[1])
        
    return fileData

In [61]:
#===== Initial Variables ===
molFolder = 'mols/'
rmin = 0
rmax = 100
pandasFile = "mols.pkl"
#===== Check that the pickle file exists ===
if os.path.isfile(pandasFile):
    masterdf = pd.read_pickle(pandasFile)
    rmin = masterdf.shape[0] + 1
else:
    datCols = ['numAtoms', 'dbindex', 'A', 'B', 'C', 'dipole', 'polar',
               'homo', 'lumo', 'gap', 'spatialSize', 'zpe',
               'U0', 'U', 'H', 'G', 'Cv',
               'atomList', 'atomCoords', 'freqs',
               'SMILES', 'InChI']
    masterdf = pd.DataFrame(columns=datCols)
    rmin = 0
    
masterdf=pd.DataFrame( columns=datCols)
#===== Iterate over file
for indx in indexes:
    xyzFile="mols/dsgdb9nsd_%06d.xyz"%(indx+1)
    masterdf.loc[indx]=readFormatedXYZ(xyzFile)
    if (indx % 10000 == 0):
        print(indx)
        masterdf.to_pickle(pandasFile)
#====== Save and
masterdf.to_pickle(pandasFile)

In [6]:
#pickle it!
pandasFile = "original.pkl"
masterdf = pd.read_pickle(pandasFile)
print(masterdf.columns)
masterdf.head()


Index([u'numAtoms', u'dbindex', u'A', u'B', u'C', u'dipole', u'polar', u'homo', u'lumo', u'gap', u'spatialSize', u'zpe', u'U0', u'U', u'H', u'G', u'Cv', u'atomList', u'atomCoords', u'freqs', u'SMILES', u'InChI'], dtype='object')
Out[6]:
numAtoms dbindex A B C dipole polar homo lumo gap ... U0 U H G Cv atomList atomCoords freqs SMILES InChI
0 5 1 157.71180 157.709970 157.706990 0.0000 13.21 -0.3877 0.1171 0.5048 ... -40.478930 -40.476062 -40.475117 -40.498597 6.469 [C, H, H, H, H] [[-0.0126981359, 1.0858041578, 0.0080009958], ... [1341.307, 1341.3284, 1341.365, 1562.6731, 156... C InChI=1S/CH4/h1H4
1 4 2 293.60975 293.541110 191.393970 1.6256 9.46 -0.2570 0.0829 0.3399 ... -56.525887 -56.523026 -56.522082 -56.544961 6.316 [N, H, H, H] [[-0.0404260543, 1.0241077531, 0.0625637998], ... [1103.8733, 1684.1158, 1684.3072, 3458.7145, 3... N InChI=1S/H3N/h1H3
2 3 3 799.58812 437.903860 282.945450 1.8511 6.31 -0.2928 0.0687 0.3615 ... -76.404702 -76.401867 -76.400922 -76.422349 6.002 [O, H, H] [[-0.0343604951, 0.9775395708, 0.0076015923], ... [1671.4222, 3803.6305, 3907.698] O InChI=1S/H2O/h1H2
3 4 4 0.00000 35.610036 35.610036 0.0000 16.28 -0.2845 0.0506 0.3351 ... -77.308427 -77.305527 -77.304583 -77.327429 8.574 [C, C, H, H] [[0.5995394918, 0.0, 1.0], [-0.5995394918, 0.0... [549.7648, 549.7648, 795.2713, 795.2713, 2078.... C#C InChI=1S/C2H2/c1-2/h1-2H
4 3 5 0.00000 44.593883 44.593883 2.8937 12.99 -0.3604 0.0191 0.3796 ... -93.411888 -93.409370 -93.408425 -93.431246 6.278 [C, N, H] [[-0.0133239314, 1.1324657151, 0.0082758861], ... [799.0101, 799.0101, 2198.4393, 3490.3686] C#N InChI=1S/CHN/c1-2/h1H

5 rows × 22 columns

1.3 Remove uncharacterized molecules

Within the article they note that 3k molecules did not obtain converged geoemtries consistant with the intial guess, i.e. radically different geometries. We filter these out.


In [76]:
#preprocessing the dataframes
df=masterdf.copy()
#Filter out the uncharacterized molecules
badmols=pd.read_csv('uncharacterized.csv')
df=df[~df['dbindex'].isin(badmols['Index'])]
print("Bad molecules: %s"%(str(badmols.shape)))
print("Before %s -> After %s "%(masterdf.shape[0],df.shape[0]))


Bad molecules: (3054, 5)
Before 133884 -> After 130830 

1.4 Add other computable descriptors

We add other descriptors to the dataframe. We create a row operating function and use pd.apply().

1.4.1 Molecular Weight (molWeight)


In [77]:
#load atomic data
atomdf=pd.read_csv('atomref.csv')
weightDict={ row['Element']:row['Mass'] for indx,row in atomdf.iterrows()}
print(weightDict)
#molecular weight function
def molecularWeight( row ):
    weightDict={'H': 1.00794, 'C': 12.0107, 'F': 18.998403, 'O': 15.9994, 'N': 14.0067}
    return np.sum([ weightDict[atom] for atom in row['atomList'] ])
#apply function to each row
df['molWeight']=df.apply(molecularWeight, axis=1)
atomdf.head()


{'H': 1.00794, 'C': 12.0107, 'F': 18.998403, 'O': 15.9994, 'N': 14.0067}
Out[77]:
Element Mass ZPVE_Hartree U (0 K)_Hartree U (298.15 K)_Hartree H (298.15 K)_Hartree G (298.15 K)_Hartree CV_Hartree_Cal/(Mol Kelvin)
0 H 1.007940 0 -0.500273 -0.498857 -0.497912 -0.510927 2.981
1 C 12.010700 0 -37.846772 -37.845355 -37.844411 -37.861317 2.981
2 N 14.006700 0 -54.583861 -54.582445 -54.581501 -54.598897 2.981
3 O 15.999400 0 -75.064579 -75.063163 -75.062219 -75.079532 2.981
4 F 18.998403 0 -99.718730 -99.717314 -99.716370 -99.733544 2.981

1.4.2 Standard Enthalphy of Atomization (Hatom)

$$ H_a=\Delta_{at}H^{\theta}= \sum_{atom} H^{\theta}_{atom} - H^{\theta}_{mol} $$

Can be understood as the energy required to break the molecule completely into its component atoms.


In [78]:
#get dict to put into function
enthalphyDict={ row['Element']:row['H (298.15 K)_Hartree'] for indx,row in atomdf.iterrows()}
print(enthalphyDict)
# define function
def enthalphyAtomization( row ):
    enthalphyDict={'H': -0.49791199999999997, 'C': -37.844411, 'F': -99.71637, 'O': -75.062219, 'N': -54.581501}
    return np.sum([ enthalphyDict[atom] for atom in row['atomList'] ])-row['H']

#apply function to each row
df['Hatom']=df.apply(enthalphyAtomization, axis=1)

df.head()


{'H': -0.49791199999999997, 'C': -37.844411, 'F': -99.71637, 'O': -75.062219, 'N': -54.581501}
Out[78]:
numAtoms dbindex A B C dipole polar homo lumo gap ... H G Cv atomList atomCoords freqs SMILES InChI molWeight Hatom
0 5 1 157.71180 157.709970 157.706990 0.0000 13.21 -0.3877 0.1171 0.5048 ... -40.475117 -40.498597 6.469 [C, H, H, H, H] [[-0.0126981359, 1.0858041578, 0.0080009958], ... [1341.307, 1341.3284, 1341.365, 1562.6731, 156... C InChI=1S/CH4/h1H4 16.04246 0.639058
1 4 2 293.60975 293.541110 191.393970 1.6256 9.46 -0.2570 0.0829 0.3399 ... -56.522082 -56.544961 6.316 [N, H, H, H] [[-0.0404260543, 1.0241077531, 0.0625637998], ... [1103.8733, 1684.1158, 1684.3072, 3458.7145, 3... N InChI=1S/H3N/h1H3 17.03052 0.446845
2 3 3 799.58812 437.903860 282.945450 1.8511 6.31 -0.2928 0.0687 0.3615 ... -76.400922 -76.422349 6.002 [O, H, H] [[-0.0343604951, 0.9775395708, 0.0076015923], ... [1671.4222, 3803.6305, 3907.698] O InChI=1S/H2O/h1H2 18.01528 0.342879
3 4 4 0.00000 35.610036 35.610036 0.0000 16.28 -0.2845 0.0506 0.3351 ... -77.304583 -77.327429 8.574 [C, C, H, H] [[0.5995394918, 0.0, 1.0], [-0.5995394918, 0.0... [549.7648, 549.7648, 795.2713, 795.2713, 2078.... C#C InChI=1S/C2H2/c1-2/h1-2H 26.03728 0.619937
4 3 5 0.00000 44.593883 44.593883 2.8937 12.99 -0.3604 0.0191 0.3796 ... -93.408425 -93.431246 6.278 [C, N, H] [[-0.0133239314, 1.1324657151, 0.0082758861], ... [799.0101, 799.0101, 2198.4393, 3490.3686] C#N InChI=1S/CHN/c1-2/h1H 27.02534 0.484601

5 rows × 24 columns


In [103]:
# save new descriptors to file
df.to_pickle('mol.pkl')

In [100]:
# load pickle file
df = pd.read_pickle('mol.pkl')

2) Data Visualization: Exploring the Dataset

2.1) Visualization

For visualization we use imolecule and rdkit which already has some ipython bindings.


In [30]:
nMols=10
indexes=np.random.randint(0,df.shape[0],nMols)
mols=[ Chem.MolFromSmiles(df.ix[i,'SMILES']) for i in indexes]
names=[ "Mol # %d"%(i) for i in indexes]
p = Draw.MolsToGridImage( mols, molsPerRow=5, subImgSize=(200, 200), legends=names)
print("=== Random selection of %d molecules ==="%nMols)
p


=== Random selection of 10 molecules ===
Out[30]:

In [99]:
randomMolecule=np.random.randint(0,df.shape[0])
smileStr=df.loc[randomMolecule,'SMILES']
print("== 2D ==")
mymol=pybel.readstring("smiles", smileStr)
mymol


== 2D ==
Out[99]:
Multiple Molecules - Open Babel Depiction C N O N

In [100]:
print("== 3D Rotatable Molecule==")
imolecule.draw(smileStr, drawing_type="ball and stick", shader="phong")


== 3D Rotatable Molecule==

2.2) Size: Number of atoms, Molecular Weight, Spatial Extensivity

spatialSize = Electronic Spatial Extent

The electronic spatial extent is a single number that attempts to describe the size of a molecule. This number is computed as the expectation value of electron density times the distance from the center of mass of a molecule. Because the information is condensed down to a single number, it does not distinguish between long chains and more globular molecules.


In [30]:
interest=['numAtoms','spatialSize','molWeight']
for i in interest:
    rangeSize=max(df[i])-min(df[i])
    sns_cycle_palette()
    # histogram or a distribution plot?
    if rangeSize < 40:
        bins = np.arange(min(df[i]),max(df[i])+1)
        plt.hist(df[i],bins)
    else:
        sns.distplot(df[i])
    plt.title('Distribution of molecules')
    plt.xlabel(i)
    plt.ylabel('Frequency')
    plt.show()
# show stats
print('=== Stats ===')
df[interest].describe().transpose()


=== Stats ===
Out[30]:
count mean std min 25% 50% 75% max
numAtoms 130830 18.032538 2.943693 3.00000 16.000000 18.00000 20.00000 29.000000
spatialSize 130830 1189.416537 280.471126 19.00020 1017.441875 1147.22465 1309.04710 3374.753200
molWeight 130830 122.721698 7.570695 16.04246 121.179640 125.12528 127.14116 152.038398

Observations:

  • Electronic Spatial Extent follows a normal distribution.
  • Molecular Weight has a bi-modal distribution, it appears to be the mix of 2 normal distributions around 110 and 130.
  • The number of atoms is heavily concentrated between 16 and 20.

2.2 Energetics: Internal energy (U), Enthalpy (H), Gibbs free energy (G), Zero Point Energy (zpe)

Internal energy (U)

Is the sum of all the microscopic energies such as translational kinetic energy, vibrational and rotational kinetic energy and potential energy from intermolecular forces.

Enthalpy $ H = U + pV $

Includes the pressure and Volume of the system. In a sense it accounts for energy transferred to the environment at constant pressure through expansion or heating.

Gibbs free energy $ G = U + pV - TS $

Measures the "usefulness" or process-initiating work obtainable from a thermodynamic system at a constant temperature and pressure.

Zero Point Energy (zpe)

Zero-point energy is the energy that remains when all other energy is removed from a system, i.e the energy of a molecule at T=0.


In [44]:
sns_cycle_palette()
sns.jointplot(df['G'], df['U'], kind="scatter",size=8);
plt.show()
sns_cycle_palette()
sns.jointplot(df['U0'], df['zpe'],kind="reg",size=8);
plt.show()

# show stats
interest=['U','H','G','zpe']
print('=== Stats ===')
df[interest].describe().transpose()


=== Stats ===
Out[44]:
count mean std min 25% 50% 75% max
U 130830 -410.812341 39.891166 -714.560153 -437.870873 -416.800560 -387.031228 -40.476062
H 130830 -410.811397 39.891166 -714.559209 -437.869929 -416.799616 -387.030284 -40.475117
G 130830 -410.854217 39.891884 -714.602138 -437.911833 -416.841337 -387.074526 -40.498597
zpe 130830 0.149090 0.033138 0.015951 0.125638 0.148630 0.171397 0.273944

Observations:

  • U0,U,H,G are all extremely correlated, this makes sense since they defined by vary similar formulas (i.e U=U0+zpe ). We can just keep one measure to reduce data size and number of variables.
  • zpe is a very small number compared to the rest, while quite similar in distribuition to the other energetics it does have some variability.

2.3 Rotational Constants (A,B,C)

Molecules are quantized (due to Quantum Mecchanics) rotating systems, so rotational energy and the angular momentum only takes certain fixed values, these are related to the rotational constants. These values are normally ordered as $A\ge B \ge C$. Based on the relationship of these constants we can classify the molecules as rotors:

  • Spherical rotors where $A = B = C$.
  • Linear molecules where $A >> B=C $ .
  • Oblate rotors where olny $A=B$
  • Prolate where olny $B=C$.
  • Asymmetric rotors where $A \neq B \neq C$.

We can count how many of these we might find in our data to make sense if we should classify molecules this way:


In [49]:
AeqB=np.abs(df['A'] - df['B'])/df['A'] < 0.005
BeqC=np.abs(df['B'] - df['C'])/df['B'] < 0.005
AeqC=np.abs(df['A'] - df['C'])/df['A'] < 0.005

largeRatio= df['A']/df['B'] < 0.25
smallRatio= df['A']/df['B'] > 100
largeA= np.logical_and(df['A'] >= 80 , df['A']/df['B'] > 4)
##clasifiers
spherical=np.logical_and(np.logical_and(AeqB,BeqC),AeqC)

linear=np.logical_or(np.logical_or(largeRatio,smallRatio),largeA)

checked=np.logical_or(linear,spherical)
oblate=np.logical_and(AeqB,np.logical_not(checked)) 

checked=np.logical_or(checked,oblate)
prolate=np.logical_and(BeqC,np.logical_not(checked)) 

checked=np.logical_or(checked,prolate)
asym=np.logical_not(checked)

print("Spherical %d mols"%( sum( spherical)) )
print("Linear %d mols"%( sum(linear) ) )
print("Oblate %d mols"%( sum(oblate) ) )
print("Prolate %d mols"%( sum(prolate) ) )
print("Asymmetric %d mols"%( sum(asym) ) )
# create new data column
df['rotType']='Asym'
for index,row in df.iterrows():
    if spherical[index]:
        df.ix[index,'rotType']='sph'
    elif linear[index]:
        df.ix[index,'rotType']='lin'
    elif oblate[index]:
        df.ix[index,'rotType']='ob'
    elif prolate[index]:
        df.ix[index,'rotType']='pro'


Spherical 7 mols
Linear 25 mols
Oblate 83 mols
Prolate 367 mols
Asymmetric 130348 mols

In [50]:
print("== Spherical Molecules ==")
subdf=df[ spherical ]
mols=[ Chem.MolFromSmiles(row['SMILES']) for indx,row in subdf.iterrows() ] 
names=[ "(# %d) %3.2f, %3.2f, %3.2f"%(row['dbindex'],row['A'],row['B'],row['C'])  for indx,row in subdf.iterrows() ]
p = Draw.MolsToGridImage( mols, molsPerRow=5, legends=names)
p


== Spherical Molecules ==
Out[50]:

In [51]:
print("== Linear ==")
subdf=df[ linear ]
mols=[ Chem.MolFromSmiles(row['SMILES']) for indx,row in subdf.iterrows() ] 
names=[ "(# %d) %3.2f, %3.2f, %3.2f"%(row['dbindex'],row['A'],row['B'],row['C'])  for indx,row in subdf.iterrows() ]
p = Draw.MolsToGridImage( mols, molsPerRow=5, legends=names)
p


== Linear ==
Out[51]:

In [52]:
print("== Some Oblate ==")
subdf=df[ oblate ]
subdf= subdf[subdf['dbindex'].isin(random.sample(subdf['dbindex'], 15))]
#subdf=subdf[:15]
mols=[ Chem.MolFromSmiles(row['SMILES']) for indx,row in subdf.iterrows() ] 
names=[ "(# %d) %3.2f, %3.2f, %3.2f"%(row['dbindex'],row['A'],row['B'],row['C'])  for indx,row in subdf.iterrows() ]
p = Draw.MolsToGridImage( mols, molsPerRow=5, legends=names)
p


== Some Oblate ==
Out[52]:

In [53]:
print("== Some Prolate ==")
subdf=df[ prolate ]
subdf= subdf[subdf['dbindex'].isin(random.sample(subdf['dbindex'], 15))]
mols=[ Chem.MolFromSmiles(row['SMILES']) for indx,row in subdf.iterrows() ] 
names=[ "(# %d) %3.2f, %3.2f, %3.2f"%(row['dbindex'],row['A'],row['B'],row['C'])  for indx,row in subdf.iterrows() ]
p = Draw.MolsToGridImage( mols, molsPerRow=5, legends=names)
p


== Some Prolate ==
Out[53]:

In [54]:
print("== Some Asymmetric ==")
subdf=df[ asym ]
subdf= subdf[subdf['dbindex'].isin(random.sample(subdf['dbindex'], 15))]
mols=[ Chem.MolFromSmiles(row['SMILES']) for indx,row in subdf.iterrows() ] 
names=[ "(# %d) %3.2f, %3.2f, %3.2f"%(row['dbindex'],row['A'],row['B'],row['C'])  for indx,row in subdf.iterrows() ]
p = Draw.MolsToGridImage( mols, molsPerRow=5, legends=names)
p


== Some Asymmetric ==
Out[54]:

In [85]:
# show stats
interest=['A','B','C']

print('=== Stats ===')
df.groupby('rotType')[interest].describe()


=== Stats ===
Out[85]:
A B C
rotType
Asym count 130348.000000 130348.000000 130348.000000
mean 3.428143 1.400986 1.122086
std 2.674430 1.296514 0.849326
min 1.405150 0.337120 0.331180
25% 2.554785 1.091817 0.911447
50% 3.088760 1.370480 1.081650
75% 3.832645 1.654280 1.281630
max 799.588120 437.903860 282.945450
lin count 25.000000 25.000000 25.000000
mean 34181.745950 8.506396 8.283852
std 130569.093627 13.280076 12.812553
min 0.000000 0.377100 0.377100
25% 0.000000 0.782570 0.782570
50% 80.462250 2.048960 2.048960
75% 159.871170 8.593230 8.593210
max 619867.683140 44.593883 44.593883
ob count 83.000000 83.000000 83.000000
mean 6.354241 6.350577 4.146749
std 32.011338 32.004122 20.864110
min 1.382500 1.381480 0.806320
25% 1.825840 1.820495 1.013025
50% 2.100840 2.096210 1.430060
75% 2.660040 2.655590 2.097040
max 293.609750 293.541110 191.393970
pro count 367.000000 367.000000 367.000000
mean 4.787481 1.378116 1.375165
std 2.701848 0.684939 0.684509
min 1.534850 0.373300 0.371660
25% 3.159955 0.907110 0.904440
50% 4.125120 1.292390 1.292330
75% 5.799905 1.728520 1.721730
max 24.204020 5.385450 5.385240
sph count 7.000000 7.000000 7.000000
mean 25.382989 25.381147 25.380613
std 58.370757 58.370755 58.369675
min 1.497270 1.497230 1.497190
25% 2.442800 2.442665 2.442625
50% 3.455190 3.445180 3.445020
75% 5.065530 5.065160 5.064920
max 157.711800 157.709970 157.706990

In [102]:
#exclude large coefficient molecules since they have very different values than skew a lot the stats
subdf=df[np.logical_and(df['A'] < 10,df['B'] < 10) ]
groupsA=[group['A'] for index,group in subdf.groupby('rotType') ]
groupsB=[group['B'] for index,group in subdf.groupby('rotType') ]
groupsC=[group['C'] for index,group in subdf.groupby('rotType') ]
names=[index for index,group in subdf.groupby('rotType')  ]
plt.title("Distribuition for different type of rotors")
sns.violinplot(groupsA,names=names,alpha=0.5,color='green',label='A')
sns.violinplot(groupsB,names=names,alpha=0.5,color='blue')
sns.violinplot(groupsC,names=names,alpha=0.5,color='red',label='A')
#custom legend
import matplotlib.patches as mpatches

green_patch = mpatches.Patch(color='green',alpha=0.5, label='A')
blue_patch = mpatches.Patch(color='blue',alpha=0.5, label='B')
red_patch = mpatches.Patch(color='red',alpha=0.5, label='C')
plt.legend(handles=[green_patch,blue_patch,red_patch])
plt.legend()
plt.show()
print("")




In [38]:
# show stats
interest=['A','B','C']
#exclude linear molecules since they have very different A values than skew a lot the stats
subdf=df[ np.logical_not(linear) ]

for i in interest:
    sns_cycle_palette()
    asubdf=subdf[ subdf[i] < 50 ]
    # histogram, seaborn has better
    sns.distplot(asubdf[i],rug=False, hist=True,hist_kws={"alpha": 0.5},label=i);
    plt.xlim([0,10])

plt.xlabel("Rot Constant value")
plt.title('Distribution of molecules')
plt.ylabel('Frequency')
plt.legend()
plt.show()


Observations:

  • B and C have normal distributions, $B \ge C$ explains the slight skew. A has around 2x bigger mean and std.
  • There are a few outlier cases, such as small molecules which have high rotational numbers ($H_2O$) or linear molecules (Huge $A$). These outliers mess-up the statistics of the distribuitions, removing them, we get ''nice'' distribuitions.
  • Would be interesting to see if each class of rotator also reflects on the distribution of other properties.
  • Might be worth considering only using the subset of asymmetric molecules since they represent the great majority of all the data.
  • How does molecular weight affect the roational constants? Intuition would dictate that more weight = less freedom of movement and this would reflect on smaller rot-constants. Linear molecules are except from this rational.

2.3 Special Properties

Electronic Band Gap (gap)

Calculated as $E_{gap}=E_{HOMO}-E_{LUMO}$, also known as "band gap", where HOMO (highest occupied molecular orbital) and LUMO (lowest unoccupied molecular orbital) are ''frontier'' orbitals for the molecule. Basically, it's how much energy you have to feed into the molecule to kick it from the ground (most stable) state into an excited state. The gap energy can tell us about what wavelengths the compound can absorb.

Molar Heat Capacity (Cv)

At the molecular level, temperature is the average kinetic energy of an ensamble of molecules. Heat capacity characterizes the amount of heat required to change a body's temperature by a given amount.

It is related to molecular complexity, more complex molecules have multiple degrees of freedom (translation, rotation and vibration)...when an ideal gas of ''simple'' molecules (low heat capacity) absorbs heat, all the energy will go to a low number of degrees of freedom...while for a gas with more complex molecules the absorbed energy is partitioned among many more kinds of motions, giving a smaller rise temperature (high heat capacity).

Dipole Moment (dipole)

Dipole moments arise from differences in electronegativity.

It measures the separation of positive and negative electrical charges in a system of electric charges, in this case electrons. The larger the difference in electronegativity, the larger the dipole moment. A zero dipole implies that the negative and positive forces are canceling each other.

Polarizability (polar)

Neutral nonpolar species have spherically symmetric arrangements of electrons in their electron clouds. When in the presence of an electric field, their electron clouds can be distorted. Polarizability is a measure of this distortion.

More electrons or electrons more far away from the nuclear charge, meaning increased polarizability.


In [47]:
interest=['Hatom','Cv','gap','dipole','polar']
# show stats
print('=== Stats ===')
df[interest].describe().transpose()


=== Stats ===
Out[47]:
count mean std min 25% 50% 75% max
Hatom 130830 2.830376 0.385466 0.342879 2.581609 2.835485 3.078835 4.211903
Cv 130830 31.620447 4.067484 6.002000 28.955250 31.578500 34.298000 46.969000
gap 130830 0.252044 0.047192 0.024600 0.217000 0.250200 0.289400 0.622100
dipole 130830 2.672963 1.503480 0.000000 1.577800 2.475300 3.596375 29.556400
polar 130830 75.281410 8.173457 6.310000 70.480000 75.600000 80.610000 196.620000

In [43]:
names=["Enthalphy of atomization","Heat Capacity","Electronic band gap","Dipole Moment","Polarizability"]
for index,i in enumerate(interest):
    sns_cycle_palette()
    # histogram, seaborn has better
    sns.distplot(df[i],rug=False, hist=True);
    plt.xlabel(i)
    plt.ylabel('Frequency')
    plt.title(names[index]+" across all molecules")
    plt.show()