第 6 回配信 : エラーと Overlay の関係性 HPC システムズ Gaussian 入門メールニュース Gaussian 計算エラー対処 虎の巻 前回の配信では Gaussian プログラムにおける Overlay と Link の概念と それらの役割について解説しました 実際のエラー対処

HPCシステムズ Gaussian入門メールニュース

「Gaussian計算エラー対処・虎の巻」

第６回配信： エラーとOverlayの関係性

前回の配信では、Gaussian プログラムにおける Overlay と Link の概念と、それらの役割について解説しま

した。実際のエラー対処において「どの Overlay でエラーが起こったのか」という視点は極めて有効であり、

それがわかっただけで問題解決できる場合もありますし、できなくとも、より高度なエラー対処の足掛かりに

なると筆者は考えています。それを踏まえ、今回の配信ではエラー時の Overlay の確認方法や、各 Overlay

でよくあるエラーとその原因の探り方について見ていきましょう。

(6-1) Overlay と Link の確認方法

まずは「どの Overlay（Link）でエラーが起こったのか」の確認方法ですが、これを知るには通常

インプット

ファイルで「#P」を指定しなければなりません

。第１回配信で「#P」を指定する癖をつけましょうと書いたのは、

主にこのためです。「#P」ではなく、通常の「#」を指定すると、以下のようなアウトプットファイル（抜粋）が得ら

れます。

Standard orientation: --- Center Atomic Atomic Coordinates (Angstroms)

Number Number Type X Y Z

--- 1 6 0 0.000000 0.000000 -0.539640 2 8 0 0.000000 0.000000 0.687677 3 1 0 0.000000 0.939404 -1.131785 4 1 0 0.000000 -0.939404 -1.131785 --- Rotational constants (GHZ): 284.1170175 37.5107865 33.1359809 Standard basis: 6-31G(d) (6D, 7F)

There are 18 symmetry adapted basis functions of A1 symmetry. There are 2 symmetry adapted basis functions of A2 symmetry. There are 6 symmetry adapted basis functions of B1 symmetry. There are 8 symmetry adapted basis functions of B2 symmetry. Integral buffers will be 131072 words long.

Raffenetti 2 integral format.

Two-electron integral symmetry is turned on.

34 basis functions, 64 primitive gaussians, 34 cartesian basis functions 8 alpha electrons 8 beta electrons

nuclear repulsion energy 30.8309745236 Hartrees.

NAtoms= 4 NActive= 4 NUniq= 3 SFac= 1.78D+00 NAtFMM= 60 NAOKFM=F Big=F One-electron integrals computed using PRISM.

NBasis= 34 RedAO= T NBF= 18 2 6 8 NBsUse= 34 1.00D-06 NBFU= 18 2 6 8

Harris functional with IExCor= 402 diagonalized for initial guess.

この出力情報を詳しく調べると、一番上の部分は分子座標情報ですから、Overlay でいうと Overlay 2 にあ

たります。その次の部分は基底関数に関する情報であり、Overlay 3 での出力となります。そして一番下の

行は「initial guess」と書かれていますので、Overlay 4 に関連した出力情報です。このように、どの出力情報

がどの Overlay に相当するのかは、#P 指定なしでも一応可能ですが、これを読み解くことは初心者にとって

事実上不可能なことがわかります。（ここでは見やすいように、Overlay 3 の出力部分に青い帯を付けました

が、当然実際のアウトプットファイルにはそのような Overlay ごとの識別サインはありません。）

しかし、#P を指定すると、どの瞬間に Link（あるいは Overlay）が切り替わるのかがアウトプットファイルに明

確に示されます。上記の出力に相当する#P 出力は、以下のようになります。

Standard orientation:

--- Center Atomic Atomic Coordinates (Angstroms)

Number Number Type X Y Z

--- 1 6 0 0.000000 0.000000 -0.539640 2 8 0 0.000000 0.000000 0.687677 3 1 0 0.000000 0.939404 -1.131785 4 1 0 0.000000 -0.939404 -1.131785 --- Rotational constants (GHZ): 284.1170175 37.5107865 33.1359809

Leave Link 202 at Wed Sep 18 10:04:18 2013, MaxMem= 134217728 cpu: 0.0 (Enter /usr/local/g09/l301.exe)

