Example of GRRM : HCOOH

Effective usage of GRRM11
1. Automated Reaction Route Mapping
2. Quantum Chemical Synthon Analysis
3. Double-end Analysis
4. L-ADDF Search
5. Cluster Structure Sampling
6. Excited-State Analysis
7. Giant-System Analysis with ONIOM
8. Exploration of A + B type Reaction Routes

1

Exploration from a randomly
generated initial structure

Exploration from HCOOH# GRRM/B3LYP/6-31G
01
H-0.452596548000
O-0.452596548000
O-0.452596548000
C-0.452596548000
H-0.452596548000
Options
GauMem=100
GauProc=4This will be copied to the Root Section of a Gaussian Input.
0.034871834140
0.643013717456-1.276412777953-0.053534532177
0.652061759532

2-1.807771705495-1.034349518878
0.242743092264
0.147136623161
0.979506709948

Optimized geometries of HCOOH at B3LYP/6-31G# GRRM/B3LYP/6-31G

This will be copied to the Root Section of a Gaussian Input

01
H
0.0
0.0
0.0
1
All atoms are randomly located, and
O
0.0
0.0
0.0
2
an optimization generates an initial
O
0.0
0.0
0.0
3
C
0.0
0.0
0.0
4
structure.
H
0.0
0.0
0.0
5
Options
NRUN=10 Start from 10 random structures
GauMem=100
GauProc=434EQ_list.log

Output
List of Equilibrium Structures• .log file : History, Amount of Compt., CPU time etc.

# Geometry of EQ 0, SYMMETRY = Cs
H
-0.452596548000
0.034871834140-1.807771705495
O-0.452596548000
0.643013717456-1.034349518878
O-0.452596548000-1.276412777953
0.242743092264
C-0.452596548000-0.053534532177
0.147136623161
H-0.452596548000
0.652061759532
0.979506709948
Energy = -189.686345910766 Electronic Energy in hartree
Spin(**2) = 0.000000000000
ZPVE= 0.033269502369
Zero Point Vibrational Energy in hartree
Normal mode eigenvalues : nmode = 90.013844532 0.017705839 0.041505728 0.043865003 0.063830393
0.075651074 0.114134257 0.376915991 0.486564365• EQ_list.log : List of Equilibrium Structures• EQn.log : History of ADDF from EQn, Structure andEnergy profiles along the ADDF path
• TS_list.log : List of Transition structures• TSm.log : Details of IRC Traces from TSm• DC_list.log : List of Dissociated Structures

Vibrational Eigenvalue fi
Vibrational frequency in cm-1 can be obtained bysqrt(fi /1822.88853006256)*219474.638170777• DCm.log : Meta IRC Traces from DCm
5

61

EQn.log

TS_list.log

PROFILE OF SHS-PATH 1

List of Transition Structures

Initial geometry (negative direction of mode 1)
H-0.273098290958
0.034871834140-1.807771705495
O-0.470810700199
0.643013717456-1.034349518878
O-0.461175866010-1.276412777953
0.242743092264
C-0.434217247418-0.053534532177
0.147136623161
H-0.425702172635
0.652061759532
0.979506709948# STEP 1 E(Harmonic) = 0.001382748031
GENERATION = 1
H-0.427536180483
0.033981498650-1.806811316248
O-0.450816718537
0.643825451874-1.035605196155
O-0.448164828929-1.276628047144
0.242624857766
C-0.470612689786-0.053172482437
0.150438437675
H-0.361723584567
0.639174894802
0.961037157834
energy : -189.685073971923
Spin(**2) : 0.000000000000# STEP 2 E(Harmonic) = 0.015363867007~~~~~~
Another EQ was reached over the TOP of SHS-PATH# Geometry of TS 0, SYMMETRY = C2v
H-0.452579896769-0.739634023903-1.194712914423
O-0.452568723828
0.622484509496-1.029441079169
O-0.452556518725-1.176022414821
0.106195948153
C-0.452547834023
0.119374517004
0.165753891233
H-0.452526289353
0.696601689310
1.079890807874
Energy = -189.621177046091
Spin(**2) = 0.000000000000
ZPVE= 0.028318278608
Normal mode eigenvalues : nmode = 9-0.134645676 0.026856941 0.040118575 0.052256091 0.054695965
0.065556445 0.091473665 0.154739450 0.409107921
CONNECTION : 0 – 0 A TS connecting EQ0 and EQ0~~~~~~
CONNECTION : 1 – 0 A TS connecting EQ1 and EQ0~~~~~~A TS connecting EQ0 and DC, whose
CONNECTION : 0 – DC
structure can be found in TSn.log~~~~~~
CONNECTION : 0 - ?? A TS connecting EQ0 and ??, which is described in

