What’s New in Chemistry in the Wolfram Language

Jason Biggs, Wolfram Research

Join the Conversation #WolframTechConf

Chemistry in the Wolfram Language

ImportExport

Import["http://www.rcsb.org/pdb/download/downloadFile.do?fileFormat=pdb&compression=NO&structureId=1tf6","PDB",ImageSize-> Medium]

Import["ExampleData/aspirin.mol"]

Import["ExampleData/fluoxetine.mol"]

Import["ExampleData/caffeine.xyz", "Rendering" -> "Spacefilling"]

ChemicalData

ChemicalData was introduced in version 6, and features curated data for over 40k molecules, drawn from authoratative sources such as NIST, PubChem,

ChemicalData["Propane","CombustionHeat"]

ChemicalData["Tetraphenoxysilane","CombustionHeat"]

ChemicalData["3Methoxy1,2Propanediol","HenryLawConstant"]

ChemicalData["IsononanoylPeroxide","VanDerWaalsConstants"]

ThermodynamicData["CisButene","MolarIsobaricHeatCapacity",{"Temperature"->Quantity[85,"DegreesCelsius"]}]

EntityValue

Text[Grid[Outer[ChemicalData,{"Glycine","LTyrosine","LLeucine","LSerine"},{"Name","ColorStructureDiagram","SpaceFillingMoleculePlot"}]]]

{(Entity["glycine", Entity["Chemical","Glycine"], "Entity[\"Chemical\", \"Glycine\"]", "chemical"]), (Entity["L-tyrosine", Entity["Chemical","LTyrosine"], "Entity[\"Chemical\", \"LTyrosine\"]", "chemical"]), (Entity["L-leucine", Entity["Chemical","LLeucine"], "Entity[\"Chemical\", \"LLeucine\"]", "chemical"]), (Entity["L-serine", Entity["Chemical","LSerine"], "Entity[\"Chemical\", \"LSerine\"]", "chemical"])}

EntityValue[{(Entity["glycine", Entity["Chemical","Glycine"], "Entity[\"Chemical\", \"Glycine\"]", "chemical"]), (Entity["L-tyrosine", Entity["Chemical","LTyrosine"], "Entity[\"Chemical\", \"LTyrosine\"]", "chemical"]), (Entity["L-leucine", Entity["Chemical","LLeucine"], "Entity[\"Chemical\", \"LLeucine\"]", "chemical"]), (Entity["L-serine", Entity["Chemical","LSerine"], "Entity[\"Chemical\", \"LSerine\"]", "chemical"])},{(EntityProperty["structure diagram", EntityProperty["Chemical","ColorStructureDiagram"], "EntityProperty[\"Chemical\", \"ColorStructureDiagram\"]"]),(EntityProperty["melting point", EntityProperty["Chemical","MeltingPoint"], "EntityProperty[\"Chemical\", \"MeltingPoint\"]"])},"Dataset"]

Wolfram DataRepository

GDB-9 Database

The Data Repository has a large EntityStores of

PrependTo[$EntityStores,ResourceData["GDB-9 Database"]]

RandomEntity["GDB9Chemical"]

EntityProperties["GDB9Chemical"]

Dataset[<|#->#["Description"]&/@{(EntityProperty["G", EntityProperty["GDB9Chemical","FreeEnergy25C"], "EntityProperty[\"GDB9Chemical\", \"FreeEnergy25C\"]"]),(EntityProperty["Cv", EntityProperty["GDB9Chemical","HeatCapacity25C"], "EntityProperty[\"GDB9Chemical\", \"HeatCapacity25C\"]"]),(EntityProperty["HOMO", EntityProperty["GDB9Chemical","HOMO"], "EntityProperty[\"GDB9Chemical\", \"HOMO\"]"]),(EntityProperty["HOMO-LUMO gap", EntityProperty["GDB9Chemical","HOMOLUMOGap"], "EntityProperty[\"GDB9Chemical\", \"HOMOLUMOGap\"]"]),(EntityProperty["InChI", EntityProperty["GDB9Chemical","InChI"], "EntityProperty[\"GDB9Chemical\", \"InChI\"]"])}|>]

RandomEntity["GDB9Chemical"]
%["Dataset"]

Pitt Quantum Repository

PrependTo[$EntityStores,ResourceData["Pitt Quantum Repository"]]

structureGraph=Graph[UndirectedEdge@@@#["EdgeRules"]]&;
nRings=Length@FindCycle[structureGraph@#,8,All]&;
ringsData=#->nRings[#]&/@EntityList["PQRChemical"]//Apply[Association];
Histogram[ringsData]

EntityValue[Keys@Sort[ringsData][[-5;;]],"MoleculePlot"]

Virtual Exploratory heterocyclic library

The VEHICLe database, short for Virtual Exploratory heterocyclic library, contains all possible aromatic molecules meeting the following constraints: 1) neutral, 2) composed only of carbon, nitrogen, oxygen, sulfur, and hydrogen, 3) all atoms in rings except hydrogen and carbonyl oxygens 4) either 1 or 2 aromatic rings. The majority of these compounds have not yet been synthesized in the laboratory.

