A parallel eigensolver using contour
integration for generalized eigenvalue
problems in molecular simulation
1

Tetsuya SAKURAI 1) , Hiroto TADANO 1) , Tsutomu IKEGAMI 2) ,and Umpei NAGASHIMA 2)1) Department of Computer Science, University of Tsukuba,
1-1-1 Tennodai, Tsukuba, Ibaraki 305-8573, Japan
2) National Institute of Advanced Industrial Science and Technology,
1-1-1 Umezono, Tsukuba, Ibaraki 305-8568, Japan

Abstract
In this paper, we consider a parallel method for computing interior eigenvalues and corresponding eigenvectors of generalized eigenvalue problems whicharise from molecular orbital computation of biochemistry applications. Matrices in such applications are sparse but often have relatively large numberof nonzero elements, and we may require some eigenpairs in a speciﬁc partof the spectrum. We use a contour integration to construct a desired subspace. Some properties of the subspace obtained by numerical integration
are discussed, and then a parallel implementation is presented. We report
the numerical aspects and parallel performance of our method with matrices
derived from molecular orbital computation.
Keywords: Interior eigenvalue problems, contour integration, Rayleigh-Ritz
procedure, biochemistry application.
MSC(2000): 65F15, 65E05, 65Y05

1

Introduction

In this paper, we consider a parallel method for computing some eigenpairs (λ, x)satisfying Ax = λBx in a speciﬁc part of the spectrum, where A, B ∈ Rn×n are symmetric and B is positive deﬁnite, which arise from molecular orbital computation.
Molecular orbital computation has been performed to study a reaction mechanismof enzymes, an electronic structure of the photosynthetic systems, etc.
The fragment molecular orbital (FMO) method [4] enables one to calculate thetotal energy of molecule without performing the calculation of whole molecule. In[1], a full electron calculation for a large molecule was performed by using the FMO
method on a massive cluster computer. The FMO-MO method [3, 12], which works
as a post-process of the FMO method, gives a good approximation for the canonicalMOs without SCF iterations. A huge Fock matrix is generated, and we need limited
number of eigenpairs corresponding to frontier orbitals. Therefore a large-scaleinterior eigenvalue problem appears in the FMO-MO method.
1

Corresponding Author Email: sakurai@cs.tsukuba.ac.jp

With growing parallel eﬃciency in the Fock matrix construction, diagonalization
process becomes bottleneck in MO calculations. Since the Fock matrix has relatively
large number of nonzero elements due to the base function of middle-range interaction of molecules, a sparse factorization of the shifted matrix for spectral transformation may not be feasible. An alternative to factorization is oﬀered by iterativesolvers. Unfortunately, it is often unattractive for parallel computation, because
iterative process is performed successively. Therefore shift-and-invert approach is
not eﬀective to treat interior eigenvalue problems in such applications.
Our method for ﬁnding eigenpairs in a given physical domain is based on contour integral presented in [6]. One of major advantages of our eigensolver is thatit does not require inner/outer loops; an inner loop to construct an approximate
subspace, and an outer loop to update approximate eigenvectors. Moreover, thevalues which are used in a numerical integration can be evaluated independently oneach integration node, it provides a variety of parallel programming models [7]. We
recently proposed a Rayleigh-Ritz type method [8] in order to improve numericalstability. The block method [2] performs well if there are multiple eigenvalues in theinterested spectrum.
The computation at each integral node involves linear system solutions where
the coeﬃcient matrices are derived from A and B. In [5], we found that a Krylov
subspace iterative method in conjunction with a preconditioning using a complete
factorization for an approximated coeﬃcient matrix is eﬀective to solve such linear
systems.
In the next section, we explain our eigensolver using a contour integration. Then
we present some numerical properties of the subspace obtained by a numerical approximation of the contour integration. In Section 3, we show an implementation of
our method. In Section 4, some numerical examples are given to verify the numericalproperties and parallel performance of the methods with the matrices derived frommolecular orbital computation.

2

An eigensolver using a contour integration

