The Molecule Function

Jason D Biggs, Wolfram Research

Join the Conversation #WolframTechConf

ABSTRACT

Since version 6, ChemicalData has provided access to thousands of curated chemical entities and their properties.  With the Molecule function, the number of molecules you can create for computations is uncountable.  Create molecules using chemical identifiers like IUPAC names or SMILES strings, or by explicit lists of atoms and bonds.  Generate structure diagrams and three-dimensional conformations for the molecules at runtime.  Search for and highlight chemical substructures using MoleculePattern.  Modify your molecules by adding and deleting atoms or bonds, or changing bond lengths and angles, or perform a geometry optimization using the MMFF force field.

Chemistry in the Wolfram Language so far

The Wolfram Language has been working on chemistry content since version 6

ChemicalData and EntityValue are portals into the Wolfram Knowledgebase, and contain curated information on over 40,000 chemical entities

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: demand for more

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

ChemTools

The MoleculeViewer package

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 Molecule

Does not rely on curated data, constructed on the fly from the input expression

Can generate 2D and 3D coordinates, compute properties, etc, based only on the structure.

The Molecule acts as both a constructor and an expression.

Molecule[...] can take many inputs and evaluates if possible to the form (RefLink[Molecule, paclet:ref/Molecule])[ atoms,bonds,…].

?Molecule

Molecule["benzene"]//InputForm

Create a Molecule from atoms and bonds

atoms = {"O","C","C","C","N","C","C","C","C",Atom["N","FormalCharge"->1],Atom["O","FormalCharge"->-1],"O","H","H","H","H","H","H","H","H"};
bonds = {Bond[{1,2},"Single"],Bond[{2,3},"Single"],Bond[{3,4},"Aromatic"],Bond[{4,5},"Single"],Bond[{4,6},"Aromatic"],Bond[{6,7},"Aromatic"],Bond[{7,8},"Aromatic"],Bond[{8,9},"Aromatic"],Bond[{9,10},"Single"],Bond[{10,11},"Single"],Bond[{10,12},"Double"],Bond[{9,3},"Aromatic"],Bond[{1,13},"Single"],Bond[{2,14},"Single"],Bond[{2,15},"Single"],Bond[{5,16},"Single"],Bond[{5,17},"Single"],Bond[{6,18},"Single"],Bond[{7,19},"Single"],Bond[{8,20},"Single"]};

mol1 =Molecule[atoms,bonds]

Hydrogen atoms may be omitted from the atom list if their presence can be inferred from the valence and bonding of the atoms present.

atoms = {"O","C","C","C","N","C","C","C","C",Atom["N","FormalCharge"->1],Atom["O","FormalCharge"->-1],"O"};
bonds = {Bond[{1,2},"Single"],Bond[{2,3},"Single"],Bond[{3,4},"Aromatic"],Bond[{4,5},"Single"],Bond[{4,6},"Aromatic"],Bond[{6,7},"Aromatic"],Bond[{7,8},"Aromatic"],Bond[{8,9},"Aromatic"],Bond[{9,10},"Single"],Bond[{10,11},"Single"],Bond[{10,12},"Double"],Bond[{9,3},"Aromatic"]};

mol2 = Molecule[atoms,bonds]

The molecules have different expressions

mol1===mol2

ByteCount[mol1]

ByteCount[mol2]

But chemically they are equivalent

MoleculeEquivalentQ[mol1,mol2]

True

Create a Molecule using an Entity

When given an Entity, Molecule will query the knowledgebase for the relevant properties

EntityValue[(Entity["propane", Entity["Chemical","Propane"], "Entity[\"Chemical\", \"Propane\"]", "chemical"]),{"VertexTypes","EdgeRules","EdgeTypes","FormalCharges","NonStandardIsotopeNumbers","VertexCoordinates","AtomPositions"},"PropertyAssociation"]

Molecule[LinguisticAssistant]

These two molecules are the same, but one is querying the knowledgebase for information, the other is creating the Molecule programmatically from the name

MoleculeEquivalentQ[Molecule[LinguisticAssistant], Molecule["propane"]]

Some chemical entities don’t have a single molecular representation - polymers for example.  These can’t make molecules,

Molecule@(Entity["poly(butyl acrylate)", Entity["Chemical","PolyButylAcrylate"], "Entity[\"Chemical\", \"PolyButylAcrylate\"]", "chemical"])

EntityValue[(Entity["poly(butyl acrylate)", Entity["Chemical","PolyButylAcrylate"], "Entity[\"Chemical\", \"PolyButylAcrylate\"]", "chemical"]),{"VertexTypes","EdgeRules","EdgeTypes","FormalCharges","NonStandardIsotopeNumbers","VertexCoordinates","AtomPositions"},"PropertyAssociation"]