PrependTo[$EntityStores,ResourceData["VEHICLe"]];

RandomEntity["HeterocyclicChemical",5]

EntityValue[%,"ColorStructureDiagram"]

WolframAlpha

Wolfram Alpha draws on ChemicalData, but also has more resources

Alpha can be queried using the natural-language input operator

WolframAlpha’s Step-by-Step solutions give not only the answer, but help you understand where the answer comes from

(*cached version*)

(* cached version *)

Chemistry in the Wolfram Language: User demand for more interactivity

The entities in ChemicalData enable a wide variety of computation, but is limited to the curated data presently in the Knowledgebase

The functions used to generate the visual representations are not exposed or documented

Stereochemistry information is absent from ChemicalData

No functionality exists for creating new chemical entities, or modifying existing entities

Attaching a legend to a MoleculePlot

How can I access the internal function that plots a molecule from a formatted XYZ file?

Changing the atom color for plotting XYZ-files

Plot 3D model of DNA in Mathematica

Retrieving coordinate positions using ChemicalData for input to quantum modelling software

Faking MoleculePlot for compounds where it isn’t available

Plotting electronic orbitals with Wolfram Language

The MoleculeViewer package

A Wolfram Language molecule data structure

Who would be interested in a computable chemical object?

Who is the target audience?

Chemistry students and educators

Constructing and manipulating molecules in 3D gives “hands on” experience

querying for structural properties as a function of 3D structure

Academic researchers and grad students

Computational chemists

Cheminformatics

Substructure screening

Database querying for similar molecules

Quantitative structure-activity relationship (QSAR)

Molecular modeling

Quantum chemistry

Conformational analysis

Molecular dynamics

Me:

Bob:

What form would a computable chemical object take?

Similar to a Graph

The vertices are atoms, the edges are bonds

Graph[UndirectedEdge@@@ChemicalData["Benzophenone","EdgeRules"]]

By adding the coordinates and labels from ChemicalData, we can get a decent representation

ChemicalData["Benzophenone", "VertexTypes"]

graph=Graph[Range @ Length @ ChemicalData["Benzophenone", "VertexTypes"],
	UndirectedEdge@@@ChemicalData["Benzophenone", "EdgeRules"],
	VertexLabels -> Thread[Range @ Length @ ChemicalData["Benzophenone", "VertexTypes"]->ChemicalData["Benzophenone", "VertexTypes"]],
	VertexCoordinates -> ChemicalData["Benzophenone", "VertexCoordinates"]
]

By calling on a specific ChemicalData property, we can get a more visually appealing Graph, with bond orders indicated

ChemicalData["Benzophenone","StructureGraph"]

The Graph is a computable object in the Wolfram Language, and so we could ask for the longest minimal path between two atoms

With[
	{
		path = Last @ SortBy[Length][
 				Union @ Flatten[
   					Table[FindShortestPath[graph, a, b],
    						{a, Length @ VertexList @ graph},
    						{b, a, Length @ VertexList @ graph}
    					],
   					1
   				]
 			]
	},
	HighlightGraph[graph, PathGraph @ path]
]

But with many distinct differences

The vertices are not distinct; in a Graph, duplicate vertices are pruned