In this section, we describe an eigensolver using a contour integration presented in[8]. This method ﬁnds eigenpairs inside a given circle. Let (λi , xi ), 1 ≤ i ≤ n be
eigenpairs of the matrix pencil (A, B). Suppose that m eigenvalues λ1 , . . . , λm arelocated inside the circle Γ with the center γ ∈ R and the radius ρ > 0 in the complex
plane. For a nonzero vector v ∈ Rn , let∫1sk =z k (zB − A)−1 Bv dz, k = 0, 1, . . . , m − 1,
(1)
2πi Γ
and let S = [s0 , . . . , sm−1 ] ∈ Rn×m .
Suppose that v is expanded by the eigenvectors {x1 , . . . , xn } as
v=

n
∑
i=1

Then we have the following theorem([8]).
2αi xi .

(2)

Theorem 1 If λ1 , . . . , λm are distinct and αj ̸= 0 for 1 ≤ j ≤ m thenspan{s0 , . . . , sm−1 } = span{x1 , . . . , xm }.
Let Q = [q1 , . . . , qm ] ∈ Rn×m be an orthonormal matrix derived from S =[s0 , . . . , sm−1 ]. Theorem 1 implies that the eigenpairs (λi , xi ), 1 ≤ i ≤ m can beextracted by using the Rayleigh-Ritz procedure with the projected matrix pencil(QT AQ, QT BQ).
We approximate the contour integration (1) via the N -point trapezoidal rule:N −1 (

1 ∑
ŝk =
N j=0

ωj − γρ)k+1(ωj B − A)−1 Bv,

k = 0, 1, . . . , m − 1

where ωj = γ + ρ exp(2πi(j + 1/2)/N ) and N is a positive integer.
In this computation, we solve the systems of linear equations(ωj B − A)yj = Bv,

j = 0, 1, . . . , N − 1(3)for y0 , . . . , yN −1 ∈ Rn . Note that the solutions yN/2 , . . . , yN −1 are obtained fromthe relation yj = ȳN −j−1 . Thus we only need to solve N/2 systems, and(()k+1 )∑
ωj − γ
1
ŝk =
2 Re
yj ,
N j=0
ρ
N/2−1

k = 0, 1, . . . .

(4)

Let θj = exp(2πi(j + 1/2)/N ), then the following relations are hold.
Lemma 2 Let η be a real number with |η| ≠ 1. For an integer k (1 ≤ k < N ), thefollowing holds.
N −1
k+1
1 ∑ θjηk=.
(5)
N j=0 θj − η
1 + ηN
Proof. If |η| < 1, we have
N −1
N −1
k+1
1 ∑ θj1 ∑ θjk=
N j=0 θj − η
N j=0 1 − θηj
∞
∑

=

p=0
∞
∑

=

η

p 1

N
−1
∑

N j=0

θk−p

(−1)q η N q+k .

(6)

q=0

