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:
- Findingthetransitionstate(TS)ofthereaction,
- ComputingtheIntrinsicReactionCoordinate(IRC)ofthereaction,
- ForeachgeometryalongtheIRC,computingthelocalizedorbitals,
- Foreachlocalizedorbital,computetheorbitalcentroid,
- 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