とても最近に、やらかした計算上のバグ、失敗 [Top][失敗知識データベース（畑村創造工学研究所←JSTから移動）]

（い）本日のバグ（レポート１０／１３、２０１１）
（Ｚ）本日のバグ（レポート１／２１、２０１１）
（Ｙ）本日のバグ（レポート３／２、２０１０）
（Ｘ）本日のバグ（レポート９／４、２００９、１２／４、２００９追記）
（Ｗ）本日のバグ（レポート８／５、２００８）
（Ｖ）本日のバグ（レポート２／８、２００８）
（Ｕ）本日のバグ（レポート１２／１７、２００６）
（Ｔ）本日のバグ（レポート８／２１、２００６）
（Ｓ）本日のバグ（レポート９／１６、２００５）
（Ｒ）本日のバグ（レポート２／１１、２００５）
（Ｑ）本日のバグ（レポート９／２６、２００５）
（Ｐ）本日のバグ（レポート１／１２、２００５）
（Ｏ）本日のバグ（レポート２／２５、２００４）
（Ｎ）本日のバグ（レポート４／９、２００３）
（Ｍ）本日のバグ（レポート１２／５、２００２）
（Ｌ）本日のバグ（レポート１２／１９、２００１）
（Ｋ）本日のバグ（レポート７／１６、２００１）
（Ｊ）本日のバグ（レポート３／２、２００１）
（Ｉ）本日のバグ（レポート９／１８、２０００）
（Ｈ）本日のバグ（レポート８／２９、２０００）
（Ｇ）本日のバグ（レポート５／２５、２０００）
（Ｆ）本日のバグ（レポート２／２４、２０００）
（Ｅ）本日のバグ（レポート１／１７、２０００）
（Ｄ）本日のバグ（レポート４／２０、１９９９）
（Ｃ）本日のバグ（レポート１／２０、１９９９）
（Ｂ）本日のバグ（レポート１１／１８、１９９８）
（Ａ）本日のバグ（レポート１１／２、１９９８）
（９）本日のバグ（レポート８／４、１９９８）
（８）本日のバグ（レポート３／１７、１９９８）
（７）本日のバグ（レポート１／２０、１９９８）
（６）本日のバグ（レポート１２／５、１９９７）
（５）本日のバグ（レポート１０／１５、１９９７）
（４）本日のバグ（レポート６／１８、１９９７）
（３）本日のバグ（レポート３／２６、１９９７）
（２）本日のバグ（レポート２／６、１９９７）
（１）本日のバグ（レポート１／２３、１９９７）

（い）本日のバグ（レポート１０／１３、２０１１） [目次]

今回もまたやってしまいました。今回の失敗も大変単純でした。
ある並列計算（OpenMPを使用）で、力の値がおかしいことに気付きました。計算対象の構造の持つ対称性から、そもそも力がゼロになるべきなのに、有限の値が得られていました。それも同じ条件で計算すると、出てくる結果がまちまちで、たまに力がゼロ（正しい結果）になることもありました（計算毎に値が違う）。全エネルギーなどの値は、どの場合でも正しい結果を与えました。明らかに、並列計算関連でおかしなことになっていると分かる状況ですが、原因究明に割と時間がかかりました。
原因は、力の計算部分での足し上げのところでした。逐次計算なら問題ないのですが、並列で計算する場合に何も対処していないとおかしな結果となる部分がありました。対処の仕方は、REDUCTION（OpenMP）によって変数を指定しておくことです。


!$OMP DO REDUCTION(+:ZFORC2)

ZFORC2が足し上げ対象となる変数。この１行がなかったために、毎回計算するたびに異なった結果（力の値）を与えていました。
そしてもっと問題なのが、この失敗は既に過去に経験していたことです （参考文献参照）。並列化されたプログラムのバージョン（版）管理にも問題がありました。筆者は、いわゆるバージョン管理ソフト等を使用していません。で、今回、古いバージョン（版）のプログラムで並列計算を行なってしまいました。筆者は、そのことに全く気付きませんでした。本当は、大元の最新プログラムと比較検討することが容易に可能であり、そうすべきだったのですが、全くそれを行ないませんでした。
REDUCTIONは、並列計算で数値の足し上げを各ＣＰＵ毎に行なう時、それぞれのＣＰＵ上で計算された値をまとめあげるための指示を与えます。そうしないと正しい足し上げ結果とならずに、実行毎に値が異なったりもする。ただ、資料によってはREDUCTIONで指定する変数はスカラー変数のみで配列変数は駄目という記述も散見されます。筆者の場合、REDUCTIONで指定している変数 (ZFORC2)は配列変数（並列計算に対しては独立）です。一応、それでも正しく動いていますが（非並列時の結果と整合している）、これに関しては更に検討・調査を行なっています（調査に結果、新しいバージョン〔FORTRAN〕では配列変数の指定が可能な模様）。

関連ページ：OpenMP導入に関するメモのメモ１参照。

参考文献：小林一昭、「徹底解剖第一原理計算：第九回”計算上の落し穴（問題点）”」、金属、第８０巻、第６号（２０１０年６月号）、５８頁（日本語）、[情報](J-GLOBAL, jst)

問題は、これまでも何度も指摘したように、同じか或いは類似した失敗を繰り返していることである。これにはより根源的（研究姿勢、環境等）な問題、欠陥の存在が示唆される（のかもしれない）。

（Ｚ）本日のバグ（レポート１／２１、２０１１） [目次]

今回の問題は、バグとは言えないものですが、ここに記しておきます。あるバンド計算プログラムを動かしている時、これまで問題なく動いていたものが、新たな系を扱ったところ突然計算途中で止まってしまうようになりました。計算が突然停止したことに対する直接的なエラーメッセージはなく、流したジョブのエラーメッセージも（筆者にとっては）意味不明なものでした。最初は、単純な入力データの間違いと思い、いろいろ調べましたがどのようにしても計算が途中で止まってしまいました。このバンド計算プログラムの場合、対称性の設定や、単純に単位胞内の原子数の設定を間違っても、計算が（何のエラーメッセージなく）止まってしまうのですが、今回はそれとは異なりました。入力等の設定間違いの場合、ジョブ全体に関してのエラー表示はないのですが、今回は意味不明ながらエラー表示（メッセージ）がありました。

ちょっと悩んだ末に、原因はメモリが足りないということに気付きました。案の定、同じ系で計算規模を少なくすると（今回の場合、ｋ点数を少なくしてみた）、計算は完結するようになりました。メモリー不足の場合、これまでは意味のあるエラー表示（見てすぐに判断出来るということ）が出ていたのですが、今回は少なくとも筆者には全く意味の分からない表示（メッセージ）しかありませんでした。このため原因に気付くのが遅れました。

今回使用したプログラムは、メモリを浪費する版だったので、少し大きな系を実行させようとすると障害（今回のように突然止まる）が起きることは、十分に気を付けていればすぐに「メモリ周りが怪しい。」と気付いたところですが、そうなりませんでした（やはりエラー表示が理解出来なかったのが痛かった。が、これを理解出来る方が筆者はおかしいと思う。残念ながら当該表示〔メッセージ〕は、ここで明らかにすることが出来ません）。
尚、当該プログラムは、よりメモリを消費しない版に改良し、当初扱った系でも問題無く動く（計算が完結する）ようになりました。

（Ｙ）本日のバグ（レポート３／２、２０１０）[目次]

今回も盛大にやってしまいました。非常に単純ながら大きな失敗をしてしまいました。
今回の失敗はプログラムのバグではありません。実は、計算としては間違った計算ではなく正しい計算を行なっています。しかし大変な失敗計算になっています。

で、どのような間違いを犯したかと言うと、ある系のバンド計算を行なう時、原子数６０個／単位胞、１個当りの電子数は４個（価電子数）という前提で、 ２４１バンドとして計算してしまったことです。ここでスピンは縮退して（非磁性の場合）います。２４１としたのは、２４０バンドまで電子が詰まっているとして（←ここ重要）、１個分多くしました。扱った系はギャップが空いていて、それもワイドギャップなので１バンド分多いだけで問題ありません。少なくともより小さい系での検証では、全く問題ありませんでした。

しかし、これは大変勘違いな計算をしていることになります。電子数２４０個の系で、スピンが縮退している場合、必要なバンド数は１２０です。実際は、金属である場合もあるので、これより多めのバンド数で設定します。広いギャップが空いている系なら、多めと言ってもほんの少しで十分です。１２１でも大丈夫なはずです（バンドが非常に狭い場合などは注意が必要）。実際は、バンド数、１２３で計算（２４１としていた場合と比べて同じ全エネルギーを返すことを確認）しています。

バンド数を２４１として計算しても間違いではありません。この場合系は金属ではなく、バンドギャップを持っているので、２４１バンドの内１２１バンド分が空のバンドとして計算されています。これらの空のバンドは、実質計算結果に寄与しません。空のバンドからの寄与が必要な、GW近似などでは重要ですが、今回の計算はごく普通のバンド計算（ギャップの空いた系の基底状態を求める電子状態計算）なので空のバンドの影響はありません。つまり正しい結果を返している訳です。その意味では間違いではない。

しかしながら、これは原子数６０個という筆者がこれまで扱ってきた系でも、最も大きな部類に入ります。勿論、昨今のバンド計算ではもっとずっと大きな系の計算も平気で行なわれていたりします。しかし、それでも１２０バンド程度で済む計算を、２４０バンド（正確には２４１バンド）で行なっていました。バンド計算の計算時間は、大体バンド数に比例します（注：グラムシュミットの直交化は、バンド数の２乗のオーダー）、メモリーもおよそ比例します。つまり２倍分の余計な計算をしてしまっていたのです。それもかなり時間が経つまで気付きませんでした（気付いたのも偶然に近い。ふと定義、設定を見ていて気付いた）。これはちょっと悲しいです。
バンド計算を始めてから２５年、未だにこのような初歩的なところで間違いを犯しています。情けない限りです。

尚、１２０バンドまで電子が詰まっている系で、２４１バンド（１２１バンドは空のバンド）という、大量の非占有バンドがある状況の計算は、これまで計算した経験がなかったが、今回試すことが出来た。計算を進めると全エネルギーは完全一致ではなくなったが、精度的には依然問題ないレベルの差しかない。
（３／５、２０１０）この間違いに気付いたのは、既に述べたように本当に偶然だった。ふとインクルードファイル” PACVPP”（参考例、今回のものとは異なります）を見ていたら、「あれっ？」と思い気付いた（その前に、別のプログラムソースを見ていて違和感を憶えたのが始まり）。ちょうどその時、別件でソフトウェアの導入で苦労していて、漸くうまくいきかけていた頃で、気分的（大変心許ない話だが）に気付き易い状況にあったのかもしれない。

（Ｘ）本日のバグ（レポート９／４、２００９）[目次]

今回の失敗は、プログラム上のバグではない。これとは別に、最近小さなレベルでのバグの発見や、検証計算があった。例えば、未定義の配列を使用していて、通常は動いていたが、ある条件の計算で突然実行が止まって（←当該部分に、IF判定によるSTOP命令があった）、気付くということがあった。

で、今回の間違いは、形成エネルギー計算における、単純で深刻な誤りだった。形成エネルギーに関しては、形成エネルギーと凝集エネルギーの関係を参照して欲しい。計算で、化合物の各構成元素１原子当たりのエネルギーに換算するためのユニットセル内の原子の数が間違っていた。少ない数で割ったため、化合物のある構成元素の原子１個当たりのエネルギーを過剰に低く（安定に）評価したため、得られる形成エネルギーがかなり損になる（不安定な）値になってしまった。

あまりにべらぼうに損（不安定）になる値であったので、再度検討をしてみた結果、ユニットセルの原子数の値が間違っていることに気付いた。ちょっと情けない話である。検算はちゃんと行なわなければならない。（追記：１２／４、２００９）
今回追加として扱うべき誤りがありました。基本的に今回の誤りは、過去の” （Ｂ）本日のバグ（レポート１１／１８、１９９８）” と同じものです。入力データとしての擬ポテンシャルデータに間違いがありました。元素A、B、Cがあったとして、ABCという順で擬ポテンシャルデータを作るべきところを、ABDと全く関係ない元素Dの擬ポテンシャルが入ったデータを入力データとして使ってしまいました。他の関連入力データ、設定は、ABCという前提で計算したので、これでは全く正しい計算になっていません。原因は、確実に詰めた訳ではありませんが、元々ACDという入力データがあり、これからABCを作成しようと、エディタ上で編集した際に間違いが起きたと考えられます。ACDからDを削除後、ACにBをAとCの間に挿入するのが正しい作業手順でした。おそらく、ACDにBを、ABCDとなるように挿入し、削除すべきはDなのにC を削除し、ABDを作ってしまったようです。
残念ながら、それでも筆者のプログラムは正しく動いてしまいます。正しく動作しますが、当然結果は正しくありません。この間違いに気付いたのは、全エネルギーから形成エネルギーを求めると、目茶苦茶不安定になる結果が出たためです。通常は、不安定になると言っても、数eV／原子以下のエネルギー差ですが、数a.u.（ハートリー原子単位）／系全体の差（原子当りでも、数十eV 以上）が出ました。いくら何でもこいつはおかしいと思って、データを調べていたら先の擬ポテンシャル並びに関する誤りに気付きました。

（Ｗ）本日のバグ（レポート８／５、２００８）[目次]

今回のバグは、計算した系の総電子数にゴミが出るというものです。それは、サブルーチンXCFFT で総電子数（ユニットセル内の系の電子の総数）を計算すると、以下のような値（総電子数８０個の場合）になることから発覚しました。


 REAL TOTAL CHARGE =    80.0026633146275827       IN XCFFT

本当の総電子数は８０個／ユニットセルですが何故か、0.0026633146275827 とゴミのような余計な部分が出てきます。もし正しく総電子数が求まるなら、以下のような値になります。


 REAL TOTAL CHARGE =    80.0000000000007       IN XCFFT

ゴミ（誤差）部分は無視出来るほど小さな値になります。この問題は、これまでとは異なる数値演算ライブラリのFFTルーチンを利用した場合に起こりました。当初は、これに全く気付きませんでした。それは総電子数にゴミがあるのに、全エネルギーなど他の値に問題が生じなかった（総電子数が正しく出る場合と値が完全に一致する）ためです。こんなに電子数の値に違いがあるのに、何で全エネルギー等が一致するのが最初不思議でした（後述。←現在でも完全に裏付けは取れていない）。で、結局原因と思われる部分として以下のルーチンが浮かび上がりました。このルーチンの一番外側のIF文の中にあるIF文、 "IF ( MOD(I,IFX2).NE.0 ) THEN"ですが、これはFFTにおけるメッシュに関係するもので、筆者が元々使用していたFTTルーチン(MFFT)では、３次元複素フーリエ変換において、メッシュ数：IFX2が2n-1でした（nは適当な正の整数）。 FFTとして使用するメッシュは、2n-2個分必要で、それより一つメッシュが多い仕様でした。従って、2n+1番目のFFT用の配列は、実際の物理量の計算には関わらない部分なので、総電子数計算の際に相当する部分の電荷密度の値（CHGB1(I ← IFX2)）をゼロとしました。ただ、以下のルーチンでは問題があります。一番外側のIF判定は、CHGB1がゼロ以下の場合という制限を与えています。従って、そうでない場合、次のIF判定である、"IF ( MOD(I,IFX2).NE.0 ) THEN"が意味を持ちません。
MFFT使用時は、以下のような判定処理をしないでも正しい値が得られるようになっていました。あくまで、2n+1番目の部分は作業用でしか使用されないし、それが結果に影響を与えないようになっていました。今回は、MFFTではなく他のFFTルーチン（数値演算ライブラリ）を使用しました。そもそも今回のFFTでは、2n+1番目の部分は必要なかったのですが、元々のプログラムがMFFTで動く仕様になっていあたため、それと整合させるため、2n+1（実際使用されるのは、 2n-2まで）として計算していました。で、筆者は以下のルーチンで、2n+1部分の処理は行なわれていると思っていました。前述のように、2n+1番目に相当するCHGB1が必ずゼロ以下である必然性はなく、実際2n+1番目にゴミ（としての数値）が残ったままで計算が行なわれてしまった訳です。そしてこれが結果として総電子数のゴミとして現れました。


      ICCCC=0
      DO 20 I=1,KSUM                                                    
        IF( CHGB1(I).LE.0.0D0 ) THEN                                      

          IF ( MOD(I,IFX2).NE.0 ) THEN

              IF ( CHGB1(I).LE.-1.0D-5 ) THEN                            
                WRITE(6,610) I,CHGB1(I)                                  
  610           FORMAT(1H ,'**WARNING CHG.DEN<0.0 AT ',                
     &                                  I5,2D15.7,'***')               
              ICCCC=ICCCC+1
              END IF                                                   
              CHGB1(I) = 1.D-40                                           
            ELSE                                                        
              CHGB1(I) = 1.D-40                                           
          END IF                                                        
        END IF                                                          
   20 CONTINUE                                                          
      S = 0.D0                                                          
      DO 1010 I=1,KSUM                                                  
        S = S+CHGB1(I)                                                    
 1010 CONTINUE                                                          
      S=S*VCEL - SCHGPC                                                 
      IF (MOD(ITER,10).EQ.0) THEN                                       
          WRITE (6,*) 'REAL TOTAL CHARGE = ',S,' IN XCFFT'              
      END IF

この問題の解決は、以下のようにルーチンを変更することによって達成されます。


      ICCCC=0
      DO 20 I=1,KSUM                                                    
        IF( CHGB1(I).LE.0.0D0 ) THEN                                      
              IF ( CHGB1(I).LE.-1.0D-5 ) THEN                            
                WRITE(6,610) I,CHGB1(I)                                  
  610           FORMAT(1H ,'**WARNING CHG.DEN<0.0 AT ',                
     &                                  I5,2D15.7,'***')               
              ICCCC=ICCCC+1
              END IF                                                   
              CHGB1(I) = 1.D-40                                           
        END IF                                                          

        IF ( MOD(I,IFX2).EQ.0 ) THEN
              CHGB1(I) = 1.D-40                                           
        END IF                                                        

   20 CONTINUE                                                          
      S = 0.D0                                                          
      DO 1010 I=1,KSUM                                                  
        S = S+CHGB1(I)                                                    
 1010 CONTINUE                                                          
      S=S*VCEL - SCHGPC
      IF (MOD(ITER,10).EQ.0) THEN                                       
          WRITE (6,*) 'REAL TOTAL CHARGE = ',S,' IN XCFFT'              
      END IF

以上にあるように、2n+1であることを判定するIF文（深紅色部分）を独立させてやれば、総電子数の値に現れるゴミ数値の問題は起こらなくなります。
最後に、今回のバグでは総電子数に問題がありましたが、全エネルギー等に問題はありませんでした。おそらくこれは、総電子数は逆格子空間の電荷密度 ρ(G)を逆フーリエ変換(Inverse FFT = IFFT)して、ρ(G) → IFFT → ρ(r) を求めました。ρ(r)は実空間での電荷密度（バンド指標、k点指標などは省略）。全エネルギー等の値は、畳み込みなどの過程を経た項が元となっています。この場合、先の逆フーリエ変換したものを更にフーリエ変換して元の逆格子空間に戻しています。これにより影響を与えた2n-1部分（実空間）が、再び2n-1部分のみに集約（逆空間）されて結果に影響を与えなかったのではないかと考えています。ここら辺は、完全には確かめ切れていません（↓追記も参照↓）。

（２／９、２００９追記）変数SCHGPCにも問題があることが判明しました。 SCHGPCは、部分内殻補正を使用する場合に必要となる変数ですが、これの初期化は部分内殻補正をするサブルーチン内のみで行なわれていました。部分内殻補正をしない場合でも上記ルーチンにあるように、総電荷数の足し上げの時に使用されています（この時初期化されていない）。状況によってはこの時、 SCHGPCに自動で初期化されず、不定の値が収まってしまう可能性があることが判明しました。プログラムの最初で、SCHGPC = 0.0D0と初期化して置くと、この問題は解消されます。これは同じプログラムでも、システムによって問題無かったり、少しプログラムが異なるだけで現れなかったりすることが分かっています（調査中）。ただ、SCHGPCのような変数は必ず初期化して置くことで、常に使用するのにある条件の時だけ初期化（或いは数値の代入）が行なわれるようにしないことです。

（Ｖ）本日のバグ（レポート２／８、２００８）[目次]

また間違った擬ポテンシャルを読み込んで計算していました。これは、これまでのバグレポート（Ｊ）、（Ｂ）と同じ失敗です。元素としては同じものながら、価電子の扱いが異なる擬ポテンシャルを誤って用いていました。これは比較的浅い内殻の電子を価電子として扱うものと、扱わない擬ポテンシャルがあり、内殻を考慮しない計算の場合に、内殻を価電子として考慮した擬ポテンシャルを内殻非考慮と思い込んで使用していました。全く過去の教訓が生かされていません。そもそも相変わらず、このような場合の誤りを事前、或いは少なくとも計算実行初期に確認（チェック）する機能が導入されていません。情けない話です。

（Ｕ）本日のバグ（レポート１２／１７、２００６） [目次]