Standard basis: 6-31G(d) (6D, 7F)

Ernie: Thresh= 0.10000D-02 Tol= 0.10000D-05 Strict=F.

There are 18 symmetry adapted basis functions of A1 symmetry. There are 2 symmetry adapted basis functions of A2 symmetry. There are 6 symmetry adapted basis functions of B1 symmetry. There are 8 symmetry adapted basis functions of B2 symmetry. Integral buffers will be 131072 words long.

Raffenetti 2 integral format.

Two-electron integral symmetry is turned on.

34 basis functions, 64 primitive gaussians, 34 cartesian basis functions 8 alpha electrons 8 beta electrons

nuclear repulsion energy 30.8309745236 Hartrees. IExCor= 402 DFT=T Ex+Corr=B3LYP ExCW=0 ScaHFX= 0.200000

ScaDFX= 0.800000 0.720000 1.000000 0.810000 ScalE2= 1.000000 1.000000 IRadAn= 0 IRanWt= -1 IRanGd= 0 ICorTp=0

NAtoms= 4 NActive= 4 NUniq= 3 SFac= 1.78D+00 NAtFMM= 60 NAOKFM=F Big=F Leave Link 301 at Wed Sep 18 10:04:18 2013, MaxMem= 134217728 cpu: 0.1 (Enter /usr/local/g09/l302.exe)

NPDir=0 NMtPBC= 1 NCelOv= 1 NCel= 1 NClECP= 1 NCelD= 1 NCelK= 1 NCelE2= 1 NClLst= 1 CellRange= 0.0.

One-electron integrals computed using PRISM. One-electron integral symmetry used in STVInt NBasis= 34 RedAO= T NBF= 18 2 6 8 NBsUse= 34 1.00D-06 NBFU= 18 2 6 8 Precomputing XC quadrature grid using

IXCGrd= 2 IRadAn= 0 IRanWt= -1 IRanGd= 0 AccXCQ= 0.00D+00. NRdTot= 244 NPtTot= 30822 NUsed= 32595 NTot= 32627

NSgBfM= 34 34 34 34 34 NAtAll= 4 4.

Leave Link 302 at Wed Sep 18 10:04:19 2013, MaxMem= 134217728 cpu: 0.4 (Enter /usr/local/g09/l303.exe)

DipDrv: MaxL=1.

Leave Link 303 at Wed Sep 18 10:04:19 2013, MaxMem= 134217728 cpu: 0.1 (Enter /usr/local/g09/l401.exe)

Harris functional with IExCor= 402 diagonalized for initial guess.;

例えば、「

Leave Link 202

」メッセージおよび「

Enter ……/l301.exe

」メッセージで L202（Overlay 2）→L301

（Overlay 3）の切り替わりが、「

Leave Link 301

」メッセージおよび「

Enter ……/l302.exe

」メッセージで L301

→L302 の切り替わりがそれぞれ明示されており、ユーザは特別な知識がなくとも、どこからどこまでがどの

Link（Overlay）なのかを知ることができます。同様に、次節で解説しますが、「どの Overlay でエラーが起こ

ったのか」も、これらのメッセージを追っていけば容易にわかります。

(6-2) Overlay からエラー原因を探る

前回配信で述べましたように、各 Overlay には固有の役割が与えられています。逆に言えば、各 Overlay

は与えられた役割以外の事はしないとも言えます。このことは、Gaussian のエラーの種類には Overlay ごと

に特徴や上限があることを意味しており、エラーが起こった Overlay からエラーの原因究明や対処が（ある

程度）可能であることを示しています。

下表に各 Overlay 内でよくあるエラー（

赤文字

は特によくありがちなエラー）をまとめました。

Overlay

大まかな役割

よくあるエラー

Overlay 0

計算の開始処理 ・Link 0 コマンドの不正指定 _{・キーワードの不正指定}

Overlay 1

分子構造の制御 ・原子座標の不正指定 ・ONIOM 計算入力の不正指定 ・構造最適化計算において更新後の座標が不良

Overlay 2

対称性の制御 ・原子座標の不正指定

Overlay 3

Fock 行列に必要な項の計算 ・基底関数の不正指定 ・電荷とスピン多重度の不整合 ・Guess=Read で、存在しない chk ファイルの読み込み

