4P085

Elongation 局所構造最適化による全系の効率的安定構造探索(九大院・総理工1, SCNU2, JST-CREST3) Liu Kai1, 折本 裕一1, Gu Feng Long2, 3, 青木 百合子1, 3

Efficient search for optimized geometry of whole system by elongation localoptimization method(Kyushu Univ.1, SCNU2, JST-CREST3) Kai Liu1，Yuuichi Orimoto1，Feng Long Gu2，Yuriko Aoki1,3

I. Introduction
One of the most successful applications of molecular quantum mechanics is the reproduction and prediction ofmolecular conformation. However, it’s still a grand challenge to perform the geometry optimization on large systemin practices. Due to a large number of atoms, the standard ab initio calculations are beyond the computational reach.
Currently, some fragment-based methods have been successfully proposed and developed to calculate large system
within a reasonable cost, for example: divide and conquer (DC)1, elongation method (ELG)2, fragment molecularorbitals (FMO)3 and systematic fragmentation method (SFM)4. Despite the complicate interactions of large system,
it’s very difficult to locate the equilibrium structure. Here we present the implementation of geometry optimization
based on the elongation method, which is famous for its high accuracy and efficiency.

II. Theoretical approach
A. Elongation method
ELG procedure generally is analogous to experimental polymer
chain synthesis as shown in Fig.1. First, a suitable size of initialmonomers (starting cluster) is chosen to initiate the ELG procedure.
The canonical molecular orbitals (CMOs) of the starting cluster,generated by conventional SCF calculation, are transformed to anorthogonal atomic basis (OAO). In the following step, the OAO-based
density matrix (DOAO) is to partition into frozen region (A) and activeregion (B). The frozen region is assumed to be far away from thechain propagation point, while the active region is consisted of theremaining part of the starting cluster. After the separate
diagonalization of the subspace DOAO (A) and DOAO (B), a set of
regional localized molecular orbitals (RLMOs) for both A and B
regions are obtained. Then, one attacking unit (C) is added to theFig.1 The flowchart of elongation method
chain propagation point for the next ELG step. Because of thenegligible interactions between A and C, only B and C will be included in the ELG SCF calculations. After the ELG
HF-SCF converges, the CMOs of B and C regions will be localized again to form a new frozen region (RLMOs B’)and a new active region (RLMOs C’). Then a new attacking unit (D) is added to repeat the above procedures until
the desired length is reached. The important feature of the ELG method is that the equations during SCF are solved
only for small subunits instead of the whole system, and the calculations of two-electron repulsion integrals (ERIs)
between A and M can be partly omitted by cutoff technique.

B. Elongation geometry optimization
In the framework of HF calculation, the first derivative (gradient) of total energy (E) with respect to the nuclear
coordinate X A at atomic orbitals basis, can be written as:𝝏𝑯𝒄𝒐𝒓𝒆𝝏𝑺𝝁𝝂𝝏(𝝁𝝂||𝝈𝝀) 𝝏𝑽𝑵𝑵
𝝏𝑬
𝟏𝝁𝝂= � 𝑫𝝂𝝁+ � 𝑫𝝂𝝁 𝑫𝝀𝝈+− � 𝑸𝝂𝝁𝝏𝑿𝑨
𝝏𝑿𝑨
𝝏𝑿𝑨
𝝏𝑿𝑨
𝝏𝑿𝑨
𝟐𝝁𝝂𝝁𝝂𝝀𝝈𝝁𝝂

where (µν||σλ) denotes two-electron integral, D is the density matrix, Hcore and S correspond to the core Hamiltonianand overlap matrices of system, respectively. The nuclear-nuclear repulsion is defined by V NN . The energy-weighted

N/2density matrix Q of conventional method is defined asQ νμ = ∑i ni εi Cμi Cνi , where n is the matrix of occupancy
number (density matrix in MO representation). In the ELG method, the equations during SCF are solved only for
active subspaces defined by B and M fragments. Therefore, the non-diagonal εLMOis employed in the ELG -OPTijAOLMOmethod. As the coefficients of the whole system CABM consist of CA (coefficients of frozen part, localized by theAOMOELG localization procedure) and CBM(coefficients of active region and attacking unit, transformed from CBM
afterthe ELG SCF calculation). Therefore, the energy-weighted density matrix Q of the ELG-HF-OPT method in atomic
N/2
basis can be rewritten asQ νμ = ∑ij ni εij Cμi Cνj . Because of the boundary effects between A and BM arisen by thetails after the ELG localization procedure, the coordinates of the one unit of the BM region, which is the closest tothe frozen region, is fixed in the gradient calculation to reduce these effects.

III. Results and discussion
A. Non-bonding model system: (HF) n=48
The linear poly-hydrogen fluoride (poly-HF) molecules are optimized by both elongation and conventional
restricted Hartree-Fork method (RHF) with different basis sets. The energy differences ∆E (∆E=Eelongtion-Econventional)of STO-3G, 6-31G and 6-31G(d,p) basis sets are -4.49×10-7, 9.09×10-9 and -3.26×10-5 Hartree/atom, respectively.
The negative values means the ELG-OPT locates an even lower ground state than conventional results. It indicatesthat for a flat energy potential surface of system, like linear poly-HF molecules, ELG-OPT may produce a morepromising candidate for the most stable geometry.

F

RF...HH

F

RF...HH

6.00E-02
4.00E-02
2.00E-02

STO-3G
6-31G
6-31G(d,p)

CONV provides ametastable structure.

0.00E+00

0-2.00E-02

0-4,800.2

Total energy

8.00E-02

Final step of optimization iteration(au)-4,800.0

RF...H bond length difference(Å)
1.00E-01

2040606-31G(d,p)-4,800.4-4,800.680100

ELG-OPT
CONV-OPTEelg lower than Econv
Telg shorter than Tconv-4,800.8-4,801.0

5 10 15 20 25 30 35 40 45 50-4,801.2

E: Energy T: TimeFig. 2 The differences of bond length and optimization iteration of poly-HF systems. ∆R=Relongation

B. Bonding system: extended ployalanine.
As a model of large biosystem, 20 units ofextended alanines, is optimized by both the ELGHF-OPT and conventional method using 6-31Gbasis set. The RMSD between the ELG andconventional optimized structure (excludinghydrogen) is 0.18. The bond length is also
compared in details and shown in Fig.3. The
difference of total energy between ELG-OPT andconventional result is 2.37×10-7 Hartree/atom.
These small differences in structure and totalenergy show the well reproduction of elongationgeometry optimization method.

-Rconventional

The difference of bond length(Å)
8.00E-04
6.00E-04

6-31G

4.00E-04
2.00E-04
5.00E-19-2.00E-04 020406080

100-4.00E-04

Fig. 3 The difference of bond length between ELG andconventional optimized structure.

Compared to conventional results, it indicates that
the ELG-OPT can well reproduce the calculations and may locate a more stable structure than conventional one.
1. W. T. Yang, Phys. Rev. Lett. 1991, 66, 1438.
2. A. Imamura, Y. Aoki, K. Maekawa, J. Chem. Phys. 1991, 95, 5419.
3. K. Kitaura, E. Ikeo, T. Asada, T. Nakano, M. Uebayasi, Chem. Phys. Lett. 1999, 313, 701.
4. V. Deev, M. A. Collins, J. Chem. Phys. 2005, 122, 154102.