今回のバグは、非常に単純なものです。バグ（と言うか問題）のあったプログラムは、バンド計算プログラム本体ではなく、バンド計算結果を使って解析を行なうごく簡単なプログラム（www.nims.go.jp/cmscへジャンプ、バンドギャップのある系の、直接、間接を判定する）です。このプログラム自身は一応（？）問題無いのですが、重大な制限があります。それはユニットセル内の総電子数が１６個（h-BN はボロン、窒素が各２個ずつで構成される）でないと正しく動作しません。勿論、系が半導体（絶縁体、つまりギャップが存在する系）であることも必要です。またスピンは縮退しているとして、１バンドに２個電子が詰まることを想定しています。

ところが筆者はこれを、ボロン、窒素がユニットセルに各４個ずつ存在する系（総電子数３２個＝１６バンド）で、そのまま使用してしまいました。プログラムを見れば分かるように、この場合完全に詰まっている８番目と９番目のバンドが計算の対象となり、これによって得られた判定結果は全く意味のないものになります。バンドギャップの値も負の値という変なものになってしまいます（←本当はここで最低限、エラーを返す〔警告する〕ようにしなければならない）。最初、筆者はこの問題に気付かず、暫く経ってから結果のおかしなところ（先の負のバンドギャップ値）に気付き、そこからプログラムの不備に辿り着きました。対策は、応急的ですが判定すべきバンドを８番目、９番目から、１６番目（ＶＢＭがあるバンド）、１７番目（ＣＢＭがあるバンド）に変えることで正しい結果が出ました。ＶＢＭ：価電子帯の頂上。ＣＢＭ：伝導帯の底。

この問題はプログラムを最初に作った段階で認識していたのですが、将来への対策を全く施して（考えて）いませんでした。で、案の定、問題を忘れてそのままでは計算していけない系に適用し、挙げ句に間違いに気付くのに暫く時間を要するという失敗の連鎖を行なってしまいました。

（１２／２７、２００６）間違いは別にまだあった。電子数が増える方ではなく、減る方もあった。つまり、ユニットセル内の全電子数が１６→８に減る場合もあった。この時、系がギャップを持つなら、電子に占有されたバンド数は４であり、非占有の一番下のバンド（空のバンド）は下から５番目ということになります。筆者は、これも見逃して占有８バンド、最低非占有９バンド目として計算していました。勿論、正しい結果は得られていませんでした（バンドギャップも負の値という変な値となっていた。←ここで気付くべきだった）。

（Ｔ）本日のバグ（レポート８／２１、２００６） [目次]

今回のバグは既に過去に言及したものと同じものです。そのバグは、（Ｋ）です。古いバージョンのプログラムを動かしたところ、計算されたストレスの値に、同じ条件での計算結果にかかわらず一致しない場合が出てきました。最初、原因不明でしたが、新しいバージョンのプログラムと、問題のる当該プログラムの差分を取り、詳しく調べたところ問題点が判明しました。その内容は（Ｋ）でのものと全く同じでした。つまり、サブルーチンFORZFBがストレス計算に対応していませんでした。これを修正することにより問題は解決しました（←完全に一致することを確認）。
プログラムのバージョン管理をちゃんと行なう必要があります。古いバージョンのプログラムでも、時としてそれを使用して計算結果を得なければならない事態もあり得るので、古いバージョンのプログラムに対する、バグ等への対応を怠らずに行なう必要があります。
いろいろな事情で、古いバージョンのプログラムを使用せざるを得ない状況は存在し、当該プログラムより後のバージョンでは対応済みのバグが、そのプログラムに残っている可能性があります。そのような場合の対応は、使用者、作成者本人もそのことを忘れる、或いは認識していない可能性もあり、悩ましい問題の一つです。出来るだけ、バグの記録を詳細に残し、どのバージョン、どの（バージョンの）プログラムまで対応したか、きちんと管理する必要があります（まあ大抵それが出来ていないことが多い）。

（Ｓ）本日のバグ（レポート９／１６、２００５） [目次]

本日の失敗（→教訓）は、結晶対称性に関してのものです。今回扱った系は、 C₆B₂系が主なものです。これは仮想物質なのですが、元々計算を始めた、Na₃As構造のC₆B₂以外の仮想構造の計算で問題が発生しました。C₆B₂ （Na₃As構造）については、第１７回化合物新磁性材料研究会[講演 ]参照。

C₆B₂（Na₃As構造）の結晶対称性は、 P6₃/mmcとhcp構造と同じ高い対称性を持ちます。しかし、この構造は完全な仮想構造で、（凝集）エネルギー的にも安定でないことが計算により分かっています。そこで、これと類似の仮想構造の候補を考え、それらの安定構造と電子状態を求めることとしました。
最近のコードは大変良く出来ていて、与えた系（ここでは C₆B₂系仮想物質）の結晶対称性を自動でも求めてくれたりします。これまで扱っていた、C₆B₂ （Na₃As構造）と類似の新たな仮想C₆B₂について、その結晶対称性を自動判定したところ、その結晶対称性では筆者のプログラムで計算出来ない（結晶対称性の判定テスト部分で止まる）ことが判明しました。当初、この原因が全く分からず、いろいろなことを試しましたが全てうまく行きませんでした。

C₆B₂（Na₃As構造）の原子座標は、筆者の計算プログラムでは、


       25.9099482945        0.0000000000        0.0000000000
        0.0000000000        4.7796946433        0.0000000000
        0.0000000000        2.3898473216        4.1393369833
           0  COORDINATES  0:NORMALIZED  1:CARTESIAN 
        0.2500000000        0.0000000000        0.0000000000
        0.7500000000        0.0000000000        0.0000000000
        0.0000557124        0.3333333333        0.3333333333
       -0.0000557124        0.6666666667        0.6666666667
        0.4999442876        0.3333333333        0.3333333333
        0.5000557124        0.6666666667        0.6666666667
        0.2500000000        0.3333333333        0.3333333333
        0.7500000000        0.6666666667        0.6666666667

となります。上記は既に構造最適化を終了しています。" 0 COORDINATES 0:NORMALIZED 1:CARTESIAN "の行以降が原子座標で、最初の６行が炭素（６個）、最後の７、８行がホウ素（２個）です（以下同じ）。" 0 COORDINATES 0:NORMALIZED 1:CARTESIAN "行の上にある数値は、格子座標です。座標は、筆者の計算では(z,x,y)となっています。従って、一番左側の数値がc軸（＝z軸）のものです。また、"0.6666666667 0.6666666667"というような表現も一般的でありません。通常は、"0.6666666667 0.3333333333"などと表現します。

まず以下の仮想構造での原子座標、


        0.2500000000        0.0000000000        0.0000000000
        0.2500000000        0.3333333333        0.3333333333
        0.7500000000        0.0000000000        0.0000000000
        0.7500000000        0.6666666667        0.6666666667
        0.0000000000        0.3333333333        0.3333333333
        0.5000000000        0.3333333333        0.3333333333
        0.0000000000        0.6666666667        0.6666666667
        0.5000000000        0.6666666667        0.6666666667

を考え、あるコード（←筆者のものでない）で結晶対称性を判定をしたところ、１８７番（P6m2←本当は数字６の上にバーが付く）という解答でした。ところが、筆者のバンド計算プログラムコードで１８７番の結晶対称性として計算させようとすると、結晶対称性判定部分で止まってしまいました。筆者のコードは結晶対称性の処理に不完全なところがあり、他の対称性では上記座標でも計算が通ってしまうのですが、この場合、計算結果が明らかにおかしなものとなります。この問題に際して、当初は結晶対称性に関してのオペレーター（演算子）の設定が筆者のコード上で間違っているのではないかと、検証、検討を試みました。とにかく、判定を通るように強引な設定で計算を走らせてみたりとか、１８７番の対称性は、P6/mmm対称性からインバージョン（反転）を除いたものなので、P6/mmmで強引に計算させてみたりしましたが、どれもうまく行きませんでした。

そこで今度は、同僚の新井さんに相談、別の計算コードで結晶対称性の判定を行なってもらいました。筆者は、この構造は１８７番の対称性ではないのではとも考えていたのですが、このコードでも判定結果は１８７番でした、ただ原子座標の設定（原点の取り方）がちょっと上記のものと異なるという結果が出てきました。つまり以下のような座標の移動（変換）が必要であることが分かりました。


      (1) 元の原子座標
        0.2500000000        0.0000000000        0.0000000000
        0.2500000000        0.3333333333        0.3333333333
        0.7500000000        0.0000000000        0.0000000000
        0.7500000000        0.6666666667        0.6666666667
        0.0000000000        0.3333333333        0.3333333333
        0.5000000000        0.3333333333        0.3333333333
        0.0000000000        0.6666666667        0.6666666667
        0.5000000000        0.6666666667        0.6666666667

     (2) x,y座標を、0.3333333333だけずらす。
        0.2500000000        0.3333333333        0.3333333333
        0.2500000000        0.6666666667        0.6666666667
        0.7500000000        0.3333333333        0.3333333333
        0.7500000000        0.0000000000        0.0000000000
        0.0000000000        0.6666666667        0.6666666667
        0.5000000000        0.6666666667        0.6666666667
        0.0000000000        0.0000000000        0.0000000000
        0.5000000000        0.0000000000        0.0000000000

     (3) z座標を、0.25だけずらす。
        0.5000000000        0.3333333333        0.3333333333
        0.5000000000        0.6666666667        0.6666666667
        0.0000000000        0.3333333333        0.3333333333
        0.0000000000        0.0000000000        0.0000000000
        0.2500000000        0.6666666667        0.6666666667
        0.7500000000        0.6666666667        0.6666666667
        0.2500000000        0.0000000000        0.0000000000
        0.7500000000        0.0000000000        0.0000000000

     (4) 単なる並べ替え。
        0.0000000000        0.3333333333        0.3333333333
        0.0000000000        0.0000000000        0.0000000000
        0.5000000000        0.3333333333        0.3333333333
        0.5000000000        0.6666666667        0.6666666667
        0.2500000000        0.6666666667        0.6666666667
        0.7500000000        0.6666666667        0.6666666667
        0.2500000000        0.0000000000        0.0000000000
        0.7500000000        0.0000000000        0.0000000000

     (5) 単なる並べ替え。
        0.0000000000        0.0000000000        0.0000000000
        0.0000000000        0.3333333333        0.3333333333
        0.5000000000        0.3333333333        0.3333333333
        0.5000000000        0.6666666667        0.6666666667
        0.2500000000        0.6666666667        0.6666666667
        0.7500000000        0.6666666667        0.6666666667
        0.2500000000        0.0000000000        0.0000000000
        0.7500000000        0.0000000000        0.0000000000

上記の最後の原子座標と、最初の原子座標は等価なものです。ただ、最後のもの(5)では、座標(0,0,0)に炭素原子が置かれているようになっています。対称性に関してのオペレーターは、対象となる原子座標の配置による影響を受けるため、正しい配置にしておかないと、判定に失敗したり、正しい計算が出来ないということがこれで分かりました（←もっと早く認識すべきだった）。
上記最後の原子座標は、新井さんが使用したコード上で示された（変換された）座標で、筆者のコードでもこの座標では、１８７番の結晶対称性で正しく計算を遂行することが出来ました。最初に筆者の使用したコード（←筆者独自のコードでない）では、おそらく座標の変換を内部で自動的に行なっているものと推定されます。コードに慣れていないため、その過程を出力結果から見い出せていない可能性があります。

今回、異なるコード（手法、手段）を使用した結果の比較によって、問題を解決することが出来ました。

【参考】対称性、群論関連の[サイト]

（Ｒ）本日のバグ（レポート２／１１、２００５） [目次]

本日の失敗は、ごくごく単純です。しかし、だからこそ重要。六方晶構造の単位胞での基本格子ベクトルを与える時、a,b,c各軸とそれらの各軸同士がなす各度α、β、γにおいてa,b軸のなす各度（γ）を、筆者は６０度としていました。ここではc軸をz方向とし、a軸、c軸及びb軸、c軸のなす各度はいずれも９０度とします。あるソフトで、γを６０度として、六方晶格子の格子データを入力して計算させると、どうも計算結果がおかしな値になり、困っていました。教科書（結晶格子について載っている固体物理関係のもの）をちゃんと見ると、γは１２０度となっていました。で、γを１２０度として当該ソフトで計算すると、確かに妥当な結果が得られました。
この場合、γを６０度として格子を作って話を進めても、それ自身は間違いではないです。少なくとも、それに対応するようにソフトウェアコードが動くようにさせていれば問題ありません。自作のものならそれでＯＫなのですが、一般のもの（フリーソフトなり市販品なり）ではそうとも行きません。γを６０度とするのは一般的ではなくほとんどの場合、慣例としてγは１２０度として格子パラメータは設定されます。筆者はこのことを十分認識していませんでした（←筆者の勉強不足が問題）。
因みに、このソフトはWyckoffパラメーターを入力すれば良いようになっています。この時、筆者は初めて、Wyckoffパラメータの存在と意味を認識しました。結晶構造データベース[Crystal Lattice Structures Database](The Naval Research Laboratory)には、ちゃんと例示されている結晶構造のWyckoffパラメーターが表示してあります。筆者はこのサイトを大変重用していて改めて良くできていると思いました。

（Ｑ）本日のバグ（レポート９／２６、２００５） [目次]

何と、本ページ自身に失敗があった。（Ｑ）レポートを飛ばして（Ｐ）の次を（Ｒ）としていた。アルファベットの順番のチェックを怠っていた。悲しい。

（Ｐ）本日のバグ（レポート１／１２、２００５） [目次]

今回のバグ（と言うより失敗）は２つあります。一つはバンド計算時の座標入力上のもので、今一つは擬ポテンシャル作成時の入力データ間違い（ミス）です。いずれも比較的（または非常に）単純なものす。１／１２以降、大変多忙な日が続き、実際のアップロードは、２／７（２００５）となりました。

（１）座標データの間違い。これは六方晶構造の物質の計算において、ユニットセルではなく、スーパーセルとして計算しようとした時に問題が発生しました。ユニットセルの規模の計算では、その系の持つ対称性から、ユニットセル内の原子が受ける力は（多くの場合）ゼロでした。←六方晶と言っても、ウルツ鉱型構造のように内部パラメーターを持つ場合、原子が力を受ける場合があります。今回の場合、ユニットセルの規模では原子は力を受けないような系だったのですが、それをスーパーセルとして計算を行なうと、想定していない力を原子が受ける事態に遭遇しました。これはユニットセルとして計算できる系をユニットセル数個分（例：１×１→ ２×２）のスーパーセルと単純に拡張しただけで、スーパーセル内の原子はユニットセルでの場合と同様、原子は力受けないはずにも関わらず、実際の計算では力が出てきてしまいました。
当初、筆者は、このようなことはキュービック（立方体。或はもっと緩めに四角→直方体）の場合ではこのような力が出てくることはなかったので、六方晶に特有の問題（例：ｋ点の問題とか）と全く根拠のない推定のもとに問題を放っていました。しかし、それは全くの事実誤認でした。最近、プログラムの出力に関して若干の改良を施し、それによるテスト計算を上記のスーパーセルからなる系に対し行なってみた結果、明らかにおかしな結果が出てきました。それは、改良した計算ではユニットセル（スーパーセル）内の各原子間の距離を計算し、大きさ順に並べる結果を出力させるようにしたのですが、それが明らかに予想した結果と異なるものでした。これは前述の力の問題以上に深刻で、あり得ないことでした。これらの事実からより詳細な検討を行なった結果、筆者のデータ入力の仕方が誤っていたことが判明しました。データの入力を正しいものにすると、あり得ない力は消え、原子間距離などに対し妥当な結果が得られました。今回、プログラム的には、変更箇所は全く無く、データの入力の仕方を変えるだけで事足りるものでした。あり得ないことやおかしな結果が出てきたら、希望的、楽観的な解釈（予想）から問題が解決した（している）と考えず、徹底的に検討するべきだと改めて痛感しました。

（２）今一つの失敗は、擬ポテンシャル作成時のもので、明らかに存在しないエネルギー準位を指定して擬ポテンシャル(Os)を作成していました。この場合、 Os（オスミウム）は、１ｓ～６ｐまで１５の準位が存在するとして計算すべきものでした。ここで、６ｐは空の準位として考えます（←空にするのは、基底状態のs、dを考える時で、ｐを扱う時は、６ｐをイオン化させて空でない準位として扱う）。ところが実際は、この１５番目の６ｐを１６番目と指定して擬ポテンシャル作成していました。この状態で作成したOsは、バンド構造としてはそう変なところはなかったのですが、実験値と比べて平衡格子定数が非常に伸びた値になってしまいました。非常に長い間、この間違いに気付かなかったのですが、ごく最近入力データを良く見ていたら１５個までしか準位が存在しないのに１６番目の準位が指定されている間違いに気付きました。
早速、間違いを修正、作成し直したOs擬ポテンシャルで、平衡格子定数を求めると正しい妥当な結果が得られました（LDAによる計算ながら若干格子定数は伸び目となる）。そもそもあり得ない準位を指定しても計算が通り、もっともらしい結果が得られたということも問題の一つであります。これについても今後検討する必要があります。

（Ｏ）本日のバグ（レポート２／２５、２００４） [目次]

久々のバグです。これはバンド計算でよく使用される高速フーリエ変換に関するものです。最近は、LiBCや、 MgB₂などのボライド系の計算が多く、あまり遷移金属の計算は行なっていませんでした。そして、ある遷移金属の計算をしてみたところ、大変計算が遅いことに気付きました。遷移金属の場合、3d、4dなどのｄ電子を価電子として考慮する必要があり、その分計算が面倒になるのですが、今回の計算は通常のバルク系で、hcp構造（つまりユニットセル内に原子は２個）が計算対象であり、計算規模も先のLiBC、MgB₂と比べて、ほとんど変わらないものでした。
しかしながら、実際に計算してみるとめったやたら計算速度が遅いことが判明しました。遷移金属だから特殊な事情があるのかとも考えましたが、使用するプログラムは、LiBC、MgB₂の時と同じもので、計算条件や設定もさして変わらないので、そのような”特殊”な事情は存在しないと判断されます。

原因は気が付けば、今回も大変単純なものでした。単に、高速フーリエ変換で使用するパラメーターの設定の仕方が妥当なものでなかったからでした。高速フーリエ変換（ＦＦＴ）では、主に電荷密度などに対して実空間⇔逆格子空間の変換を行なうのですが、その時に使用する格子の数（空間のメッシュ数）は、２の倍数である場合が最も高速となります。そうでない場合でも２、３、５などからなる倍数で、１７とか２３のような値の倍数の選択は計算速度を非常に遅くするので推奨されません。ＦＦＴのプログラムによっては２の倍数或いは、これに準じる２、３、５の倍数でないものはそもそも計算できない設定となっているものもあります。
今回筆者が使用したＦＦＴルーチンは、２（２、３、５）以外でも計算可能なもので、３４（１７×２）を使っていました。これはＦＦＴ計算部分を大変遅くすることが判明しています。

（どうやって気付いたか）
たまたま、上記のＦＦＴルーチンが使用できないマシン上で、昔から使っているMFFT（２、３、５の倍数しか使えないＦＦＴルーチン）を使用したため、３４ではエラー表示が出たため気付きました。このことがなければ、今でも気付かないままでした。通常は、メッシュ数が２、３、５程度の倍数でない場合は、それに近いより大きい２、３、５位までの倍数になる値で設定していたのですが、今回はそのまま気付かずに３４（１７×２）のままで計算していました。この計算では、結局３２（２の５乗）の値も使えました（←本来なら３６だが、筆者は更に勘違いしていて３２でよいものをわざわざ３４で計算していた）。３４と３２でのバンド計算全体での速度差は、粗い評価ですが１０倍以上異なることが判明しています。

もう既に何度も書いていますが、計算がどこかおかしい時、たとえば（１）計算結果は正しいが動作がおかしい（例：計算速度がやたら遅いなど）、（２）動作速度などは問題ないが計算結果が微妙に（或いは何となく）おかしい、（３）計算結果も動作もどこか（微妙に）おかしい、などの場合、計算に何か問題があるはずで、放置せずに原因究明に努めることが必要と考えます。

（Ｎ）本日のバグ（レポート４／９、２００３）[目次]

プログラムというものは、これまでやったことがない新しいことをしようとすると、トラブルに見舞われることが良くあります（筆者の場合特に）。
今回のバグは、筆者の使っているバンド計算プログラムにおいて、新しい系を計算した時に起こりました。勿論、プログラムは、この初めて扱う系にも仕様上は正しく計算できるはずでした。
しかし、実際計算を行なうと、ストレスなどが正しく計算できていないことが判明しました。