Reached another EQ different from the started EQ 7

TSm.log

8

Quantum Chemical Synthon Analysis

TSm.log(Synthon prediction from decomposition analysis)

INITIAL STRUCTURE
H
0.299046126297
0.824914186341-1.543088134585
O-0.527191800919
0.602385717426-1.071167510316
O-0.483843412292-1.274397964889
0.264284928724
C-0.391162088155-0.064922267476
0.144548168047
H-0.276089959676
0.619520579434
0.994380478018~~~~~~IRC FOLLOWING (FORWARD) STARTING FROM FIRST-ORDER SADDLE~~~~~~IRC FOLLOWING ALONG BACKWARD DIRECTION~~~~~~•

Forward and Backward correspond to either reactant or product•

The last structure of Forward or Backward is EQ, DC, or ??

•DC or ?? in TS_list.log should be confirmed in the last structure of Forward
or Backward in TSm.log•?? Indicates an unsuccessful optimization or an imaginary frequency
9

TS_list.log

Input# Geometry of TS 2, SYMMETRY = CsC
0.229022877916-0.923253337269
0.358038127892
H
0.726068922113-0.942691597054
1.319554214820
H
0.644845702569-1.554721845183-0.416706366822
C-1.053219478319-0.229836209480
0.173955228531
O-0.513955621118
0.894468104516-0.245975228151
O-2.220201824375-0.590109061205
0.354703370666
N
2.264782055458
0.403351849702-0.207377902944
H
2.785759310501
0.105522928710-1.023821667847
H
1.508553545021
1.051792813467-0.402612865724
H
2.857103270848
0.669565839676
0.570429335491
Energy = -284.231895209385
Spin(**2) = 0.000000000000
ZPVE= 0.073957342609
Normal mode eigenvalues : nmode = 24
-0.004846823 0.000187008 0.000440696 0.001744631 0.003440002
0.005361809 0.010240801 0.012068208 0.013136015 0.015585156
0.030332068 0.033688846 0.036637877 0.039287499 0.059837651
0.082224122 0.106827611 0.112198995 0.114951553 0.391070332
0.419509691 0.455191436 0.507963113 0.513310523
CONNECTION : 0 - DC Structure of this DC is found in TS2.log# GRRM/B3LYP/6-31G
01
C
0.692268578114-0.859602740672
H
0.677740993312-1.514459597380
H
0.678682736478-1.510652355957
C-0.611168011571-0.040816289021
O-0.394870440670
1.303201201109
O-1.727513327852-0.556298092642
N
1.863759281644
0.033038180758
H
2.433071741324
0.008847406074
H
0.603583762900
1.414439955155
H
2.429316521329
0.010021320926
Options
GauMem=800
GauProc=8FirstOnly
First step ADDF only

Maeda, S.; Ohno, K. Chem. Lett. 2004, 33,
1372.; Maeda, S.; Ohno, K. Chem. Phys.
Lett. 2004, 398, 240.
10-0.001060332460
0.876536586179-0.881491724884-0.000086983445-0.000453147261
0.000681594335
0.000826841830-0.836395739236-0.000747436209
0.840645157227

Glycine Structure

11

12

2

EQ0.log

EQ0.log

PROFILE OF SHS-PATH 12

PROFILE OF SHS-PATH 22

