Development of Divide-and-Conquer Quantum Chemical Codefor Biomolecules and Nano Materials
Masato Kobayashi 1, Patchreenart Saparpakorn 2,3, and Hiromi Nakai 2,4,51

Waseda Institute for Advanced Study (WIAS), Waseda University, Japan
School of Advanced Science and Engineering, Waseda University, Japan
3
Faculty of Science, Kasetsart University, Thailand
4Research Institute for Science and Engineering, Waseda University, Japan
5CREST, Japan Science and Technology Agency, Japan
2Introduction

Straightforward applications of standard quantum chemical
methods to biomolecules and nano materials are essentially
impossible because of the undesired computational timerequired for the calculations, scaling at least cubically with thesystem size N. Up to now, many types of accelerating
techniques for the quantum chemical calculations, which reducethe complexity to O(N), have been developed. The divide-andconquer (DC) quantum chemical method firstly proposed by
Yang and Lee [1] is a pioneering method that diminishes thecomputational cost for the Hartree–Fock (HF) and density
functional theory (DFT) calculations. The authors andcoworkers have extended the DC method to the post-HF
electron correlation calculations (e.g., the second-order Møller–Plesset (MP2) theory [2]) with the help of the energy density
analysis (EDA) [3]. In 2009, they released the DC quantum
chemistry program as a part of GAMESS package [4]. The DCbased O(N) methods are summarized in the review articles [5,6]In this Paper, we present recent developments in the DC
method and show their applications to biomolecules and nanomaterials.
2DC-MP2 Energy and Its Gradient

The DC-MP2 method [2,5] utilizes the subsystem molecularorbitals (MOs) obtained by solving the DC-HF equation:
F α C αp = ε αp S α C αp .
(1)Here, the superscript α denotes the spatially partitionedsubsystem. The subsystem MOs are constructed by the linear
combination of atomic orbitals corresponding to the disjoint
subsystem (referred to as the central region) α [S(α)] and its
neighboring region called the buffer region [B(α)]. Fα and Sα are
the subsystem Fock and overlap matrices, respectively, which
are the submatrices of the entire Fock and overlap matrices Fand S.
The total DC-MP2 correlation energy is estimated bysumming up the correlation energies corresponding to individualsubsystems as follows

ΔEDC-MP2 =subsystem∑αα.
ΔEMP2(2)Here, the correlation energy of subsystem α is evaluated withsubsystem MOs with the assistance of EDA as follows:occ(α ) vir(α )αΔEMP2= ∑ ∑i, ja ,b∑α Cμα μ jα aα bα ( 2tα − tα ) , (3)μ∈S ( )

*iij , abij ,bawhere {iα, jα, …} and {aα, bα, …} represent the occupied andα
virtual MOs of subsystem α, respectively. tij ,ab is the effective
two-electron excitation amplitude for subsystem α, which isgiven in the MP2 case bytijα, ab = − aα bα iα jα(ε + ε − ε − ε ) .
ααααabij(4)Recently, we have derived the gradient of the DC-MP2energy. The conventional MP2 energy gradient is expressed as

∂H μν
∂S
∂ΔEMP2(2)= ∑ Dμν+ ∑ Wμν(2) μν∂Q∂Q∂Qμνμνcore+ ∑ Γ (2)μνλσ
μνλσ(2)∂ μλ νσ.
∂Q(5)Here, D(2) and W are standard and energy-weighted perturbed
one-body density matrices, respectively, and Γ(2) is theperturbed two-body density matrix. In the present DC-MP2
gradient method, the one-body density matrices, D(2) and W(2),are evaluated in the same way as the unperturbed density matrix:(2)(2),αDμν= ∑ Pμνα Dμν,
(6)αwhere P is the partition matrix that is introduced to avoid
multiple counting of electronic contribution in the buffer region.
When constructing subsystem perturbed density matrices, weapplied the Z-vector method [7] to solve the coupled perturbed
HF (CPHF) equations of subsystems. The two-body matrix, Γ(2),can be divided into the separable (S) and non-separable (NS)terms. The S term is expressed as the product of one-bodymatrices, which can be evaluated in the DC manner as Eq. (6).
In the present method, we evaluated the NS term as follows:
NS,DCNS,αΓ μνλσ
= ∑ Pμνα Γ μνλσ,(7)αNS,αΓ μνλσ
= ∑∑ tijα, abCμαi Cναa Cλα j Cσαb .

(8)i , j a ,b-2315.78-2315.80
Energy [hartree]

1Conv. MP2DC-MP2 (n bcorr = 3)-2315.82

DC-MP2 (n bcorr = 4)
DC-MP2 (n bcorr = 5)-2315.84-2315.86-2315.88
246810Cycle

Fig. 1. Energy convergence process in the geometry
optimization of C60H62 in DC and conventional MP2
calculations with the 6-31G** basis set.

Table 1. Errors of optimized geometrical parameters of C60H62 at theDC-MP2/6-31G** level. The energy errors are given at the bottom.

nbcorr = 3

nbcorr = 4

nbcorr = 5

Bond length

MAE

0.202

0.022

0.015

[pm]

MaxE

0.522

0.071

0.059

Bond angle

MAE

1.003

0.062

0.020

[°]

MaxE

2.506

0.186

0.043

11.142

3.354

2.023

Energy [mEh]