では今回扱ったのはどのような系だったのか、物質名は具体的にはまだ挙げられないのですが（正直、そんなに秘密にするほどの系ではない）、ある３元系で、アルカリ金属（擬ポテンシャルにおいて、ｓ、ｐのみ非局所）、遷移金属（ｓ、ｐ、ｄとも非局所）、酸素（ｓ、ｐのみ非局所）からなる化合物（酸化物）でした。内殻補正は、アルカリ金属と遷移金属の擬ポテンシャルに対して施されていました（今回は、これは影響なし）。
問題となったのは、非局所の部分に対応する、局所ポテンシャル部分の扱いのところでした。これは、計算結果がおかしいことに気付いた後、検証のため、元素の座標の入力の順序を、アルカリ金属、遷移金属、酸素の順番を、遷移金属、アルカリ金属、酸素と変えると、計算結果のある部分に違いが存在することから判明しました（当然、本当は元素の座標の順序の入れ換えで結果に影響はないはず）。違いがあったのは、全エネルギーの値が順序入れ換えで異なり、更に細かく調べると、局所擬ポテンシャル部分に関わる値がおかいしことが判明しました。以下が、当該する部分のプログラムルーチンです、


      DO 500 II=1,KTYP                                                  
      ETOT1 = ETOT1 + FLOAT(IATOM(II))*PAI*ACHG(II)                     
     &      * ( AC(II,1)/BC(II,1) + AC(II,2)/BC(II,2) )/UNIVOL          
      TOTCH=TOTCH + FLOAT(IATOM(II))*ACHG(II)                           
C                                                                       
      IF (ITER.EQ.0) THEN                                               
        IF (NLSPD(II).EQ.1) THEN
          READ(16,*) MESHR,NMES,DX,RAD,VD,VDNL
          DO 1210 N=1,MESHR                                             
          VDD(N,II)=VD(N)                                               
          VDDNL(N,II)=VDNL(N)
 1210 CONTINUE                                                          
        ELSE
          DO 1212 N=1,MESHR                                             
          VDD(N,II)  =0.0D0
          VDDNL(N,II)=0.0D0
 1212 CONTINUE                                                          
        END IF          
      ELSE                                                              
          DO 1211 N=1,MESHR                                             
          VD(N)=VDD(N,II)                                               
          VDNL(N)=VDDNL(N,II)
 1211 CONTINUE                                                          
      END IF                                                            
C                                                                       
      S=0.0D0                                                           
      DO 1200 N=1,MESHR                                                 
      S=S + OMO(N)*VDNL(N)*(RAD(N)**3)
 1200 CONTINUE

色（深紅）のついた配列VDNLですが、実はこれが正しくゼロ初期化できないことがバグの原因でした。
３元系なので、DO 500ループ繰り返し数は３で、その都度、各元素に対応した局所ポテンシャル（VDNL）を読み込むのですが、VDNLは、ｓ、ｐ、ｄが非局所の場合は読み込まずゼロ初期化するようにしています（IF文による分岐）。しかし、実際に初期化しているのは、VDDNLで、VDNL自身ではありません。ところが、後の部分（S=S + OMO(N)*VDNL(N)*(RAD(N)**3)）でVDNLが使われています。この部分は、まだDO 500ループ内です。つまり、もし３元系の元素の順番が上記アルカリ金属、遷移金属、酸素なら、最初、アルカリ金属（ｓ、ｐ非局所）で、 VDNLが読み込まれ、次の遷移金属の時は、VDNLはゼロとして扱われるべきですが、上記ルーチンでは（上記で示していない部分を含めて）、VDNLは初期化されておらず、前のアルカリ金属の値がそのまま残っています。このため、Sが余計に計算され、局所部分の値がおかしくなります。

（解決法）
非常に簡単で、上記ループで、



          DO 1212 N=1,MESHR                                             
          VDD(N,II)  =0.0D0
          VDDNL(N,II)=0.0D0
          VDNL(N)    =0.0D0
 1212 CONTINUE

と初期化をきちんとすれば問題は解決します。

（何故、今まで気付かなかったか？）
例えば、LiBC（LiBCの研究ページ）や、MgB₂（MgB₂の研究ページ）の計算を筆者は行なっていて、LiBCも同じ３元系で、上記と同様な問題が生じてもおかしくなかったのですが、”幸運”にも、偶然このバグは発生しないようになっていました。それは、LiBCでは、Liのみが、ｓ、ｐ、ｄ非局所であり、この最初の元素がｓ、ｐ、ｄ非局所なら、VDNLの読み込みは行なわれず、局所ポテンシャル部分の計算結果はおかしくなりませんでした（MgB₂ も同様の理由）。
但し、これは、偶然筆者の使用するシステムが、暗黙の内に未代入、未使用の配列の値をゼロ初期化しておいてくれたためで、計算システムの仕様によっては、ゼロ初期化されずに、おかしな値となる可能性があります。

今回の教訓は、今までちゃんと正しく動いたプログラムでも、新しい系、条件の計算では、いきなりおかしな結果となることがあること、おかしな結果を与える原因（バグ）も決して複雑なものであるとは限らず、今回のように大変単純なものであることの方が、むしろ多いということです。しかし、これによりプログラムはより強力（Robust）で信頼性の高いものになっていく訳です。
構成元素の原子座標の順番を変えることで、バグを発見できたのは、バグ発見の手法、手順としては典型的と言えます（参考ページ）。

（Ｍ）本日のバグ（レポート１２／５、２００２） [目次]

実に１年ぶりのレポートです（ここ最近あんましコードを書いていなかった）。

今回のバグ・失敗も大変単純なものでした。

（バグ・失敗詳細）
ｋ点の座標を読み込むルーチンを作成して、そのｋ点座標を基にしてバンド計算を行なうようにプログラムを変更し、実行させたのですが、テスト段階で既に結果が分かっている系での、ｋ点座標を読み込ませて計算した結果が、正しい結果と一致しないという事態に遭遇しました。

（原因）
これは割と簡単に発見できました。ｋ点座標を読み込むためのプログラムでは、まず大元のデータから必要なｋ点座標を抽出し、それを掃き出すプログラムを作り、その掃き出したデータを、実際にバンド計算プログラムから読み込むようにさせました（この２段階の操作に問題があった←何でこのような面倒なことになるかについては、いろいろ事情があったりする）。
筆者の扱うバンド計算プログラムは、ベクトルマシン用のチューンの影響がまだ生きていて、なるべく３次元配列（例：Ａ（１０００、１００、１００））の最初の次元での配列の数を大きめにとるようになっています。この影響で、六方晶系を扱う場合、ｃ軸（普通はz座標〔ｃ軸方向の格子定数は大抵a,b軸より長い〕で、３次元配列なら、最後の次元〔つまりｚ座標〕）をｘ座標として扱うようにしていました。

このため第一段階のｋ点座標抽出・掃き出しプログラムにおいて、以下の様に元データ中から座標を読む（READ文）段階で、ｚ座標→ｘ座標に変換することとしました（変換する手としては他にも沢山ある：今回のREDA文での方法を推奨している訳ではない）。


      READ(20,*) KVZ(NNN),KVY(NNN),KVX(NNN),KQWGT(NNN)

本来は、KVX(NNN),KVY(NNN),KVZ(NNN)というｋ点座標の並びを、READ文で読み込む時の配列の並びを、KVZ(NNN),KVY(NNN),KVX(NNN)として、KVX(NNN)→ KVZ(NNN)、KVZ(NNN)→KVX(NNN)となるようにしました。
これ一回で座標の変換は終了なのですが、筆者はやってしまいました。この掃き出したデータを、バンド計算プログラム側で読み込むルーチンのREAD文でも、再び、


      READ(61,*) KVZ(NNN),KVY(NNN),KVX(NNN),KQWGT(NNN)

としてしまいました。これは変換を２度行なっている（つまり、ｘ→ｚ→ｘ）ことであり、元に戻していることに他なりません。
このため計算結果が、一致しなかった（元の正しい結果は、上記第一段階での変換した座標での計算結果）のでした。上記のバンド計算側のREAD文を、


      READ(61,*) KVX(NNN),KVY(NNN),KVZ(NNN),KQWGT(NNN)

とすると、比較すべき正しい結果と、計算結果の一致を見ることができました。
（どうやって気付いたか）エラーの原因は、バンド計算プログラム側で読み込んだｋ点座標をWRITE文で確認してみると、最初の元データの座標と同じ並びに戻っていることから気付きました。

（補足）
尚、上記計算において、六方晶の格子座標の取り方の違いによる座標変換が必要でした。上記バグとは関係ないのですが、参考までにその手続きを以下に示します。
筆者の計算プログラムでは、実空間での六方晶の格子を、


       13.1338124721        0.0000000000        0.0000000000
        0.0000000000        5.1692358601        0.0000000000
        0.0000000000        2.5846179299        4.4766895732

のように与えています。これに対応する逆空間での格子座標は、


   0.478398       0.000000       0.000000
   0.000000       1.215496      -0.701767
   0.000000       0.000000       1.403534

です。バグの説明のところでも述べましたが、x座標がc軸となっています。ここでy、z座標の取り方に、上記とは異なる取り方があります。それは実空間格子座標では、


       13.1338124721        0.0000000000        0.0000000000
        0.0000000000        2.5846179299       -4.4766895732
        0.0000000000        2.5846179299        4.4766895732

という取り方です。これに対応する逆空間格子座標は、


   0.478398       0.000000       0.000000
   0.000000       1.215496      -0.701767
   0.000000       1.215496       0.701767

となります。基のｋ点の取り方が、整数表示になっていたので、上記逆空間格子座標を使って、目的のｋ点座標を求めることが出来ます。以下に、座標を変換している部分を示します。


C
      READ(60,*) KVX(NNN),KVY(NNN),KVZ(NNN),KQWGT(NNN)  ←読み込み（整数）
C
      KQSUM = KQSUM + KQWGT(NNN)      
C
      VX(NNN) = DFLOAT(KVX(NNN))/24.0D0  ←実数化
      VY(NNN) = DFLOAT(KVY(NNN))/24.0D0
      VZ(NNN) = DFLOAT(KVZ(NNN))/24.0D0
      QWGT(NNN) = DFLOAT(KQWGT(NNN))/DFLOAT(KKX*KKY*KKZ)
C
C           R22 = RLTV(2,2)
C           R23 = RLTV(2,3)
C           R32 = RLTV(3,2)
C           R33 = RLTV(3,3)
           R22 =  RLTV(2,2)  ←座標変換１
           R23 =  RLTV(2,2)
           R32 =  RLTV(3,2)
           R33 = -RLTV(3,2)
           WRITE (6,*) "R22 = ",RLTV(2,2)
           WRITE (6,*) "R23 = ",RLTV(2,3)
           WRITE (6,*) "R32 = ",RLTV(3,2)
           WRITE (6,*) "R33 = ",RLTV(3,3)
C
            Q(1) = VX(NNN) ←座標変換２
            Q(2) = VY(NNN)                                          
            Q(3) = VZ(NNN)                                          
            QX=RLTV(1,1)*Q(1)+RLTV(1,2)*Q(2)+RLTV(1,3)*Q(3)      
            QY=RLTV(2,1)*Q(1)+R22*Q(2)+R23*Q(3)      
            QZ=RLTV(3,1)*Q(1)+R32*Q(2)+R33*Q(3)      
            VX(NNN) = QX
            VY(NNN) = QY
            VZ(NNN) = QZ
C
          QSUM = QSUM + QWGT(NNN)      
C
C      WRITE(6,*) KNUM,VX(NNN),VY(NNN),VZ(NNN),QWGT(NNN)
C     WRITE(6,402) KNUM+1,VX(NNN),VY(NNN),VZ(NNN),QWGT(NNN)
      WRITE(6,*) KVX(NNN),KVY(NNN),KVZ(NNN),KQWGT(NNN)
      WRITE(6,402) KNUM,VX(NNN),VY(NNN),VZ(NNN),QWGT(NNN)
  402 FORMAT((I4,4(F15.9,2X)))
  400 CONTINUE

上記ルーチンは、あまりというかほとんど洗練されていません（２／１３、２００３）。

（Ｌ）本日のバグ（レポート１２／１９、２００１） [目次]

久々に間違ってしまいました。相変わらず今回のバグも単純なものです。

プログラムに新しい対称性（Wurtzite型構造）に対応するよう改良を施す過程で生じた間違いです。筆者の計算では、対称性は次のような番号で指定されていました。


C-----KBZTYP=1 : WHOLE    B.Z.  ,   2 : SIMPLE CUBIC IR.B.Z.            
C-----       3 : BCC   IR.B.Z.  ,   4 : FCC          IR.B.Z.            
C-----       5 : DIAMOND STRUCTURE IR.B.Z.                              
C-----       6 : HEXAGONAL STRUCTURE IR.B.Z.                            
C-----       7 : P6/MMM(ALB2) STRUCTURE IR.B.Z.
C-----       8 : TETRAGONAL STRUCTURE IR.B.Z FOR BETA-TIN               
C-----       9 : A-PBO2 STRUCTURE IR.B.Z FOR SIO2                       
C-----      10 : TETRAGONAL STRUCTURE IR.B.Z FOR RUTILE                 
C-----      11 : ORTHOROMBIC, CACL2 TYPE                                
C-----      12 : FCO   IR.B.Z.
C-----      13 : PERPVSKITE(ORTHO)  IR.B.Z.

これに新たに、Wurtzite型構造（六方晶型、P6₃mc）を加えるにあたってちょっと困りました。もう既に0以上の番号（番号0は対称性を考慮しない表面系用）が埋まっていて、HCP構造のところに、この新たな対称性のための番号を配することができなかったのです。そこで仕方なく以下のような番号付けをしました。


C-----KBZTYP=1 : WHOLE    B.Z.  ,   2 : SIMPLE CUBIC IR.B.Z.            
C-----       3 : BCC   IR.B.Z.  ,   4 : FCC          IR.B.Z.            
C-----       5 : DIAMOND STRUCTURE IR.B.Z.                              
C-----       6 : HEXAGONAL STRUCTURE IR.B.Z.                            
C-----      -6 : WURTZITE STRUCTURE IR.B.Z.
C-----       7 : P6/MMM(ALB2) STRUCTURE IR.B.Z.
C-----       8 : TETRAGONAL STRUCTURE IR.B.Z FOR BETA-TIN               
C-----       9 : A-PBO2 STRUCTURE IR.B.Z FOR SIO2                       
C-----      10 : TETRAGONAL STRUCTURE IR.B.Z FOR RUTILE                 
C-----      11 : ORTHOROMBIC, CACL2 TYPE                                
C-----      12 : FCO   IR.B.Z.
C-----      13 : PERPVSKITE(ORTHO)  IR.B.Z.

何と、”-6”という負の番号を強引に定義したのです。間違い（バグ）の原因は、番号が負の整数値だったことです。プログラム内で、IF判定で各対称性に対応した計算を行なっている時、これまでは対称性を考慮しない場合（番号、 0,1）と、対称性を考慮する場合（番号2以上）というようになっていたので、 IF判定で、IF (KBZTYP.GE.2)或いは、IF (KBZTYP.LE.1)という部分がプログラム内に何箇所かありました。筆者は、この部分への対応を忘れてしまったのです。
また、この番号を担うプログラム内の変数が、”KBZTYP”だけでなく、もう一つ別に、”NBZTYP”というものがあり、筆者はこの”NBZTYP”に関する部分の変更にのみ注目してしまったことが、”KBZTYP”絡みの部分への注意を失わせることとなりました。

当然の如く、計算結果に問題が生じました。計算して得られる力や、ストレスの値で、明らかに対称性からゼロにならねばならないものに有限の（無視できない）値が出てきてしまったのです。そして、その原因が暫くの間全く分かりませんでした。
いろいろ試行している間、与えた結晶内原子の座標の設定に問題があるのではとも考えたのですが、このWurtzite型構造による系を対称性を考慮しない条件で計算した場合では、正しく力、ストレスの値が求まる（この場合の”正しい”は、ゼロになるところは、ちゃんとゼロになるという意味）ことが判明し、結晶内の原子位置の設定の問題ではないということとなりました。
バグを見つけた発端は、対称性が正しく計算されている場合と今回の計算の結果との比較からでした。それは、対称性が正しく計算されている場合のストレスの結果、


 TOTAL SUMM   1 =  -3.391163378869085E-004
              2 =   1.273323054860103E-003
              3 =   1.716355131468628E-004
              4 =  -4.301059916385391E-004
              5 =   1.596389500226893E-004
              6 =   4.345718045303151E-004
 TOTAL SUMM OP1 =  -3.391163378869083E-004
              2 =   1.129377263005734E-021
              3 =  -1.411721578757167E-020
              4 =   2.232906445887955E-006
              5 =  -1.129377263005734E-021
              6 =   2.232906445888011E-006

と、今回の場合の同様の結果、


 TOTAL SUMM   1 =   1.523900247301507E-004
              2 =   1.641143043416956E-005
              3 =   1.089318435141048E-005
              4 =   2.612009135650878E-004
              5 =   1.295694267740399E-004
              6 =   2.076013932175158E-004

との比較でした。見てお分かりの通り、後者では、”TOTAL SUMM OP”の計算結果部分が抜けています。この最後の部分の結果は、それまで得られたストレスの値に対称オペレーターをかけて、対称性をきちんと考慮した結果を与えるもので、この結果の部分が表示されていないということは、最後に行なわれる対称オペレーターをかける計算が行なわれていないことを意味します。

調べてみると、当該部分の対称性考慮、非考慮の判定が、”IF (NBZTYP.LE.1) GO TO 9000”となっていました。これはこれまでの計算なら、番号0,1の場合（対称性を考慮しない）は、対称オペレータをかけずに、文番号９０００に飛ぶというのでよかったのですが、今回のWurtzite型構造に対応する番号は、-6で、これもこのままでは対称性を考慮しない計算として話が進んで（文番号９０００に飛んで）しまっていた訳です。前述のように、変数” KBZTYP”に関しての判定部分のチェックを全く見過ごしていたことが間違いのもとでありました。他の箇所も含め、変数”KBZTYP”の関わる部分（電荷密度、力の計算部分）の修正を行なうことで、力、ストレスとも妥当な結果を与えるようになりました。

今回の問題は、

（１）事前のプログラム設計の不備。番号が割り付けられず、負の整数値を導入したところで駄目駄目と言えます。
（２）NBZTYP、KBZTYPという事実上同じ機能を持った変数が独立に存在して、プログラム内で使用されている（プログラム改良時に一方が忘れられる可能性があり、バグの元になる）。

と言えます。

（Ｋ）本日のバグ（レポート７／１６、２００１） [目次]

今回のは単純で大きな間違い（バグ）でありました。

問題はストレスの計算に関しての部分で、症状は基本的に同じ計算を行なっているはずなのに計算値が異なるというものです。まず、その結果を以下に示してみます。

（ケース１）


 TOTAL SUMM OP1 =   8.363462325643757E-004 <-- ここ
              2 =   1.129377263005734E-021
              3 =  -1.129377263005734E-021
              4 =  -5.152418845876215E-004 <-- ここ
              5 =   1.044673968280304E-020
              6 =  -5.152418845876215E-004 <-- ここ
（中略）
 TOTAL ENERGY FOR 401-TH ITERATION= -51.5028059      -0.5150281D+02
 MIXING RATE =   0.500000000000000     
  ITER=401 ET(H)=  0.1493166D-06 ET(M)=  0.004061 DC=  0.3162807D-06

（ケース２）


 TOTAL SUMM OP1 =  -3.459023787458713E-003 <-- ここ
              2 =   0.000000000000000E+000
              3 =   1.807003620809174E-020
              4 =   1.641161644709054E-003 <-- ここ
              5 =  -2.258754526011467E-020
              6 =   1.641161644709055E-003 <-- ここ
（中略）
 TOTAL ENERGY FOR 401-TH ITERATION= -51.5028059      -0.5150281D+02
 MIXING RATE =   0.500000000000000     
  ITER=401 ET(H)=  0.1493166D-06 ET(M)=  0.004061 DC=  0.3162807D-06

ケース１とケース２で、”ここ”で示した値（ユニットセルにかかる圧力〔ストレス〕）の値が異なっています。一方、全エネルギー（TOTAL ENERGY）の値は一致しています。上記には示されていませんが、力の値も一致します。値が異なるのはストレス部分のみです。上記２つの場合での計算条件の違いは、ケース１では電子状態のみの計算を行ない、最後にストレスの値を計算し、それを表示させる。ケース２は、セルに関しての構造最適化を行なう過程で、最初に電子状態計算を行なわせ、セルの構造最適化を始める直前でのストレスの値を計算表示させるものでした。いずれも上記の場合では、４０１回目（イタレーション）での値です。

この状況から見ると、原因追求は簡単かとも思えたのですが、バグ発見当日には原因を突き止められませんでした。
状況として、（１）値が異なるのはストレスのみ、（２）全エネルギーなど他の値に違いが無い、（３）計算条件における違いは、セルに関しての構造最適化をさせるか、させないかだけの差、などが挙げられます。特に（３）から、バグが潜んでいそうなのは、計算条件に関わる部分の可能性が高く、比較のための（小規模な）テスト計算を繰り返したのですが、ストレスに関しての値が異なることがますます明白になる一方で、どうして値が異なってしまうのかさっぱり分からないという状態が翌々日まで続きました。

原因究明のきっかけとなったのは、ケース１、ケース２の計算出力結果を、 diffにかけて違いを比較検討したことでした。このこと自体は、バグ発見当日にも行なったのですが、その時点では問題点を見い出せませんでした。翌々日になって改めて見直してみると、違いが存在することに気付きました。それは、以下の部分です。


