Technical Report, University of Tsukuba, CS-TR-08-14 (2008)

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 the present paper, we consider a parallel method for computing interior
eigenvalues and corresponding eigenvectors of generalized eigenvalue problems
that arise in molecular orbital computation of biochemistry applications. Matrices in such applications are sparse but often have a relatively large numberof nonzero elements, and we may require some eigenpairs in a specific part ofthe spectrum. We use a contour integration to construct a desired subspace.
Properties of the subspace obtained by numerical integration are discussed,
and a parallel implementation is then presented. We report the numerical aspects and parallel performance of the proposed 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 the present paper, we consider a parallel method for computing eigenpairs (λ, x)satisfying Ax = λBx in a specific part of the spectrum that arise in molecularorbital computation, where A, B ∈ Rn×n are symmetric and B is positive definite.
Molecular orbital computation is performed in order to investigate phenomena suchas the reaction mechanisms of enzymes and the electronic structure of photosynthetic
systems.
The fragment molecular orbital (FMO) method [4] enables the total energy ofa molecule to be calculated without performing molecular orbital computation for
the entire molecule. In [1], a full electron calculation for a large molecule was
performed 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 canonical MOs without SCF iterations. A huge Fock matrix is
generated, and we require a limited number of eigenpairs corresponding to frontier
1

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

orbitals. Therefore, a large-scale interior eigenvalue problem occurs in the FMO-MO
method.
With growing parallel efficiency in the Fock matrix construction, the diagonalization process becomes a bottleneck in MO calculations. Since the Fock matrix
has a relatively large number of nonzero elements due to the base function of the
middle-range interaction of molecules, a sparse factorization of the shifted matrix
for spectral transformation may not be feasible. An alternative to factorization isoffered by iterative solvers. Unfortunately, this is usually not an attractive option
for parallel computation, because the iterative process is performed successively.
Therefore, the shift-and-invert approach is not effective for treating interior eigenvalue problems in such applications.
The proposed method for finding eigenpairs in a given physical domain is basedon the contour integral presented in [6]. One of the major advantages of the proposed
eigensolver is that it does not require an inner loop to construct an approximate
subspace or an outer loop to update approximate eigenvectors. Moreover, the values
that are used in a numerical integration can be evaluated independently on each
integration node, which 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 in which
the coefficient 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 coefficient matrix, is effective for solving suchlinear 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 present an implementationof the proposed method. In Section 4, numerical examples are presented in orderto verify the numerical properties and the parallel performance of the method withthe matrices derived from molecular orbital computation.

2

An eigensolver using contour integration

In this section, we describe an eigensolver using the contour integration presented
in [8]. This method finds 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 center γ ∈ R and 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 .

2

Suppose that v is expanded by the eigenvectors {x1 , . . . , xn } asv=n"αi xi .

(2)

i=1

Then, we have the following theorem([8]).
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 canbe extracted 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 following systems of linear equations(ωj B − A)yj = Bv,

j = 0, 1, . . . , N − 1(3)for y0 , . . . , yN −1 ∈ Cn . The solutions yN/2 , . . . , yN −1 are obtained from the relation
yj = ȳN −j−1 . Thus, we only need to solve N/2 systems, and%#$k+1 &
N/2−1
1 "ωj − γ
ŝk =
2 Re
yj ,
N j=0
ρ

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

(4)Note that ŝk are real vectors.
Letting θj = exp(2πi(j + 1/2)/N ), the following relations 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

3

η

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
N −1

1 " p
θ =
N j=0 j'(−1)q if p = qN for q ∈ Z.
0otherwise(7)

Similarly, for the case in which |η| > 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

!

Thus, (5) is obtained from (6) and (8).
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 =it follows thatŝk

N −1 #

1 "
=
N j=0

N −1 #

ωj − γ
ρ

1 "
=
N j=0

ωj − γρ1
N j=0

ωj − γ
ρ

=

=

N
−1 #
"

n
"
αi

ρ
i=1

%

$k+1

1
xi ,
ωj − λi(ωj B − A)−1 Bv$k+1 "ni=1αi (ωj B − A)−1 Bxi$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

ρ

·

4

ηik
xi .
1 + ηiN(9)!

This proves the theorem.
For vectors a, b ∈ Rn , we define 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 .
Theorem 4 Assume that the eigenvectors x1 , . . . , xn are normalized to be B-orthogonal,i.e., %xi , xj &B = δij . Then,|%xi , ŝk &B | =( αi

η 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 eigensolver is proposed in [2], which enables us to obtaindegenerate eigenvalues. In the eigensolver, a matrix V = [v1 , . . . , vL ] ∈ Rn×L is usedinstead of v in (3), where v1 , . . . , vL are linearly independent, and positive integer
L is the 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 the following 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.
In choosing 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 from Ŝ(:, 1 : m̂). The
algorithm is shown below.
Algorithm (block method using contour integration):
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* k+1 +)
2. Compute Ŝk = N/2−1
2Re
θj Yj , k = 0, . . . , M − 1
j=0

3. Perform singular value decomposition U T ΣW = [Ŝ0 , . . . , ŜM −1 ], and findm̂ such that |σj |/|σ1 | ≥ δ for 1 ≤ j ≤ m̂

4. Construct an orthonormal basis Q̃ from Ŝ(:, 1 : m̂)5

5. Form Ã = Q̃T AQ̃ and B̃ = Q̃T B Q̂6. Compute the eigenpairs (λ̂i , ŵl ), 1 ≤ l ≤ m̂ of the projected pencil (Ã, B̃)
7. Set x̂i = Q̃ŵi , 1 ≤ l ≤ m̂

Note that we can use U (:, 1 : m̂) instead of Q̃ so as to avoid the orthogonalizationof Ŝ(:, 1 : m̂).