Initial geometry (positive direction of mode 18)~~~~~~
GENERATION = 1
C
0.704548498814-0.741116935273-0.126321896943
H
0.555038468712-2.319903451614
1.976816850989
H
0.683372436623-1.683165680060-0.658432283213
C-0.578949905598-0.031018435066-0.001789867057
O-0.476352511700
1.359110857350
0.059088142868
O-1.685407307918-0.570704842910-0.013677878776
N
1.897754907198-0.073465480394-0.036504634400
H
2.702162293155-0.467051156845-0.520529374586
H
0.374633193317
1.696818067805-0.277459800688
H
2.129938359518
0.475165331871
0.791523766035
energy : -284.166219569826
This is a DC of NH2CHCOOH + H
Spin(**2) : 0.000000000000

Initial geometry (positive direction of mode 20)~~~~~~
C
2.085568310772-1.555922826638
H
2.025909252773-2.099175762747
H
3.042414055749-1.830570201822
C-1.500160566931
0.453787226214
O-0.689606968874
1.324781529310
O-2.445631699592-0.293171875315
N
2.220059624468-0.023207551327
H
2.693061080022
0.294589735679
H
1.316607948607
0.436826833210
H
2.863884763055
0.271329193564
energy : -284.203751045137
Spin(**2) : 0.000000000000# STEP 16 E(Harmonic) = 4.157448647915

Old DC was reachedThis is not listed in DC_list.log

Minimum search on the final SH was failed
13

Double-end-ADDF (d-ADDF)
Sphere-Expansion-Mode-0.227698296854
0.735223768196-0.704432982106
0.008435038057-0.089551191140
0.117873087918
0.120360416317-0.786965605219
0.225132094619
1.019953946341

This is CH2NH3 + CO2

SHS optimization failed becauseof a flat PES in the DC region

14

d-ADDF-TS/EQ (2PSHS/SCW)

Sphere-Contraction-Mode•d-ADDF-TS (2PSHS): Find a TS directly connecting two
points•d-ADDF-EQ (SCW): Find EQs between two points

Reactant or Product

ADDF passes through a TS region as a sphere expansion

A minimum on PES is a minimum on a hypersphere of a certain radius

Product or Reactant

Procedures finding a multi-step reaction path :1.

Exploration of intermediate EQ (EQs) by SCW

2.

Exploration of TS between neighboring EQs by 2PSHS

Quick search connecting a reactant and a product15

Ring ⇔ Cage (C10H10)

16

Ring ⇔ Cage (C10H10)

1. Optimization of Reactant and Product2. Exploration of EQs by SCW

Top view of the Cage
(Pentaprismane)

17

18

3

Ring ⇔ Cage (C10H10)

Input# SCW/B3LYP/6-31G

3. Exploration of TSs by successive SCWs

One Structure of Reactant or Product01
C
2.445962211022-5.078046387022
1.291611782326
C
2.313863867720-3.896710330530
2.316180910401
C
0.785812904391-3.622434879896
2.546590288136~~~~~~
Reactant
C
2.200136081408-6.095456140104
1.697094495718
C
2.071261719400-5.039832830138
2.723796656501
C
1.136478821867-4.070037184363
2.828097321401~~~~~~
Options
The Other Structure of Reactant or Product
GauMem=800
GauProc=8•

Structures of intermediates are listed in EQ_list.log•

EQs are listed in the reverse order in EQ_list.log

EQ3

EQ2

EQ1

EQ0

19

Ring ⇔ Cage (C10H10)

Input# 2PSHS/B3LYP/6-31G
01
C
2.445962211022
C
2.313863867720
C
0.785812904391~~~~~~
Reactant
C
2.071278856174
C
2.104963849879
C
0.764704277745~~~~~~
Options
GauMem=800
GauProc=8•20

4. IRC traces from each TS

One Structure of Reactant or Product-5.078046387022-3.896710330530-3.622434879896

1.291611782326
2.316180910401
2.546590288136-5.573029589518-4.121231952680-3.478859068434

1.732444377651
2.114651850808
2.679681180328

EQ0

EQ1

Neighboring Structure

TS is listed at the end of the .log file21

Important point 2

Important point 1• All intermediate are not necessarily explored by one SCW.

• Result depends on the choice of the center of the hyperspherecentercenter

1.
2.
3.

22

SCW between the fullerene and the ring.
SCWs between neighboring intermediates
2PSHS between neighboring EQs with no intermediates by SCW
23

