Message Boards

GROUPS:

Staff Picks Chemistry Data Science Material Sciences Curated Data Wolfram Language Machine Learning Wolfram Summer School Neural Networks

Sathvik Ajay Iyengar

[WSS20] Neural Networks for Software-Interface Molecular Vibration Data

Sathvik Ajay Iyengar, Rice University

Posted 1 year ago

2937 Views

2 Replies

13 Total Likes

Follow this post

enter image description here

Deep learning is becoming commonplace in many fields of science and technology. The ability to identify and predict patterns in large data sets is a powerful tool, and can in fact, replace many brute-force and calculation/simulation intensive methods. Traditionally, identifying molecular vibrational modes, and hence the corresponding IR-spectra for organic molecules is a non-trivial process.Theoretically, it is performed by specialized ab initio quantum chemistry software, and is a time-consuming process. Experimentally, it is determined by specialized spectrometers. This project is aimed at predicting such inherent patterns between parametrized molecular physical/chemical properties, and their vibrational modes.
Here, we focus on two major aspects. First, we introduce a method to write, execute and animate input and data files to a third party quantum chemistry software (Orca) through the Wolfram Language (WL). Next we use this freshly generated vibrational spectra data to try and predict patterns using a variety of machine learning functionalities.

Introduction

An infrared or vibrational spectrometer is an instrument that experimental chemists, materials scientists and other researchers use to identify or “characterize” chemical substances. When the molecules of an unknown chemical substance are made to interact with infrared radiation, the molecules absorb the radiation, shift to an excited state, and dissipate energy in the form of molecular vibrations. These vibrations are either symmetric or antisymmetric, giving rise to many sub-varieties (stretching, scissoring, rocking, wagging, twisting, out-of-plane). While the near-IR is responsible for molecular vibrations, the far-IR corresponds to more rotation dominant behavior.

In[]:=

Style/@Text/@TextSentences[WikipediaData["Infrared spectroscopy"]][[;;3]]//Column

Out[]=

Infrared spectroscopy (IR spectroscopy or vibrational spectroscopy) involves the interaction of infrared radiation with matter.

It covers a range of techniques, mostly based on absorption spectroscopy.

As with all spectroscopic techniques, it can be used to identify and study chemical substances.

Now, the question is, given enough descriptor parameters for a molecule, or as it is commonly known as a "molecular fingerprint" (discussed further below), can we predict the vibrational modes that it would exhibit? Vibrations are the result of chemical bonds executing harmonic oscillations, and are highly specific to the type of atoms involved in the bond, and the overall structural and energy landscape. However, they are also not trivial, and are commonly calculated using quantum chemistry software that rely on Density Functional Theory (DFT) principles such as the Born - Oppenheimer (BO) approximations, among others. Hence, this is a perfect example of hard-to-compute data with a rich set of identifiable patterns. This project deals with the generation of such data, and the construction of neural network models to best predict vibrational modes, given a unique molecule.

This is what an IR spectra plot for an organic molecule “aniline” looks like. The frequency, commonly measured in wavenumber (1/cm), has a corresponding value of transmittance.

In[]:=

WebImageSearch["Infrared spectroscopy"];

In[]:=

%[[3]]//Magnify[#,0.75]&

Out[]=

The lower the transmittance, the greater the absorption and hence the vibrational mode. Thus, we can characterize the “most significant” frequencies as the ones with the lowest transmittance values.

Highlighted in yellow are the bonds that exhibit the various stretching modes as shown above for aniline, in decreasing order of frequency:

In[]:=