Vertex and edge properties, such as charge, isotope, bond-order, are very important

hyperedges - connections between multiple vertices.

stereo elements - collections of 4 or 5 vertices with a specific 3D spatial relationship

Coordinates are not just for visualization, but directly influence molecular properties

The Molecule

What’s in a name?

Technically a molecule must be a connected structure with no net charge.  This is too restrictive for our purposes.

MolecularEntity is more precise, but would be far too confusing as a WL Head.  Another precise name would be an AtomContainer, as used by the CDK toolkit

For now, let’s use Molecule.  It’s a simple, like Graph or Region, and lends itself to compound function names like MoleculePlot, MoleculeNearest, MoleculeReplace, etc

InputForm for a Molecule

?Molecule

Manual entry

Create a simple molecule

Molecule[{"C","C","C"}, {Bond[{1,2}],Bond[{2,3}]}]

Then you can query the Knowledgebase to find a matching ChemicalData Entity

ToChemicalEntity@%

Adding atom properties is done the same way as with Graph

m = Molecule[{"C",Property["C",{"FormalCharge"->1}],"C"}, {Bond[{1,2}],Bond[{2,3}]}]

MoleculePlot@%

After input, the properties are moved to the options.

InputForm@m

From an Entity

Create a Molecule from an Entity

Molecule@(Entity["propane-1,1,1,2,3,3,3-d 7", Entity["Chemical","Propane(1,1,1,2,3,3,3)D7"], "Entity[\"Chemical\", \"Propane(1,1,1,2,3,3,3)D7\"]", "chemical"])

PrettyForm@%

From a chemical name

A molecule can be created from a valid chemical name, using the OPSIN parser

namesList = PersistentValue["ChemicalNames"];

Molecule@Echo@RandomChoice@namesList

PrettyForm@%

ToChemicalEntity@Molecule[]

From a SMILES string

smilesList = PersistentValue["PQRSmiles"];

smiles = RandomChoice @ smilesList

m = Molecule @"CCCCOc1nc(N)nc(n1)N"

m = Molecule @"CC(CCc1ccccc1)COc1nc(N)nc(n1)[NH3+]"

PrettyForm@%

Molecule visualization

ChemicalData

ChemicalData has had the properties “ColorStructureDiagram” and “MoleculePlot”

EntityValue[(Entity["adenosine triphosphate", Entity["Chemical","AdenosineTriphosphate"], "Entity[\"Chemical\", \"AdenosineTriphosphate\"]", "chemical"]),{"MoleculePlot","ColorStructureDiagram"}]

The lack of double bonds in the 3D plot has long been a source of complaints.  For this, the “MultiBondMoleculePlot” was created

EntityValue[(Entity["adenosine triphosphate", Entity["Chemical","AdenosineTriphosphate"], "Entity[\"Chemical\", \"AdenosineTriphosphate\"]", "chemical"]),{"MultiBondMoleculePlot","ColorStructureDiagram"}]

MoleculePlot - 2D depictions