Overlay 4

initial guess の生成 ・initial guess の不正指定 _{・initial guess となる軌道入れ替えの不正指定}

Overlay 5

SCF 収束手続き ・SCF エネルギーの収束不良 _{（現象は１種類だが、根本の原因はさまざま）}

Overlay 6

各種プロパティ計算 ・NBO 計算のための追加入力の不正指定 _{・電子密度フィッティンググリッドの不正指定}

Overlay 7

AO積分の解析的微分 ・原子に働く力が異常に大きすぎる _{・詳細振動解析のための追加入力の不正指定}

Overlay 8

２電子積分変換 ・分子軌道凍結（積分変換領域）の不正指定 ・積分変換データによるディスク容量オーバー ・メモリ不足

Overlay 9

電子相関計算および励起状態計算 ・電子相関／励起状態計算パートの計算条件不正指定 ・電子相関／励起状態計算パートにおける収束異常 ・計算中間データによるディスク容量オーバー ・メモリ不足

Overlay 10

軌道係数の解析的微分 ・CPHF 方程式解法の収束不良 ・電磁場摂動の不正指定 ・メモリ不足

Overlay 11

軌道係数の微分に必要な項の計算 ・SAC-CI 計算で、CPMOD 方程式解法の収束不良

Overlay 99

計算の終了処理 ・構造最適化計算において分子構造が未収束

これを踏まえて、以下に例を示します。

GradGradGradGradGradGradGradGradGradGradGradGradGradGradGradGradGradGrad Leave Link 103 at Wed Feb 20 14:37:58 2013, MaxMem= 65536000 cpu: 0.0 (Enter /usr/local/g09/l202.exe)

Standard orientation:

--- Center Atomic Atomic Coordinates (Angstroms)

Number Number Type X Y Z

--- 1 6 0 0.000000 0.000000 -0.539640 2 8 0 0.000000 0.000000 0.687677 3 1 0 0.000000 0.939404 -1.131785 4 1 0 0.000000 0.939404 -1.131785 --- Distance matrix (angstroms):

1 2 3 4 1 C 0.000000

2 O 1.227317 0.000000

3 H 1.110457 2.047663 0.000000

4 H 1.110457 2.047663 0.000000 0.000000 Small interatomic distances encountered: 4 3 Problem with the distance matrix.

Error termination via Lnk1e in /usr/local/g09/l202.exe at Wed Feb 20 14:37:58 2013.

前節の Link 切り替わりメッセージや「Error termination」のメッセージから、この計算は Overlay 2（Link202）

でエラー停止していることがわかります。上記の表より、Overlay 2 の役割は「対称性の制御」であり、そこで

起こり得るエラーは「原子座標の不定指定」であると推測されます。事実、この計算では 3 番目と 4 番目の

原子座標が重なってしまっており、それがエラーの原因となっていることがわかります。（実際の場面では、

4 番目の原子の y 座標に、本来あるべきマイナス符号を付け忘れることによって起こり得ます。）このように、

エラーで停止した Overlay が把握できれば、それだけで問題解決に繋がることがあります。

しかし逆に、

結果的にはある Overlay で停止しているが、その根本の原因はその Overlay ではなく、もっ

と上流の Overlay にあるという可能性

も常に意識しておかねばなりません。これは、

ある Overlay でエラー

が起きても、停止せず計算を続行することがある

からです。ただその場合も、根本の原因となる Overlay で

何らかの警告的なメッセージが出力されていることが多いので、そのようなメッセージの存在に気を付けな

がら根気よくアウトプットファイルを見ていけば、問題解決に繋がります。

例を示しましょう。

#P B3LYP/Gen Opt

AuH geometry optimization 0 1 Au 0.00000000 0.00000000 0.00000000 H 0.00000000 0.00000000 1.50000000 Au 0 SDD **** H 0 aug-cc-pVDZ ****

これは、Au 原子に SDD 基底関数、H 原子に aug-cc-pVDZ 基底関数を用いて、AuH 分子の構造最適化を

行うことを意図したインプットファイルです。重要なポイントは、SDD が重原子に対する有効内殻ポテンシャ

