Computing the arrows of chemical reactions

Pietro Vidossich and AgustiLledós

Departament de Química, UniversitatAutònoma de Barcelona, 08193 Cerdanyola del Vallés, Spain

Supporting Information

Thestepsrequiredtocomputethearrowsofachemicalreactioninclude:

  1. Findingthetransitionstate(TS)ofthereaction,
  2. ComputingtheIntrinsicReactionCoordinate(IRC)ofthereaction,
  3. ForeachgeometryalongtheIRC,computingthelocalizedorbitals,
  4. Foreachlocalizedorbital,computetheorbitalcentroid,
  5. Visualizingresults

Inthefollowing,weprovide,stepbystep,therequiredinputfilesforthecalculationsandthescriptstoperformtheanalysis.Wewillusethewell-knownaldolreactionasanillustrativeexample.WepresentcalculationsforthepopularcomputationalquantumchemistrycodeGaussian(versionGaussian09[1]). Calculationsare performed at the B3LYP/6-31G* level.WiththehelpoftheGaussianuser'sreferenceguidethestudentmayeasilyunderstandandincasemodifythekeywordsoftheinputfiles appearing below.

Step1.Findthetransitionstate(TS)

Thealdolreactionhasbeenthesubjectofseveralcomputationalstudies.Inarecentwork,ZhangandHouk[2] investigatedthepossibilityofpurewaterfacilitatingtheketo-enolisomerization. ThestudentmaydownloadthepaperandtheassociatedSupportingInformation,wheretheauthorsreportedthecoordinatesofthetransitionstate.Zhangetal.basedthecalculationontheB3LYP/6-311++G-(3d,3p)levelandtheCPCMsolvationmodel. Here, for computational convenience, we use a smaller basis set (6-31G*) and perform calculations in the gas-phase. Thus, firstwe computetheTSattheB3LYP/6-31G*level of theory.

Gaussianinputfile:

%nproc=1

%mem=2GB

#pB3LYP/6-31G*opt=(ts,noeigentest,modredundant,calcfc)Freq

ts

01

81.2127151.206352-0.341895

61.250570-0.091791-0.255667

60.172155-0.851059-0.719697

1-0.358873-0.494516-1.595360

10.205959-1.930187-0.598447

6-1.2442700.0416920.476195

8-1.1399451.2849380.175412

10.1458601.469589-0.187330

6-2.490310-0.7070470.041660

1-2.743598-0.466695-0.994527

1-2.385199-1.7895590.161183

1-3.321193-0.3712900.676849

62.408003-0.6691690.510120

12.455147-0.2252081.511521

12.348539-1.7572470.588947

13.342780-0.4042030.001571

1-0.848475-0.2967621.451795

Executiontimewas6 min 27 secseconasingle(1)coreofourlocalclusterbasedonAMD64nodes(16corespernode,64GBRAMpernode).