3

An implementation

When matrices A and B are large, the computational costs for solving systems oflinear equations (10) are dominant. In the present implementation, the matrices
have a relatively large number of nonzero elements, which could make sparse direct
methods impractical due to the large number of fill-ins required. Instead, we apply
a Krylov subspace iterative method with the complete factorization preconditioner
presented in [5].
In the preconditioning method of the present study, a complete factorizationof the approximate matrix C̃j = (c̃ij ) of the coefficient matrix Cj := ωj B − A is
performed. The approximate matrix C̃j is obtained from the drop-thresholding of
the original coefficient matrix Cj . Since there are fewer nonzero elements in C̃j thanin Cj , we expect fewer nonzero elements in the preconditioner than the matrix factor
obtained from complete factorization of Cj . Drop-thresholding is defined 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 multiplications are performed in
the iterative process.
We can easily extend the method for the case in which several circular regions
are given. Suppose that Nc circles Γ(0) , . . . , Γ(NC −1) with center (γ (l) , ρ(l) ) and radius
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.
First, 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 the results of the6

linear solver, and then perform the Rayleigh-Ritz procedure to extract eigenpairs ineach circle. Since the computational costs to solve systems of linear equations dominate the computations, we can expect high performance in parallel computations.
The algorithm for parallel implementation is summarized as follows:
Parallel implementation (circles Γ0 , . . . , ΓNC −1 ):1. Broadcast sparse matrix data of A and B
2. do l = 0, 1, . . . , NC − 1
3.
4.
5.

(l)

(l)

Solve (ωj B − A)Yj
(l)

(l)

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

(l)

(l)

(l)

Compute singular values σ1 , . . . , σM of Ŝ (l) = [Ŝ0 , . . . , ŜM −1 ],
(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

4

Numerical examples

In this section, we present numerical examples with the matrices derived from the
FMO-MO method[10].
Example 1 The test matrices A and B are derived from the computation of themolecular orbitals of an eight-base-pair model DNA. The size of the matrices isn = 1, 980, and the number of nonzero elements is 728, 080.
The CPU was a Core2Duo (2.2 GHz) with 2 GB of memory. Computation was
performed using MATLAB 7.5 with double-precision arithmetic. The systems oflinear 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 1, 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 Γ. The eigencomponents
associated with eigenvalues that are located outside Γ decay exponentially to zeroin ŝ0 . For vectors ŝk , k ≥ 1, we can obtain similar results. In Figure 2, we show7

|%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.
A total of 16 eigenvalues are located in the interval [−0.21, −0.19], and all of theeigenvalues in this interval were obtained with sufficient accuracy.

Figure 1: |%xi , ŝ0 &B | (γ = −0.2, ρ = 0.01, N = 32)

Figure 2: |%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 a target molecule for anticancer agents. The Fockmatrix was constructed by the method described in [10]. The dimension of the Fockmatrix was 96, 234, and the number of non-zero elements was 457 million.
The calculations were performed on a P32 subsystem of the AIST Super Cluster(ASC) at National Institute of Advanced Industrial Science and Technology Agency(AIST). Each computing node is an Opteron Dual processor (2.0 GHz) with 6 GBof memory 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 placed eight circles around the highest occupied molecular orbital (HOMO) andthe lowest unoccupied molecular orbital (LUMO) using the approximate eigenvalues
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

obtained from the results of the FMO method. The total number of systems to besolved was 8 × (32/2) = 128 with 16 right-hand-side vectors. Since each system of
linear equations was solved on a single node of a dual processor, 16 nodes (32 PUs)
were used for one circle. Thus, a total of 256 PUs were used.
The preconditioned COCG method[11] was used as the iterative linear solver.
The stopping criterion for the relative residual norm was 10−10 . The preconditioner
was constructed by applying a complete factorization for an approximate coefficientmatrix, which was obtained from drop-thresholding of the original coefficient 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 are shown in Figure 3. 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 4 shows the amount of time of each step in the algorithm for the circle(7)
Γ , where gray, light gray, white, and black segments represent the factorization
times for the approximate matrix, forward and backward substitution in iteration,
matrix-vector multiplication, and the Rayleigh-Ritz procedure, respectively.

9Figure 3: Timing results of the method with 256 PUs in seconds. Each bar represents
the wall-clock time for one circle. For each circle, 32 PUs were used, for a total of
256 PUs. (": Broadcast, ": Linear system solver ": Rayleigh-Ritz procedure)

Figure 4: Time for each process for one circle. Each bar represents the time on
one integration node. A total of 16 systems of linear equations were solved. (":Factorization, ": Forward/backward Substitution, ": Matrix-vector multiply, ":Rayleigh-Ritz procedure)105Conclusions

A parallel method for computing interior eigenvalues and corresponding eigenvectorsof generalized eigenvalue problems that arise from molecular orbital computation
was considered. In the present application, matrices often have a relatively large
number of nonzero elements, and eigenpairs are required in a specific part of thespectrum.
We applied a Rayleigh-Ritz procedure with numerical integration. In this method,a contour integration is used to construct a desired subspace. The proposed eigensolver does not require inner and outer loops to construct an approximate subspace
and update approximate eigenvectors. For the computation of the contour integration, we solve a 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 by which 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.
The numerical experiments indicate that the proposed method gives good parallel
performance and that desired eigenpairs can be obtained with sufficient accuracy. In
the future, we intend to estimate appropriate parameters with theoretical analysis
and compare the proposed method with other methods.

Acknowledgements
The present study was supported in part by CREST of the Japan Science andTechnology 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 filter 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, Definition of molecularorbitals in fragment molecular orbital method, Chem. Phys. Lett. 364 (2002)
139–143.

11

[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.
[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