ル（ECP）を含んだ基底関数だということです。ECP とは擬（pseudo）ポテンシャル法ともいい、化学現象に重

要な価電子（外殻の電子）のみを顕に計算し、内殻の電子の影響は外側の電子が内側から感じるポテンシ

ャルに置き換えてしまうという近似法です。

さて、このインプットで計算を行うと、以下のようなメッセージと共に計算が停止してしまいます。

Cycle 128 Pass 1 IDiag 1:

E= -3125.56918337291 Delta-E= -0.000337617204 Rises=F Damp=F DIIS: error= 1.62D-03 at cycle 128 NSaved= 7.

NSaved= 7 IEnMin= 1 EnMin= -3125.56932595477 IErMin= 1 ErrMin= 1.28D-03 ErrMax= 1.62D-03 EMaxC= 1.00D+00 BMatC= 1.13D-04 BMatP= 8.14D-05

IDIUse=1 WtCom= 1.00D+00 WtEn= 0.00D+00

Coeff-Com: -0.478D+00 0.563D+00-0.177D+01 0.160D+01-0.397D+01 0.664D+00 Coeff-Com: 0.439D+01

Coeff: -0.478D+00 0.563D+00-0.177D+01 0.160D+01-0.397D+01 0.664D+00 Coeff: 0.439D+01

Gap= 0.000 Goal= None Shift= 0.000

RMSDP=7.51D-02 MaxDP=1.53D+00 DE=-3.38D-04 OVMax= 0.00D+00 >>>>>>>>>> Convergence criterion not met.

SCF Done: E(RB3LYP) = -3125.56918337 A.U. after 129 cycles Convg = 0.7505D-01 -V/T = 13.4059

KE= 2.519425965616D+02 PE=-4.492678436684D+03 EE= 1.083997897717D+03 Convergence failure -- run terminated.

Error termination via Lnk1e in /usr/local/g09/l502.exe at Thu Sep 19 20:36:02 2013. Job cpu time: 0 days 0 hours 0 minutes 22.1 seconds.

エラーメッセージを見る限り、計算は L502、すなわち Overlay 5 で停止していることがわかります。Overlay 5

でのエラー原因は「SCF エネルギーの収束不良」ですから、一般論としては、繰り返し回数を増やしたり、収

束条件を変えたりすれば、問題解決する可能性があります。しかしこの場合は、一番上の行に「Cycle 128」

と出力されているように、128 回もの SCF 手続きを行ってもなお収束しておらず（SCF 手続きの詳細に関し

ては第４回配信をご参照下さい）、AuH という特段に異端でもない分子で、これだけ SCF 手続きを行っても

収束しないというのは少し変です。こういう場合は、根本の原因が Overlay 5 ではなく、もっと上流の Overlay

にあるため、単純に繰り返し回数上限を増やしても結果は好転しないことがほとんどです。

停止した Overlay が根本的な原因でないと予想される場合、次の方策としては２通り考えられます。１つは、

アウトプットファイルを先頭から見ていき、最初に奇妙な挙動が見られる Overlay をつきとめる方法で、もう１

つは、逆にアウトプットファイルの末尾から上の方へたどり、上流の Overlay での出力情報からエラー原因を

探る方法です。どちらが良いかはケースバイケースですが、筆者の場合、比較的短いアウトプットでは先頭

から、長いアウトプットでは末尾からまず調べていき。最終的には両方の調査結果を照らし合わせてエラー

原因を推測することが多いです。

いずれの方法にせよ、この計算例のアウトプットファイルを調べていくと、以下のような出力情報を見つけ

ることができます。

Leave Link 202 at Thu Sep 19 20:35:54 2013, MaxMem= 39321600 cpu: 0.0 (Enter /usr/local/g09/l301.exe)

Basis read from rwf: (5D, 7F)

No pseudopotential information found on rwf file.

There are 22 symmetry adapted basis functions of A1 symmetry. There are 3 symmetry adapted basis functions of A2 symmetry. There are 10 symmetry adapted basis functions of B1 symmetry. There are 10 symmetry adapted basis functions of B2 symmetry. Integral buffers will be 131072 words long.

Raffenetti 2 integral format.

Two-electron integral symmetry is turned on.

45 basis functions, 76 primitive gaussians, 48 cartesian basis functions 40 alpha electrons 40 beta electrons