The last step follows from the fact that
1 ∑ p
θ =
N j=0 jN −1{(−1)q if p = qN for q ∈ Z.
0otherwise

3(7)

Similarly, for the case |η| > 1, we have
) k+1
N −1
N −1 (k+1
1 ∑ θj1 ∑ −1 θj=θ
N j=0 θj − η
N j=0η
1 − ηj
)
)( N
−1
∞ (
∑
∑
−1
1
θp+k+1
=
p+1
η
N j=0
p=0
∞
∑

=

(−1)q−1 η −N q+k .

(8)

q=1

If follows from (6) and (8), we have the equation (5).

¤Theorem 3 Let ηi = (λi − γ)/ρ, 1 ≤ i ≤ n, then
ŝk =

n
∑
αi
i=1

ρ

·

ηik
xi .
1 + ηiN

Proof. Since(ωj B − A)−1 Bxi =

(9)

1
xi ,
ωj − λiit follows thatN −1 (

ŝk

1 ∑
=
N j=0

ωj − γ
ρ

1 ∑
=
N j=0

ωj − γρ1
N j=0

ωj − γρN −1 (

=

=

N
−1 (
∑

n
∑
αi

ρ
i=1

(

)k+1
(ωj B − A)−1 Bv)k+1 ∑nαi (ωj B − A)−1 Bxi

i=1)k+1 ∑nαixi
ωj − λi
i=1
)

N −1
k+1
1 ∑ θjN j=0 θj − ηixi .

Therefore, from Lemma 2, we have
ŝk =

n
∑
αi
i=1

ρ

·

ηik
xi .
1 + ηiN¤

This proves the theorem.
For vectors a, b ∈ Rn , we deﬁne the B inner product 〈·, ·〉B by〈a, b〉B = aT Bb.

From Theorem 3, we have the following estimation for eigencomponents that areincluded in ŝk .
4Input A, B ∈ Rn×n , V ∈ Rn×L , N , M , δ
Output (λ̂i , x̂i ), 1 ≤ i ≤ m̂
1. Solve (ωj B − A)Yj = BV for Yj , j = 0, . . . , N/2 − 1()∑
2. Compute Ŝk = N/2−1
2Re θjk+1 Yj , k = 0, . . . , M − 1
j=0
3. Perform singular value decomposition U T ΣW = [Ŝ0 , . . . , ŜM −1 ], and ﬁndm̂ such that |σj |/|σ1 | ≥ δ for 1 ≤ j ≤ m̂6. Construct an orthonormal basis Q̃ from Ŝ(:, 1 : m̂)7. Form Ã = Q̃T AQ̃ and B̃ = Q̃T B Q̂
8. Compute the eigenpairs (λ̂i , ŵl ), 1 ≤ l ≤ m̂ of the projected pencil (Ã, B̃)
9. Set x̂i = Q̃ŵi , 1 ≤ l ≤ m̂
Figure 1: Algorithm of the block method using the contour integration
Theorem 4 Suppose that the eivenvectors x1 , . . . , xn are normalized to be B-orthogonal,i.e., 〈xi , xj 〉B = δij . Then
{ αi|〈xi , ŝk 〉B | =η k + O(|ηi |N +k )
ρ i
O(|ηi |−N +k )(1 ≤ i ≤ m)
(m + 1 ≤ i ≤ n).

From this theorem, we can see that the eigencomponent in ŝk with respect to λi
decays exponentially with the ratio |ηi | = |λi − γ|/ρ, if λi is located outside Γ.
A block variant of the method is proposed in [2], which enable us to obtaindegenerate eigenvalues. In this method, a matrix V = [v1 , . . . , vL ] ∈ Rn×L is usedinstead of v in (3), where v1 , . . . , vL are linearly independent, and positive integer
L is a block size. Then the numerical integration (4) is represented as(()k+1 )
N/2−1
1 ∑ωj − γ
Ŝk =
2 Re
Yj
N j=0
ρ

k = 0, 1, . . . , M − 1,with systems of linear equations(ωj B − A)Yj = BV,

j = 0, 1, . . . , N/2 − 1,

(10)where M is a positive integer chosen so that M ≥ m/L.
For the choice of the size of the subspace, we use the singular value of Ŝ =[Ŝ0 , . . . , ŜM −1 ]. Let U ΣW = Ŝ be a singular value decomposition of Ŝ, and let
Σ = diag(σ1 , . . . , σM ). Let m̂ be an integer such that σj /σ1 ≥ δ for 1 ≤ j ≤ m̂ andσj /σ1 < δ for j > m̂. Then the orthonormal matrix Q̃ ∈ Rn×m̂ is constructed fromŜ(:, 1 : m̂). The algorithm is shown in Figure 1. Note that we can use U (:, 1 : m̂)instead of Q̃ to avoid the orthogonalize process of Ŝ(:, 1 : m̂).

5

3

An implementation

When matrices A and B are large, the computational costs for solving systems oflinear equations (10) are dominant in our method. In our application, matrices
have relatively large number of nonzero elements, which could make sparse direct
methods impractical due to large number of ﬁll-ins. Instead, we apply a Krylov
subspace iterative method with a complete factorization preconditioner presented in[5].
In our preconditioning method, a complete factorization of the approximate matrix C̃j = (c̃ij ) of the coeﬃcient matrix Cj := ωj B − A is performed. The approximate matrix C̃j is obtained from drop-thresholding of the original coeﬃcient matrix
Cj . Due to less nonzero entries in C̃j than Cj , we expect fewer nonzeroes in thepreconditioner than the matrix factor obtained from complete factorization of Cj .
The drop-thresholding is deﬁned as follows:{cij(|cij | > ε)c̃ij =,
0(|cij | ≤ ε)where ε is a small positive number.
Since systems (10) can be solved independently for j, we solve N/2 systems
Cj Yj = BV,

j = 0, 1, . . . , N/2 − 1(11)

on each computing node of clusters. Before starting the iterative process of the
Krylov subspace method, the approximate matrix C̃j is factorized, and then forward/backward substitutions and matrix-vector multiplies are performed in the iterative process.
We can easily extend the method for the case that several number of circularregions are given. Suppose that Nc circles Γ(0) , . . . , Γ(NC −1) with center and radius(γ (l) , ρ(l) ), 0 ≤ l ≤ NC − 1 are given. Then we solve NC × (N/2) linear systems
(l)

(l)

(ωj B − A)Yj

= BV,

j = 0, . . . , N/2 − 1,(l)l = 0, . . . , NC − 1,

where ωj = γ (l) + ρ(l) exp(2πi(j + 1/2)/N ), j = 0, . . . , N/2 − 1 are equidistributed
points on the l-th circle Γ(l) . We solve these linear systems simultaneously on distributed computing nodes.
In the beginning, we broadcast the sparse matrix data of A and B to all computing nodes. Then we start to solve systems of linear equations on the computing
nodes simultaneously. After that, we construct a subspace for each circle using theresults of the linear solver, and then perform the Rayleigh-Ritz procedure to extract
eigenpairs in each circle. The computational costs to solve systems of linear equations are dominant in whole computations, thus we can expect high performance inparallel computations. The algorithm is shown in Figure 2.

6

1. Broadcast sparse matrix data of A and B
2. do l = 0, 1, . . . , NC − 1
(l)

(l)

(l)

3.

Solve (ωj B − A)Yj

4.

AllReduce 2Re(θj )k+1 Yj , j = 0, . . . , N/2 − 1, and compute()∑N/2−1(l)(l)(l)
Ŝk = j=0
2Re (θj )k+1 Yj , k = 0, . . . , M − 1

5.

Compute singular values σ1 , . . . , σM of Ŝ (l) = [Ŝ0 , . . . , ŜM −1 ],(l)= BV for Yj , j = 0, . . . , N/2 − 1
(l)

(l)

(l)

(l)

(l)

(l)

(l)

and set m̂(l) such that |σj |/|σ1 | ≥ δ for 1 ≤ j ≤ m̂(l)6.

Construct the orthnormal basis Q̃(l) from Ŝ (l) (:, 1 : m̂(l) )7.

Form Ã(l) = (Q̃(l) )T AQ̃(l) and B̃ (l) = (Q̃(l) )T BQ(l)8.

Compute the eigenpairs (λ̂i , ŵi ), 1 ≤ i ≤ m̂(l) of the projected

(l)

(l)

pencil (Ã(l) , B̃ (l) )
9.

(l)

(l)

Set x̂i = Q̃(l) ŵi , 1 ≤ i ≤ m̂(l)10. end do

Figure 2: A parallel implementation with circles Γ0 , . . . , ΓNC −1

4

Numerical examples

In this section, we show some numerical examples with the matrices that are derived
from the FMO-MO method[10].
Example 1 The test matrices A and B are derived from computation of the molecular orbitals of an eight-base-pair model DNA. The size of matrices is n = 1, 980,
and the number of nonzero elements is 728, 080.
The CPU was Core2Duo(2.2GHz) with 2GB of memory. Computation was performed in MATLAB 7.5 with double-precision arithmetic. The systems of linear
equations were solved by the sparse direct solver UMFPACK in MATLAB. The elements of V were distributed randomly on the interval [0, 1] by a random number
generator. The exact eigenpairs (λi , xi ) 1 ≤ i ≤ n were evaluated by the MATLAB
command eig.
In Figure 3, we show |〈xi , ŝ0 〉B | with respect to eigenvalues λi ∈ [−0.25, −0.16].
The parameters were γ = −0.2, ρ = 0.01, N = 32, L = 12, M = 12, and δ =10−12 . In this case, the interval [−0.21, −0.19] is inside and on Γ. We can see thatthe eigencomponents associated with eigenvalues that are located outside Γ decay
exponentially to zero in ŝ0 . For vectors ŝk , k ≥ 1, we can obtain similar results. In
Figure 4, we show |〈x̂i , Q̃〉B |. We can see that Q̃ includes only the eigencomponents
around the interval [−0.21, −0.19].
In Table 1, the error and residual of approximate eigenvalues in Γ are shown. 16
eigenvalues are located in the interval [−0.21, −0.19], and all the eigenvalues in thisinterval were obtained with suﬃcient accuracy.
7

Figure 3: |〈xi , ŝ0 〉B | (γ = −0.2, ρ = 0.01, N = 32)

Figure 4: |〈x̃i , Q̃〉B | (γ = −0.2, ρ = 0.01, N = 32)
Example 2 The test matrices were obtained from an epidermal growth factor receptor (EGFR) protein which is known as a target molecule for anticancer agent.
The Fock matrix was constructed by the method described in [10]. The dimension of
the Fock matrix was 96, 234, and the number of non-zero elements was 457 million.
The calculations were performed on P32 subsystem of the AIST Super Cluster(ASC) in National Institute of Advanced Industrial Science and Technology Agency(AIST). Each computing node is Opteron Dual processor (2.0GHz) with 6 GB ofmemory, and is interconnected by both Myrinet and gigabit ethernet. Intel C and
Fortran compiler 9.1 were used to compile the codes with Intel Math Kernel Library
10.0.
The number of integration points was N = 32, and the block size was L = 16.
We put 8 circles around the highest occupied molecular orbital (HOMO) and thelowest unoccupied molecular orbital (LUMO) using the approximate eigenvalues
obtained from the results of the FMO method. Then the total number of systems
to solve was 8 × (32/2) = 128 with 16 right-hand side vectors. Each system of linear
equations was solved on single node of dual processor, thus 16 nodes (32 PUs) wereused for one circle. Therefore total 256 PUs were employed.

8

Table 1: Errors and residuals in Example 1.
i
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16

|λ̂i − λi |λ̂i−0.192423480937650−0.193053272393176−0.193867518039814−0.194749418147652−0.195768218806505−0.196048479728466−0.196721896535061−0.197042777317239−0.200591957877910−0.200835717839500−0.203655629362533−0.203800310807316−0.203932047603002−0.204805988179016−0.204824251443765−0.205054724756940−15

2.9 × 10
3.7 × 10−15
1.1 × 10−15
2.0 × 10−14
2.5 × 10−15
7.2 × 10−16
1.2 × 10−15
4.7 × 10−16
1.2 × 10−15
5.8 × 10−16
2.5 × 10−16
1.4 × 10−16
6.1 × 10−16
3.4 × 10−15
2.6 × 10−15
1.3 × 10−14∥Ax̂i − λ̂i B x̂i ∥2
9.4 × 10−15
1.6 × 10−14
1.0 × 10−14
2.0 × 10−14
2.5 × 10−14
2.8 × 10−14
4.0 × 10−14
5.5 × 10−14
5.5 × 10−13
1.9 × 10−13
5.6 × 10−14
4.4 × 10−14
4.4 × 10−14
1.6 × 10−14
3.4 × 10−14
2.8 × 10−14

The preconditioned COCG method[11] was used for iterative linear solver. The
stopping criterion for the relative residual norm was 10−10 . The preconditioner wasconstructed by applying a complete factorization for an approximate coeﬃcient matrix which was obtained from drop-thresholding of the original coeﬃcient matrix.
The drop-thresholding parameter was 2 × 10−4 . The complete factorization is performed by a sparse direct solver in the PARDISO library[9]. The elements of V were
distributed randomly on the interval [0, 1] by a random number generator.
The timing results were shown in Figure 5. The wall-clock time was 602 seconds,
and 94 eigenpairs were obtained. The time to broadcast the matrix data to allcomputing nodes was 36.4 seconds, and the time to solve linear systems was 563.8
seconds. The maximum residual of these eigenpairs was 3.4 × 10−10 .
Figure 6 shows the amount time of each step in the algorithm for the circle Γ(7) ,where gray, light gray, white, and black parts represent the time for factorization of
approximate matrix, forward and backward substitution in iteration, matrix-matrix
multiply, and others, respectively.

5

Conclusions

A parallel method for computing interior eigenvalues and corresponding eigenvectorsof generalized eigenvalue problems which arise from molecular orbital computation
was considered. In our application, matrices often have relatively large number of9Figure 5: Timing results of the method with 256 PUs in seconds. Each bar represents
the wall-clock time for one circle. 32 PUs were used for each circle, and total 256 PUs
were used. (¥: Broadcast, ¥: Linear system solver ¥: Rayleigh-Ritz procedure)

Figure 6: Amount time of each process for one circle. Each bar represents the time on
one integration node. 16 systems of linear equations were solved. (¥: Factorization,¥: Forward/backward Substitution, ¥: Matrix-vector multiply, ¥: Rayleigh-Ritz
procedure)10

nonzero elements, and we require some eigenpairs in a speciﬁc part of the spectrum.
We applied a Rayleigh-Ritz type method with numerical integration. In thismethod, a contour integration is used to construct a desired subspace. Our eigensolver does not require inner/outer loops to construct an approximate subspace andupdate approximate eigenvectors. For the computation of the contour integration,
we solve a certain number of systems of linear equations. When A and B are large,
the computational costs for solving these systems of linear equations are dominant.
Since these linear systems can be solved independently for each integration node,
the process to derive the subspace is performed in parallel.
We showed that the subspace obtained by a numerical integration includes eigenvectors corresponding to the eigenvalues around the interested interval. The i-th
eigencomponent in the subspace decays exponentially to zero with a rate |λi − γ|/ρ,
if the corresponding eigenvalue λi is located outside the interval.
Through the numerical experiments, we can see that the presented method gives
good parallel performance, and we can obtain desired eigenpairs with suﬃcient accuracy. Estimation of appropriate parameters with theoretical analysis and comparison
with other methods are part of our future work.

Acknowledgements
This work was supported in part by CREST of the Japan Science and Technology
Agency (JST).

References[1] T. Ikegami, T. Ishida, D. G. Fedorov, K. Kitaura, Y. Inadomi, H. Umeda,
M. Yokokawa and S. Sekiguchi, Full electron calculation beyond 20,000 atoms:
ground electronic state of photosynthetic proteins, Proceedings of Supercomputing 2005, Seattle, Washington, USA. 2005.
[2] T. Ikegami, T. Sakurai, U. Nagashima, A ﬁlter diagonalization for generalized
eigenvalue problems based on the Sakurai-Sugiura projection method, Technical
Report CS-TR-08-13, Tsukuba, 2008.
[3] Y. Inadomi, T. Nakano, K. Kitaura, and U. Nagashima, Deﬁnition of molecularorbitals in fragment molecular orbital method, Chem. Phys. Lett. 364 (2002)
139–143.
[4] K. Kitaura, T. Sawai, T. Asada, T. Nakano and M. Uebayasi, Pair interaction
molecular orbital method: an approximate computational method for molecularinteractions. Chem. Phys. Lett. 312 (1999) 319–324.
[5] M. Okada, T. Sakurai, K. Teranishi, A preconditioner for krylov subspacemethod using a sparse direct solver. Abstract of Precond’07, Toulouse, 2007.
[6] T. Sakurai and H. Sugiura, A projection method for generalized eigenvalue problems, J. Comput. Appl. Math. 159 (2003) 119–128.
11

[7] T. Sakurai, Y. Kodaki, H. Tadano, H. Umeda, Y. Inadomi, T. Watanabe, U. Nagashima, A master-worker type eigensolver for molecular orbital computations.
Kagstrom, B., Elmroth, E., Dongarra, Wasniewski, J. (eds.) PARA’06. LNCS
No. 4699, 617–625, Springer, Heidelberg (2007).
[8] T. Sakurai, T. Tadano, CIRR: a Rayleigh-Ritz type method with contour integral for generalized eigenvalue problems, In: The First China-Japan-Korea Joint
Conference on Numerical Mathematics, Special Issue of Hokkaido Math. J. 36(2007) 745–757.
[9] O. Schenk and K. Gärtner, Solving unsymmetric sparse systems of linear equations with PARDISO, Future Generation Comput. Sys. 20 (2004) 475–487.
[10] H. Umeda, Y. Inadomi, T. Watanabe, T. Ishimoto, U. Nagashima, Gridenabled parallel Fock matrix construction, Abstract of MATH/CHEM/COMP
2007 (MCC’07), Dubrovnik, Croatia, 2007.
[11] H. A. van der Vorst and J. B. M. Melissen, A Petrov-Galerkin type method
for solving Ax = b, where A is a symmetric complex matrix, IEEE Trans. on
Magnetics 26 (1990) 706–708.
[12] T. Watanabe, Y. Inadomi, T. Ishimoto, H. Umeda, T. Sakurai, and U. Nagashima, Molecular orbital calculation for large molecule, J. Comput. Chem.
Japan, 6 (2007) 217–226.

12