3c3
<   0.500000000000000       0.000000000000000E+000  0.000000000000000E+000
---
>   0.500000000000000       0.000000000000000E+000   10.0000000000000     
15c15
<          400         400         400           1         401           0
---
>          400        2000          40           1        2000           0
45,46c45,46
<  IOVE      =          400
<  IMDI      =          400
---
>  IOVE      =           40
>  IMDI      =         2000
50,51c50,51
<  ISTOP     =          401
<  ISTMD     =          401
---
>  ISTOP     =         2000
>  ISTMD     =         2000
72c72
<  DTIMUC    =  0.000000000000000E+000
---
>  DTIMUC    =   10.0000000000000     
2855c2855
<  CALL FORCE
---
>  CALL FORZFB
2901,2906c2901,2906 （↓ストレスの計算途中の値）
<  SIGNL     4  1 =  -0.275069507467294     
<               2 =   1.189155365786973E-003
<               3 =   2.122295586945811E-004
<               4 =  -0.274503255358969     
<               5 =   1.073767340287859E-003
<               6 =  -0.275494497853064     
---
>  SIGNL     4  1 =  -0.279364877487318     
>               2 =   3.965559338907192E-003
>               3 =   2.529826782671025E-003
>               4 =  -0.268815066930155     
>               5 =  -3.305458596529087E-004
>               6 =  -0.276869879223285     
（以下略）

上記（diffコマンドによる差の表示結果）で、（主に前半部分の）計算条件の違いによる違いは別として、二つのことが指摘できます。
（い）”CALL FORCE”と”CALL FORZFB”の違い。（ろ）ストレス計算途中の値の違い。まず、（い）からサブルーチンの呼び出しの仕方が異なることが判明します。（ろ）に関しては、diffをかける前から、ストレスの計算途中の各段階での値のどの部分に違いがあるか調べた段階でこのことは判明はしていました。ストレスに関しては、その計算途中の値全てが異なっていた訳ではなく、ある特定のルーチンでの結果（値）がおかしいことは分かっていたのですが、どうしてこれが間違いに繋がるのか不明でした。それを解決に導いたのが、（い）のサブルーチン呼び出しの違いでした。

サブルーチンFORCE は、非局所擬ポテンシャル部分の計算を行なうルーチンで、このルーチンでは原子に働く力において、非局所擬ポテンシャル部分が寄与する部分の計算を行なっています。一方、サブルーチンFORZFB は、サブルーチンFORCEで行なう計算の内、力に関しての部分を除いた計算を行ないます。これは電子状態のみを求める計算では、力の計算が必要ではないため、この場合、より簡便で計算時間を消費しないルーチンに切替えて計算した方が良いためです。

ここで問題だったのが、このサブルーチンFORZFBは”電子状態”計算のみ用でしかなく、”電子状態＋ストレス”（原子は動かないとする場合）用になっていませんでした。以下が、サブルーチンFORZFBの該当部分です。


      DO 1000 IA=1,KATM                                                 
      IF (NLSPD(KFTYPE(IA)).EQ.1) THEN
          LNUM = 4
      ELSE
          LNUM = 9
      END IF
      DO 1003 L=1,LNUM
!XOCL SPREAD NOBARRIER DO /IP
      DO 1001 IK=1,KV3                                                  
      DO 1002 IBAN=NBD1,NBD2                                            
      ZFC(IBAN,IK,IA,L) =DCMPLX(0.0D0,0.0D0)                            
 1002 CONTINUE                                                          
 1001 CONTINUE                                                          
!XOCL END SPREAD
 1003 CONTINUE
 1000 CONTINUE                                                          
C!XOCL END PARALLEL
C     VPP-PARALLEL START
C!XOCL PARALLEL REGION
!XOCL SPREAD NOBARRIER DO /IP
      DO 2100 IK=1,KV3                                                  
      DO 2110 IA=1,KATM                                                 
      CS=1.0D0/(WS(KFTYPE(IA))*UNIVOL)                                  
      CP=1.0D0/(WP(KFTYPE(IA))*UNIVOL)                                  
      CWL(1)=CS
      CWL(2)=CP
      CWL(3)=CP
      CWL(4)=CP
C
      IF (NLSPD(KFTYPE(IA)).EQ.1) THEN
          LNUM = 4
      ELSE
          LNUM = 9
          CD=1.0D0/(WD(KFTYPE(IA))*UNIVOL)
          CWL(5)=CD
          CWL(6)=CD
          CWL(7)=CD
          CWL(8)=CD
          CWL(9)=CD
        （正）CWL(10)=CDを付加
      END IF
      DO 3222 L=1,LNUM
      DO 3200 IBAN=NBD1,NBD2                                            
C     1991 11/28  I ---> I1                                             
      DO 3510 I=1,IBA(IK)                                               
      I1  = NBASE(I,IK)                                                 
      L1  = IG1(I1)+KX1                                                 
      L2  = IG2(I1)+KY1                                                 
      L3  = IG3(I1)+KZ1                                                 
      ZTMP=ZZZ(I,IBAN,IK)*DCONJG( ZFM2( I1 ,IA ) )                     
      ZFC(IBAN,IK,IA,L)=ZFC(IBAN,IK,IA,L)+
     &                  ZTMP*SSS(I,IK,KFTYPE(IA),L)
 3510 CONTINUE                                                          
      ZFC(IBAN,IK,IA,L)=CWL(L)*ZFC(IBAN,IK,IA,L)                            
 3200 CONTINUE                                                          
 3222 CONTINUE
C
 2110 CONTINUE                                                          
 2100 CONTINUE

配列SSS()が非局所擬ポテンシャル項に関するもので、通常の電子状態の計算では、ｌ＝２（ｄ軌道）まで計算する場合、ｓ、ｐ、ｄ合わせて９つ（LNUM = 9）までの計算が必要となります。実は、本プログラムではストレスの計算でのｄ軌道部分からの寄与として、LNUM = 10に相当する部分の計算が必要でした（深紅色で示した部分）。これを忘れたためストレス部分の計算結果が一致しない事態となった訳です。勿論、LNUM = 10として計算をし直すと、結果は一致するようになりました。

（最近行なった、MgB₂の計算に関して）
一軸圧縮による、MgB₂の計算を最近行なって、学会発表や、投稿を行なっていたのですが、上記問題の影響はあるかどうかについて言及する必要があります。
結論から言うと、確かにこの影響はあったが、不幸中の幸いにも最終結果や結論に影響を与えるほどの誤差（0.1 GPaのオーダー以下）になりませんでした。MgB₂においては、このLNUM = 10に相当する部分の寄与が結果として小さかったのが幸いしました。また、電子状態のみの計算（サブルーチンFORCEでは、LNUM = 10に相当する計算は正しく行なっていた）でのチェック計算でも、MgB₂は十分正しく収束（含むストレス）していることが確認できていました。

今回、問題に気付いたのは、MgB₂以外に、最近他の理論計算の予測（プレプリントサーバーで発見）で、CuB₂、AgB₂、 AuB₂などが高い超伝導転移温度を示すというものがあったため、筆者もこの系の計算を試みてみようと思い、平衡格子定数をセルの最適化から求めようとしたら、予想と大分異なった値になったことでした。いろいろチェックする内に今回のバグにまで行き当たった次第です。実際、当該ルーチン部分をコードしてから数年経つのに、これまで全く気付かなかったのが問題です（あまりストレス計算を必要とする機会がなかったというのもあるが）。

今回の教訓としては、これまでと異なる未経験の系を計算してみておかしな結果が出てきた時は、プログラムのバグも疑えということです（常にバグであるという訳ではない。←実はこれも問題：バグかバグでないかを確認するのも困難が伴うことが多い）。

（厳密に言うと）
上記の修正は、ストレスの計算を行なわないで、電子状態だけを計算するだけなら、LNUM = 10の部分だけ余計に計算することになるので、若干無駄な計算をしていることになります。ただ、今回のバグの遠因として、いろいろな計算条件毎に計算を場合分けすることにより、計算が複雑化し、それがバグ発生の要因となったとも言えます。従って、今回は上記処方で、”よし”としました（勿論、もっとエレガントな方法もあるかと思うが、、、）。

（更にお恥ずかしい話）
何と、このバグは既に報告（[バグレポート３]）されていた！。筆者に学習能力や、プログラム管理能力が無いことが歴然！。
ただ、現在と過去とで、本バグとその対策（対処法）に対する認識の違いは興味深い。

（Ｊ）本日のバグ（レポート３／２、２００１）[目次]

新世紀発の失敗です。今回は、現時点で原因の不明な錯誤によるものです。

今回のは、プログラムのバグではなく、入力としてのデータの間違いによるものでした。問題は二つあり、現行の筆者のシステム（プログラム実行に至る、過程として）では、今回の失敗を完全に回避することができないことと、長期間この失敗に気付けなかったことです。

失敗は、バンド計算で間違った擬ポテンシャルのデータを入力として使い続けてしまったことです。まず、バルクの計算を行ない、平衡格子定数、バンド構造等を求めました。この段階で、この物質は仮想的な物質だったため、他に比較すべき実験値、理論値の資料がありませんでした。バンド構造に、特に問題はなかったのですが、全エネルギーの出方にはちょっと変なところがありました。
それは、平衡格子定数を求めるため、格子定数-全エネルギーの計算を行なっていると、全エネルギーの最小値近傍で、全エネルギーの値がふらつくことでした。つまり、非常に細かく格子定数を変化させると、全エネルギーが、大きくなったり、小さくなったりで、ぎざぎざとなることです。これは、格子定数の変化による平面波数の急な変化によるものではないです（この場合の格子定数の変化よりずっと小さい）。また他の系では、このようなことはほとんどありませんでした。
で、この時点では、格子定数の変化（0.002 a.u.）が小さいものによるもの（誤差）で、単に他の系より全エネルギーが敏感なだけと思い、そのまま計算を次の段階（バルクから表面系へ）に進めてしまいました。

ちょっとまだ曇行きが怪しくなってきているので、現在記述を凍結中です（７／１６、２００１）。

（Ｉ）本日のバグ（レポート９／１８、２０００） [目次][関連]

今回の過ちは、前回のバグ[レポート]と深く関連しています。

（９／１９、２０００）表面の双極子モーメントを求めるために、筆者は３次元の実空間に関しての積分を行なっていました。これに必要な実空間電荷密度は、フーリエ変換（注：高速ではない）により求めていました。従って、それに要する計算時間は、４重のループ（実空間用に３重、逆空間用に１重）になり、大変膨大なものになります。Alphaチップマシン（21264使用）でさえ、計算に２、３日もかかってしまいました。これは甚だ効率の悪いものであると言えます。

偶然、数年前に書かれたコードサンプル（筆者のものではない）を見ている内に、この実空間での積分とは異なる手法が存在することに気が付きました。
それは、実空間での電子の電荷密度に関する双極子モーメントの計算式（Ｚ方向を表面垂直方向とする）は、

μ = ∫zρ(r)dr

を計算することによって求められます。積分はスラブモデルの場合、スラブの中心から、真空層の中心（スーパーセル全体の半分に相当）までの体積積分です。
筆者は、この積分計算を、まず、バンド計算から求めた、逆空間での電子の電荷密度ρ(G)を、フーリエ変換（注：高速ＦＦＴではない）してρ(r)を求め、それをもとに上式に右辺の空間に関して（３次元）の積分を行ない、双極子モーメントを求めました。これは、３次元空間に関してのメッシュの数にもよりますが、大変時間を要する計算となります。実際、Alphaチップ（Alpha21264）搭載のワークステーションを使用しても、２、３日以上かかってしまいました。これは甚だ計算時間を消費（浪費）した計算だと言えます。

ところが、よく考えてみると、これはもっとずっと簡単に計算可能であることが分かりました。きっかけは、昔ＪＲＣＡＴの森川先生から貰ったプログラムコードを見ていて、筆者と同じように双極子モーメントの計算をするルーチンがあり、その内容が筆者のものと全然異なることに気付いたことでした。

上記の、μを求める式で、ρ(r)はフーリエ逆変換（ρ(G)→ρ(r)を逆変換とする）を使って、

ρ(r) = Σρ(G)e^iGr

と表現できます。これを、先のμを求める式に代入すると、

μ = ∫z(Σρ(G)e^iGr)dr

となり、更に、

μ = Σρ(G)∫z(e^iGr)dr

となります。上の式では、（定）積分部分は最早、数値積分に頼らずに解析的に解けてしまいます（もう忘れかけているが、部分積分を使う）。つまり、逆格子Ｇに関しての和のループだけで計算可能になってしまう訳です。解析的の解く上で、上記の積分は、更に、Gr（ベクトル積）の部分が、 G_xx + G_yy + G_zzとなり、

∫z(e^iGr)dr = ∫z(e^{i(G_xx + G_yy
+ G_zz)})dxdydz

= ∫(e^iG_xx)dx∫(e^iG_yy)dy ∫z(e^iG_zz)dz となります。

正に、解析的な手法に勝る、数値計算なしと言えます。

実際の計算、検証等は、学会（日本物理学会第５５回年次大会）後に行なう予定です。

（Ｈ）本日のバグ（レポート８／２９、２０００） [目次]

今回もごく卑近な間違いです。バンド計算で得られた電子の電荷密度（逆空間）を読み込み、それを実空間へフーリエ変換（注：高速フーリエではない）して実空間の電荷密度にして、これを積分して全電荷数（電子数）を求める計算を行なう計算で、間違っていました。
扱った系は表面系（スーパーセル＋スラブモデル[*]）で、最初はスーパーセル全体の電荷数（電子数）を計算し、この段階では総電子数は、その系の真の電子数と十分な精度に一致していました。ところが、スーパーセルの半分、スラブの中心から、真空層の中心までを積分範囲として、電荷密度の和をとると、セル全体の総電子数の半分にならず若干誤差が出ました。

系（総電子８１個からなる）全体の総電子数


 total chrage        =    81.0000000025603

スラブの中心から、真空層の中心（系の半分に相当）までの総電子数


 TOTAL CHRGE(FIX)   =    40.6109851973210

上記値を見ると、総電子数が、40.5になっていませんでした。最初、この程度は誤差として許容できるのではと思ったのですが、系全体で積分すると正しく、総電子数は８１（誤差：0.0000000025603）になるのに、積分領域を半分にすると、0.1109851973210もの誤差が生じるのは、やはりおかしいということで、原因を追求することにしました。

最初、系が完全に対称になっていないため、スラブの中心の右側、左側（上側、下側とも言えるが、左右で表現することにする）で電荷密度分布が対称になっていないのではないかと考えましたが、スラブでの原子の配置は左右で完全に対象になっており、この配置で電荷密度分布が非対称になっていることはあり得ないと判断しました（もしそれで非対称になる場合は、元のバンド計算そのものがおかしいことになる）。
筆者の計算では、スラブの底を座標ゼロとして、スラブ＋真空層からなる系（スーパーセル）でバンド計算を行なっていましたが、これを真空層＋スラブ＋真空層とし、スーパーセルの中心が、スラブの中心と一致するようにしてバンド計算を行ない、それで得られた電荷密度で、上記積分計算を行なってみましたが、結果は変わりませんでした。当然、バンド計算においても結果は、いずれの場合でも同じ結果（全エネルギー等）を与えます。

ここで、万策尽きかけたのですが、問題は積分の刻みの取り方にあるのではないかと思い、今度は、積分の仕方に着目してみることにしました。
積分計算では、仮にＺ方向を表面垂直方向とすると、フーリエ変換（注：高速フーリエではない）した実空間電荷密度ρ(x,y,z)を、まずx,y方向で積分し、その後、z方向で積分します（z方向の積分刻み数５００）。実際の計算ループでは、１～５００までの繰り返しを行ないますが、この場合、積分範囲をAとすると、刻み幅を、hz=A/500とし、実際に積分を行なう座標は、z=(k-1)*hz、としていました（kは繰り返しループでのループ変数）。
筆者は、z=k*hz、の場合の計算をしてみることにしました。そうすると積分の結果は、


 TOTAL CHRGE(FIX)   =    40.3890148052776

となりました。これは、先の結果、


 TOTAL CHRGE(FIX)   =    40.6109851973210

と比べると、丁度両方を足すと、ほぼ正確な総電子数８１になることが判明しました。
系全体を積分した場合は、正しい総電子数が得られる。積分が１～５００と０～４９９で値が異なり、両方を足すとほぼ正確な総電子数になる。この段階で、筆者は積分手法に重大な誤りがあることに、気付いてしまいました。これまでの積分は、単に刻み幅分の長方形を足しているだけでした。つまり、積分を、


         TOTC   = TOTC   + CHG(I,J,K)

としていました。これは数値積分における、台形公式になっていません。台形公式として積分するには、


         TOTC   = TOTC   + (CHG(I,J,K) +CHG(I+1,J,K) )*0.5D0

としないといけません。このため積分計算の精度が悪いものになっていました。この場合、系全体での積分では、電荷密度分布がスラブの中心から左右で完全に対称になっているため、ちょうど右半分の区間の積分と左半分の区間の積分の和が、上記の台形公式と等価となっているため正確な値が得られていました。

（予告）
（９／１２、２０００）この計算は、非常に効率の悪いものであることが判明しました。わざわざ普通のフーリエ変換を行なって実空間での積分をする必要がありませんでした。近々、これに関してのレポートをする予定です。現在、筆者は大変多忙のため、暫くお待ち下さい。

（Ｇ）本日のバグ（レポート５／２５、２０００） [目次]

また、久々にやらかしてしまいました（この”久々”にも[意味]あり）。今回の誤りはOpenMPによるＳＭＰ下での並列計算で起こりました。

この前のOpenMP絡みの[バグ]以来、並列計算は行なっていなかったのですが、計算速度向上のため、しばらくぶりに計算を行なってみました。そうしたところ、計算結果がオリジナルの結果（OpenMPなしの単一ＣＰＵ動作による結果）と一致しないことが判明しました。

まず、OpenMPによる並列動作結果を以下に示します。


 SUB-NAME =     KBMSD   ----------------- TIME =     0.000000
  OCCUP2 :    WIDTH=  5.000000000000000E-003
 ---------- THE FERMI ENERGY =-0.3609612356673470D-01   81.000000
 SUB-NAME =     FERMI   ----------------- TIME =     0.000000
 CALL FORCE
 SUB-NAME =     FORCE   ----------------- TIME =     0.000000
 SUB-NAME =     CHAVER  ----------------- TIME =     0.000000
 SUB-NAME =     XCFFT   ----------------- TIME =     0.000000
 SUB-NAME =     FORLOC  ----------------- TIME =     0.000000
 MINUS CHG ICCCC =            0
 TOTAL ENERGY FOR   1-TH ITERATION= -95.1142476      -0.9511425D+02
 SUB-NAME =     ENERGY  ----------------- TIME =     0.000000
  ITER=  1 ET(H)=  0.1000000D+06 ET(M)=********** DC=  0.5504872D-01
 SUB-NAME =     CONV2   ----------------- TIME =     0.000000
 READ 2 --- REINIT OR NOT
 READ 2 --- REINIT OR NOT           0
 EVOUT  ITER =            1
 MINUS CHG ICCCC =            7
 DTIME =    2.00000000000000     
 SUB-NAME =     C3FFT   ----------------- TIME =     0.000000
 SUB-NAME =     MSD     ----------------- TIME =     0.000000
 ---------- THE FERMI ENERGY =-0.8970905953935529D-01   81.000000
 SUB-NAME =     FERMI   ----------------- TIME =     0.000000
 CALL FORZFB
 SUB-NAME =     FORCE   ----------------- TIME =     0.000000
 SUB-NAME =     CHAVER  ----------------- TIME =     0.000000
 SUB-NAME =     XCFFT   ----------------- TIME =     0.000000
 SUB-NAME =     FORLOC  ----------------- TIME =     0.000000
 MINUS CHG ICCCC =           13
 TOTAL ENERGY FOR   2-TH ITERATION=-103.3155412      -0.1033155D+03
 SUB-NAME =     ENERGY  ----------------- TIME =     0.000000
  ITER=  2 ET(H)=  0.8201294D+01 ET(M)=********** DC=  0.5515755D-01
 SUB-NAME =     CONV2   ----------------- TIME =     0.000000
 MINUS CHG ICCCC =            7
 ---------- THE FERMI ENERGY =-0.4685106501394076D-01   81.000000
 CALL FORZFB
 MINUS CHG ICCCC =           23
 TOTAL ENERGY FOR   3-TH ITERATION= -99.3989027      -0.9939890D+02
 >> ETOOLD.LT.ETONEW <<
 ETONEW-ETTOOLD=   3.91663855637013     
  ITER=  3 ET(H)=  0.1000000D+01 ET(M)=********** DC=  0.5920008D-01
 MINUS CHG ICCCC =            7
 ---------- THE FERMI ENERGY =-0.2932751999252538D-01   81.000000
 CALL FORZFB