nuclear repulsion energy 31.1687590324 Hartrees.

Warning! Au atom 1 has 79 valence electrons but only 36 basis functions. This is less than a minimal basis set!

IExCor= 402 DFT=T Ex+Corr=B3LYP ExCW=0 ScaHFX= 0.200000

ScaDFX= 0.800000 0.720000 1.000000 0.810000 ScalE2= 1.000000 1.000000 IRadAn= 0 IRanWt= -1 IRanGd= 0 ICorTp=0

NAtoms= 2 NActive= 2 NUniq= 2 SFac= 1.00D+00 NAtFMM= 60 NAOKFM=F Big=F No density basis found on file 724.

Leave Link 301 at Thu Sep 19 20:35:55 2013, MaxMem= 39321600 cpu: 0.2 (Enter /usr/local/g09/l302.exe)

NPDir=0 NMtPBC= 1 NCelOv= 1 NCel= 1 NClECP= 1 NCelD= 1 NCelK= 1 NCelE2= 1 NClLst= 1 CellRange= 0.0.

この出力情報から、L301（Overlay 3）で、基底関数に関する警告メッセージが発せられていることがわかり

ます。警告メッセージでは「Atom 1（Au 原子）は 79 個の電子を持っているにも関わらず、基底関数が 36 個

しかありません。これは最小の基底関数より小さいです！」とありますが、そもそも Au 原子は ECP を用いて

価電子のみを顕に計算するつもりですから、本来ならば 79 個の電子全てを顕に考慮する必要はありませ

ん。つまり、この計算のエラーの根本的な原因は、Overlay 3 での基底関数の読み込みに際し、ユーザが意

図した ECP は採用されず、外側の価電子用の基底関数のみを用いて全電子計算しようとしていることであ

るということがわかります。本来計算する意図のない内殻電子を、基底関数なしで計算しようとしているわけ

ですから、これでは計算が上手く走るはずがありません。

そこで、このメッセージを踏まえてインプットを以下のように書き換えてみます。

#P B3LYP/Gen Opt Pseudo=Read

AuH geometry optimization 0 1 Au 0.00000000 0.00000000 0.00000000 H 0.00000000 0.00000000 1.50000000 Au 0 SDD **** H 0 aug-cc-pVDZ **** Au 0 SDD

すなわち、「Pseudo=Read」キーワードで、追加入力セクションから顕に「擬（pseudo）ポテンシャルを読み込

む」ように指定し、追加入力セクションで「Au に SDD 擬ポテンシャルを採用する」ことを指定します。すると、

Leave Link 202 at Fri Sep 20 10:05:16 2013, MaxMem= 39321600 cpu: 0.1 (Enter /usr/local/g09/l301.exe)

General basis read from cards: (5D, 7F) Centers: 1 SDD **** Centers: 2 aug-cc-pVDZ **** =============================================================================== Pseudopotential Parameters =============================================================================== Center Atomic Valence Angular Power

Number Number Electrons Momentum of R Exponent Coefficient SO-Coeffient =============================================================================== 1 79 19 G and up 2 1.0000000 0.00000000 0.00000000 S - G 2 13.2051000 426.84667900 0.00000000 2 6.6025500 37.00708300 0.00000000 P - G 2 10.4520200 261.19958000 0.00000000 2 5.2260100 26.96249600 0.00000000 D - G 2 7.8511000 124.79066600 0.00000000 2 3.9255500 16.30072600 0.00000000 F - G 2 4.7898200 30.49008900 0.00000000 2 2.3949100 5.17107400 0.00000000 2 1

No pseudopotential on this center.

===============================================================================

（中略）

Two-electron integral symmetry is turned on.

45 basis functions, 76 primitive gaussians, 48 cartesian basis functions 10 alpha electrons 10 beta electrons

というアウトプットが得られますが、これは、Au の ECP が正しく読み込まれており、かつ、顕に計算するべき

電子数が計 20 個と、正しく内殻の電子が除外されていることを示しています。

この例では警告メッセージに顕に基底関数のことが言及されているので、そのメッセージさえ見つかれば、

Overlay のことを知らなくてもエラー対処できるかもしれません。しかし「何か警告らしきメッセージがあるが、