Create a Molecule from a Chemical Name

Version 12 of the Wolfram Language contains the OPSIN parser for chemical names, which can convert a systematic chemical name into a chemical graph - from there we can create a Molecule

SeedRandom[42]
RandomChoice@names

Here is the internal function we can use to parse this to a SMILES string

nts=Chemistry`OPSINLink`NameToSMILES

Convert the chemical name to a SMILES string

nts["2-methylpent-1-en-3-yne"]

You can make a Molecule from this SMILES string, or from the name directly

Molecule@nts["2-methylpent-1-en-3-yne"]
Molecule["2-methylpent-1-en-3-yne"]

SameQ[%,%%]

Further neat examples

Molecule@"(2S)-N-[(1S)-2-carbamoyl-1-[[(2S)-1-[(2S)-2-[[(1S)-1-[[(1S)-1-carbamoyl-2-(4-methoxyphenyl)ethyl]carbamoyl]-4-(diaminomethylideneamino)butyl]carbamoyl]pyrrolidin-1-yl]-1-oxo-3-sulfanylidene-propan-2-yl]carbamoyl]ethyl]-2-[[(2S)-2-[[(2S)-3-(4-hydroxyphenyl)-2-[[2-(1-sulfanylcyclohexyl)acetyl]amino]propanoyl]amino]-3-phenyl-propanoyl]amino]pentanediamide"

mol = Molecule@"[(2R,3S,4R,6R)-6-[(1R,2'R,3R,3'S,4'R,5R,6R)-3-[(2S,3R,4R,5S,6R)-2-[(2R,3S,4R,5S,6S)-6-[(1R,1'R,2S,2'R,4R,5R,6R,6'S)-5-acetyl-5-hydroxy-4-methyl-2'-(2-methylpropanoyloxy)spiro[3,7,9-trioxabicyclo[4.3.0]nonane-2,8'-4,7,9-trioxabicyclo[4.3.0]nonane]-3'-yl]oxy-4-hydroxy-5-methoxy-2-(methoxymethyl)oxan-3-yl]oxy-3-hydroxy-5-methoxy-6-methyl-oxan-4-yl]oxy-4'-hydroxy-1,2',5-trimethyl-spiro[4,7,9-trioxabicyclo[4.3.0]nonane-8,6'-oxane]-3'-yl]oxy-4-hydroxy-2-methyl-oxan-3-yl] 3,5-dichloro-4-hydroxy-2-methoxy-6-methyl-benzoate"

MoleculePlot@mol

MoleculePlot3D@mol

Create a Molecule from a SMILES string

({{({{▪ SMILES chemical format.}, {▪ Commonly used to describe the structure of chemical molecules. }, {▪ SMILES is an acronym for Simplified Molecular Input Line Entry Specification.}, {▪ Used in cheminformatics applications and in chemistry databases to represent chemical formulas.
▪ Developed in the 1980s by Arthur Weininger and David Weininger.}}), ({{}, {}})}})

Import has been able to read a SMILES string for some time, to create a Graph

ImportString["CN1C=NC2=C1C(=O)N(C(=O)N2C)C", {"SMILES","EdgeRules",1}]//Graph[#,GraphLayout->"SpringEmbedding"]&

Now we can create a Molecule from the same input:

Molecule["CN1C=NC2=C1C(=O)N(C(=O)N2)C"]

Molecule["C(CCCCCCCCCCCCCCCCCC3)N1C=NC2=C1C(=O)N(C(=O)N2)C3"]

MoleculePlot3D@%

These SMILES strings can be quite fun to mix and match, rearrange to create novel molecules which may or may not be synthesizable

Take a SMILES string like “O[C@@H](C[C@H](CC(=O)O)O)CCn1c(c2ccc(cc2)F)c(c(c1C(C)C)S(=O)(=O)N1CCOCC1)c1ccccc1” and start adding on to it

Molecule["O[C@@H](C[C@H](CC(=O)O)O)CCn1c(c2ccc(cc2)F)c(c(c1C(C)C)S(=O)(=O)N1CCOCC1)c1ccccc1"]

Now take another random SMILES and add it anywhere

"O=C[C@@H](Cc1ccc2c(c1)OCO2)C"

Molecule["O[C@@H](C[C@H](CC(=O)O)O)CCn1c(c2ccc(cc2)F)c(c(c1C(C)C)S(=O)(OOC[C@@H](Cc1ccc2c(c1)OCO2)C)N1CCOCC1)c1ccccc1"]

MoleculePlot@%

MoleculePlot3D@%%

StereochemistryElements