Function[bond,{Image@MoleculePlot3D[Molecule@"Aniline",{Bond[First@bond,"Single"|"Aromatic"]}],Style[Text@bond/.{{a_,b_},c_}a<>"-"<>b<>" "<>c,FontSize20]}]/@{{{"N","H"},"~3400

-1

"},{{"C","H"},"~3000

-1

"},{{"C","C"},"~1500

-1

"},{{"C","N"},"~1300

-1

"}}//Transpose//Grid[#,FrameAll]&

Out[]=


N-H ~3400 -1 cm	C-H ~3000 -1 cm	C-C ~1500 -1 cm	C-N ~1300 -1 cm

Even molecules have fingerprints~

We can identify each molecule with a unique fingerprint in the form of the “radial distribution function”(RDF), a model that describes the variation of atom density as a function of distance from a reference atom:

In[]:=

WebImageSearch[SearchQueryString@"Radial distribution function dr"][[{4,7}]]//ImageCollage[#,Automatic,400]&

Out[]=

RDF values in the WL are a set of discrete values in a {210} numeric array, which is based a list of molecular descriptors calculated by Dragon[1, 2]

In[]:=

Molecule[#]["RDF"]&/@{"ethanol","acetic acid","acetone"}

Out[]=

NumericArray

Type: Real64

Dimensions: {210}

,NumericArray

Type: Real64

Dimensions: {210}

,NumericArray

Type: Real64

Dimensions: {210}



Although the differences seem subtle, they are based on highly sensitive parameters.

In[]:=

MatrixPlot[Partition[Normal@#,9],FrameFalse]&/@%//Grid[{#,{"ethanol","acetic acid","acetone"}},FrameAll]&

Out[]=


ethanol	acetic acid	acetone

The RDF in the WL has 7 crucial parts (1 unweighted and 6 weighted) that make it ideal for fingerprinting and “deploy-worthy” in NNs:

Pure radial distribution parameter

Molecular mass

Van der Waals volume

Sanderson electronegativity

Polarizability

Ionization potential

I-state (intrinsic-state from local vertex variants)

Data sets- Let’s automate things!

Now, we need to find a suitable software that will generate our dataset for us. Among a wide range of software that can do this, I found Orca[3] to be the most suitable. It has a very simple installation procedure, and is windows-compatible. An easy to follow documentation is a huge plus.

In[]:=

Style/@Text/@TextSentences[WikipediaData["ORCA (quantum chemistry program)"]][[;;2]]//Column

Out[]=

"ORCA is an ab initio quantum chemistry program package that contains modern electronic structure methods including density functional theory, many-body perturbation, coupled cluster, multireference methods, and semi-empirical quantum chemistry methods."

"Its main field of application is larger molecules, transition metal complexes, and their spectroscopic properties."

An exciting bonus here is that not only can Orca generate data sets to do IR spectra (vibrational; absorption) analysis, but it can also be extended to study and predict patterns in Raman (vibrational; scattering), UV-vis (electro-optical), and NMR (magnetic) spectra.

General notes:

Computation time: ~15 min for 20 atom system [geometry optimization + spectra calculations]

Command line run software

From molecule name to spectra, all on our dear Mathematica

We will perform execution and acquisition of data from Orca all through the WL. This will help internalize our entire process of data set generation. Orca requires 2 types of files to calculate the Hessian values for various organic molecules. They are the input and the data files. The data file is fairly straightforward; using the built-in function Molecule[], we can export .xyz data. The input file is very short and describes the method of energy minimization etc. which can be written as a string expression in the WL. First we download the software.

From here, code can be run only if Orca and its executable files are installed at "C:\\ Orca" from the link above.

(*CreateDirectory[FileNameJoin[{NotebookDirectory,"orca data"}]]*)

We can acquire ordered molecule name lists from the “GDB-9” database. We will leave it uncompressed until necessary, to save time.

In[]:=

dataGDB=CloudGet@CloudObject["https://www.wolframcloud.com/obj/jasonb/Compressed_GDB-9_molecules"];

Next, we shall consider smaller molecules (<=20 atom systems) so that the computation time does not overflow, while at the same time ensuring that we still have a diverse representation of functional groups (and hence a diverse representation of vibrational frequencies). Additionally, some uncommon molecules do not have IUPAC name in the PubChem service, and hence those are removed.

newMoldata[n_,m_]:={#,Molecule@#}&/@Last/@First/@(Select[(ServiceExecute["PubChem","CompoundProperties",{"Molecule"#,"Property"{"IUPACName"}}]&/@Select[Uncompress/@dataGDB[[n;;m]],BondCount@#≤20&]),!StringContainsQ[ToString@FullForm@#[[1]],"Missing"]&])

Each time we wish to add data for new molecules, we take a part of the GDB database in serial order (serial 81 through 120 in this example, meaning the first 80 have been completed and saved in the “orca data” folder in the NotebookDirectory):

In[]:=

new=newMoldata[121,140];

As the quantum chemistry simulations are time-consuming, we should perform a sanity check to see if some molecules have already been computed, then define a list of fresh molecules to compute:

In[]:=

currentMoldata={#,Molecule@#}&/@StringReplace[(ToLowerCase/@StringCases[#,"data"~~_~~x__~~"_out_"->x]&/@FileNames[___~~"_out_IR.txt",FileNameJoin[{NotebookDirectory[],"orca data"}]]//Flatten),StartOfString~~"n-""N-"];

A visual check if we have diversity in the types of organic bonds:

In[]:=

Image/@(MoleculePlot/@RandomSample[%[[All,2]],40])//ImageCollage[#,Automatic,{1000,1000},BackgroundWhite]&//Magnify[#,0.25]&

Out[]=

List of molecules to compute (example run):

In[]:=

toDo=Select[new,!IntersectingQ[List[#[[2]]],currentMoldata[[All,2]]]&][[All,1]];

In[]:=

%//Partition[#,5]&//Grid

Out[]=

pent-1-yne	butanenitrile	2-(methylamino)acetonitrile	3-methoxyprop-1-yne	2-methoxyacetonitrile
but-3-yn-1-ol	3-hydroxypropanenitrile	butanal	N-ethylformamide	ethyl formate
2-methoxyacetaldehyde	3-hydroxypropanal	pentane	butan-1-ol	1-methoxypropane

And here we automate the process:

In[]:=

molToIR[mol_String]:=Module[{name=StringReplace[mol," """]},(*Writinggeometryoptimizationinputfileandxyzdatafile*)MapThread[Export[FileNameJoin[{NotebookDirectory[],"orca data",name<>#1}],#2]&,{{"_in_geo.txt","_guess.txt"},{"!B3LYP DEF2-SVP OPT D4"<>"\n"<>"* xyzfile 0 1 "<>name<>"_guess.txt"<>"\n",ExportString[Molecule[mol],"XYZ"]}}];(*Runninggeometryoptimization;forCaffeine-TOTALRUNTIME:0days0hours19minutes44seconds333msec*)RunProcess[{"C:\\Orca\\orca.exe",name<>"_in_geo.txt"},"StandardOutput",ProcessDirectory->FileNameJoin[{NotebookDirectory[],"orca data"}]];(*Geometryoptimizedfileisnowthenewdatafile*)Export[FileNameJoin[{NotebookDirectory[],"orca data",name<>"_in_IR.txt"}],"!BP86 DEF2-SVP FREQ"<>"\n"<>"* xyzfile 0 1 "<>name<>"_in_geo.xyz"<>"\n"];(*RunningandexportingIRspectraprediction;forCaffeine-TOTALRUNTIME:0days0hours14minutes7seconds203msec*)MapThread[Export[FileNameJoin[{NotebookDirectory[],"orca data",name<>#1}],#2]&,{{"_out_full.txt","_out_IR.txt"},Function[s,{s,StringTake[s,First/@(First/@StringPosition[s,#]&/@{"IR SPECTRUM","The first"})-1]}]@RunProcess[{"C:\\Orca\\orca.exe",name<>"_in_IR.txt"},"StandardOutput",ProcessDirectory->FileNameJoin[{NotebookDirectory[],"orca data"}]]}]]

In[]:=

molToIR["caffeine"]

Out[]=

{C:\Users\Sathvik Iyengar\Desktop\WSS\Sathvik\orca data\caffeine_out_full.txt,C:\Users\Sathvik Iyengar\Desktop\WSS\Sathvik\orca data\caffeine_out_IR.txt}

Our Data are on Cloud (nine)

We will need to access this data often. It will be saved in the form:

{ { {molecule name1, RDF1}, { {frequency (1, 1), relative transmittance (1, 1)}, {}, ... }, { {molecule name2, RDF2}, ... }, ...}

In[]:=

irData=Function[name,(*moleculename*){{#,Molecule[First@#]["RDF"]}&@StringReplace[ToLowerCase@StringCases[name,"data"~~_~~x__~~"_out_"->x],StartOfString~~"n-""N-"],(*{Frequency,Relativetransmittance}*)ToExpression/@StringSplit[#,WhitespaceCharacter..]&/@StringCases[#,_~~_~~":"~~WhitespaceCharacter..~~(x:(NumberString~~WhitespaceCharacter..~~NumberString))x]&@Import[name]}]/@FileNames[___~~"_out_IR.txt",FileNameJoin[{NotebookDirectory[],"orca data"}]]//Quiet;

However, this data is not ready yet. We need to do some cleaning, and make it variable by choice by introducing a parameter:

Some peaks are very insignificant and can be discarded (very high relative transmittance values)

Some values fall outside the usual IR frequency range of 400-4000 cm-1, we can specify a new range to restrict our search

Finally, some peaks are extremely close, separated by <10 cm-1. A parameter is used to set that frequency width and return the most prominent peak within it

The molecules are also in name order, better to RandomSample to be safe.

In[]:=

irDataFilter[transSpec_,freqSpec_:30,freqRange_List:{400,4000}]:=Module[{newData=irData},(*settingatransmittancevaluethreshold*)newData[[All,2]]=Function[data,Select[data[[2]],#[[2]]<transSpec&]]/@irData;(*texttopreventsettingparameterstooextreme,atleast1frequencymustremain*)If[MemberQ[Length/@newData[[All,2]],0],Return@Print["Transmittance threhold too low as some molecules have no peaks at this value"],Continue];(*settingafrequencyrange*)newData[[All,2]]=Function[data,Select[data[[2]],Between[#[[1]],freqRange]&]]/@newData;(*settingaparametertoexcludeoverlapswithinthefrequencywidth,eg:freqSpec=30willgivethepairwithlowertransmittancebetween1310cm-1and1320cm-1astheyare10cm-1apart*)newData[[All,2]]=Table[First/@ResourceFunction["MaximalSubsets"][SortBy[#,Last]&/@(Function[list,Nearest[list,#,{Infinity,freqSpec}]&/@list]@newData[[All,2]][[n]])],{n,Length@newData[[All,2]]}];newData[[All,2]]=Function[i,SortBy[i[[2]],#[[1]]&]]/@newData;newData[[All,2]]=(#+RandomReal[{0.0001,0.001},Dimensions@#])&/@irData[[All,2]];(*texttopreventsettingparameterstooextreme,atleast1frequencymustremain*)If[MemberQ[Length/@newData[[All,2]],0],Return@Print["Either the frequency exclusion width too high or the frequency range is too narrow as some molecules have no peaks at this value"],Continue];RandomSample@newData]

In[]:=

filteredData=irDataFilter[30];

In[]:=

CloudPut[filteredData,"IR_data_function",Permissions"Public"]

Out[]=

CloudObject[

https://www.wolframcloud.com/obj/sathvik/IR_data_function

]

Alternatively, it can be exported as a “.mx”:

In[]:=

Export[FileNameJoin[{NotebookDirectory[],"irDataSet.mx"}],irData]

Out[]=

C:\Users\Sathvik Iyengar\Desktop\WSS\Sathvik\irDataSet.mx

So you have the numbers, but can you animate it?

Numbers are hard to understand, let us try to animate these molecular vibrations.

In[]:=

irData=CloudGet["https://www.wolframcloud.com/obj/sathvik/IR_data_function"];

And here it goes :

In[]:=

molToGIF[mol_String/;MemberQ[currentMoldata[[All,1]],mol]True]:=Module[{name=StringReplace[mol," """],atomNum,freqN,freq,files,vibData},(*InputfiletograbHessians*)RunProcess[{"C:\\Orca\\orca_pltvib.exe",name<>"_in_IR.hess","all"},"StandardOutput",ProcessDirectory->FileNameJoin[{NotebookDirectory[],"orca data"}]];(*cleaningupthedata*)freqN=Length@irData[[First@@Position[irData,name],2]];freq=irData[[First@@Position[irData,name],2,All,1]];files=Import[#,"Text"]&/@FileNames[name~~"_in_IR.hess.v"~~___,FileNameJoin[{NotebookDirectory[],"orca data"}]][[-freqN;;]];atomNum=First@StringCases[First@files,StartOfString~~x:NumberStringx];(*TabletomakeGIFsforALLvibrationalmodes*)Table[vibData=(StringRiffle[#,"\n"]&/@(Prepend[#,atomNum]&/@Partition[Rest@StringCases[files[[i]],LetterCharacter~~Repeated[Whitespace~~NumberString,3]],ToExpression@atomNum]));(*Tabletoputtogethervariousxyztrajectoriesforasinglevibrationanimation*)Table[Export[FileNameJoin[{NotebookDirectory[],"animations",name<>"vib_"<>ToString@freq[[i]]<>"cm-1_"<>ToString@n<>".xyz"}],vibData[[n]],"Text"],{n,Length@vibData}];(*PlotusingMoleculePlot3Dandexportitasalistofimages*)Export[FileNameJoin[{NotebookDirectory[],"animations",name<>"vib_"<>ToString@freq[[i]]<>"cm-1_"<>".gif"}],Image/@(MoleculePlot3D[#]&/@(Import[First@#,"XYZ"]&/@(FileNames[name~~"vib_"~~ToString@freq[[i]]~~"cm-1_"~~ToString@#~~".xyz",FileNameJoin[{NotebookDirectory[],"animations"}]]&/@Range[20]))),"AnimationRepetitions"∞],{i,Length@files}];]

Let's try this for a random molecule :

molToGIF[random=First@RandomChoice@currentMoldata]

Now we have the data and coordinate points! We need to specify the frame rates and adjust the animation so it can run smoothly [5]

In[]:=

makeAnimation[list_,delayList_:{.03}]:=DynamicModule[{l=Length[list],delays=Abs@Flatten[{delayList}],times,totalTime,delta=.03,frames},times=Round[.5+PadRight[delays,l,delays]/delta];frames=Flatten@Table[Table[list[[i]],{times[[i]]}],{i,Length[times]}];totalTime=Length[frames];EventHandler[Dynamic[frames[[Clock[{1,totalTime,1},totalTimedelta]]],TrackedSymbols{}],{"MouseUp",2}Null]]

There are quite a few vibrational modes for this molecule ~20, so we will just animate 3 random ones. These vibrations are highly accurate as they retain all trajectory data that were calculated :

In[]:=

Column[makeAnimation[Import[#]]&/@RandomSample[FileNames[random~~__~~".gif",FileNameJoin[{NotebookDirectory[],"animations"}]],3]]

Out[]=

Predicting frequencies using self normalizing neural nets (SNNs)

NNs are nit - picky towards data format

First, we try to work on this problem statement- Given a trained molecule (and hence its RDF), and an arbitrary transmission intensity, can we assign the intensity value to a suitable vibrational frequency?

In[]:=

irData=CloudGet["https://www.wolframcloud.com/obj/sathvik/IR_data_function"];

We include copies of random variate data for each molecule to account for variability in intensity

In[]:=

irDataRandom=Table[Table[Function[sd,(Flatten@(Append[Quiet@Normal[irData[[All,1,2]][[k]]],RandomVariate[NormalDistribution[#,sd]]]&@irData[[All,2,All,2]][[k,j]])irData[[All,2,All,1]][[k,j]])]/@Range[0.01(*variatestartrange*),10(*stop*),0.2(*stepsize*)],{j,Length@irData[[All,2,All,2]][[k]]}],{k,Length@irData}];

We split the dataset twice, such that we have “forTrain” which accounts for a “traininset” and a “validationset”, and we have “forTest” which is exclusively kept for testing new molecules.

In[]:=

{forTrain,forTest}=ResourceFunction["TrainTestSplit"][irDataRandom];

In[]:=

Length@forTrain

Out[]=

161

In[]:=

Length@forTest

Out[]=

In[]:=

{trainingset,validationset}=ResourceFunction["TrainTestSplit"][forTrain];{trainingset,validationset}={RandomSample@Flatten@trainingset,RandomSample@Flatten@validationset};

SNN time!

We will use a self-normalizing NN:

In[]:=

SNNnet[dropoutrate_,nhidden_,nlayers_]:=NetChain[Join[Table[NetChain[{LinearLayer[nhidden],ElementwiseLayer["SELU"],DropoutLayer[dropoutrate,"Method""AlphaDropout"]}],nlayers-1],{LinearLayer[1]}]]

In[]:=

test=SNNnet[0.001,200,2]

Out[]=

NetChain



uniniti

alized

Input port:	array
Output port:	vector (size: 1)



Standardize our data:

In[]:=

extractor=FeatureExtraction[N@Keys[trainingset],"StandardizedVector"];trainStandardized=extractor[Keys[trainingset]]Values[trainingset];testStandardized=extractor[Keys[validationset]]Values[validationset];

experiment=NetTrain[test,trainingset,All,TargetDevice"CPU",TimeGoal60,ValidationSetvalidationset];experiment["EvolutionPlots"]

Out[]=

Loss

	rounds
loss

	validation
	training



Clearly, the minimization of loss is not very efficient, let us test this-

In[]:=

trainedNet=experiment["TrainedN