{m1, m2} = Molecule[#,"RemoveHydrogens"->True]&/@{"adenosine","adenosine triphosphate"}

MoleculePlot@m1

ContainsStructureQ[m2,m1]

MoleculePlot[m2,"Highlights"->{m1}]

MoleculePlot3D

fullerene =Molecule@ (Entity["buckminsterfullerene", Entity["Chemical","FullereneC60"], "Entity[\"Chemical\", \"FullereneC60\"]", "chemical"])

Using SMARTS strings for substructure detection, highlight all 5-membered rings in fullerene

MoleculePlot3D[
fullerene,
"Highlights" ->{"C1@C@C@C@C1"},
"Rendering"->"MultiBond"]

Stereochemistry in 3D

m = Molecule@"L-alanine"

MoleculePlot[m,"ShowHydrogens"->True]

PrettyForm@m

MoleculePlot[m,"LabelIndices"->True]

MoleculePlot3D[m,"Highlights" ->{{1},{3},{4},{9}},"Rendering"->"MultiBond"]

m2=(DefinitionSymbol[Chemistry`Molecule, Molecule])[
	{"N", "C", "C", "C", "O", "O", "H", "H", "H", "H", "H", "H", "H"},
	{
		(DefinitionSymbol[Chemistry`Bond, Bond]) @ {1, 2},
		(DefinitionSymbol[Chemistry`Bond, Bond]) @ {2, 3},
		(DefinitionSymbol[Chemistry`Bond, Bond]) @ {2, 4},
		(DefinitionSymbol[Chemistry`Bond, Bond])[{4, 5}, "Double"],
		(DefinitionSymbol[Chemistry`Bond, Bond]) @ {4, 6},
		(DefinitionSymbol[Chemistry`Bond, Bond]) @ {1, 7},
		(DefinitionSymbol[Chemistry`Bond, Bond]) @ {1, 8},
		(DefinitionSymbol[Chemistry`Bond, Bond]) @ {2, 9},
		(DefinitionSymbol[Chemistry`Bond, Bond]) @ {3, 10},
		(DefinitionSymbol[Chemistry`Bond, Bond]) @ {3, 11},
		(DefinitionSymbol[Chemistry`Bond, Bond]) @ {3, 12},
		(DefinitionSymbol[Chemistry`Bond, Bond]) @ {6, 13},
		(DefinitionSymbol[Chemistry`StereoElement, StereoElement])[2, "Counterclockwise", {1, 3, 4, 9}]
	}
]

MoleculePlot3D[#,"Highlights" ->{{1},{3},{4},{9}},"Rendering"->"MultiBond"]&/@{m,m2}

We can modify the s

Structure-based property prediction

Joback Method - predicting properties based on substructure matching

A group contribution method developed in the 1980’s, attempts to find thermodynamic properties by counting the presence of certain substructures.

m = Molecule@(Entity["hydrastine", Entity["Chemical","Hydrastine"], "Entity[\"Chemical\", \"Hydrastine\"]", "chemical"])

MoleculePlot[m, "Highlights"->{"JobackContributions"}]

QSAR properties

These vectors of numeric properties can be used to predict molecular properties and biological activity.  Applications to machine learning in the WL will be forthcoming.

m = Molecule["tetrahydroimidazo[4,5,1-jk][1,4]benzodiazepinone"]

m = Molecule["tetrahydroimidazo[4,5,1-jk][1,4]benzodiazepinone"]

m["QSARDataset"]

m2=Molecule["cytidine monophosphate"]

m2["QSARDataset"]

Available 2D descriptors:

Available 3D descriptors:

Fingerprint similarity searching

Available fingerprint methods

SeedRandom[42];
smiles = RandomSample[smilesList,10000];
mols = Table[
Molecule[smiles[[n]],"RemoveHydrogens"->True],
{n,10000}];~Monitor~n

mnf = MoleculeNearest[mols]

m = Molecule[Echo @ RandomChoice @ Complement[smilesList, smiles],"RemoveHydrogens"->True]

Find the 5 molecules most similar to a query

mnf[m,5]

Grid[Partition[
MoleculePlot/@Prepend[%,m],2]]

3D structure manipulation

m = Molecule@"Oc1ccc(cc1)C1=CC(=O)c2c([OH+]1)cc(c(c2O)[C@@H]1O[C@H](C)[C@H](C(=O)[C@H]1O[C@@H]1O[C@H](C)[C@@H]([C@H]([C@@H]1O)O)O)O)O"

MoleculeManipulate@m

Constrained geometry optimization

We can use the Merck Molecular Force Field to find local energy minima for generated 3D structures.  Here we can generate a potential energy surface for the torsional coordinate of ethane

m = Molecule["Ethane"]

MoleculePlot[m,"LabelIndices"->True]

{#,MoleculeEnergyMinimize[
m,
"GeometryConstraints" -> <|
"TorsionConstraints" -> {{{6,2,1,4},#}}
|>
];m["ForceFieldEnergy"]} &/@Range[-180.,180.,5.]

ListLinePlot@%

<<Chemistry`

Future directions

Similarity searching through molecules in the Knowledgebase

Substructure screening through molecules in the Knowledgebase

Enhanced chirality detection - cumulene stereochemistry, chiral axes

Vibrational frequency determination

Molecular orbital calculations

? - Join us for the Chemistry roundtable discussion at lunch on Friday