?StereochemistryElements

The option value describes the local geometry around an atom or bond, and does not specify the absolute stereochemistry(e.g. stereocenters "R" or "S").

The value of (RefLink[StereochemistryElements, paclet:ref/StereochemistryElements]) is used for generating two-dimensional structure diagrams and for embedding the molecule in three dimensions.

The value of the (RefLink[StereochemistryElements, paclet:ref/StereochemistryElements]) option can be obtained from a molecule via (RefLink[MoleculeValue, paclet:ref/MoleculeValue]) using the property name "StereochemistryElements".

Elements are entered as a list of associations, with a "StereoType" key.

Chiral centers

ccw=Molecule[{Atom["C"],Atom["C"],Atom["C"],Atom["C"],Atom["O"],Atom["H"],Atom["H"],Atom["H"],Atom["H"],Atom["H"],Atom["H"],Atom["H"],Atom["H"],Atom["H"],Atom["H"]},{Bond[{1,2},"Single"],Bond[{2,3},"Single"],Bond[{3,4},"Single"],Bond[{2,5},"Single"],Bond[{1,6},"Single"],Bond[{1,7},"Single"],Bond[{1,8},"Single"],Bond[{2,9},"Single"],Bond[{3,10},"Single"],Bond[{3,11},"Single"],Bond[{4,12},"Single"],Bond[{4,13},"Single"],Bond[{4,14},"Single"],Bond[{5,15},"Single"]},StereochemistryElements->{\[LeftAssociation]"StereoType"->"Tetrahedral","ChiralCenter"->2,"Direction"->"Counterclockwise","FiducialAtom"->1,"Ligands"->{3,5,9}\[RightAssociation]}]

cw = Molecule[{Atom["C"],Atom["C"],Atom["C"],Atom["C"],Atom["O"],Atom["H"],Atom["H"],Atom["H"],Atom["H"],Atom["H"],Atom["H"],Atom["H"],Atom["H"],Atom["H"],Atom["H"]},{Bond[{1,2},"Single"],Bond[{2,3},"Single"],Bond[{3,4},"Single"],Bond[{2,5},"Single"],Bond[{1,6},"Single"],Bond[{1,7},"Single"],Bond[{1,8},"Single"],Bond[{2,9},"Single"],Bond[{3,10},"Single"],Bond[{3,11},"Single"],Bond[{4,12},"Single"],Bond[{4,13},"Single"],Bond[{4,14},"Single"],Bond[{5,15},"Single"]},StereochemistryElements->{\[LeftAssociation]"StereoType"->"Tetrahedral","ChiralCenter"->2,"Direction"->"Clockwise","FiducialAtom"->1,"Ligands"->{3,5,9}\[RightAssociation]}]

The absolute chirality is computed based on the local assignments

cw["CIPChirality",2]

ccw["CIPChirality",2]

MoleculeEquivalentQ[cw,ccw]

MoleculeEquivalentQ[cw,ccw,SameTest->"Enantiomers"]

If the absolute chirality is entered, the

cw2 = Molecule[{Atom["C"],Atom["C"],Atom["C"],Atom["C"],Atom["O"],Atom["H"],Atom["H"],Atom["H"],Atom["H"],Atom["H"],Atom["H"],Atom["H"],Atom["H"],Atom["H"],Atom["H"]},{Bond[{1,2},"Single"],Bond[{2,3},"Single"],Bond[{3,4},"Single"],Bond[{2,5},"Single"],Bond[{1,6},"Single"],Bond[{1,7},"Single"],Bond[{1,8},"Single"],Bond[{2,9},"Single"],Bond[{3,10},"Single"],Bond[{3,11},"Single"],Bond[{4,12},"Single"],Bond[{4,13},"Single"],Bond[{4,14},"Single"],Bond[{5,15},"Single"]},StereochemistryElements->{\[LeftAssociation]"StereoType"->"Tetrahedral","ChiralCenter"->2,"Direction"->"S"\[RightAssociation]}]

Options@%