それがどういう意味なのかわかりにくい」こともしばしばあり、その場合は警告メッセージの Overlay と前節の

表を参考にしながら、エラーを類推して対処していくことになるでしょう。

(6-付録) HF および DFT 振動計算の理論的背景

ここでは HF/DFT 計算の分子振動計算の理論的背景をごく簡単に概説します。これまでの配信と同様、

これにより第２部の配信内容の理解をより深めることができますが、あくまで付録ですので、興味のない方は

読み飛ばしていただいて構いません。特に今回の内容は、難解と感じる方も多いと思いますので、

理論の

結果だけ知りたい方は、最後の２段落だけ読んで下さい

。また、分子振動はエネルギーの原子核座標に関

する２次微分で与えられますが、以下の議論はエネルギーの２次微分量であれば一般的に成り立ちますの

で、振動計算だけでなく、分極率、分子磁化率、核磁気共鳴などにも適用可能です。

前回の配信の通り、HF/DFT 計算での分子の全エネルギーは以下のように与えられます。

nuc

2

1

V

J

C

C

C

C

n

n

h

C

C

n

E

ij j j i i j i i i i i







 

 

         

(6.1)

また式(6.1)の原子核座標 X に関する（１次）微分表現は以下の式で与えられます。

X

V

X

S

C

C

n

X

J

C

C

C

C

n

n

X

h

C

C

n

X

E

i i i i i ij j j i i j i i i i i





























 

 

 

nuc

2

1































(6.2)

各項の詳細は前回・前々回配信の付録に譲りますが、ここで重要な点は、式(6.1)を微分するにあたっては、

厄介な軌道係数行列 C の微分は特別な方法によって回避することができ、最終的な表現が式(6.2)で表さ

れるように C の微分を含まない形で書き表されるということです。したがって、新たに計算しなければならな

い微分量は、比較的容易に計算可能な h, J, S, V

nuc

のみです。



さて、分子の振動を計算するには、式(6.2)を原子核座標 Y でさらに微分してやる必要があり、その表現は、

以下のようになります。（なお、Y=X の場合も、以下の議論は成立します。）

XY XY

Q

P

X

Y

E

_

_







2

(6.3)

ただし

X

Y

V

X

Y

S

C

C

n

X

Y

J

C

C

C

C

n

n

X

Y

h

C

C

n

P

i i i i i ij j j i i j i i i i i XY

































 

 

 

nuc 2 2 2 2

2

1

             



(6.4)













 

 

 































i i i i i ij j j i i j i i i i i XY

X

S

Y

C

C

n

X

J

C

C

Y

C

C

n

n

X

h

Y

C

C

n

Q

             



(6.5)

ここで、P

XY

_{は軌道係数 C の微分を含まない項、Q}

XY

_{は C の微分を含む項を表しています。エネルギー２次}

微分では１次微分の時と異なり、Q

XY

_{において C の微分を回避することはできず、まともに計算せねばなりま}

せん。このため、分子振動の計算は力の計算に比べて数倍程度の時間を要します。

さて、原子核座標（分子構造）をわずかに変位させると、分子軌道の形もわずかに歪みます。この歪みは、

元の分子軌道が分子構造の歪みによって他の分子軌道と少しずつ混ざることによって生じたと解釈すると、

実用的な計算が可能となります。そうすると、軌道係数の微分は、他の分子軌道係数の線形結合の形で書

けると考えることができます。









all a Y ai a i

U

C

Y

C





(6.6)

ここで行列 U

Y

_は

coupled perturbed Hartree-Fock (CPHF)係数

行列と呼ばれ、摂動 Y によって分子軌道

がどのくらい混ざり合うかを表す行列です。式(6.6)を式(6.5)に代入すると、





































i i a a Y ai i i ij j j i a a Y ai j i i i a a Y ai i XY

X

S

C

C

U

n

X

J

C

C

C

C

U

n

n

X

h

C

C

U

n

Q































2

2

2

(6.7)

となり、この中で未知の量は CPHF 係数行列 U

Y

_{だけとなります。}

残りの問題は U

Y

をどう決めるかですが、これは「Fock 行列が満たすべき方程式はいかなる摂動 Y によっ

ても不変となる」ことを利用して決定します。具体的には、（MO 表現の）Fock 行列が満たすべき関係式

)