Reverse choice of the center should also be tested.

24

4

Large-ADDF (l-ADDF)

Important point 3••

Only two sets of multi-step pathways can be found between thereactant and the product.
An intermediate not located linearly between reactant and product isdifficult to be searched （ex. A-B + C + D  B-C-D-A  A-B-C-D）
Reactant

Intermediate

To find many routes,
several intermediates
need to be assumed.

• Full-ADDF (f-ADDF): All ADDs are traced.
• l-ADDF: Only larger ADDs are traced.
Our postulate: The larger ADD
path has a TS at the shorter
distance.
Hammond postulate: A TS at
the shorter distance has themore stable product.
Bell-Evans-Polanyi principle:
The more stable product has a
TS at the shorter distance.

The larger ADD has the lower TS and the more stable product.

Product

Maeda, S.; Ohno, K. J. Phys. Chem. A 2007, 111, 4527.
25

Example of l-ADDF : HCOOH

26

Exploration from a HCOOH molecule# GRRM/B3LYP/6-31G

This will be copied to the Root Section of a Gaussian Input.

01
H-0.452596548000
0.034871834140-1.807771705495
O-0.452596548000
0.643013717456-1.034349518878
O-0.452596548000-1.276412777953
0.242743092264
C-0.452596548000-0.053534532177
0.147136623161
H-0.452596548000
0.652061759532
0.979506709948
Options
Optimized structure of HCOOH at B3LYP/6-31G
GauMem=100
GauProc=4
LADD=5 Trace five largest ADDs

f-ADDF: 80,926 gradients and 2,672 Hessian (38 h on 4 CPU-cores)
l-ADDF (N = 5): 7,617 gradients and 314 Hessian (3.7 h on 4 CPU-cores)

27

28

Input

l-ADDF Exploration of (H2O)8# GRRM/RHF/6-31G

Water Monomer

Lowest 20 minima at 0 K

Lowest 20 minima at 400 K

29

01
O
-0.000000000000
0.010468354749
0.754819362610 1
H0.000000000000
0.771670165130
1.319188188547 1
H-0.000000000000-0.747077150089
1.324086214974 1
O-0.000000000000
0.010468354749
0.754819362610 2
H0.000000000000
0.771670165130
1.319188188547 2
H-0.000000000000-0.747077150089
1.324086214974 2~~~~~~
O-0.000000000000
0.010468354749
0.754819362610 8
H0.000000000000
0.771670165130
1.319188188547 8
H-0.000000000000-0.747077150089
1.324086214974 8
Options
Start from 24 random structures
NRUN=24
LADD=5
Apply ADDF to 24 lowest EQs in the EQ list for
NLowest=24
Temperature=500.0 the Free energies at 0~500 K.
EQOnlyNo search for TS
UpDC=12DownDC=1230

5Definition of DCs•

Thermodynamic simulation

Upper limit of the bond length is assumed as• Superposition Approach
Eg. Wales, D. J.; et al. Adv. Chem. Phys. 2000, 115, 1.

0.1×N×(RA + RB)
1.

For an upward search by ADDF: UpDC = N (Default = 10)2.

For a downward search by IRC: DownDC = N (Default = 8)•

Recommended value

1.

Covalent bonds: UpDC = 10, DownDC = 8 (Default)

2.

Hydrogen bonds: UpDC = 12, DownDC = 12

3.

Dispersion forces : UpDC = 15, DownDC = 15

• Population at Temperature T

Pi (T ) =

( )
( )
∑
Zi T
X
Z T
k =1 k• Zi(T) : Partition function at T
(Harmonic approx.)• Denominator: Sum of partition
functions of all structures
Thermodynamic simulation can be made on the basis of the data obtained by GRRM
31

32

Thermal distributions of H+(H2O)8

Thermodynamic analysis