ThestudentmayopentheGaussianoutputfilewithavisualizationprograme.g.MOLDEN(freelyavailableat

Step2.ComputetheIntrinsicReactionCoordinate(IRC)

CopyandpastetheoptimizedTSgeometryinanewinputfilewhichrequeststhecalculationoftheIRC.

Gaussianinputfile:

%nproc=1

%mem=2GB

#pB3LYP/6-31G*irc=(calcAll,maxpoints=40)

irc

01

8 1.290100 1.151980 -0.388876

6 1.257622 -0.149111 -0.206514

6 0.143858 -0.904902 -0.543520

1 -0.384038 -0.709482 -1.473993

1 0.091286 -1.950187 -0.286397

6 -1.294168 0.117313 0.421290

8 -1.153401 1.263502 -0.114205

1 0.219123 1.584706 -0.460692

6 -2.541036 -0.659734 0.097009

1 -2.853650 -0.505712 -0.946387

1 -2.429225 -1.736351 0.260269

1 -3.375158 -0.306922 0.722275

6 2.473421 -0.682482 0.466746

1 2.229567 -1.152379 1.431289

1 2.967656 -1.450087 -0.148459

1 3.215167 0.109631 0.663701

1 -0.847258 -0.071945 1.408795

Executiontimewas8 h32 minonasinglecoreofourlocalcluster. The execution time may be reduced in a parallel run, or reducing the number of IRC point to compute.

AttheendoftheGaussianoutputfile,asummaryofenergyofthesystemasafunctionoftheIRCisgiven(searchthefilefor “Summaryofreactionpathfollowing”).ThestudentmayplottheenergyprofileandvisualizethechangeofgeometryfromreactantstoproductsalongtheIRC.

Step3.ComputelocalizedorbitalsforeachgeometryalongtheIRC

TheGaussianinputfilefororbitallocalizationisstructuredasacompoundjob,inwhichfirstthewavefunctionisoptimizedforagivengeometryandthentheorbitalsarelocalized.TheIOp(4/9)specifiesoptionsforthelocalizationalgorithm,specificallyIOp(4/9=10212)asksfortheBoyslocalizationcriterion.

Gaussianinputfile:

%nprocshared=1

%Chk=aldol

%mem=2GB

#pB3LYP/6-31G*Symmetry=NoneNoSymm

scf

01

81.2737361.099808-0.348871

61.244837-0.193737-0.166557

60.120823-0.950417-0.498671

1-0.389216-0.763746-1.441830

10.073798-1.996753-0.242804

6-1.3001880.0588480.456834

8-1.1656951.216876-0.076485

10.1597571.532050-0.403814

6-2.554781-0.7085380.136533

1-2.867981-0.553785-0.906597

1-2.445984-1.7857130.299799

1-3.386922-0.3531900.762649

62.458942-0.7327650.505798

12.214101-1.2026431.470075

12.952626-1.500712-0.109025

13.2018010.0583750.703631

1-0.865091-0.1174141.451820

--Link1--

%nprocshared=1

%Chk=aldol

%mem=2GB

#pB3LYP/6-31G*Geom=CheckGuess(NoSymm,Read,Local,Save,Only)Pop=FullSymmetry=NoneNoSymm

IOp(4/9=10212)

localization

01

AcalculationsuchastheaboveistobeperformedforeachgeometryalongtheIRC, including the TS (the execution time is <1 min on a single core).Actually,inourexperience,lessthanavailablegeometriesmaybeenoughtovisualizearrows.OnaLinuxworkstationandusingthebashscriptinglanguage,geometriesmaybeextractedusingacommandsuchas:

grep-A23"PointNumber: 1 PathNumber: 1"$1|tail-n17|awk'{print $2,$3,$4,$5}'tmp.xyz

Wheretmp.xyzisatemporaryfileinwhichthecoordinatesofPointNumber1alongPathNumber1arestored.Thenumbers “23” and “17” intheabovecommandlinearetailoredforthenumberofatomsinthepresentexample.Allothergeometriesmaybeextractedchangingthe “PointNumber” and “PathNumber”.Wesuggestperformingcalculationsinseparatedirectories,oneforeachgeometry.

4.Computetheorbitalcentroids

ForeachgeometryalongtheIRC,thecentroidsofthelocalizedorbitalshavetobecomputed.Thelocalizedorbitalsarestoredinthe.chkfile.Thisfilehastobeprocessedinordertoextracttheorbitalcentroids.UsingthebashscriptingenvironmentunderLinuxandtheg09utilitiesformchk,cubegenandcubman,theanalysiswouldproceedasfollows(assumingpoint1_path1.chkisthe.chkfilefromtheg09calculation,and 21 isthenumberofvalenceorbitalsinthesystem).

4.1.Formatthe.chkfileusingformchk:

formchkpoint1_path1.chkpoint1_path1.fchk

4.2.Generate.cubefilesforeachlocalizedorbitalusingcubegen(thesemaybevisualizedbysome

visualizationsoftware):

foriin`seq1 21`;doecho$i;cubegen0MO=$ipoint1_path1.fchk$i.cube0h;done

4.3.Computetheorbitaldensityusingthecubmanutility:

foriin`seq1 21`;doprintf'%s\n'SQ$i.cubey$i\squared.cubey|

cubman;done

4.4.Usethecubmanutilitytocomputetheorbitalcentroids(DipAEintheoutput):

foriin`seq1 21`;doprintf'%s\n'P$i\squared.cubey|cubman;done

|grep"DipAE="|awk'{print"X",-$2*0.529177,-$3*0.529177,-$4*0.529177}'

Thelocalizedorbitalcentroidsareprintedonthescreenandmaybecopied

5.Visualizethearrows

Avisualizationsoftwarecapableofvisualizingseveralgeometriesinthesamewindowisrequiredtodisplaythearrows.OnesuchsoftwareisVMD(freelyavailableat few structures around the TS, the VMDvisualizationwindowshoulddisplaythefollowingimage.

Notes: This exercise was prepared with the popular Gaussian (version 09) program [1]. Other programs, such as CPMD ( and CP2K (ww.cp2k.org), perform orbital localization and directly print the localized orbitals centroids. In computing localized orbitals with Gaussian09, the following comments may be of help:

- it may happen g09 fails to localize the orbitals. Changing the orbital localization criterion may be of help.

- we first tried to set up this exercise with the semiempirical method PM6. The method fairly compares with B3LYP/6-31G* results in terms of structures and energetics, but unfortunately g09 failed to localize the orbitals.

- localized orbitals are provided by g09 in the molecule’s “standard orientation”. It may happen, when performing orbital localization, that the code reorients the molecule.

References

[1] Gaussian09,RevisionA.1,Frisch,M.J.;Trucks,G.W.;Schlegel,H.B.;Scuseria,G.E.;Robb,M.A.;Cheeseman,J.R.;Scalmani,G.;Barone,V.;Mennucci,B.;Petersson,G.A.;Nakatsuji,H.;Caricato,M.;Li,X.;Hratchian,H.P.;Izmaylov,A.F.;Bloino,J.;Zheng,G.;Sonnenberg,J.L.;Hada,M.;Ehara,M.;Toyota,K.;Fukuda,R.;Hasegawa,J.;Ishida,M.;Nakajima,T.;Honda,Y.;Kitao,O.;Nakai,H.;Vreven,T.;Montgomery,Jr.,J.A.;Peralta,J.E.;Ogliaro,F.;Bearpark,M.;Heyd,J.J.;Brothers,E.;Kudin,K.N.;Staroverov,V.N.;Kobayashi,R.;Normand,J.;Raghavachari,K.;Rendell,A.;Burant,J.C.;Iyengar,S.S.;Tomasi,J.;Cossi,M.;Rega,N.;Millam,J.M.;Klene,M.;Knox,J.E.;Cross,J.B.;Bakken,V.;Adamo,C.;Jaramillo,J.;Gomperts,R.;Stratmann,R.E.;Yazyev,O.;Austin,A.J.;Cammi,R.;Pomelli,C.;Ochterski,J.W.;Martin,R.L.;Morokuma,K.;Zakrzewski,V.G.;Voth,G.A.;Salvador,P.;Dannenberg,J.J.;Dapprich,S.;Daniels,A.D.;Farkas,Ö.;Foresman,J.B.;Ortiz,J.V.;Cioslowski,J.;Fox,D.J.Gaussian,Inc.,WallingfordCT,2009.

[2] Zhang X, Houk KN (2005) J. Org. Chem. 70:9712