(

0

m

n

J

n

h

f

k mnkk k mn mn











(6.8)

の両辺を微分した、

0















k mnkk k mn

Y

J

n

Y

h

(6.9)

という関係式を利用します。式(6.9)は「分子構造を Y によっていかに変位させようと、式(6.8)は常に満たされ

ねばならない」ということを意味しています。式(6.9)を少し変形すると、









0

) ( ) (

_





































k k k n m Y mnkk k n m Y mn

J

Y

C

C

C

C

J

n

h

Y

C

C

h





















(6.10)

ただし









_



























Y

J

C

C

C

C

J

Y

h

C

C

h

_mn(Y) _{m n}

,

_mnkk(Y) _{m n k k}

(6.11)

となりますから、詳細は省きますが、これに式(6.6)を代入して整理すると、最終的に以下の

CPHF 方程式

と

呼ばれる方程式が導かれます。

Y mn unocc a occ i Y ai mnai

U

B

A



 

(6.12)

ただし、



_{mnai nmai}



_{ma ni}



_{m n}



mnai

J

J

A



2















(6.13)















ai nmai Y ai Y mn n k Y mnkk Y mn Y mn

h

J

S

S

J

B

( ) ( )



( )

2

( )

(6.14)

A は既知の２電子積分より、B

Y

は既知の１電子積分・２電子積分およびそれらの微分によりそれぞれ計算さ

れる行列であり、いずれも容易に計算することができる量です。CPHF 方程式(6.12)は連立一次方程式の

形をしており、これを解くことによって U

Y

_{が得られます。これで、最終的に得たい量(6.3)を計算するために}

必要な量全てが揃ったことになります。

ここで、エネルギー２次微分に関する以上の流れをまとめ、Gaussian での Overlay との関係性について考

えてみたいと思います（実際の Gaussian プログラムでは、効率的な計算を行うため、必ずしも以下の流れに

沿っているわけではありませんが、本質的には差し支えありません）。振動計算を行うためには、式(6.3)、す

なわち、P

XY

_{と Q}

XY

_{の計算が必要ですが、このうち Q}

XY

_{を求めるには CPHF 方程式(6.12)を解かねばなりませ}

ん。さらに方程式(6.12)を解くためには、式(6.13)、(6.14)を計算して行列 A と B

Y

_{を用意してやる必要があり}

ます。この A と B

Y

の計算を Overlay 11 で行います。次に、A と B

Y

が求まったら、CPHF 方程式(6.12)を解く

作業を Overlay 10 で行います。そして、この方程式を解いて CPHF 係数 U

Y

_{が求まったら、P}

XY

_{と Q}

XY

_の計

算が可能となりますので、それを Overlay 7 で行います。

前回配信では、話が難しくなることを避けるために、微分計算パート（Overlay7、10、11）の役割について

ほとんど触れませんでしたが、実際の振動計算アウトプットファイルから Overlay の流れを読み取ると、大抵

Overlay 11→10→7 という流れを含んでいることが読み取れます。これは上記の内容と一致していることが

わかります。また、MP2、CCSD、SAC/SAC-CI などの電子相関計算の構造最適化計算では、（１次微分で

あっても）顕に軌道係数の微分を行わねばならないため、同様な Overlay の流れとなります。一方、単なる

HF/DFT 構造最適化計算の場合は、軌道係数の微分が回避できますから、CPHF 方程式を解く必要はあり

ません。そのため、（前回配信の付録でも触れたように）１電子積分と２電子積分の微分を行う Overlay 7 ま

でが呼び出され、CPHF 方程式に関係する Overlay 10、11 は呼び出されません。

今回の内容は以上です。また「第２部 Gaussian でエラーが起こる仕組みを知る」もこれにて終了し、次回

からは「第３部 Gaussian エラー対処各論」に入ります。次回（第７回配信）では、「Link 0 コマンド指定に関

するエラー」を解説致します。

☆本メールニュースの内容は、過去配信分も含め、弊社ホームページ上に掲載されております。配信内容

のご感想やご希望に関する簡単なアンケートも行っておりますので、よろしければご意見をお寄せ下さい。

http://www.hpc.co.jp/gaussian_nyumon.html