次に、正しい（オリジナル）結果を示します。


 SUB-NAME =     KBMSD   ----------------- TIME =     0.000000
  OCCUP2 :    WIDTH=  5.000000000000000E-003
 ---------- THE FERMI ENERGY =-0.3609612356673101D-01   81.000000
 SUB-NAME =     FERMI   ----------------- TIME =     0.000000
 CALL FORCE
 SUB-NAME =     FORCE   ----------------- TIME =     0.000000
 SUB-NAME =     CHAVER  ----------------- TIME =     0.000000
 SUB-NAME =     XCFFT   ----------------- TIME =     0.000000
 SUB-NAME =     FORLOC  ----------------- TIME =     0.000000
 MINUS CHG ICCCC =            0
 TOTAL ENERGY FOR   1-TH ITERATION= -83.4446440      -0.8344464D+02
 SUB-NAME =     ENERGY  ----------------- TIME =     0.000000
  ITER=  1 ET(H)=  0.1000000D+06 ET(M)=********** DC=  0.5504872D-01
 SUB-NAME =     CONV2   ----------------- TIME =     0.000000
 READ 2 --- REINIT OR NOT
 READ 2 --- REINIT OR NOT           0
 EVOUT  ITER =            1
 MINUS CHG ICCCC =            7
 DTIME =    2.00000000000000     
 SUB-NAME =     C3FFT   ----------------- TIME =     0.000000
 SUB-NAME =     MSD     ----------------- TIME =     0.000000
 ---------- THE FERMI ENERGY =-0.8970905953933028D-01   81.000000
 SUB-NAME =     FERMI   ----------------- TIME =     0.000000
 CALL FORZFB
 SUB-NAME =     FORCE   ----------------- TIME =     0.000000
 SUB-NAME =     CHAVER  ----------------- TIME =     0.000000
 SUB-NAME =     XCFFT   ----------------- TIME =     0.000000
 SUB-NAME =     FORLOC  ----------------- TIME =     0.000000
 MINUS CHG ICCCC =           13
 TOTAL ENERGY FOR   2-TH ITERATION= -91.6459376      -0.9164594D+02
 SUB-NAME =     ENERGY  ----------------- TIME =     0.000000
  ITER=  2 ET(H)=  0.8201294D+01 ET(M)=********** DC=  0.5515755D-01
 SUB-NAME =     CONV2   ----------------- TIME =     0.000000
 MINUS CHG ICCCC =            7
 ---------- THE FERMI ENERGY =-0.4685106501392522D-01   81.000000
 CALL FORZFB
 MINUS CHG ICCCC =           23
 TOTAL ENERGY FOR   3-TH ITERATION= -87.7292991      -0.8772930D+02
 >> ETOOLD.LT.ETONEW <<
 ETONEW-ETTOOLD=   3.91663855636882     
  ITER=  3 ET(H)=  0.1000000D+01 ET(M)=********** DC=  0.5920008D-01
 MINUS CHG ICCCC =            7
 ---------- THE FERMI ENERGY =-0.2932751999250589D-01   81.000000
 CALL FORZFB

（説明）
深紅色と深桃色で、それぞれ示されているのが、全エネルギーの値です。お分かりのように、全エネルギーの値は、それぞれ全く一致していません。

以前の[バグ]が解決した時点では、オリジナルの結果との比較を怠っていました。並列計算で、１ＣＰＵのみの実行結果と、多ＣＰＵ（２ないし４）での実行結果との比較（これらの結果は一致）したことで、” よし”としていました。この点で詰めが甘かったと言えます。

不一致の原因究明の作業過程で浮かび上がった事実は、

OpenMPによる計算では、１ＣＰＵのみの計算結果と、多ＣＰＵの計算結果は一致。
OpenMPによる計算結果と、オリジナルの結果は全エネルギーが不一致。
全エネルギーの結果は第一回目のイタレーションの結果から不一致。
オリジナル計算のソースプログラムを、OpenMP並列計算を行なっているディレクトリに複写し、OpenMP計算と全く同じ環境（違いは、並列動作しないだけ）で（単一ＣＰＵで）計算しても、計算結果は一致しない（当然、この結果はオリジナルの結果とは一致した）。
OpenMP用のソースプログラムを、OpenMP用のコンパイルオプション（- omp）なしで、コンパイルしても結果は変わらなかった。
OpenMPによる並列計算結果と、オリジナルの結果を、よくよく比べてみると、全エネルギーは一致していないが、フェルミエネルギーの値（深橙色で表示）は、互いに一致していることが分かった（重要）。

特に、最後の項目は非常に重要です。初期の段階から、第一回目のイタレーションでの全エネルギーが一致しないことは判明していましたが、第一回目のフェルミエネルギーの値は一致していました（非常に小さな値の部分での不一致は除く）。筆者は、この段階で、最初のフェルミエネルギーの計算をするサブルーチンまでは正しく計算が行なわれていると判断しましたが、それ以降については、全エネルギーが一致しないので、フェルミエネルギーも、これ以降は一致していないだろうと思い込んでしまいました。この計算では、全エネルギーの計算の前にフェルミエネルギーの計算が行なわれます（そうでないと、全エネルギーは計算できない）。

この判断は完全な誤りで、よくよく並列計算の結果と、オリジナルの結果を比べてみると、その後のイタレーションでもフェルミエネルギーの値は誤差を除いて完全に一致していました。これは、非常に重要なことで、普通数イタレーション計算を経た後で、全エネルギーが一致しないのに、フェルミエネルギーだけが一致することはありえません。そこで、各イタレーションでの全エネルギーの値とフェルミエネルギーの値を以下に示したいと思います。まず、OpenMPによる並列計算の場合を示します。


 TOTAL ENERGY FOR   1-TH ITERATION= -95.1142476      -0.9511425D+02
 TOTAL ENERGY FOR   2-TH ITERATION=-103.3155412      -0.1033155D+03
 TOTAL ENERGY FOR   3-TH ITERATION= -99.3989027      -0.9939890D+02
 TOTAL ENERGY FOR   4-TH ITERATION= -97.7885536      -0.9778855D+02
 TOTAL ENERGY FOR   5-TH ITERATION= -97.6782727      -0.9767827D+02

 1-TH ----- THE FERMI ENERGY =-0.3609612356673470D-01   81.000000
 2-TH ----- THE FERMI ENERGY =-0.8970905953935529D-01   81.000000
 3-TH ----- THE FERMI ENERGY =-0.4685106501394076D-01   81.000000
 4-TH ----- THE FERMI ENERGY =-0.2932751999252538D-01   81.000000
 5-TH ----- THE FERMI ENERGY =-0.1644041912340696D-01   81.000000
 6-TH ----- THE FERMI ENERGY =-0.6217050465534285D-02   81.000000

次に、オリジナルの場合を示します。


 TOTAL ENERGY FOR   1-TH ITERATION= -83.4446440      -0.8344464D+02
 TOTAL ENERGY FOR   2-TH ITERATION= -91.6459376      -0.9164594D+02
 TOTAL ENERGY FOR   3-TH ITERATION= -87.7292991      -0.8772930D+02
 TOTAL ENERGY FOR   4-TH ITERATION= -86.1189500      -0.8611895D+02
 TOTAL ENERGY FOR   5-TH ITERATION= -86.0086691      -0.8600867D+02

 1-TH ----- THE FERMI ENERGY =-0.3609612356673101D-01   81.000000
 2-TH ----- THE FERMI ENERGY =-0.8970905953933028D-01   81.000000
 3-TH ----- THE FERMI ENERGY =-0.4685106501392522D-01   81.000000
 4-TH ----- THE FERMI ENERGY =-0.2932751999250589D-01   81.000000
 5-TH ----- THE FERMI ENERGY =-0.1644041912338591D-01   81.000000
 6-TH ----- THE FERMI ENERGY =-0.6217050465513631D-02   81.000000

以上を比較してみると、確かに、小数点以下の非常に小さな部分を除けば、フェルミエネルギーの値は一致し、全エネルギーの値は異なっています。前述のように全エネルギーが異なる場合、バンド構造そのものも異なるはずで、フェルミエネルギーは一致しないはずなのに、実際は一致しています。この場合考えられるのが、バンド計算そのもの（セルフコンシステントな計算過程：つまりイタレーション）には影響しない部分で不一致が起こっていることが考えられます。そこで、OpenMPによる並列計算とオリジナルな計算での各イタレーションでの全エネルギーの値の差をとってみると、それはイタレーション回数に関係なく、11.6696036（絶対値）という一定値であることが分かりました。

（原因究明）
この、11.6696036という値は、調べてみると、部分内殻補正で使用していた補正項の値であることが判明しました。この補正項は、バンド計算の最初で設定されるもので、値そのものは全エネルギーの計算の時に、全エネルギーの値に足し込まれるだけで、バンド計算そのものの結果には全く影響を与えないものです。前回の[バグ]の時にデバッグ中に、この計算部分がバグの原因ではと削除しておいたのを完全に忘れていました。元に戻さずに計算を行なったため、全エネルギーの値が一致しないように見え、並列計算に深刻なバグがあるのではと慌てた次第です。

実際は、バグとしてはそう深刻なものではなく、ほんの数行の修正（今回の場合は削除した部分の追加）で、問題を解決することができました。問題は、むしろバグを直したと思って、２、３ヶ月もほったらかしておくと、何をしたかを忘れてしまうことです。作業記録をきちんと残し、やり残したことや、無用な変更、削除を元に戻しておくことを怠らないよう（うっかり忘れないよう）にし、その時できることは、その時に完結させておく（記録も残す）ようにしましょう。

（Ｆ）本日のバグ（レポート２／２４、２０００） [目次]

今回のバグは、部分内殻補正絡みものです。

筆者のプログラムで、部分内殻補正の計算を行なう場合、扱う物質や系によってストレスの値が発散する（異常に大きな数になる）ことがありました。ただ、扱っていた系が表面系だったり（普通、表面系ではユニットセルにかかる圧力〔ストレス〕は使わない。表面ストレスを求めようという場合は別だが）、その時点でストレスは必要なかったりで、さして気に止めていませんでした。

しかし、やはりしっくりこないこともあり、原因を突き止めることにしました。経験上バルク系では、ストレスが異常に起こることはほとんどなかったのですが、最近テストとして計算していたバルク系で、ストレス値が10¹⁸もの値になる事態に遭遇していたので、これをデバッグしてみることにしました。以下に、異常なストレス値が出る場合の出力結果を示します。


 TOTAL STR 1  1 =  -0.106726751049972     
              2 =   1.114843552513939E-003
              3 =   3.786085519029700E-004
              4 =  -0.107180229756837     
              5 =  -6.187637920119319E-006
              6 =  -0.107625292219891     
 REAL TOTAL CHARGE =    8.00229133076466      IN XSTPC
 ICOUNT =         2244 XSTPC END
 STMQ,SS1,SSX,ZXXX(1) =   4.499630388065659E-003  0.000000000000000E+000
  7.977901970862270E+020 (-1.118441627010112E+022,379875.028700487)
 ZCHG(1),ZRHPC(1),ZVXC(1),ZRRPC(1),ZEXC(1) = 
 (6.713436944043497E-002,0.000000000000000E+000)
 (4.196147755972463E-003,0.000000000000000E+000)
 (-0.279550593402226,0.000000000000000E+000)
 (0.000000000000000E+000,0.000000000000000E+000)
 (-0.216469172365056,0.000000000000000E+000)
 TOTAL STR 2  1 =   7.977901970862270E+020
              2 =   0.000000000000000E+000
              3 =   0.000000000000000E+000
              4 =   7.977901970862274E+020
              5 =   0.000000000000000E+000
              6 =   7.977901970862264E+020
 ZVX,ZEX =  -8.052755104079717E+020 -3.244955581708985E-002
 STM1(XY)=   1.128607285201954E-002 -0.163484318525950     
 STU,V,W =   0.126066049507907      -5.284113835626290E-004
 -8.052755104079717E+020
 STOT    =  -8.052755104079717E+020
 TOTAL STR 3  1 =  -8.052755104079578E+020
              2 =   5.204170427930421E-018
              3 =   2.981555974335137E-018
              4 =  -8.052755104079581E+020
              5 =   5.175371307723775E-019
              6 =  -8.052755104079573E+020
 SIGNL     1  1 =   0.221485143478849     
              2 =   0.000000000000000E+000
              3 =   0.000000000000000E+000
              4 =   0.221485143478849     
              5 =   0.000000000000000E+000
              6 =   0.221485143478849     
 SIGNL     2  1 =  -3.008891331718425E-002
              2 =  -4.277013915213738E-005
              3 =  -1.448585533509181E-005
              4 =  -3.008246681184892E-002
              5 =   6.375125928747571E-007
              6 =  -3.007457410262290E-002
 SIGNL     3  1 =  -6.993430363272542E-002
              2 =   9.378111455720368E-004
              3 =   3.224221320190277E-004
              4 =  -7.021242131084532E-002
              5 =   3.474492615092482E-006
              6 =  -7.048630456989730E-002
 SIGNL     4  1 =  -0.135626705862546     
              2 =   9.203622245516506E-004
              3 =   3.093997338064650E-004
              4 =  -0.136006849184124     
              5 =  -1.467699662472126E-005
              6 =  -0.136396710849138     
 SIGNL     5  1 =   0.135041992087479     
              2 =  -2.647062024157973E-003
              3 =  -9.026027297115934E-004
              4 =   0.136016126412455     
              5 =   9.309930437597757E-006
              6 =   0.136972147395874     
 SIGEWA       1 =   5.254142857619389E-002
 SIGEWA       2 =   7.470500746495105E-018
 SIGEWA       3 =  -4.548484389417145E-018
 SIGEWA       4 =   5.254142857619389E-002
 SIGEWA       5 =   1.892189202705232E-018
 SIGEWA       6 =   5.254142857619388E-002
 SIGNL        1 =   0.120877212753873     
 SIGNL        2 =  -8.316587931864228E-004
 SIGNL        3 =  -2.852667192211925E-004
 SIGNL        4 =   0.121199532584486     
 SIGNL        5 =  -1.255060979156260E-006
 SIGNL        6 =   0.121499701353065     
 TOTAL STR 4  1 =   0.000000000000000E+000
              2 =  -8.316587931864153E-004
              3 =  -2.852667192211971E-004
              4 =   0.000000000000000E+000
              5 =  -1.255060979154368E-006
              6 =   0.000000000000000E+000
 ETOT1,TOTCH    =   0.316363262756181        8.00000000000000     
 TOTAL STR 5  1 =   0.000000000000000E+000
              2 =   0.000000000000000E+000
              3 =   0.000000000000000E+000
              4 =   0.000000000000000E+000
              5 =   0.000000000000000E+000
              6 =   0.000000000000000E+000
 TOTAL SUMM   1 =  -7.485313321730875E+018
              2 =   2.831847593275287E-004
              3 =   9.334183268177589E-005
              4 =  -7.485313321730745E+018
              5 =  -7.442698899273169E-006
              6 =  -7.485313321731027E+018
 TOTAL SUMM OP1 =  -7.485313321730884E+018
              2 =   7.058607893785836E-023
              3 =   0.000000000000000E+000
              4 =  -7.485313321730884E+018
              5 =  -2.258754526011467E-021
              6 =  -7.485313321730885E+018
 CALL FORLOC AND MD
 NOW MD SET EPP          11
 EPP(1,1) =    1.00000000000000       0.000000000000000E+000
 EPP(1,2) =   0.000000000000000E+000  0.000000000000000E+000
 EPP(1,3) =   0.000000000000000E+000  0.000000000000000E+000
 EPP(2,1) =   0.000000000000000E+000  0.000000000000000E+000
 EPP(2,2) =    1.00000000000000       0.000000000000000E+000
 EPP(2,3) =   0.000000000000000E+000  0.000000000000000E+000
 EPP(3,1) =   0.000000000000000E+000  0.000000000000000E+000
 EPP(3,2) =   0.000000000000000E+000  0.000000000000000E+000
 EPP(3,3) =    1.00000000000000       0.000000000000000E+000
 TOTAL ENERGY FOR  11-TH ITERATION= -32.4511839      -0.3245118D+02
 MIXING RATE =   0.800000000000000     
  ITER= 11 ET(H)=  0.1239419D-01 ET(M)=337.122064 DC=  0.1087433D-01

もうデバッグ用の計算になっているので、１１回イタレーションを進めた段階でプログラムを止めています（計算は収束していない）。深紅色（CRIMSON）で表示されている部分が異常に大きな値をしめしているストレス値です。もともとデバッグも考慮して、各計算途中段階でのストレスの値を出力するようにしていました（上記ではデバッグ用に新たに出力させたものもある）。この結果から、かなり早い段階（”TOTAL STR 2 1 = ”と表示されたところ）で、値が異常になっていることが分かります。

筆者がまず最初に注目したのは、この”TOTAL STR 2 1 = ”の段階の表示部分で、これらの変数（配列変数）は最初の値（つまりA(1)とかB(1)）しか表示していないことでした。筆者のバンド計算プログラムでは、逆格子ベクトルは大きさ順に並べるので、最初の配列変数の値（例 ZCHG(1)）は、Ｇ＝０での値となります。Ｇ＝０では発散項などが生じる場合があり、扱いに特別な考慮が必要な場合があり、ストレスの計算でも変数の１番目は別途計算する場合があり、それがTOTAL STR 2 1の部分だった訳です。

この段階で値がおかしいので、さらに詳しく値を表示させることにし、


 STMQ,SS1,SSX,ZXXX(1) =   4.499630388065659E-003  0.000000000000000E+000
  7.977901970862270E+020 (-1.118441627010112E+022,379875.028700487)
 ZCHG(1),ZRHPC(1),ZVXC(1),ZRRPC(1),ZEXC(1) = 
 (6.713436944043497E-002,0.000000000000000E+000)
 (4.196147755972463E-003,0.000000000000000E+000)
 (-0.279550593402226,0.000000000000000E+000)
 (0.000000000000000E+000,0.000000000000000E+000)
 (-0.216469172365056,0.000000000000000E+000)

という結果を得ました。この計算結果から逆格子空間での電荷密度、部分内殻電荷密度、交換相関項の値（ZCHG(1),ZRHPC(1),ZVXC(1),ZRRPC(1),ZEXC(1)）には異常はなく、SSX、そしてZXXX(1)の値に問題があることが分かります。

特に問題となるのが、ZXXX(1)で、これと連動してSSXという値もおかしくなります。ストレス値の異常の根本原因は、変数ZXXXにありそうだということで、今度はZXXX(1)だけでなく、ZXXXの全ての値を表示させてみました。結果としてZXXXは全ての値が異常（発散している）であることが分かりました。

ZXXXは、部分内殻補正導入によって生じる、新たなストレス項（交換相関項が関わる）を計算するための変数で、それは以下の実空間での交換相関項に関しての計算、


C  -- X ------------                                                    
          DO 1020 I=1,KSUM                                              
            RS=CO2*( 1.0D0/CHGB1(I) )**THIRD
            RH1   = -0.458D0*CO3/RS-(0.44D0*CO3*RS + 3.432D0)/          
     &                                          (RS+7.8D0)**2           
            RH2   = -0.458D0/RS-0.44D0/(RS+7.8D0)                       
C                                                                       
            RHV(I)=RH1 - RH2 
            ZV1D(I)=ZRX(I)*(RHV(I))/CHGB1(I)
 1020     CONTINUE

の、変数ZV1Dをフーリエ変換（実空間→逆空間をフーリエ変換と定義）したものがZXXXとなります。問題なのは、深紅色で示した、1/CHGB1(I)の部分で、CHGB1(I)は実空間での電荷密度です。実空間の電荷密度は、逆格子空間で計算された電荷密度ZCHGを逆フーリエ変換して求めます。部分内殻補正を考慮した場合は、ZCHG+ZRHPCと部分内殻補正電荷密度がZCHGに加わります。普通、実空間電荷密度が負の値になることはありえないのですが、場合によっては負の値になってしまうことがあります。それは系が表面系で、真空領域がある場合（どうしてもフーリエ変換等の過程で誤差が入り込む）や、今回の部分内殻補正を考慮した場合です。部分内殻補正用の電荷密度を加えて、（逆）フーリエ変換する間に、何らかの誤差が入り込んで電荷密度が負になる場合があります（どういう兼ね合いでそうなるかは更に原因調査中）。

そして実空間での電荷密度（部分内殻補正電荷密度込み）CHGB1が負の値になった場合、CHGB1の値を10^-40にする処理をしていました。このため今回の場合、1/CHGB1が発散してしまう事態となった訳です。分母の項は、この発散を打ち消す程の値にならないことも判明しました。分母の項は RS=CO2*( 1.0D0/CHGB1(I) )**THIRDの逆数のオーダーなのですが、このTHIRDが１／３のため、分母はCHGB1の１／３乗のオーダーで、1/CHGB1をかけることにより発散してしまいます。

対処の仕方は、RCHGB(I) = 1.0D0/CHGB1(I)という変数を定義し、CHGB1が負の値になった場合、RCHGBはゼロになるようにしておきます。そして、


C            ZV1D(I)=ZRX(I)*(RHV(I))/CHGB1(I) <-- 2/18, 2000
            ZV1D(I)=ZRX(I)*(RHV(I))*RCHGB(I)

とすれば、ZXXX（ZV1Dをフーリエ変換）の発散を防ぐことができます。

（Ｅ）本日のバグ（レポート１／１７、２０００） [目次]

今回のバグは比較的軽微なものです（場合によっては深刻）。

問題はサブルーチンPCC 内の配列変数CHGPCです。この変数は、サブルーチンＰＣＣ内で定義されている局所的な配列変数で、READ文により読み込まれた部分内殻補正用電荷密度の総電荷数が格納されます。問題は、この読み込みは最初のイタレーションでしか行なわれないのですが、CHGPCはその後のイタレーションでも使用されます。
普通、その場合は、サブルーチンの引数としてCHGPCはメインプログラム側へ引き継がれるようになっていなければならないのに、実際は、SUBROUTINE PCC(KTPCC,RHPC,SCHGPC,EPC,PCM)と、CHGPCは（）内の引数の中に存在しません。