We assessed the performance of the present DC-MP2 energy
gradient in the geometry optimization of an all-trans polyene
chain, C60H62. A C2H2 (or C2H3 for the edges) unit was adopted
as the central region and the adjacent left and right nb units wereadopted as the corresponding buffer region. Figure 1 shows theenergy convergence process in DC and conventional MP2calculations. Here, the buffer size for the preceding DC-HF
calculation was fixed at nbHF = 6 , and that for the MP2
calculation was varied between n = 3 − 5 . As the buffer size
increases in the DC-MP2 method, the optimization process
becomes closer to the conventional MP2 calculation and yields
lower converged energy. Table 1 gives the mean absolute error(MAE) and the maximum error (MaxE) of the optimized
geometrical parameters and the energy error of C60H62 from theconventional results. The geometrical errors also decrease as thebuffer size increases. Note that the computational time for one
optimization cycle by DC-MP2 method with nbcorr = 5 wasaround 15 times as fast as that by conventional MP2.
corrb

3 Application to the Interaction between HIV-1
Reverse Transcriptase and MK-4965 Inhibitor
Reverse transcriptase (RT) of HIV-1 is an essential enzyme
converting the single-stranded viral RNA genome into doublestranded proviral DNA prior to its integration into the hostgenomic DNA. MK-4965 (Fig. 2) is a novel non-nucleosidereverse transcriptase inhibitor (NNRTI) containing diaryl ether
and indazole moieties, which reveals high levels of potency
against wild-type (WT) and some important HIV-1 RT mutants,that is, 0.23 and 0.39 nM in IC50 for WT and Y181C HIV-1 RTs,
respectively.
Recently, the authors applied the spin-component scaled(SCS) version of DC-MP2 to elucidate the binding of MK-4965
in the WT and Y181C HIV-1 RT binding pockets [8]. The DCMP2 method not only reduces the computational cost but alsoenables to obtain the partial energy distributions by EDA, which
was applied to estimate the partial interaction energies between
HIV-1 RTs and MK-4965 subsystems. Figure 3 illustrates theinteraction energy components of each amino acid residue in
WT and Y181C HIV-1 RTs. It clearly shows the vitalinteraction between Lys102 and SL3 in both WT and Y181C.
The second important interaction is seen between Lys103 andSL3. The more negative interaction energy about 3 kcal/mol isfound in Y181C. The aforementioned distribution of interactionenergies in each amino acid was clarified in detail by using theanalysis of the geometry and atomic charges.

Fig. 3. Interaction energies (in kcal/mol) between SLs andeach amino acid obtained by the DC-SCS-MP2/6-31G**
calculation.
4After the first release of our DC program in GAMESS, wehave continuously developed the DC method and its code for
effectively treating biomolecules and nano materials. Especiallyfor MP2 calculations, we have proposed the two-level
parallelization method [9] that enables to treat huge systems
with high accuracy by using massively parallel computers.
Acknowledgements
Some of the present calculations were performed at theResearch Center for Computational Science (RCCS), Okazaki
Research Facilities, National Institutes of Natural Sciences(NINS). This study was supported in part by Grants-in-Aid forYoung Scientists (B) “KAKENHI 22750016” and for
Challenging Exploratory Research “KAKENHI 22655008”
from the Ministry of Education, Culture, Sports, Science andTechnology (MEXT), Japan; by the Strategic Program for
Innovative Research (SPIRE), MEXT; and by theComputational Materials Science Initiative (CMSI), Japan.
References
[1]

[2]

[3]
[4]
[5]

[6]

[7]

[8]

[9]

Fig. 2. Structure and subsystems of MK-4965.

Concluding Remarks

W. Yang and T.-S. Lee, A density-matrix divide-and-conquer
approach for electronic structure calculations of large molecules, J.
Chem. Phys., 103, 5674, 1995.
M. Kobayashi, Y. Imamura, and H. Nakai, Alternative linearscaling methodology for the second-order Møller–Plesset
perturbation calculation based on the divide-and-conquer method,
J. Chem. Phys., 127, 074103, 2007.
H. Nakai, Energy density analysis with Kohn–Sham orbitals,
Chem. Phys. Lett., 363, 73, 2002.
http://www.msg.chem.iastate.edu/gamess/
M. Kobayashi and H. Nakai, Divide-and-conquer approaches toquantum chemistry: Theory and implementation, Linear-scaling
techniques in computational chemistry and physics: Methods andapplications, Springer, Dordrecht, 97–127, 2011.
M. Kobayashi and H. Nakai, How does it become possible to treat
delocalized and/or open-shell systems in fragmentation-based
linear-scaling electronic structure calculations? The case of thedivide-and-conquer method, Phys. Chem. Chem. Phys., in press(DOI: 10.1039/c2cp40153c).
N.C. Handy and H.F. Schaefer, On the evaluation of analytic
energy derivatives for correlated wave functions, J. Chem. Phys.,81, 5031, 1984.
P. Saparpakorn, M. Kobayashi, S. Hannongbua, and H. Nakai,Divide-and-conquer based quantum chemical study for interaction
between HIV-1 reverse transcriptase and MK-4965 inhibitor, Int. J.
Quantum Chem., in press.
M. Katouda, M. Kobayshi, H. Nakai, and S. Nagase, Two-level
hierarchical parallelization of second-order Møller–Plesset
perturbation calculations in divide-and-conquer method, J. Comput.
Chem., 32, 2756, 2011.