Located 168 local minima (MP2/6-311++G(3df,2p)//B3LYP/6-311+G**)

List of Equilibrium Structures

Rotation# Geometry of EQ 0, SYMMETRY = C1
O
-3.990749673348
0.812671330496
0.759596654549
H-4.623998138945
1.477984504789
0.989567688782
H-4.102553621535-0.016833813104
1.233805951983~~~~~~
Energy = -608.028031295279 Electronic energy
Spin(**2) = 0.000000000000
ZPVE= 0.215589990623
Vibration
Normal mode eigenvalues : nmode = 66
0.000010181 0.000032843 0.000087596 0.000141739 0.000189141
0.000289300 0.000410092 0.000619427 0.000760112 0.000883149
0.001138962 0.001194830 0.001290467 0.001517249 0.001737303~~~~~~•

Population (by the superposition approach) on the basis of harmonic vibrationalfrequencies reproduced three (IR) observations at different temperatures
1.
2.
3.

Partition functions can be computable at arbitrary temperature.

33

JACS 2000, 122, 1398: a mixture of chains and a Zundel single-ring, at ~170 K
Science 2004, 304, 1134: dominated by chains, at >> 170 K
Science 2005, 308, 1765: dominated by a Zundel single-ring, at <<170 K
Luo, Y.; Maeda, S.; Ohno, K. J. Phys. Chem. A 2007, 111, 10732.

34

Asymmetric synthesis by Ru-BINAP catalyst

Input# GRRM/PM6 EmpiricalDispersion=GD3

90 atoms outside the reaction center01
C-0.128588
1.878776
3.239214~~~~~~
H
1.027143-4.484840
3.683857
Frozen Atoms
P
1.563807-0.793743-0.751263
P-1.403001-1.238551-0.057742
Ru
0.135998-2.368924-1.238020
H0.578221-3.350460
0.067996Cl-0.285807-1.341860-3.384823
C
1.155131-4.429632-0.860073
O
0.985945-4.096942-2.110175
Options
GauMem = 400
GauProc = 8
Fix the reaction center at the form of TS
LADD=5
NLowest=5 Explore lower energy conformations

68 TSs leading to asymmetric selection
were searched automatically. Among theman R type was found to be the lowest inagreement with the experiment. Right
figures show five lowest TSs.
Maeda, S.; Ohno, K. J. Am. Chem. Soc.
2008, 130, 17228.
35

Reaction center

366

Application to crossings on
potential surfaces

H2CO (S0, S1, and T1 PESs)
New

Exploration of minima on seams of crossingsNew

Automated algorithm for MSX（Seam Model Function Method）
1. Approximate automated search for
MSX based on penalty function
2. Further optimization yields more
precise MSX structures.
A penalty function for automated MSX search:

) (

New

)