このため、最初のイタレーション以降、CHGPC内の値は不定となってしまい、正しい値が保持される保証は全くありませんでした。実際に調べてみると CHGPCの値は、ＤＥＣ（現ＣＯＭＰＡＱ）のAlphaマシン、富士通のＶＰＰ（並列計算）では最初のイタレーションの値が正しく保持されていました。但しこれは非常に好運だったと言えます。また、CHGPC自身も他のサブルーチンで使用されておらず（されていたら引数を受け渡していたはず）、たとえ間違った値になったとしても深刻な問題を起こす変数ではありませんでした。逆にこのため、この問題点に気付くのが遅れたとも言えます。

問題発見の端諸は、旧ＤＥＣのＳＭＰマシン上でOpenMPプログラムのテスト中に、並列計算が途中で止まってしまい、止まる箇所がサブルーチンＰＣＣであったことです。詳しく問題原因を調査した結果、上記のバグが判明しました。並列計算が止まる原因もCHGPCの値が、OpenMPによる並列計算では保持されず（ＶＰＰでの並列計算では保持【単なる偶然？】、但し全ての場合はチェックしていない）、これが元でプログラムが異常終了してしまったようです。CHGPCをメインルーチンに引き渡すようにすると並列計算は問題なく行なえるようになりました。

今回のバグは結果として軽微（OpenMPによる並列計算では障害があったが）でしたが、まかり間違えば大きなバグに繋がる可能性があります。局所変数、大域変数の定義、動作等のチェックはちゃんとするようにしましょう。コンパイラによってはオプション等で引数等のチェックは可能ですが、サブルーチン内で定義された（配列）変数が局所的か大域的かまでは普通判断できません（人が見るしかない）。

（４／２７、２００９追記）今年度、スーパーコンピュータ更新に伴い、本格的なOpenMP化作業を始めた。情けないことに再び、CHGPCに関する同じ間違いを犯してしまう。CHGPCに関して上記と同じ修正を行なうことで問題は解決する。９年経っても進歩なし。そして過去の教訓が全く生かされていない（←強く反省すべき点）。

（１２／１、２０１１追追記）再び同じ間違いをやらかす。悲しい。上記（４／２７、２００９）に関する修正が行なわれていいないプログラム（コード）がまだ残っており、実行中にエラーで計算が止まってしまった。早速、修正を施すが、全くもってなってない話（←改めて強く反省）。

（Ｄ）本日のバグ（レポート４／２０、１９９９） [目次]

３カ月ぶりのバグレポートですが、今回は、これまでのバグの中でも、最大級に深刻なものです。それは以下のようなものです。サブルーチンFORLOC （既に修正済み）において、これまでは、


C     FORCE2 AND PRESSURE ADDITION SSN(KNG1)                            
      DO 161 IA=1,KATM                                                  
      DO 162 I=1,KG                                                     
        FX = GX(I)                                                      
        FY = GY(I)                                                      
        FZ = GZ(I)                                                      
        ZFORT     =CDEXP(-ZI*(CATX(IA)*FX+CATY(IA)*FY+CATZ(IA)*FZ))    
        ZFTMP     =ZFORT  *PSC(I,KFTYPE(IA))*DCONJG(ZCHG(I))        
      ZFORC2(IA,1)=ZFORC2(IA,1)+GX(I)*ZFTMP                             
      ZFORC2(IA,2)=ZFORC2(IA,2)+GY(I)*ZFTMP                             
      ZFORC2(IA,3)=ZFORC2(IA,3)+GZ(I)*ZFTMP                             
  162 CONTINUE                                                          
      IF (KTPCC(KFTYPE(IA)).EQ.1) THEN                                  
          DO 1000 N=1,6                                                 
            ZFPC(N)=DCMPLX(0.0D0,0.0D0)                                 
 1000     CONTINUE                                                      
          DO 200 I=1,KG                                                 
            ZFORT =CDEXP(-ZI*(CATX(IA)*FX+CATY(IA)*FY+CATZ(IA)*FZ))    
            ZFTMP =ZFORT  *RHPCG(I,KFTYPE(IA))*DCONJG(ZVXC(I))      
            ZFPC(1)=ZFPC(1)+GX(I)*ZFTMP                                 
            ZFPC(2)=ZFPC(2)+GY(I)*ZFTMP                                 
            ZFPC(3)=ZFPC(3)+GZ(I)*ZFTMP                                 
            ZFTMP =ZFORT  *RHPCG(I,KFTYPE(IA))*DCONJG(ZVXCPC(I))    
            ZFPC(4)=ZFPC(4)+GX(I)*ZFTMP                                 
            ZFPC(5)=ZFPC(5)+GY(I)*ZFTMP                                 
            ZFPC(6)=ZFPC(6)+GZ(I)*ZFTMP                                 
  200     CONTINUE

としていました。問題は、ＤＯ１６２ループ内のＦＸ、ＦＹ、ＦＺです。これはフォームファクターを計算するための変数なのですが、これがＤＯ２００ループ内で定義されていません。ＤＯ１６１とＤＯ２００は完全に独立したループで、ＤＯ２００ループではＦＸ、ＦＹ、ＦＺを再定義しなければならないのに、それを怠っていました。これでは正しい計算ができるはずがありません。

修正は、


C     FORCE2 AND PRESSURE ADDITION SSN(KNG1)                            
      DO 161 IA=1,KATM                                                  
      DO 162 I=1,KG                                                     
        FX = GX(I)                                                      
        FY = GY(I)                                                      
        FZ = GZ(I)                                                      
        ZFORT     =CDEXP(-ZI*(CATX(IA)*FX+CATY(IA)*FY+CATZ(IA)*FZ))    
        ZFTMP     =ZFORT  *PSC(I,KFTYPE(IA))*DCONJG(ZCHG(I))        
      ZFORC2(IA,1)=ZFORC2(IA,1)+GX(I)*ZFTMP                             
      ZFORC2(IA,2)=ZFORC2(IA,2)+GY(I)*ZFTMP                             
      ZFORC2(IA,3)=ZFORC2(IA,3)+GZ(I)*ZFTMP                             
  162 CONTINUE                                                          
      IF (KTPCC(KFTYPE(IA)).EQ.1) THEN                                  
          DO 1000 N=1,6                                                 
            ZFPC(N)=DCMPLX(0.0D0,0.0D0)                                 
 1000     CONTINUE                                                      
          DO 200 I=1,KG                                                 
            FX = GX(I)
            FY = GY(I)
            FZ = GZ(I)
            ZFORT =CDEXP(-ZI*(CATX(IA)*FX+CATY(IA)*FY+CATZ(IA)*FZ))    
            ZFTMP =ZFORT  *RHPCG(I,KFTYPE(IA))*DCONJG(ZVXC(I))      
            ZFPC(1)=ZFPC(1)+GX(I)*ZFTMP                                 
            ZFPC(2)=ZFPC(2)+GY(I)*ZFTMP                                 
            ZFPC(3)=ZFPC(3)+GZ(I)*ZFTMP                                 
            ZFTMP =ZFORT  *RHPCG(I,KFTYPE(IA))*DCONJG(ZVXCPC(I))    
            ZFPC(4)=ZFPC(4)+GX(I)*ZFTMP                                 
            ZFPC(5)=ZFPC(5)+GY(I)*ZFTMP                                 
            ZFPC(6)=ZFPC(6)+GZ(I)*ZFTMP                                 
  200     CONTINUE

です。ＤＯ２００ループ内でＦＸ、ＦＹ、ＦＺを再定義します。このループはＰＣＣ（部分内殻補正）に関しての力の計算を行なっている部分です。

もし、ZFORTの計算で、


            ZFORT =CDEXP(-ZI*(CATX(IA)*GX(I)+CATY(IA)*GY(I)+CATZ(IA)*GZ(I)))

としていたら、この間違いはなかったのですが、わざわざＦＸ、ＦＹ、ＦＺに代入してZFORTを計算するのには、確か意味がありました（多分、ベクトル計算の上でこのようにした方が、計算上有利だったように憶えているのですが確かではありません。現在調査中）。

調べてみると、古いバージョンのプログラムでは変数ZFORTが、配列変数ZFORT(I)になっていたことが分かりました。以下にそれを示します。


C     FORCE2 AND PRESSURE ADDITION SSN(KNG1) Old version
      DO 161 IA=1,KATM                                                  
      DO 162 I=1,KG                                                     
        FX = GX(I)                                                      
        FY = GY(I)                                                      
        FZ = GZ(I)                                                      
        ZFORT(I)   =CDEXP(-ZI*(CATX(IA)*FX+CATY(IA)*FY+CATZ(IA)*FZ))
      ZFTMP       =ZFORT(I)   *PSC(I,KFTYPE(IA))*DCONJG(ZCHG(I))        
C*****ZFTMP       =ZFORM(I,IA)*PSC(I,KFTYPE(IA))*DCONJG(ZCHG(I))        
      ZFORC2(IA,1)=ZFORC2(IA,1)+GX(I)*ZFTMP                             
      ZFORC2(IA,2)=ZFORC2(IA,2)+GY(I)*ZFTMP                             
      ZFORC2(IA,3)=ZFORC2(IA,3)+GZ(I)*ZFTMP                             
  162 CONTINUE                                                          
      IF (KTPCC(KFTYPE(IA)).EQ.1) THEN                                  
          DO 1000 N=1,6                                                 
            ZFPC(N)=DCMPLX(0.0D0,0.0D0)                                 
 1000     CONTINUE                                                      
          DO 200 I=1,KG                                                 
            ZFTMP =ZFORT(I)   *RHPCG(I,KFTYPE(IA))*DCONJG(ZVXC(I))
C*****      ZFTMP =ZFORM(I,IA)*RHPCG(I,KFTYPE(IA))*DCONJG(ZVXC(I))      
            ZFPC(1)=ZFPC(1)+GX(I)*ZFTMP                                 
            ZFPC(2)=ZFPC(2)+GY(I)*ZFTMP                                 
            ZFPC(3)=ZFPC(3)+GZ(I)*ZFTMP                                 
            ZFTMP =ZFORT(I)   *RHPCG(I,KFTYPE(IA))*DCONJG(ZVXCPC(I))    
C*****      ZFTMP =ZFORM(I,IA)*RHPCG(I,KFTYPE(IA))*DCONJG(ZVXCPC(I))    
            ZFPC(4)=ZFPC(4)+GX(I)*ZFTMP                                 
            ZFPC(5)=ZFPC(5)+GY(I)*ZFTMP                                 
            ZFPC(6)=ZFPC(6)+GZ(I)*ZFTMP                                 
  200     CONTINUE

ＤＯ１６２もＤＯ２００ループも変数IA（ユニットセル内の原子数に関して）のＤＯ１６１ループの内側にあるので、ZFORT(I)（で十分）となります。注釈でのZFORT(I,IA)から、更に古いプログラムでは変数I（逆格子数）及びIA （原子数）の配列であることが分かります。これから、ZFORT(I)からZFORTへの変更は、ベクトル化云々より、メモリを節約するためだったと考えられます。そのための作業過程（ZFORT(I)→ZFORT）で、ＦＸ、ＦＹ、ＦＺの再定義を忘れたため、このようなバグが生じてしまいました。

以前（数年前）から、ＰＣＣに関しての力には、虚数の力が出てきたり、対称性から考えてありえない力が出てきたりしていました。筆者は、それを負の電荷密度（ＰＣＣを導入すると動径方向の実空間電荷をフーリエ変換し、再度実空間電荷に逆フーリエする過程で出てくる）によるものだと解釈していました。実際それによる影響もあるのですが、もっと重大で深刻な問題が潜んでいました。ＰＣＣ版を導入してから数年以上経つ今の今まで全く気付きんませんでした。

このバグ発見の経緯は、ＧＧＡ導入（まだ内容はできていません。→一応出来ていますがまだ未完成〔４／２１、２０１０〕）に向けて、どうも力の挙動がおかしいと思い、以前からＰＣＣを考慮していた計算でも力に賦に落ちないところがあったので、まずＰＣＣのみ考慮したバージョンで、簡単な系でテストをしていたら、バグに気付いた次第です。
テストは、Alを仮想的なＢＣＣ構造として、立方体のユニットセルを考え、その中に二つのAl原子を置いて、１１１方向にそれを近付けたり、遠ざけた場合（対称性からは近付けたことと同じ）での力の値を調べました。この場合、 Al原子は１１１方向上にあれば、どのような位置同士であっても、１１１方向に同じ大きさの力（方向は互いに反対）を受けます。ＰＣＣを考慮しない場合は、そうなるのですが、ＰＣＣ版ではそうなりませんでした。

これは何か、ＧＧＡ導入以前にＰＣＣによる力の計算に深刻な問題があると判断し、更にしろいろテストしていく内にバグを発見しました。バグそのものはさして複雑でも、見つけ難いものでもなく、非常に単純なものでした。しかしながら、これによる影響は大変深刻です。少なくとも、これまで行なってきたＰＣＣ考慮版の計算で、構造最適化を施したものの信頼性が（完全に駄目ではないにしろ）大分落ちてしまいます。更に、現在保有するプログラムを全て、修正しなければなりません。この作業は見落としなく、全てのプログラム（いろいろなマシンに分散して存在）を修正する作業は大変なものです。

（４／２２、１９９９）現在、修正作業を行なっています。プログラムはいろいろなディレクトリに散らばって存在するので、作業はやはり大変です。

今回の教訓としては、以前のリポートでも書いているのですが、”おかしい”、” 何か釈然としない”、”計算結果がどうもしっくりいかない”と思ったら甘い解釈はしないで、徹底的に調べてみるというものです。

（Ｃ）本日のバグ（レポート１／２０、１９９９） [目次]

今回はよくやる典型的な間違い（バグ）です。

以下のルーチンの部分に間違いがありました。


      DO 2000 JBAN=1,KEG
       DO 2010 IBAN=1,KEG
       IF (IBAN.EQ.JBAN) THEN
        DO 2100 J=1,IIBA
         DO 2110 I=1,IIBA
          I1=NBASE(I,NNN)
          J1=NBASE(J,NNN)
          CALL HPDF(HP,HD,HF,VP,VD,VF
     &             ,IA,NSEK,I,J)
          Z(IBAN,JBAN)=Z(IBAN,JBAN) + EKO(IBAN,NNN) 
     &  + DCONJG(ZZZ(I,IBAN,NNN))*RAA(I,NNN)*QC(I,J)*RAA(J,NNN)
     &   *HP*ZZZ(J,IBAN,NNN)*PC
     &   *ZFM2(I1,IA)*DCONJG(ZFM2(J1,IA))
 2110 CONTINUE
 2100 CONTINUE
 2010 CONTINUE
 2000 CONTINUE

上記ルーチンで、IBAN,JBANはバンド（インデックス）に関してのループ、 I,Jは平面波数に関してのループで、EKOはバンド計算におけるエネルギー固有値です。そしてこのEKOを配列変数Zに足し込んでいるのですが、既に分った人もいると思いますが、このEKOの位置に問題があります。このルーチンでは、バンド数×バンド数の行列Zを作り、それを最終的に対角化（固有値問題、Zに対する固有値〔先のバンド計算による固有値とは別物〕を求める）することを目的としていました。対角項にはバンド計算によるエネルギー固有値も和の一部として含んでいます。

ところが、実際対角化して得られた固有値は、とんでもなく大きな値になってしまいました。原因を突きとめるため、上記ルーチンのHPの値をゼロにしてみましたが、状況は変わりませんでした。更に、以下のようにルーチンを変更しても状況は依然変わりませんでした。


      DO 2000 JBAN=1,KEG
       DO 2010 IBAN=1,KEG
       IF (IBAN.EQ.JBAN) THEN
        DO 2100 J=1,IIBA
         DO 2110 I=1,IIBA
          I1=NBASE(I,NNN)
          J1=NBASE(J,NNN)
C          CALL HPDF(HP,HD,HF,VP,VD,VF
C     &             ,IA,NSEK,I,J)
          Z(IBAN,JBAN)=Z(IBAN,JBAN) + EKO(IBAN,NNN) 
 2110 CONTINUE
 2100 CONTINUE
 2010 CONTINUE
 2000 CONTINUE

この段階でどうしてだろうと頭を抱えそうになったのですが、他のループで ZにEKOの値のみを代入して（上記ルーチンは動かなくする）対角化すると正しい結果が出ました。そのルーチンは、


      DO 272 I=1,KEG2
        DO 273 J=1,KEG2
          Z(J,I)=CMPLX(0.0,0.0)
          IF (I.EQ.J) THEN
            Z(J,I)=EKO(I,NNN)
          END IF
  273   CONTINUE                                                        
  272 CONTINUE

としていました。ようやくこのルーチンと、先のルーチンを比較して筆者は間違いに気付きました。問題は、


        DO 2100 J=1,IIBA
         DO 2110 I=1,IIBA
          I1=NBASE(I,NNN)
          J1=NBASE(J,NNN)
C          CALL HPDF(HP,HD,HF,VP,VD,VF
C     &             ,IA,NSEK,I,J)
          Z(IBAN,JBAN)=Z(IBAN,JBAN) + EKO(IBAN,NNN) 
 2110 CONTINUE
 2100 CONTINUE

で、平面波数の２重ループ内において、エネルギー固有値EKOを足していました。つまりバンド×バンドに関しての行列の１つの対角要素に平面波数×平面波数分のEKOの和が代入されていた訳です。これでは対角化による固有値がとんでもなく大きな値（発散してもおかしくなかった）になるのも当然です。この段階になるまで気付けなかったのがちょっと情けないです。

ただ、この手の和をとるとき、とるべき変数の位置を間違えるのは、よく犯す間違いの典型例と言えます。間違い方としては、上記のように和をとる必要のないループにその変数を置いてしまう場合と、和ととるべきループ以外の場所にその変数を置いてしまう場合などです。前者では、大抵その変数は、その和のループ内では定数和と同等となり、非常に大きな絶対値（大き過ぎれば発散してエラーになる。今回は発散するまでには至らなかった）になって後々悪さをします。

このように、和をとる時の変数位置の間違いは筆者も過去何度も経験しているのですが、またやってしまいました。変数の位置のチェックは怠らないようにしましょう。

（Ｂ）本日のバグ（レポート１１／１８、１９９８） [目次]

今回も単純なミスです。

新しい系の計算をしようとして、擬ポテンシャル、初期設定用座標データなどを用意して、テスト計算していたら間違っちゃいました。

最初、テスト計算で計算が４、５回イタレーションを回した段階で破綻（全エネルギー発散）しました。この時点では、よくあることだと思い、初期設定（波動関数更新の時間刻み、電荷密度混合比、バンド幅など）を変えて再度試してみましたが、症状が一向に変わりませんでした。

この計算では、２種類の原子から成る系を計算していたのですが、これを１種類のより分かりやすいものに替えて、比較のための計算を行なってみました。これはデバッグの仕方ページの項目４に該当します。結果として、１種類での計算は破綻せずに、うまくいく（少なくとも４、５回程度での全エネルギー発散はない）ことが判明しました。

そこで、再度入力データ（擬ポテンシャル、入力座標、配列宣言部分、ファイル入出力指定箇所など）を調べ直しました。そして、見つけてしまいました。擬ポテンシャルデータに問題があったことが判明しました。

この擬ポテンシャルデータは、前の計算で使い回していた擬ポテンシャルデータに新たに必要な擬ポテンシャルのデータを付加して再使用していました。筆者は元のデータには１種類の擬ポテンシャルデータしか入っていないと思っていました。ところが実際は、元のデータには２種類のデータが格納されていました。それでも前の計算では、最初の１種類目のデータしか読み込まないようになっていたので（１種類の系の計算をしている、そのまた前の計算で２種類の系での計算をしていました）、その段階では２種類分あることを認識していたのですが、そのままにしていました。そして今回、新たに擬ポテンシャルデータを付加する時には、すっかりそのことを忘れていました。このため、新たな計算では、２種類めの擬ポテンシャルデータは新たに付加した正しいものではなく、元からあった全然異なる擬ポテンシャルデータが読み込まれてバンド計算が行なわれてしまった訳です。

これでは、計算がうまく行くはずがありません。

本日の教訓は、”入力データのチェックはしつこいほどちゃんとしよう。” です。

（Ａ）本日のバグ（レポート１１／２、１９９８） [目次]

今回のバグ（と言うより失敗、間違いの類）も単純です。

表面系の計算で、格子振動の計算をする必要があり、初めて格子振動（フォノン）の計算をしてみたところ、いきなり失敗してしまいました。
原因は、格子振動を求めるために与えた、原子位置のずれのオーダーを間違えたためでした。