MoleculePlot3D[#,Thread[{{3,5,9}}]]&/@{cw,ccw}

Double bond stereochemistry

Create a molecule with defined geometry around a double bond:

m1=Molecule[{Atom["Cl"],Atom["C"],Atom["C"],Atom["C"],Atom["Cl"],Atom["C"],Atom["Cl"]},{Bond[{1,2},"Single"],Bond[{2,3},"Single"],Bond[{3,4},"Double"],Bond[{4,5},"Single"],Bond[{4,6},"Single"],Bond[{3,7},"Single"]},StereochemistryElements->{\[LeftAssociation]"StereoType"->"DoubleBond","StereoBond"->{3,4},"Ligands"->{7,5},"Value"->"Opposite"\[RightAssociation]}]

Next create its geometric isomer:

m2=Molecule[{Atom["Cl"],Atom["C"],Atom["C"],Atom["C"],Atom["Cl"],Atom["C"],Atom["Cl"]},{Bond[{1,2},"Single"],Bond[{2,3},"Single"],Bond[{3,4},"Double"],Bond[{4,5},"Single"],Bond[{4,6},"Single"],Bond[{3,7},"Single"]},StereochemistryElements->{\[LeftAssociation]"StereoType"->"DoubleBond","StereoBond"->{3,4},"Ligands"->{7,5},"Value"->"Together"\[RightAssociation]}]

The geometry is reflected in the SMILES and InChI identifiers:

MoleculeValue[{m1,m2},"SMILES"]

MoleculeValue[{m1,m2},"InChI"]

MoleculePlot3D/@{m1,m2}

MoleculeEquivalentQ[m1,m2]

MoleculeEquivalentQ[m1,m2,SameTest->"Stereoisomers"]

Molecule properties

Computed molecule properties can be accessed by sub-values

m = Molecule@"Clc1ccc(cc1)C(c1ccc(cc1)Cl)c1ccc(cc1)Cl"

m["SMILES"]

m["InChI"]

Partition[m["Properties"],UpTo[4]]//Grid

MoleculeValue is an accessor function for properties

MoleculeValue[m,"SMILES"]

MoleculeValue[m,{"LabuteApproximateSurfaceArea","RingCount","SaturatedRingCount","AromaticRingCount"}]

?MoleculeValue

Use MoleculeValue to get properties for atoms:

MoleculeValue[{m,5},"AtomicNumber"]

Compute the Gasteiger partial charges

MoleculeValue[{m,{5,6,7}},"GasteigerPartialCharge"]

AtomList, AtomCount

?AtomList

mol=Molecule["benzene"]

m=Molecule["sodium benzoate"]

AtomList[m,Atom["FormalCharge"->Except[0]]]
AtomCount[m,Atom["FormalCharge"->Except[0]]]

AtomList[m,Atom["FormalCharge"->GreaterThan[0]]]
AtomCount[m,Atom["FormalCharge"->GreaterThan[0]]]
AtomCount[m,Atom["FormalCharge"->GreaterEqualThan[0]]]

AtomList[m,Atom["AtomicNumber"->GreaterThan[6]]]
AtomCount[m,Atom["AtomicNumber"->GreaterThan[6]]]
AtomCount[m,Atom["AtomicNumber"->GreaterEqualThan[6]]]

BondList, BondCount

mol = Molecule[];

MoleculePlot@mol

BondList[mol,"Double"]

BondList[mol,Bond[_,"Double"]]

Get all bonds where only one atom is in a ring

bl=BondList[mol,Bond[{Atom["RingMemberQ"->True],Atom["RingMemberQ"->False]}],IncludeHydrogens->False]

MoleculePlot[mol,Style[bl,]]

Get all bonds where both atoms are in a ring

bl=BondList[mol,Bond[{Atom["RingMemberQ"->True],Atom["RingMemberQ"->True]}],IncludeHydrogens->False]

MoleculePlot[mol,Style[bl,]]

MoleculePlot[mol,Bond[{Atom["RingMemberQ"->True],Atom["RingMemberQ"->True]}]]

mol=Molecule[]

MoleculePlot@mol

BondList[mol,Bond[{Atom["FormalCharge"->GreaterThan[0]],_}]]

BondList[mol,Bond[{Atom["FormalCharge"->GreaterThan[0]],Atom["FormalCharge"->LessThan[0]]}]]

BondQ

?BondQ

The bond may be a (RefLink[Bond, paclet:ref/Bond]) or a list of two atom indices.

(RefLink[BondQ, paclet:ref/BondQ]) checks for the existence of a supported bond type in the molecule expression.

An unspecified bond order will match any bond type between the given atoms.

The list of atoms in bond is treated as orderless for the purposes of matching.

m=Molecule[];
MoleculePlot[m]

Atoms 1 and 2 are bonded, 1 and 3 are not:

BondQ[m,{1,2}]
BondQ[m,{2,1}]

BondQ[m,{1,3}]

Atoms 2 and 4 have a double bond:

BondQ[m,Bond[{4,2},"Double"]]
BondQ[m,Bond[{4,2},"Single"]]

MoleculePattern, FindMoleculeSubstructure

?FindMoleculeSubstructure

m=Molecule["CSc1ccc2c(c1)N(CC[C@H]1CCCCN1C)c1c([S@]2=O)cccc1"]

Define a pattern to match sulfur bonds

A single SO bond

This pattern can be written in two equivalent ways

patt1=MoleculePattern[{"S","O"},{Bond[{1,2}]}]
patt1b=Bond[{"S","O"}]

Find the pattern inside a molecule using FindMoleculeSubstructure

m=Molecule[];
FindMoleculeSubstructure[m,patt1,All]
FindMoleculeSubstructure[m,patt1b,All]

If we use the SMARTS language, we can simply write this pattern as “S~O” where ~ signifies ‘any’ kind of bond

FindMoleculeSubstructure[m,MoleculePattern["S~O"],All]

A three-atom substructure

MoleculePattern is needed to find substructures containing more than 2 atoms:

patt2=MoleculePattern[{"S","C","C"},{Bond[{1,2},"Single"],Bond[{1,3},"Single"]}]

FindMoleculeSubstructure[m,patt2,All]

MoleculePlot[m,Keys/@%]

Define a pattern for an N-membered ring:

(ringPattern[n_] := MoleculePattern[ConstantArray[_, n],
	Map[Bond, Partition[Range @ n, 2, 1, 1]]
];
ringPattern[list_List] := MoleculePattern[
	list, 
	Bond /@ Partition[Range @ Length @ list, 2, 1, 1]
];)

ringPattern[6]

MoleculePlot[m,ringPattern[6]]

Highlight all six-membered rings with at least one nitrogen

MoleculePlot[m,ringPattern[{"N",_,_,_,_,_}]]

m4=Molecule[];

MoleculePlot[m4,{Style[ringPattern[4],RGBColor[0.91,0.5,0.26]],ringPattern[8]}]

MoleculePlot3D[m4,{Style[ringPattern[4],RGBColor[0.91,0.5,0.26]],ringPattern[8]}]

Define a pattern for a polarized bond:

bondPattern = Bond[{Atom["GasteigerPartialCharge"->LessThan[-0.1]],Atom["GasteigerPartialCharge"->GreaterThan[0.1]]}]

m=Molecule@"CN(S(=O)(=O)c1ccccc1)C1CCN(CC1)C"

There are three bonds that match this pattern:

FindMoleculeSubstructure[m,bondPattern,All]

MoleculePlot[m,bondPattern]

m["GasteigerPartialCharge",{2,3,4,5}]

MoleculeModify

Manual modification

The Molecule is not AtomQ, so it can be modified in place

m = Molecule@"CCO[P@@](=O)(Oc1ccc(cc1)N(=O)=O)C"

ReplaceAll can be used to switch bond orders

m2 = m /."Aromatic"->"Single"

Part assignment works

m2[[1,15]]=LinguisticAssistant

Evaluation converts this to an Atom

m2[[1,15]]

Hydrogens are added to achieve the right valence

m2

Part assignment can make a valid molecule into an invalid one

m = Molecule@"cyclohexane"

ReplaceAll[
m, "Single"->"Double"]

AddAtom

m = Molecule@"caffeine"

MoleculeModify[m,"AddAtom",{"C","O",(Entity["carbon", Entity["Element","Carbon"], "Entity[\"Element\", \"Carbon\"]", "element"]),(Entity["carbon-14", Entity["Isotope","Carbon14"], "Entity[\"Isotope\", \"Carbon14\"]", "isotope"])}]

DeleteAtom

MoleculeModify[m,"DeleteAtom",{3,4,5}]

ConnectedMoleculeComponents@%

TransformAtomCoordinates

m2=MoleculeModify[m,"TransformAtomCoordinates",TranslationTransform[{10,10,10}]]

Show[MoleculePlot3D/@{m,m2}]

MoleculeModify[m,"TransformAtomCoordinates",ScalingTransform[{3,2,4}]]//MoleculePlot3D

ReplaceAtom

MoleculeModify[m,"ReplaceAtom",1->"P"]

MoleculePlot3D/@{m,%}

AddBond

A bond can be added between existing atoms in the molecule, but it’s very likely to lead to a broken structure

MoleculeModify[m,"AddBond",Bond[{12,16},"Single"]]

AddBond can also be used to add an atom

MoleculeModify[m,"AddBond",Bond[{6,"O"},"Single"]]

RenumberAtoms

m=Molecule[];
m2=MoleculeModify[m,"RenumberAtoms",RandomSample@Range@AtomCount@m]

AtomList[#,_->"AtomicNumber"]&/@{m,m2}

MoleculeEquivalentQ[m,m2]

ReplaceSubstructure

m=Molecule@"CCC[C@H](OC(=O)CCC(C)C)C"

MoleculeModify[m,"ReplaceSubstructure",MoleculePattern["C(=O)O"]->"C(F)C"]

MoleculeModify[
Molecule[],"ReplaceSubstructure",MoleculePattern["c1ccccc1"]->"C1CCCCC1"]

ReplaceSubstructureList

m = Molecule@"trinitrotoluene"

MoleculeModify[m,"ReplaceSubstructureList",MoleculePattern["[#7](~[#8])~[#8]"]->"C"]

DeleteDuplicates[%,MoleculeEquivalentQ]

RemoveHydrogens

m = Molecule[];

AtomCount[m]

AtomCount[m,Atom[Except["H"]]]

Length@First@m

m2=MoleculeModify[m,"RemoveHydrogens"]

Length@First@m2

MoleculePlot/@{m,m2}//Column

MoleculeEquivalentQ[m,m2]

AddHydrogens

m3 = MoleculeModify[m2,"AddHydrogens"]

m3===m

MoleculeEquivalentQ[m,m2,m3]

MoleculeModify - Norbornane to Camphor

{norbornane,camphor}=Molecule/@{"norbornane","camphor"};
Row[{MoleculePlot@norbornane,Style[" -> ",20],MoleculePlot@camphor}]

MoleculePlot3D/@%

res=norbornane

res = MoleculeModify[res,
	"AddBond",
	{
		Bond[{"O", 2}, "Double"],
		Bond[{"C", 7}, "Single"],
		Bond[{"C", 7}, "Single"],
		Bond[{"C", 1}, "Single"]
	}
]

MoleculeEquivalentQ[res,camphor]

MoleculeModify - fun example

Take a simple starting point - cyclopentane

m = Molecule["cyclopentane"]

SetOptions[MoleculePlot3D,PlotTheme->"StickPlot"];

MoleculePlot3D[m]

Now replace every hydrogen atom with a benzene ring

m2=MoleculeModify[m,"ReplaceSubstructure","H"->Molecule@"benzene"]

MoleculePlot3D[m2]

Now say we want to replace the central 5-membered ring with a 10-membered ring.  This requires setting multiple attachment points.

m3 = MoleculeModify[m2,
	"ReplaceSubstructure",
	MoleculePattern["C1CCCC1"] -> Molecule["cyclodecane"],
	"AttachmentPoints" -> <|
		1 -> 1, 2 -> 3, 3 -> 5, 4 -> 7, 5 -> 9
	|>
]

Row @ {
	MoleculePlot[m2, MoleculePattern @ "C1CCCC1"],
	MoleculePlot[m3, MoleculePattern @ "C1CCCCCCCCC1"]
}

Here we can just highlight the fact that you can extract part of a molecule to get a new molecule:

Molecule[m2,Range[5]]

MoleculePlot3D/@{m3,m2}//Column

m3["SMILES"]

Embedding molecules in 3D

Automatic coordinates

A molecule expression can contain coordinates directly,

m=Molecule[{Atom["C"],Atom["C"],Atom["C"],Atom["H"],Atom["H"],Atom["H"],Atom["H"],Atom["H"],Atom["H"],Atom["H"],Atom["H"]},{Bond[{1,2},"Single"],Bond[{2,3},"Single"],Bond[{1,4},"Single"],Bond[{1,5},"Single"],Bond[{1,6},"Single"],Bond[{2,7},"Single"],Bond[{2,8},"Single"],Bond[{3,9},"Single"],Bond[{3,10},"Single"],Bond[{3,11},"Single"]},AtomCoordinates->QuantityArray[]]

Or the coordinates can be automatically generated

m2=Molecule[{Atom["C"],Atom["C"],Atom["C"],Atom["H"],Atom["H"],Atom["H"],Atom["H"],Atom["H"],Atom["H"],Atom["H"],Atom["H"]},{Bond[{1,2},"Single"],Bond[{2,3},"Single"],Bond[{1,4},"Single"],Bond[{1,5},"Single"],Bond[{1,6},"Single"],Bond[{2,7},"Single"],Bond[{2,8},"Single"],Bond[{3,9},"Single"],Bond[{3,10},"Single"],Bond[{3,11},"Single"]}]

MoleculePlot3D/@{m,m2}//Show

ComputeAtomCoordinates

Automatic embedding is performed using the method described in J. Chem. Inf. Comp. Sci. 55:2562-74 (2015)

Compute a distance bounds matrix, based on the chemical graph

Find a random distance matrix that satisfies these bounds.

Embed the molecule - find 3D coordinates to match the random distance matrix

Correct the structure using torsion angle preferences trained on data in the Cambridge Structural Database

m = Molecule[];

m2=MoleculeModify[m,"ComputeAtomCoordinates"]

Options@m2

The embedded coordinates are the same as the automatically generated coordinates

m["AtomCoordinates"]===m2["AtomCoordinates"]

ComputeAtomCoordinates - random conformer generation

By providing a RandomSeeding option, we can generate

mcopies =Table[ MoleculeModify[m,"ComputeAtomCoordinates",RandomSeeding->Automatic],100];//AbsoluteTiming

MoleculePlot3D/@mcopies[[;;10]]//Show

#["MMFFEnergy"]&/@mcopies//Histogram

mcopies2=MoleculeModify[#,"EnergyMinimizeAtomCoordinates"]&/@mcopies;//AbsoluteTiming

#["MMFFEnergy"]&/@mcopies2//Histogram

Select[mcopies2,QuantityMagnitude@#["MMFFEnergy"]<.2&]

Map[
#["MMFFEnergy"]&,
Select[mcopies2,QuantityMagnitude@#["MMFFEnergy"]<82&]
]//Histogram

Map[
#["MMFFEnergy"]&,
Select[mcopies2,QuantityMagnitude@#["MMFFEnergy"]>82&]
]//Histogram

AssignStereochemistryFromCoordinates

If a molecule has coordinates but no defined StereochemistryElements, they can be computed

m1=Molecule["C[C@@H]1C[C@H](C1)C(=O)OC(=O)CC1(C)CCC1"]//MoleculeModify["ComputeAtomCoordinates"]

m2=MoleculeModify[m1,"RemoveStereochemistry"]

m3=MoleculeModify[m2,"AssignStereochemistryFromCoordinates"]

MoleculeEquivalentQ[m1,m3]

EnergyMinimizeAtomCoordinates

Using an implementation of the Merck Molecular Force Field,

mol = Molecule["water"]

mol["BondAngle",{2,1,3}]

mol2=MoleculeModify[mol,"EnergyMinimizeAtomCoordinates"]

mol2["BondAngle",{2,1,3}]

mol["MMFFEnergy"]

mol2["MMFFEnergy"]

waters=Table[MoleculeModify[mol,"ComputeAtomCoordinates",RandomSeeding->Automatic],1000];

Histogram[#["BondAngle",{2,1,3}]&/@waters]

waters2=MoleculeModify[#,"EnergyMinimizeAtomCoordinates"]&/@waters;

Histogram[#["BondAngle",{2,1,3}]&/@waters2]

Geometry modifications - setting bond lengths, bond angles, and torsion angles

SetBondLength

m=Molecule[];

m["AtomCoordinates"]

MoleculePlot3D[MoleculeModify[Molecule[],"SetBondLength",{19,20} ->Quantity[#,"Angstroms"]]]&/@{1,5,20}

SetBondAngle

mols=MoleculeModify[Molecule[],"SetBondAngle",{7,1,8} ->Quantity[#,"AngularDegrees"]]&/@{15,120,175};

MoleculePlot3D[#,Style[{7,1,8},Green]]&/@mols

Geometry example - ethane torsional potential

Define ethane and get the atom indices for the torsion angle

m =MoleculeModify["ComputeAtomCoordinates"]@ Molecule["ethane"];
FindMoleculeSubstructure[m,MoleculePattern["[H]CC[H]"]]
indices = Keys@First@%

A simple torsional potential is found just by adjusting the torsion angle holding all other degrees of freedom constant

rigidRotation=Table[
mol=MoleculeModify[m,"SetTorsionAngle",indices->ang];
mol["MMFFEnergy"],
{ang,0,360,5}];
ListLinePlot[rigidRotation,DataRange->{0,360}]

A better potential can be found by doing a constrained energy minimization, allowing all other degrees of freedom to relax while holding fixing the torsion angle.

We do this by adding an additional term to the force field, with a large force constant to be sure to find the right state

relaxedRotation = Table[
mol=MoleculeModify[m,"EnergyMinimizeAtomCoordinates","GeometryConstraints"-><|"TorsionConstraints"->{<|"Atoms"->indices,"Value"->ang,"ForceConstant"->10^5|>}|>];
mol["MMFFEnergy"],
{ang,0,360,5}];
ListLinePlot[{rigidRotation,relaxedRotation},DataRange->{0,360}]

File Import - partial access to multiframe SD files

Chemical datasets often are released as large text files in the SDF format

For example, the ChEMBL 23 database contains 1.7 million molecules, with 2D coordinates for each, and can be download from https://www.ebi.ac.uk/chembl/downloads

FileByteCount@"/Users/jasonb/Downloads/chembl_23.sdf"

mols =Import["/Users/jasonb/Downloads/chembl_23.sdf",{"Molecule",10;;1000000;;100}];//AbsoluteTiming

Length@mols

In just 60 seconds, we’ve imported 10,000 molecules from a 3.6GB file

SeedRandom[401];
m=RandomChoice@mols

MoleculePlot@m

Substructure filtering - carboxylic acid

We can search for carboxylic acids within the dataset by defining a pattern for the COOH substructure

carboxylicAcidPattern = MoleculePattern[{"C","O",Atom["O","HydrogenCount"->1]},{Bond[{1,2},"Double"],Bond[{1,3},"Single"]}]

MoleculeContainsQ[m,carboxylicAcidPattern]

acids = Select[mols,MoleculeContainsQ[#,carboxylicAcidPattern]&];

Length@acids

Just over ten percent have this substructure

m = RandomChoice @ acids

MoleculePlot[RandomChoice @ acids,carboxylicAcidPattern]

MoleculePlot[RandomChoice @ acids,carboxylicAcidPattern]

Substructure filtering - using a custom pattern

m=RandomChoice@mols

MoleculePlot[m]

We can define a substructure by a set of indices - here the fused ring system has indices 1 through 9 (read off via the tooltip).

fusedRingPattern = MoleculePattern[m,{1,2,3,4,5,6,7,8,9}]

MoleculePlot[m,fusedRingPattern]

Now we search for this pattern in the larger set

filtered = Select[mols,MoleculeContainsQ[#,fusedRingPattern]&];

Length@filtered

MoleculePlot3D[RandomChoice@filtered,fusedRingPattern]

Future Topics - Chemical Reactions and Electronic Orbitals

Chemical Reactions

Need for a symbolic representation for a chemical reaction

Pattern-based reactions, which can be applied to multiple groups of reactants

Specific reactions

Storing metainformation - reaction kinetics, spectator compounds

Quantum chemistry

Computing and plotting electronic orbitals

Reading output files for common quantum chemistry packages

Stereoisomer enumeration

With an option to attempt embedding, to weed out implausible structures

Tautomer standardization

Tautomer enumeration

Experimental functions

These functions will need a bit more polishing before being promoted to the System context

?Chemistry`Experimental`*

ResonanceStructures

mol = Molecule@"phenanthrene"

Chemistry`Experimental`ResonanceStructures[mol]

DeleteDuplicates[%,MoleculeEquivalentQ]

MaximumCommonSubstructure

mols = Molecule /@{"C[C@@H]1C[C@H]1OC(=O)[C@H]1CCC1(C)C","CCC(C(=S)C1CC(C1)(C)C)(C)C"}

patt=Chemistry`Experimental`MaximumCommonSubstructure@@mols

MoleculePlot[#,patt]&/@mols

SymmetryInformation

This uses the method described in this paper to find molecular point groups.  This may be exposed as a Function Repository resource, or as an internal function

mol = Molecule["tetramethylAmmonium"]//MoleculeModify["EnergyMinimizeAtomCoordinates"]

sym=Chemistry`Experimental`SymmetryInformation[mol]

Cases[sym,assoc:KeyValuePattern["SymmetryPlane"->_]:>Show[MoleculePlot3D[mol],Graphics3D@assoc["SymmetryPlane"]],Infinity]

Fingerprint similarity

Molecule fingerprints can be used for similarity comparison, substructure filtering, and will soon be hooked up as named NetEncoders

m = Molecule[];

m2=Molecule[];

Chemistry`Experimental`MoleculeFingerprint[m,"FingerprintType"->#]&/@{"RDKit","AtomPairs","TopologicalTorsions","MACCSKeys","MorganConnectivity","MorganFeature"}

fp1=Chemistry`Experimental`MoleculeFingerprint[m];
fp2=Chemistry`Experimental`MoleculeFingerprint[m2];

1.-JaccardDissimilarity[fp1,fp2]

Chemistry`Experimental`MoleculeSimilarity[m,m2]

Using this as a distance function, we can define a Nearest function:

SeedRandom[42];
smiles = ;

Take 1000 random molecules, generate and store the fingerprints to create a nearest function:

mols = Table[
Molecule[smiles[[n]],IncludeHydrogens->False],
{n,Length@smiles}];~Monitor~n

mnf = Chemistry`Experimental`MoleculeNearest[mols,"FingerprintType"->"MACCSKeys", "SimilarityMeasure"->"Dice"]

mnf[m2,3]

code

This is a text cell.

Quit

Molecule
With[{
		smilesList = ,
		names = 
	},
	randomSmiles[] := RandomChoice[smilesList];
	randomName[] := RandomChoice[names];
]
(* symbol should have been removed *)
Chemistry`Experimental`SMARTSPattern = MoleculePattern;