1
E State1 (Q ) − E State2 (Q )X (Q ) = E State1 (Q ) + E State2 (Q ) +2α(Mean potential energy

New

New

Initial search: (4e,3o)-CASSCF/6-31G and (2e,3o)-CASSCF/6-31G(85,000 (4e,3o)-CASSCF gradients and 80,000 (2e,3o)-CASSCF gradients)
Optimization 1: (12e,10o)-CASSCF/aug-cc-pVDZ
Optimization 2: (12e,10o)-CASPT2/aug-cc-pVDZ (shown in parentheses)

A penalty function

Maeda, S.; Ohno, K.; Morokuma, K. J. Phys. Chem. A 2009, 113, 1704.

37

38

Automated search of approximate
MSX structures

Optimization of MSX structures
Geometry optimization# MIN/UB3LYP/6-31G

New

2# GRRM/UB3LYP/6-31G

Unify the calculation levels

01
Pre-optimized structure with an input using MIN
C-0.341655575913-0.061830207285-0.629143745161
O0.051090096597-0.153079538593
0.678261012364
H-0.148111884773
0.634881653256
1.224415641392
H
0.344846324087-0.329784658376-1.435884114596
Options
ModelF(Seam)
Use of model function placing EQ near MSX
Second Input
UB3LYP/6-31G
03
END
MaxStepSize = 0.1
Stable=Opt

0 1 Singlet stateC-0.341655575913-0.061830207285-0.629143745161
O0.051090096597-0.153079538593
0.678261012364
H-0.148111884773
0.634881653256
1.224415641392
H
0.344846324087-0.329784658376-1.435884114596
Options
OptX(Seam) Within crossing regions
Second Input
UB3LYP/6-31G
0 3 Triplet state
END
MaxStepSize = 0.1 Stepsize should be smaller than the default value.
Stable=Opt•••

Without this option, the closed-shell singlet surface is explored with DFT.

GRRM options such as LADD and Nlowest are available.
EQOnly is automatically assumed.
Crude MSX is obtained as EQ, which needs to be refined by using OptX.

39

Geometrical Microiteration

An Application to the Si6-Prism

・ Divide reaction center from the others in the Giant system.
・ Always optimize for variables other than the reaction center

NO

Initial Structure

Initial Structure

Calculate force on all
atoms

Optimize positions of nonreaction center atoms

Determine displacementsof all atoms

Calculate force onreaction center atoms

Converged?

Determine displacementsof reaction center atoms

YES
Exit

Algorithm for General Optimization

NO

40

Real System

Model
SystemReactant
0 kJ/mol

Reactant

Converged?

0 kJ/mol

YES
Exit

Algorithm for Microiteration

New channels open in the real system
41

42

7

Application to complex reactions

Input# GRRM/ONIOM(B3LYP/6-31G*:UFF)

Reaction center including link atoms
01
Si-Si3-0.10182 -1.1881857-0.9102208-1.0465917
H
Si-Si3-0.09963 -1.2116703-0.4296025
1.3210086
H~~~~~~
C-C_R--0.09616 2.3999742-0.7349068
2.7547234
L
H-H_--0.08955 5
C-C_R--0.09221 2.3697321
2.7146470-0.8151322
L
H-H_--0.10701 6
External Atoms
C-C_R--0.09800 -4.5129210-2.4126440-3.1593396
L~~~~~~
Options
Atoms of every-step-optimization in Microiteration
GauMem=800
GauProc=8FirstOnly
GauInpB
1
2
1.0 3
1.0 4
1.0 7
1.0 Connectivity matrix for UFF calc.,
2
3
1.0 5
1.0 8
1.0
to be copied at the end of the~~~~~~
Gaussian input file .
180
END
43

•
•

Exploration of reaction precursors by l-ADDF (cf. water cluster)Hand-making of product geometries and connecting reaction
routes by d-ADDF

44

Luo, Y.; Maeda, S.; Ohno, K. Chem. Phys. Lett. 2009, 469, 57.

Connecting reaction route by
d-ADDF

Precursor search by l-ADDF# GRRM/M062X/6-31+G*# 2PSHS/M062X/6-31+G*

02
H
0.279704309-2.419095131
0.087118953 1
O
0.279704309-1.767230166
0.810586808 1
O0.279704309-3.669699011
2.040431695 1
C0.279704309-2.467398098
1.960665441 1
H0.279704309-1.772053247
2.812917361 1
O0.349449893-2.139473528
0.851411182 2H
0.349449893-1.380136528
1.447454182 2H
0.349449893-2.898810528
1.447454182 2O
0.349449893-2.139473528
0.851411182 3H
0.349449893-1.380136528
1.447454182 3
Options
DownDC=12
UpDC=12LADD=3
Similar to the search for water clusters
NLowest=3
NRUN=16
EQOnlyStarting structure obtained by l-ADDF (Left Fig. bellow)
02H-0.462142232965-1.064413231066
0.102466212778
O-1.206875810840-0.667073981701
0.642204190946
O-0.096500753587
1.280103269769
0.399425908015~~~~~~
Reactant
H-0.084528637-1.459791639-0.466700719
O-1.206875811-0.667073982
0.642204191 Hand-made structure from aO-0.096500754
1.280103270
0.399425908 structure obtained by l-ADDF(Right Fig. below)~~~~~~
Options
GauMem=800
Shift H-atoms
GauProc=8
DownDC=12
UpDC=12

HCOOH

H 2O
OH

45

46

Ru-BINAP-catalyzed hydrogenation
Mechanism by d-ADDF(Only QM-atoms are shown)
ONIOM(B3LYP/6-311++G**:UFF)

The whole cycle obtainedby d-ADDF
Ohno, K.; Maeda, S. J. Mol. Catal.
A: Chem. 2010, 324, 133.

47

8