筆者は、以下の入力データの座標部分の元データを、


   1000  0.9500  0.1000D-07  6.5000   0
       65.8112000000        0.0000000000        0.0000000000
        0.0000000000        5.8169400000        0.0000000000
        0.0000000000        0.0000000000        5.8169400000
           0 COORDINATES  0:NORMALIZED  1:CARTESIAN 
        0.0061633466        0.0000000000        0.0000000000
        0.0663263417        0.5000000000        0.5000000000
        0.1273274737        0.0000000000        0.0000000000
        0.1885776509        0.5000000000        0.5000000000
        0.2499996946        0.0000000000        0.0000000000
        0.3114217229        0.5000000000        0.5000000000
        0.3726719932        0.0000000000        0.0000000000
        0.4336731733        0.5000000000        0.5000000000
        0.4938364117        0.0000000000        0.0000000000
        0.0022391446        0.5000000000        0.5000000000
        0.0652380151        0.0000000000        0.0000000000
        0.1266380629        0.5000000000        0.5000000000
        0.1883738854        0.0000000000        0.0000000000
        0.2499996696        0.5000000000        0.5000000000
        0.3116255713        0.0000000000        0.0000000000
        0.3733614129        0.5000000000        0.5000000000
        0.4347616726        0.0000000000        0.0000000000
        0.4977604284        0.5000000000        0.5000000000

次のように変更しました。


   1000  0.9500  0.1000D-07  6.5000   0
       65.8112000000        0.0000000000        0.0000000000
        0.0000000000        5.8169400000        0.0000000000
        0.0000000000        0.0000000000        5.8169400000
           0 COORDINATES  0:NORMALIZED  1:CARTESIAN 
        0.0061633466        0.0000000000        0.0000000000
        0.0663263417        0.5000000000        0.5000000000
        0.1273274737        0.0000000000        0.0000000000
        0.1885776509        0.5000000000        0.5000000000
        0.2499996946        0.0000000000        0.0000000000
        0.3114217229        0.5000000000        0.5000000000
        0.3726719932        0.0000000000        0.0000000000
        0.4336731733        0.5000000000        0.5000000000
        0.4938364117        0.0000000000        0.0000000000
        0.0122391446        0.5000000000        0.5000000000
        0.0652380151        0.0000000000        0.0000000000
        0.1266380629        0.5000000000        0.5000000000
        0.1883738854        0.0000000000        0.0000000000
        0.2499996696        0.5000000000        0.5000000000
        0.3116255713        0.0000000000        0.0000000000
        0.3733614129        0.5000000000        0.5000000000
        0.4347616726        0.0000000000        0.0000000000
        0.4877604284        0.5000000000        0.5000000000

赤色で示したデータの値を0.01ずらし（た値は深紅色で表示）ました。これはスラブモデルにおいて、表、裏の各表面のトップレイヤーの原子を構造最適化で得られた平衡位置からずらしていることを意味しています。そして、この段階で大きな間違いを犯していました。ここでの0.01を筆者は十分に小さなずれと思っていました。しかし実際は全然小さなものではありませんでした。

因みに格子振動の計算のための参考論文として、M. T. Yin and M. L. Cohen, Phys. Rev B26, 3259(1982)を挙げておきます（特に３２６４頁）。そこには原子位置のずれの大きさは、0.01から0.1Aとありました。筆者は、先の0.01のずれを、0.01 a.u.（大体0.005A）と誤認してしまいました。これは重大な過ちで、ここでの0.01は、表面のスーパーセルの長さ65.8112 a.u.に対する0.01だったのです。つまり実際のずれは、0.658112 a.u.（大体 0.35A）にもなっていました。更に筆者は、同じようにして0.02ずらした場合（0.7Aものずれ）の計算を行なっていました。いくつかの表面系の計算をしていて、系によっては0.02ずらすと、計算が破綻（収束しない）している場合もありました（この破綻をみて、どうもおかいしと思って、よくよく考えると気付いた次第です）。

今日の教訓は、”データ入力は慎重に行ない、十分その意味するところも把握し、何度も検討（検証）ておくべきである。”です。

格子振動の計算もまだ初めたばかりで、正直まだ本当にこれで正しく求められかも、今後検証していかなければならない。

（９）本日のバグ（レポート８／４、１９９８）[目次]

今回はバグではなく、単なるミス（失敗）です。

筆者は単純な入力データの間違いで、これまで計算した大事な結果の１／３を台無しにしてしまいました。

（経緯）
遷移金属の計算で、ＰＣＣ版の擬ポテンシャルの多くにゴーストバンドがあることが判明、これを修正すべく新たにＰＣＣを考慮した遷移金属擬ポテンしゃルを作成しました。ほとんどのものは、局所部分（Vcore(r)）の修正だけで済みました。これは遷移金属ではｓ、ｐ、ｄを非局所とし、Vcore(r)（ＢＨＳ論文で定義されている通りのもの）を記述するパラメーターを、例えば 1.0,0.5,0.5だったのを0.5,0.5,0.5にするものでした。

これを遷移金属炭化物表面の計算をしている途中で気が付き、擬ポテンシャルを修正、バルクの計算までは正しく修正したデータを使用していたのですが、肝心な表面計算用の入力データを修正するのを忘れてしまいました。実は明日（８／５）セミナーがあり、この遷移金属炭化物表面の話をするのですが、完全に真っ青な状態です。

間違いは、0.5とすべき数値を、元の1.0（旧ポテンシャル用）のままで計算してしまったことです。普通、入力データを間違えると、計算が破綻したり、収束しなかったり、明らかに間違った結果が得られるのですが、今回はもっともらしい結果が出てきたので、かなりの時間問題に気付きませんでした（間違いもあまりに単純だったこともある）。気付いたときは正に「後の祭り」状態でした（放心状態）。ただ、後になって結果を良く見てみると確かに若干変なところがありました（事前に気付けるほどではない）。

今日の教訓は、「データはくどいほどチェックすべきである。」「誤り、間違いは単純なほど気付きにくく、そして致命的。」です。

（補足、追加、８／１７、１９９８）
この一件に関して、住友電工の澤村先生から指摘（感謝）がありました。それは、擬ポテンシャル内に、その元素名、価電子数、作成時のパラメータ（切断半径、Vcore(r)のパラメータ）及び作成方法、条件の簡単な説明、使っている単位系、などなどの情報を埋め込んでおけば良いというものです。

こうすれば、バンド計算時に、バンド計算として設定した計算条件と、擬ポテンシャル内の内部情報を比較してチェックするようにしておけば、もし互いの情報に矛盾があれば、すぐに修正して計算をやり直せます（計算の正確性が増す）。これは単純な間違いや、一体何の計算をしているか分からなくなる事態を回避するために、大変有効だと考えられます。

（更に補足、８／１８、１９９８）
何と、まだ間違いがありました。それは価電子数の指定のところで、例えば、 Nbは５価（４ｄに４電子、５ｓに１電子）として擬ポテンシャルを作ったのに、バンド計算では４価として計算しているというものです（それでも計算がちゃんと収束しているから怖い）。これも大変単純な間違いであり、改めて、上記擬ポテンシャルの内部情報埋め込みの必要性を認識する次第です。

（またまた補足、９／３、１９９８）
上記の「更に補足」に関連して、各原子の価電子数が分かると、計算すべき系の総電子数が分かります。スピンを考えなければ、１バンドに２個の電子が入り（詰められ）、バンド計算において必要なバンド数が大体分かります。この大体の意味ですが、バンド数が必要な数より少ないと、計算は普通破綻してしまいますが、ではどのくらい余計にとっておく必要があるかが問題になります。

ダイヤモンド（バルク、ダイヤモンド構造でユニットセル内部に４個の場合）のようにギャップの大きな非金属系では、扱うバンド数＝系の全電子数／２としても良いのですが、これは非常に限定された話で、大抵の（特に金属的な）系では、扱うバンド数＝全電子数／２という訳にはいかないです。

系が金属的なバンド構造になる場合、ある程度多めにバンドをとっておく必要があります。例えば、Ta（５価）、C（４価）がそれぞれ９個づつから成る系（表面系）では、合計で８１個の電子が存在し、最低でも４１バンドが必要となります。
この場合、全電子数が８１なので、考えるまでもなく系は金属的にならざるおえなくなっています。では、それぞれ１０個ずつで９０電子としたらどうなるのでしょうか？。これは即答はできません。計算してみないことには分かりません。これは難しい問題で、本当にこの表面系が金属的になっているかどうかは、Ta、Cを９個づつでは判定できません。検討の必要あり。ただし、TiとCの同じ表面系では、Tiは４価なので、全電子数は７２個となり、全てのバンドは詰まっていることになりますが、構造最適化した場合のこの系の電子状態は金属的となります。

実際の計算では、４１バンドでは不足で、より上のバンドまで考慮する必要があります。ではどこまでとる必要があるかですが。これはあまりはっきりしません。以前、扱う系によっては大分上の方のバンドまで考慮に入れないといけないという論文（詳細は失念、いわゆるカー・パリネロ的な電子状態計算過程で必要なのか、部分対角化の時に特に必要であるのかも不明）を目にしたことがあります。筆者も、経験的にバンドは、十分に余裕をもって多めに取る必要があると判断しています。どのくらいかというと、これも場合場合によるのですが、先のTaとCから成る系（表面系）では、どうも４４バンド程度では駄目（どうも力の計算がおかしくなる）で、５１バンド位は必要なようです。但し、４５、４６、４７、、、と全ての場合でチェックしていないので、どこからが安全、不安の境界であるのかは不明です。また系が金属的か絶縁体か、バルクか表面かなどによっても、バンド数の取り方は変わってしまいます（どこが第一原理なのだろう^^;）。

従って、事前に擬ポテンシャルに価電子数の情報を盛り込んでおけば、実際のバンド計算で（最低）必要なバンド数がわかる訳で、バンド計算時にバンド数に問題があれば（少ない過ぎる、足りないかもしれないなどと）警告するようにできます（更に、実際は状況によりますが上記のような理由により、ある程度多めのバンド数で計算を行なう必要があります）。

（またまた補足、８／２０、１９９８）
擬ポテンシャルデータ内に埋め込む情報について列挙してみましょう。

元素名
切断半径
原子での電子配置
価電子数
擬ポテンシャルの種類（BHS, TM, Vanderbilt, etc）
作成日
作成者
作成場所
ＬＤＡの型（Wigner, Perdew-Zunger, MJW, 相対論〔考慮、非考慮〕、 etc）
ＧＧＡを考慮しているか（考慮している場合の型）
ＰＣＣ（部分内殻補正）を考慮しているか
どこまでを価電子と考えているか（価電子数と連動）
参照エネルギーの数（普通は一つ）
バンド計算時に（最低）必要なエネルギーカットオフ（※）
バンド計算時の計算状況（収束しやすいか、時間刻みの上限等）（※）
可能ならば、原子の計算での全エネルギー
原子の計算でのスピンポラリゼーションエネルギー
スピン起動相互作用対応かどうか
その他（注意事項など）

（※）この情報は、少なくとも一度はバンド計算を行なって、確かめておく必要があります。（おそらく一度の計算くらいでは、確実なデータは出てこないので、ある程度の計算経験を経てから付加できるようにしておくと良いと思われます。）

（８）本日のバグ（レポート３／１７、１９９８） [目次]

今回のバグも単純です。

プログラムの拡張を行なっていて、ヨウ素のfcc構造でのテスト計算を行なっていて、予想された結果が得られず（格子定数がおかしい）思案していました。行なった拡張部分は、最初、


C     PBE                             
      RGGAEX(I)= (EXPBE + ECPBE)
      VGGAEX(I)= (VXUPPBE + VCUPPBE) + (VXDNPBE + VCDNPBE)

というものでした。もうこれを見て何を拡張しようとしているかお分かりの方々もいるかと思います。

そうです、J. P. Perdew先生のＰＢＥコードの導入を行なおうとしていた訳です。ＰＢＥの詳細は参考文献を見てください。また、Rutgers大学のKieron Burke先生のサイト[ページ]からＰＢＥソースコードの入手が可能です（取り扱いについては、付随するドキュメント等の使用上の注意、指示に従って下さい）。

で、上のソースリスト部分は、ＰＢＥで計算されたＬＤＡ＋ＧＧＡ（ＰＢＥ）の結果（EXPBE：交換エネルギー、ECPBE：相関エネルギー、VXUPPBE：スピンアップ交換ポテンシャル、VCUPPBE：スピンアップ相関ポテンシャル、VXDNPBE：スピンダウン交換ポテンシャル、VCDNPBE：スピンダウン相関ポテンシャル）を、こちらのプログラム側で対応する変数（RGGAEX：交換相関エネルギー、 VGGAEX：交換相関ポテンシャル）に代入しています。

この場合、ＰＢＥは最初からスピンがアップ、ダウンの別を考慮した形になっていたのですが、こちらの計算では取り敢えず最初はパラの状態で考えることにしました。つまりスピンアップ側の電荷密度＝スピンダウン側の電荷密度＝全体の電荷密度／２として、計算を行ないました。いろいろ曲折はあったものの、どうにかうまくいっているのではないかと思えるところまできたのですが、どうしてもＰＢＥ（ＬＤＡにしても）の計算での格子定数が正しく出ませんでした。

どうしようもないと思っていたのですが、同僚の新井さんに相談してみるといとも簡単に問題点を指摘されてしまいました。
問題は、ポテンシャル代入部分、

VGGAEX(I)= (VXUPPBE + VCUPPBE) + (VXDNPBE + VCDNPBE)

で、アップとダウン両方足し込んでは駄目だと言われました。つまり、

VGGAEX(I)= (VXUPPBE + VCUPPBE)

または、

VGGAEX(I)= (VXDNPBE + VCDNPBE)

だけ（どちらか一方）で良いということです。これは、全電荷密度をP、アップ側をPU、ダウン側をPDとして、パラ（常磁性）では、

P = PU + PD

PU = PD = P/2

V(up) = V(down) = V

V(up)*PU + V(down)*PD = [V(up)*P + V(down)*P)]/2

= 0.5*2*V*P = V*P = V(up)*P = V(down)*P

となるため、スピンを考えない時の交換相関ポテンシャルは、アップ側、ダウン側の一方だけで良いということになります。

ここでの問題は、筆者の知識（勉強）不足であり、教訓はデバッグの常套手段にもあるように、人に相談（議論）することはデバッグにとって非常に重要であるということです（新井さんに感謝）。

ＧＧＡ（ＰＢＥ）に関しては、[拡]張編にておいおい話ていく予定です。

（７）本日のバグ（レポート１／２０、１９９８） [目次]

今回はプログラム上のバグではなく、研究上の失敗です。

Pdの擬ポテンシャルは現在、格子定数が実験値より2.6%も長くなってしまっています。これを改善しようとして、切断半径や原子での電子配置などを変更して、新たにPdの擬ポテンシャルを構築してみました。

そして格子定数が3.89Aと実験値とどんピシャリな擬ポテンシャルを作成することに成功しました。バンド構造（pngファイル、9KB）も表示してみてゴーストバンドは存在しないように見えました。

ただちょっと一番上のバンドが平たいなとは思い、ちょっと引っかかるものが残りました。しかし、この時格子定数が思った以上に良く出ていたことに気を取られ、これのもつもっと深い問題を見逃してしまいました。

そして新しいPdの擬ポテンシャルができ、これを今度配布するCD-Rにも収録しようとしてましたが、どうもバンド構造をみる毎に、一番上の部分的に平たくなっているバンドが気になり、一応もう少し調べてみることにしました。

そして、より上のバンド（pngファイル、 9KB）まで描いてみました。困ったことに出てきてしまったのです。ゴーストバンドが、、、

呆れるほどはっきりとした、（奇麗な？）ゴーストバンドがフェルミエネルギーの上5eV付近にあるじゃないですか（ショック！）。前のバンドは、６バンド分しか表示しなかったため、一番上の６バンド目が丁度ゴーストになっているのですが、これは上の表示していない７番目のバンドと重なっていて、平たい明らかにゴーストと分かる部分を、その７番目のバンドが担っていました。そのため、６番目の一部のみが平たくなっているだけで、最初このバンドの持っている致命的な問題に気付きませんでした。

更に、ゴーストバンドがフェルミエネルギーより上だったのが発見を遅らせました。フェルミエネルギーより上だったため、あまり（と言うよりほとんど）格子定数に悪い影響を及ぼしませんでした。ゴーストさえなければこの擬ポテンシャルは非常に良いものであったかもしれません。

ゴーストバンドがフェルミエネルギーよりずっと上（数Ｒｙ以上）に出てくるのなら全く問題無しとしてしまうこともできますが、5eVでは小さ過ぎます。単独ならば良いかもしれませんが、化合物や表面（とそれに原子が吸着する問題なら尚更）系ではゴーストの存在による悪影響が顕在化する可能性が高いです。従って、この新しい擬ポテンシャルは使えません。

今日の教訓は、非常に良い（都合の良い）結果だけに着目して、悪い部分や、より厳密なチェックを怠る（ポテンシャル作成時にログデル位は見ておくべきだった。ただみなくても７バンド分表示させればゴーストの存在は明白だった）ことは、研究上非常に危険であるということです。

因みにゴーストの無い正しいバンド構造（pngファイル、9KB）を示します。これも６バンド分しか表示していませんが、表示している中で一番上のバンドは、ゴーストのある場合と異なり、両端が平たくなっていません。当然、その上にもゴーストバンドは（ずっと上は不明ですが）ありません（但し、この場合は格子定数が良くない）。

（６）本日のバグ（レポート１２／５、１９９７） [目次]

今回のバグも教訓的です。

あるプログラム（自作ではなく、公開されている他の人のもの）をDEC上で動かそうとして、コンパイルは全く問題なく終了したのに、いきなり実行させると、

forrtl: error: floating divide by zero IOT trap

で止まってしまいました。メッセージから筆者は単純なゼロ割りで止まっているのだと判断しました（実は違った）。

但し、これは自分が作ったプログラムではなかったので、一体どこの箇所で止まっているのか即座には見当できませんでした。更に、DECのマシン上では、この時エラーで止まった箇所が示されていません。デバッガーにかければエラー箇所が分かるだろうと思ったのですが、f77 -g でコンパイルし、デバッガーにかけると、何か訳の分からないメッセージが沢山出るのですが、肝心なエラー箇所の表示がなく、この段階で筆者は大変困ってしまいました。デバッガー（gdb）が表示したメッセージは、


Program received signal SIGFPE, Arithmetic exception.
0x3ff808404f4 in __dpml_signal () at ref_cp/alpha_osf_exception.c:317
ref_cp/alpha_osf_exception.c:317: No such file or directory.
(gdb)

ですが、筆者には何を言っているのかさっぱり分かりません。そこで、デバッガーによる問題解決の道を諦め、別の手で行くことにしました。まずデバッグに関する常套手段を基本路線とします。但し、これは人の作ったプログラムで、この場合１、の正しい結果を与える標準的な状況が存在しません（これはテスト用のプログラムで、入力データはプログラム内部で全て設定するようになっており、その場合の正しい出力結果に関しての情報は分かっていました）。

デバッガーが駄目で、コンパイラのオプションを変えても駄目なので、やはり、プログラム内に多数のチェックポイント（WRITE文）を置いて、しらみ潰しに探すしかないと思っていたのですが、その前に６、の他のシステムで試してみようと思い、他のワークステーション（SUN、HPなど）で試してみたところ、SUNやHPでは正しく動作することが分かりました。ただ正しく動く場合、何のエラーメッセージも出てこないので、これではDECの場合のエラーは、どうして、そしてどこで起こっているのかわかりません。ただ、これは言語仕様、コンパイラの仕様、実行時の判断の細かいところの差によるのではないかと推定できました。

更に思案する内にSUNのマシンにはSUN標準のフォートランコンパイラと、富士通のフォートラン９０コンパイラの二つがインストールしてあり、まだ富士通の方は試していないので、これを試してみることにしました。そして、SUN の標準のコンパイラでは正しく稼働したのですが、富士通のコンパイラでは正しく稼働しないことが判明しました。更に幸運なことにDECの場合とは異なるエラーメッセージが出てきました。

そのエラーメッセージは、


jwe0266i-e In dx1**dx2 (dx1,dx2:real*8),dx1.eq.0.0 .and. dx2.le.0.0 (dx2=-0.3333333333330000d+00).

と言うもので、エラーの発生したサブルーチンとその発生箇所（何行目か）もはっきり示されていました。このエラーの意味は、dx1**dx2（dx1をdx2乗）する時、dx1の値はゼロで、dx2の値が負になっているというものです（富士通製のコンパイラで、日本語のマニュアルも助けになった）。これは、（もう数学というより算数は忘れたの自信がないです）1/(dx1**|dx2|)になっていて、確かにDECのエラーメッセージの言う通りゼロ割りなのですが、正直もっと親切な表現にして欲しいです（エラー箇所が出てこないのは、やはりおかいしい）。

結局この富士通フォートランでコンパイル、実行した時のエラーメッセージで原因はすぐに分かりました。それは、

A = ((1.0D0 + Z1)**Z2 - (1.0D0 - Z1)**Z2) ：Z2は負の定数

なっている部分で、Z1が1.0D0または-1.0D0の場合、0**Z2, Z2 < 0という状況に陥り、エラーになっていることが判明しました。そこで、場合分けして、 Z1が-1.0D0または1.0D0の時にゼロ割りにならないように書き換えると問題なく稼働するようになりました。

ここでの問題は、SUN、HPでは動いたものが、DEC、富士通のフォートランでは動かないことです。多分、ゼロ割りの時のエラー処理等の制限が甘い、厳しいの差と思われるのですが、何とかして欲しいものです。

（５）本日のバグ（レポート１０／１５、１９９７） [目次]

久々に愚かな間違いを犯してしまいました。

前のバグ（４）に関しては、まだ納得していない部分もあり、現在もパラメーター文に関しての調査を不定期的に続行中です。

さて、今回のバグは、サブルーチンDIAGON において、新しい対角化のサブルーチンを導入しようとして遭遇したバグです。
この新しいサブルーチンは、文献「Fortran77による数値計算ソフトウェア」、渡辺力、名取亮、小国力監修、丸善刊の、エルミート行列対角化ルーチン EIGCH（当然、これは公開できません）です。

まず筆者が、この書籍に添付されていたフロッピーから、EIGCHを引っ張ってきて、計算プログラム本体revpe_d.fにくっつけ、DIAGON 側でCALLするようにしました。
ところが、計算されて出てくる値（全エネルギー）は、全く正しくないものでした。

原因は、サブルーチンのCALL EIGCH(----)で使っていた、引き渡し用の変数の設定を間違ったためでした。

このサブルーチンの正しいCALLは以下のようになります。


      CALL EIGCH(ZZZ,IIBA,KNG11,-KGB,KGB,EPS,WWE,LW,EG,ZVN,ICON)   正

ところが最初、筆者は、次のようにして実行していました。


      CALL EIGCH(ZZZ,IIBA,IIBA,-IIBA,IIBA,EPS,WWE,LW,EG,ZVN,ICON)  誤

ここで、ZZZが対角化すべき配列（エルミート行列）、IIBAは平面波数（ｋ点毎に異なる値になる）、EPSは打ちきり誤差、WWE、LWは作業用配列、EGがエネルギー固有値用配列、ZVNが固有ベクトル（波動関数）用の配列、ICONがエラーコードとなります。

間違いは、３番目の変数IIBAです。文献の中での説明には、ちゃんと配列 ZZZの行数とあります。そこを筆者は何を間違えたのか、可変な変数IIBA（ｋ点に関するループ内で値が変わる）を引き渡してしまっています。このため、配列のメモリーか何かが壊れた（？）ため、全エネルギーがおかしくなったと考えられます。

この３番目の変数をIIBAから、正しい定義による変数KNG11にすると、ちゃんと正しい結果を与えることがわかりました（４、５番目の変数はIIBAでも KGBでも大勢に影響ありません。但し、配列の定義部分も連動して変えておく必要があります）。

今日の教訓としては、サブルーチンで引き渡している変数の定義、説明はちゃんと読んで、理解しておかなければならないことです。変な思い込みで、適当なことをやるとろくなことになりません。

筆者は、このバグに気付くのに３日もかかりました^^;;)。

（４）本日のバグ（レポート６／１８、１９９７） [目次]

これはこれまでとは、ちょっと趣を異にしたバグです。

最近、新しいＤＥＣのサーバーマシンが導入され、最新のＯＳと最新のフォートランコンパイラがインストールされていました。バージョンはそれぞれ、ＯＳが4.0A、フォートランコンパイラ（FORTRAN77仕様）が4.0でした。
ところが、自分のプログラムが、このバージョン4.0のフォートランコンパイラでは、正常に動作しないことが分かりました。コンパイルは無事に成功し、エラー表示も警告以上のものは出てきません（エラーの出方そのものも、従来正しく動いておた場合と全く同じものでした）。

そして、実行用オブジェクトを実行させると、思いもしないところで止まってしまいました。エラーはsegmentation faultというもので、最も何を言っているのかさっぱりわからないエラー表示の代表格と言えるものです。このエラー表示が出たときは、メモリ領域を壊して（何らかの不正を行なっている）いることが経験的に多いです。メモリーを壊すのは、大抵配列の領域外を参照しようとした場合など、配列に関わる扱いのところで壊し易いです。

しかし、この場合は従来のバージョン3.8以前では、このようなエラーでプログラムが停止することは全くありませんでした。更にＤＥＣ以外の機種でも一応問題なく動いています（富士通ＶＰＰ、富士通ＰＣ用フォートラン、ＨＰ、ＮＥＣのＳＸ４など）。では原因は一体何なのでしょうか？、ＤＥＣ側のコンパイラのバグも疑ってみましたが、問題はやはり筆者のプログラムにあることが（ＤＥＣ側の指摘により）わかりました。

問題は、パラメーター文で、設定（定義）した変数を、サブルーチンの引数として渡しており、引き渡されたサブルーチン側では、変数の名前が変わっていて、更に、それに値を代入して演算を行なっていました。ＤＥＣ側によると、新しいコンパイラーでは、チェックが厳しくなっており、このように実質的にパラメーター文に代入を行なっているような場合には、エラーとして検出されるようになり、segmentation faultエラーで、プログラムが止まったと言われました。

この状況を分かりやすくプログラム的に示すと以下のようになります。


  MAIN PROGRAM
   IMPLICIT REAL(A-H,O-Z)
   PARAMETER (KNVX=8)
   -------
   CALL SUB1(KNVX)
   -------
  STOP
  END
C
  SUBROUTINE SUB1(NKX)
   IMPLICIT REAL(A-H,O-Z)
   -------
   CALL SUB2(NKX)
   -------
  END
C
  SUBROUTINE SUB2(NKX)
   IMPLICIT REAL(A-H,O-Z)
   -------
   NKX = NKXH + NKXH  <--- エラー 
   -------
  END

確かに、この場合おお元はパラメーター文で定義されています。但し、この変数はサブルーチンに渡される段階で、変数名が変えられています（そして、そのサブルーチン上でパラメーター文で定義されている訳ではない）。本当に厳密なFORTRAN77の仕様が、このような扱いを禁止しているのか、少々疑問もあるので、現在これに関して調査中です。
（８／１４）確かに、パラメーター文自体は、一度定義した後に、定義し直すことができません。つまり、ルーチン内で、定義されたパラメーター変数に値を代入したりすることはできません。しかし、サブルーチンの仮引数として、サブルーチンに渡されたあとまでの、厳密な定義はあるのでしょうか？（--->識者の人々）
ここでの場合、サブルーチンに引き渡されたパラメーター変数は、そのサブルーチン内で、普通の（整数）変数として再定義されていると解釈できないでしょうか。そうなると、これはエラーではないような気がするのですが？（尚も、調査中）

ＤＥＣのフォートランコンパイラも古いバージョンでは、このエラーは引っかからなかったのですが、新しいコンパイラではエラーになるようになってしまいました。ＤＥＣ側の言い分としては、コンパイラの解釈がより厳密になったとしていますが、ユーザー側からみると（この解釈が正しいとしても）「ありがた迷惑」のような気がします。

このエラーに対する回避手段として、以下のものが考えられます。


  MAIN PROGRAM
   IMPLICIT REAL(A-H,O-Z)
   PARAMETER (KNVX=8)
   -------
   CALL SUB1(KNVX)
   -------
  STOP
  END
C
  SUBROUTINE SUB1(NKX)
   IMPLICIT REAL(A-H,O-Z)
   -------
   CALL SUB2(NKX)
   -------
  END
C
  SUBROUTINE SUB2(NNKX)
   IMPLICIT REAL(A-H,O-Z)
   -------
   NKX = NNKX <--- NNKXには代入などの演算は行なわない
   NKX = NKXH + NKXH
   -------
  END

これは、ＤＥＣの最新のフォートランコンパイラでも正しく動作することを確認しました。
（注意）但し、現在（１／２６、１９９９）のＤＥＣ（今はＣＯＭＰＡＱ）フォートランコンパイラでは、上記変更では正しく動かないです。

或いは、


  MAIN PROGRAM
   IMPLICIT REAL(A-H,O-Z)
   PARAMETER (KNVX=8)
   NKX = KNVX <--- 変更点
   -------
   CALL SUB1(NKX)
   -------
  STOP
  END
C
  SUBROUTINE SUB1(NKX)
   IMPLICIT REAL(A-H,O-Z)
   -------
   CALL SUB2(NKX)
   -------
  END
C
  SUBROUTINE SUB2(NKX)
   IMPLICIT REAL(A-H,O-Z)
   -------
   NKX = NKXH + NKXH
   -------
  END

でも問題ないはずです（これはまだ動作確認をしていません）。
現在（１０／３１）この方法を試してみたところ、テスト用の小規模な計算ですが、うまくいくことを確認しました。
現在（１／２６、１９９９）でも正しく動くことを確認。

ただし新たな問題が出てきました（１０／３１）。まず、このバグに対する最初の処方では、ＤＥＣ上のコンパイルとして、

f77 -r8 -O5 sample.f

では正しい結果を与えますが、

f77 -r8 -O4 sample.f

では、実行中segmentation faultエラーで止まってしまうことが判明しました。コンパイル時にはこれに関係しそうなエラーメッセージは出てきません。また、最適化オプション-O2,-O3,-O4の場合はsegmentation faultエラーで止まりますが、オプション-O0,-O1,-O5の場合は正しく動くことが確認できました。ここで、注目されるのは、最も最適化した場合の、-O5では正しく動いてしまうことです（普通、最適化の度合いを上げれば、エラーで止まることが多くなるのに、今回はその逆になっています）。

この場合、今後も-O5のオプションによる最適化状態で、より実際的な計算を継続していって良いものかどうか不安です（まだ試していませんが、より大規模な計算になると-O5オプションでは動かなくなる可能性があります）。当然、より古いバージョンのＤＥＣフォートランコンパイラではこのようなことは起きません。

同じく、-g0,-g1,-g2（デバッグモードに関してのオプション）では正しく動きますが、-g3オプションでは動きません。動かない場合に対して、デバッガーを使ってみましたが、訳のわからない（どうみてもエラーでない）箇所で止まっています。どうやら、メモリー領域か何かを壊していて、見当違いなところでエラーになっている（バグとしては最悪のパターン）ようです。
後の解決方法（処方）でも事情は全く同じで、-O5では動きますが、-O4では動きません。

どうもやっぱり、ＤＥＣの新しいバージョンのコンパイラの方にも何か問題がありそうで、本当の原因究明と根本的な解決のための調査が必要（でも時間がない）と考えられます。

正直な話（２／９、１９９８現在）、次のコンパイラのバージョンアップでの改善を期待していた（今までうっちゃってあった^^;）のですが、ＤＥＣが何とコンパックに身売りすることになってしまいました。Degital UNIXの将来に一抹の不安を憶えます。
（６／３０、１９９８）現在、ＤＥＣのフォートランコンパイラのver4.0で駄目でしたが、最近、最新のバージョンver5.0でもやはり駄目であることが判明しました。どうもコンパイラの仕様そのものが、上記のような記述を認めないという方針になってしまったようです。せめてオプションで回避できるようにしてくれたらと思う次第です。

実は、ＤＥＣ側に別に最新である必要もなかったので、ＯＳはそのままで、フォートランコンパイラのバージョンを従来のプログラムで問題無く動く ver.3.8にしてもらいました（これも前向きではないですが解決策の一つ）。

このバグの問題は、もともとパラメーター文で定義されていた変数に対して、設定変更（代入）を行なおうとしたためなのですが、これにはサブルーチン（SUB1、SUB2）が後から付け足されたもので、もともとメインプログラムと整合性があまりなかったことにあります。それをプログラムの拡張のために無理にくっつけた感があります。改めて、プログラムの設計や拡張は十分に練って行なう必要があることを認識しました（そうでないと、後々面倒になる場合があることを思い知る）。

（３）本日のバグ（レポート３／２６、１９９７） [目次]

最近、プログラム開発がなかったので、あまりバグも出てこなかったのですが、とうとうバグ（というよりプログラムの構造、仕様上の問題）を見つけてしまいました。

まず問題となった部分（計算結果）を以下に示します。


 TOTAL SUMM OP1 = 5.930788159663123e-04 
              2 = 0.0000000000000000e+00 
              3 = 0.0000000000000000e+00 
              4 = -1.251707515445920e-03 
              5 = 0.0000000000000000e+00 
              6 = 8.736301364726182e-04

これは正方格子におけるGaのストレスの値です。１番目、４番目、６番目の値がそれぞれｘ方向、ｙ方向、ｚ方向に対応しています。この値を見てすぐにおかしいと思ったら、貴方はするどいです。

扱った系が正方格子（正確には体心正方格子、ｂｃｔ）なので、明らかに上の３つのストレスの値の２つは値が同じでなければなりません。ところが３つとも異なる値になっています。明らかにこの計算結果は間違っています。

この間違いは、ある条件設定の場合に起こります。そうでない場合は、以下のように正しい値が出てきます。


 TOTAL SUMM OP1 =  -3.386523218593673E-003
              2 =   0.000000000000000E+000
              3 =   0.000000000000000E+000
              4 =  -3.430875185404469E-003
              5 =   0.000000000000000E+000
              6 =  -3.430875185401056E-003

対称性から、ｙ方向の値とｚ方向の値が（計算誤差を除いて）一致しています。では、どのような条件の時に値がおかしくなるかと言うと、ストレスを使ってユニットセルの形についての構造最適化を行なう場合であることがわかりました。単に、電子状態のみの計算で、最後にストレスの計算を行なう場合は上のように正しい値を返しました。

最初、筆者はこれをＶＰＰ上で計算していたので、間違った先入観に因われ、原因は並列化部分ではないかという方針で、デバッグを行ないましたが、全く原因を見つけられませんでした。

しかし、土、日（３／２２、３／２３）と休んでいる時に、家でもう一回正しく出ている場合と正しく出ていない場合の計算出力結果を比較検討してみることを思い立ちました。そして、昨日（火曜日、３／２５）に以下の違いに気付きました。

間違っている方


 ---------- THE FERMI ENERGY = 0.9234210201451082d-01   26.000000
 CALL FORZFB
 REAL TOTAL CHARGE = 25.99999999997447  IN XCFFT
 TOTAL ENERGY FOR 100-TH ITERATION=-131.9255525      -0.1319256d+03
  ITER=100 ET(H)=  0.2962622d-07 ET(M)=  0.000806 DC=  0.3412940d-06
 ---------- THE FERMI ENERGY = 0.9234210813638645d-01   26.000000
 CALL FORZFB 
 TOTAL STR 1  1 = -0.3181689275022663 
              2 = 6.268243606208947e-04 
              3 = 6.268243606208876e-04 
              4 = -0.3180186954074477 
              5 = -6.008948961606797e-05 
              6 = -0.3180186954074479

正しい方


 ---------- THE FERMI ENERGY = 0.2661819904309356D+00   26.000000
 CALL FORZFB
 REAL TOTAL CHARGE =    25.9999999999665      IN XCFFT
 TOTAL ENERGY FOR 100-TH ITERATION=-131.8150778      -0.1318151D+03
  ITER=100 ET(H)=  0.4234662D-06 ET(M)=  0.011518 DC=  0.5330057D-06
 ---------- THE FERMI ENERGY = 0.2661826508252645D+00   26.000000
 CALL FORCE 
 TOTAL STR 1  1 =  -0.526384186291203     
              2 =   1.794706243030863E-003
              3 =   1.794706243030874E-003
              4 =  -0.526410035928503     
              5 =  -4.039471827445321E-005
              6 =  -0.526410035928507

１００回目まで電子状態の計算のみを行ない、１０１回目にストレスの計算を行なうのですが、良く見ると、赤色（色の出ない場合は御容赦下さい）で示した部分、つまりサブルーチンＦＯＲＣＥ（またはＦＯＲＺＦＢ）を呼び出している部分が正しい方と間違っている方で異なることがわかります。

この場合、ＦＯＲＣＥを呼び出すのが正しいのです。ＦＯＲＺＦＢでは正しいストレスの計算が出来ません。間違った方はこのＦＯＲＺＦＢを呼び出し、ストレスの値が正しくなくなってしまているのです（ＦＯＲＣＥでは後で計算するストレスのサブルーチンで必要な配列変数の計算を行なっているが、ＦＯＲＺＦＢではそれを行なっていないため）。

そして、この違いを生ずる根本的な原因はメインプログラム部分にあります。


      IF (MOD(ITER,IMDI).EQ.IRES) THEN
        IF (ITER.EQ.1 .OR. ITER.GT.IMD) THEN                            
          WRITE (6,*) 'CALL FORCE'
          CALL FORCE(IREC8)                                             
        ELSE
          CALL FORZFB
        END IF
      ELSE
        WRITE (6,*) 'CALL FORZFB'
        CALL FORZFB
      END IF

上の赤い部分（赤くなってない人ごめんなさい^^;）の制御変数、IMDIによる条件判定に重大な問題があります。この変数は、原子間に働く力によるユニットセル内の原子の構造の最適化用のもので、ストレス用ではありません。従って２０回おきにストレスの計算をするように設定（CONTL.DATで設定）しても、その制御変数はIOVEであるため、ストレスの計算のために必要なサブルーチンＦＯＲＣＥはいつまでたっても呼び出されません。これでは正しいストレスの計算はできません。

解決策としては、いまのところ２つ考えられます。（案１）上の条件判定部分を以下のように変更する。


      IF (MOD(ITER,IMDI).EQ.IRES .OR. MOD(ITER,IOVE).EQ.IRES) THEN

この方法はまだ試していませんが、プログラムそのものを修正するので、根本的な解決策と言えます。一方、（案２）はCONTL.DATを以下のようにします。


  0
  1.0D0  0.0D0  5.0D0
  0.7D0  0.005D0  300.0D0
 -0.0034D0
  0.0D0
  0.0D0
 -0.0034D0
  0.0D0
 -0.0034D0
  0  0  2  0  12
  100  20  20  1  200  0  200  0
  0

実は、このbct-Gaの系は対称性からユニットセル内部の原子は動きません。従って原子間の働く力を計算する必要がありません。そこで上のデータの最後から２番目の行を 100 200 20 1 200 0 200 0としていたのです。これでサブルーチンは２０１回目（実際の計算は２００回で終るようになっている）にならないと呼び出されません。そこで、上のように原子間に働く力とストレスによるユニットセルの最適化も２０回毎に計算するようにすれば良い訳です。但し、この系では原子は動かないので、原子のための動力学用の時間刻みをゼロに設定してあります。

この方法は、プログラムを修正しないで、入力データの変更だけで対応できるので、作業が簡単ではありますが、ユニットセル内の原子の最適化操作（刻み時間をゼロにしたので、原子は動かないし、この系では対称性から力そのものが出てこない）を行なうので、若干計算に無駄が生じます。

ここでの重大な問題は、これはバグと言うより、プログラム設計、仕様上の重大な欠陥であり、これをこのプログラムを作った当の本人が認識できなかったことにあります（バグが仕様になっているとも、仕様がバグになっているとも言える）。
まさにしようもないプログラムと言えます。

（２）本日のバグ（レポート２／６、１９９７）[目次]

今回のバグはそう深刻なものではありませんが、よくやるバグです。

プログラムを改良中に、あるＤＯループ内の配列変数を、他のＤＯループ内に移動（複写でも同じようにバグるはず）して、コンパイルは何の問題もなく終了したのに、いざ実行してみると計算結果が明らかにおかしくなってしまいました。
原因は簡単で、移動（または複写）した配列ＡがＤＯ１００Ｋ＝１、１００というループ内にあったのに、移動（複写）先のＤＯループの変数がＫではなくＮだったためです。


      DO 100 K = 1,100
        A(K)   = A(K) + 1.0DO
  100 CONTINUE

を、


      DO 200 N = 1,100
        A(K)   = A(K) + 1.0D0
  200 CONTINUE

としてしまいました。これではちゃんと動きません。ここでの問題は、同じ繰り返し数で扱うＤＯループの制御変数（ＫやＮ）が統一されていないことです。行き当たりばったり的に作ったプログラムは、特にこのような単純な間違いで足元をすくわれてしまいます。
（自戒の意味を込めて）プログラムを作る時は、同じことをするＤＯループは同じ制御変数を使いましょう。

（１）本日のバグ（レポート１／２３、１９９７） [目次]

新しい計算ルーチン作成中に、今日の作業はここまでと、敢えて注釈を付けずに------------------------------を入れておいたのを完全に忘れ、翌日そのままコーディング作業を続けてコンパイル、実行してみるとバグの嵐だった。このプログラムはSUNのＯＥＭマシンだったが、コンパイルエラーメッセージだけでは原因が良く分からなかった。ところがＤＥＣのマシンでコンパイルすると同じ箇所でコンパイルエラーとなるが、エラーメッセージからすぐに、注釈を付けずに、----------------------------をコード内に挿入していたことがわかった。

この問題の重要な点は、筆者がＤＥＣが主に使用するマシンで、SUNのマシンはごくたまにしか使わないため、あまりSUN上でのデバッグの経験が乏しいことにあります（その前に、ちゃんとマニュアルのコンパイルエラーメッセージの説明を読んでいないのが最大の問題かもしれません^^;）。

既に昨日のことになってしまいましたが、ここで指摘したいのは、コンパイラーのエラーメッセージの意味を良く理解しておくこと（エラーメッセージの出し方や、メッセージの内容の貧弱さ、難解さにも問題がありますが）、わからないメッセージに出合ったら、ちゃんとマニュアルを読むことです。偶然、異なるシステム同士の結果の比較からエラーの原因を突き止めましたが、これはエラーメッセージの意味をちゃんと把握しておけば、そのようなことをする必要のないバグでした。

[先頭][総目次 ][最初に戻る][Top][失敗知識データベース(JST)]

とても最近に、やらかした計算上のバグ、失敗 [Top][失敗知識データベース （畑村創造工学研究所←JSTから移動）]