たまに書きます。

気になって調べたことを書いていきます。

大昔のgccのプリプロセッサに関する、極めてどうでもよい話

先日、酒を飲みながらgccだかなんだかのコンパイラの話をしていたときに、ふと、「そういえば昔のgccプリプロセッサ#pragmaディレクティブを見つけると、コンパイルを中止してゲームを起動するとかいうのがあったよな」と思い出した。 この元ネタは、もう絶版だけど、「エキスパートCプログラミング」という本に書いてある。

エキスパートCプログラミング―知られざるCの深層 (Ascii books)

エキスパートCプログラミング―知られざるCの深層 (Ascii books)

C/C++において、#pragmaは、ヘッダの二重インクルードを防止するために、伝統的な

    #ifndef GUARD
    #define GUARD
    ...
    #endif

を使用したインクルードガードに変わって#pragma once のように用いられることが多い。しかし、pragmaに続くキーワードやその動作は下記のページにもあるように、処理系ごとに定義が委ねられている。

C言語/前処理指令 - Wikibooks

当時のgccはどうも#pragmaをいけ好かないものとしていたようで、「処理系ごとに決められた動作」として、ゲームを起動するというようなことをやっていたらしい。

せっかくなので、gccの古いソースコードを掘り返してみる。 アーカイブへの直リンクは

http://ftp.tsukuba.wide.ad.jp/software/gcc/old-releases/gcc-1/

さて、まずは一番古そうなのをということで、verion 0.9 (gcc-0.9.tar.bz2, 1987/03/22リリース)を覗いてみる。 このなかの、cccp.cというファイルを見てみると、2401行目付近に、

/*
 * the behavior of the #pragma directive is implementation defined.
 * this implementation defines it as follows.
 */
do_pragma ()
{
  close (0);
  if (open ("/dev/tty", O_RDONLY) != 0)
    goto nope;
  close (1);
  if (open ("/dev/tty", O_WRONLY) != 1)
    goto nope;
  execl ("/usr/games/hack", "#pragma", 0);
  execl ("/usr/games/rogue", "#pragma", 0);
  execl ("/usr/new/emacs", "-f", "hanoi", "9", "-kill", 0);
  execl ("/usr/local/emacs", "-f", "hanoi", "9", "-kill", 0);
nope:
  fatal ("You are in a maze of twisty compiler features, all different");
}

というように、このころのコードには確かにゲームを起動するコードが仕組まれている。ちなみに、私はhackもrogueも名前だけは知っていたが、起動して遊んだことはないので何するゲームなのかよく知らない*1。 なるほどexeclはそのプロセスをそのまま呼び出すコマンドで置き換えるから、より新しいものから起動を試してみて、見つからなかったらrogueというように順々に探していくわけだ。 ちなみに、コンパイルができなかったために、本当にpragmaを見つけるとゲームが起動するのか、手元で確かめてはいない。というのも、do_pragma()という関数自体が、関数ポインタを介して呼ばれているようで、必ず呼ばれるものか、いまいち確信は持てていないためだ。

つぎに、gcc-1.21(gcc-1.21.tar.bz2, 1988/05/01リリース)を見てみる。 すると、同じくcccp.cの2993行目付近

#if 0
/* This was a fun hack, but #pragma seems to start to be useful.
   By failing to recognize it, we pass it through unchanged to cc1.  */

/*
 * the behavior of the #pragma directive is implementation defined.
 * this implementation defines it as follows.
 */
do_pragma ()
{
  close (0);
  if (open ("/dev/tty", O_RDONLY) != 0)
    goto nope;
  close (1);
  if (open ("/dev/tty", O_WRONLY) != 1)
    goto nope;
  execl ("/usr/games/hack", "#pragma", 0);
  execl ("/usr/games/rogue", "#pragma", 0);
  execl ("/usr/new/emacs", "-f", "hanoi", "9", "-kill", 0);
  execl ("/usr/local/emacs", "-f", "hanoi", "9", "-kill", 0);
nope:
  fatal ("You are in a maze of twisty compiler features, all different");
}
#endif

というように、コメントアウトされてしまっている。しかも「なんだか使い道が見えてきた(意訳)」という言い訳つきである。上記のwikibooksによると、「処理系が認識できないプラグマは無視する」と書いてあるが、この段階ではコメントにもあるように、プリプロセッサでは無視しているようだ。

エキスパートCプログラミングという本自体、1996年に出版されていることを考えると、この本が「昔」というのだから「昔のさらに昔」ということで本当に初期の一瞬だけの話だけだったようだ。ちなみに、私は1990年生まれなので、生まれる前にはもうpragmaは「ちゃんと使われていた」ことを今更知った次第である。

なお、gcc-0.9のビルド周辺に関しては、下記のページでビルドの試行錯誤などがされている。

virtuallyfun.superglobalmegacorp.com

終わり

DJに入門した話

実は3年くらい前から、DJに憧れていた。 きっかけは、知り合いのエンジニアの方が、いつぞやのプログラミング言語のカンファレンスの懇親会でDJをやっているのを見たから。その後、今の研究室に移ってからはなんとプロのDJがいて、ちょっと教えてもらうにも非常に良い環境でもあるのだが、 なかなかお金もかかるもので結局やらないままだった。(ずっと「いつか始めたいんですけどね。。。」と言うだけだった)

そうこうすること3年弱、この前我が研究室のDJの方が、こいつ↓を持って来て、プレイするところを見せてくれた。

jp.monstergodj.com

そもそも自分がなかなか始められなかったのは、DJを始める上での初期投資となる

  • タンテ×2(もしくはCD×2)

  • ミキサー

  • ヘッドホン

のコストが自分の経済力ではやれなかったからなのだけど、今の時代(ここ数年?)こんな小さなマシンであんなカッコいいことが出来るとかマジでヤバい。鞄からさっと取り出していきなりバンバン音を操りまくって。。。ていうのがカッチョよすぎる。
就職活動でなんとなく落ち着かない生活も終わったということで、迷わず再入荷のタイミングで買いました。

こんな感じ。やばい。かっこいい。 f:id:salondunord:20170711002759j:plain

今はまだ操作の仕方も全然感覚的でなくて(まだどのボタンを押したらどうなるのかもよくわかってない)、曲を繋ぐとかもってのほかなのだけど、それでも自分がPCに入れていた曲がちょっとボタンを押したりつまみをひねるだけで音が変わったりするのは本当に面白い。 それ以来、毎晩帰宅してからちょこちょこ弄って遊んでるんだけど、しっかり練習していつか僕も何かの懇親会とか結婚式の二次会とかそういうところで回させてもらえるようになりたいと思う。

以上、やめない為の決意表明でした。

MONSTER GODJ Portable, Stand-Alone DJ System and Production Studio【正規代理店品】

MONSTER GODJ Portable, Stand-Alone DJ System and Production Studio【正規代理店品】

量子化学計算のプログラムを自作してみた

普段、量子化学の計算を行う際は、GaussianやGAMESS他のソフトが用いられることが殆どで、自分でコードを書くのは計算手法の研究をされている一部の方に限られるのではないかと思います。

近年では、B3LYP汎関数の使用を始めとする、密度汎関数理論(Density Functional Theory; DFT)が用いられるケースが殆どですが、それに先立って確立された手法として、Hartree Fock法があります。 Hartree Fock (HF)法は、計算結果のエネルギーの見積もりの不正確さ*1から、今ではこの手法が単体で使用される場面は多くはありません。 しかし、近似の仕方が明確なため*2量子化学の基礎として多くの方が一度は勉強される理論かと思います。

私も、数年前に教科書でHFの理論を勉強してから、一度は自分でコードを書いておきたいという気持ちがずっとあったのですが、毎回途中やめになっており、このたび漸くそのminimalな実装をC++で書ききりました。実際に書いてみると、式を読んで分かった気になっていたものの、実際の理解は曖昧だったりといったことがよくわかりました。

コードは、下記のリポジトリにおいていますので、気になる方、これから同じことを勉強しようと言う方はどうぞ。

github.com

感想としては、結構大変でした。Gaussianなどの開発チームだったり、もしくは日々新たな計算手法を開発されている方たちが神に見えてきます。 以下に気づき、感想、苦労した点を列挙します。

1. Gauss型の短縮基底関数の正規化

これは、自分で書いてみるまで気がつかなかったことです。 (といっても書き始めたのが2年くらいは前なので、その頃に気がついたのですが。。。) 例えば、STO-3Gという最も基本的な基底関数系で、水素の1S軌道は、軌道指数(Orbital exponent, よくαと書かれます)と計数を強調して、

 H_{1s} (r) = 0.44635 g(-0.168856) + 0.535328 g(-0.623913) + 0.154329 g(-3.42525)

という感じでよく表されます。ただ、これを  g(\alpha) = \exp(-\alpha r^{2}) と考えて素直に計算すると、一番簡単なはずの重なり積分から数字が合いません。

実際には、(後で気づけば当たり前のことではあるのですが...)それぞれのプリミティブガウス関数に、それぞれを全空間で積分した際の値が1になるように、規格化因子がかけられています。また、さらに、それらを重ね合わせた関数(短縮ガウス関数)にも、規格化因子が掛けられます。したがって、より正しい表記は、

 H_{1s} (r) = \frac{1}{0.999999 } (0.44635 \frac{\exp(-0.168856 r^{2})}{0.187736} + 0.535328 \frac{\exp(-0.623913 r^{2})}{0.500326} + 0.154329\frac{\exp(-3.42525  r^{2})}{1.79444} )

となります。これらの規格化因子は、プリミティブガウス関数のそれぞれについてまず求め、その後、短縮ガウス関数について求める必要があります。なお、p軌道以上の角運動量を持つ軌道には、 x^{l}y^{m}z^{n}という因子がかかるので、これも考慮に入れる必要があります。

2.積分の実装

積分は、1966年に竹田、藤永、大旗の3名が書いた論文を参考にしました。これの実装がとにかく難しかったです。慣れれば、重なり積分と、運動エネルギー積分はまだなんとかなります。難しいのは、核引力積分(Nuclear Attraction Integral)と電子間反発(4中心)積分の部分です。

これらは、とにかく展開方法もなかなか理解できませんし、非常に入れ子の深いループになってしまいます。

現状は上記の論文をそのまま参考にしており、STO-3G基底での水素分子やHeH+イオン程度は一瞬で計算は終わりますが、試しにp軌道を分極関数として追加してみると、一気に遅くなりました。

近年のGaussianを始めとする量子化学計算プログラムは、PRISMアルゴリズムなど、積分ごとに20種類くらいのパスを使い分けるような非常に複雑な実装*3がなされているそうですが、その重要性がよーく分かりました。 今後は、できればDKRの積分方法や、Obara-Saikaの展開公式なども実装してみたいと思います(できれば、です。2電子積分の展開など、まだ全然理解できていません。)

3. 何を使って行列演算を行うか

王道はBLASなどのライブラリを使うことなのだろうと思いますが、恥ずかしながら私は他のプログラムからリンクさせる為にコンパイルした程度の経験しかありませんでした。

以前にnumpyを使って今回と同じようなものを書こうとしたことがありました*4ので、今回は比較的それと似たような感覚で使えそうなEigenというテンプレートライブラリを使用しました。

おそらくBLASなどのような高度なチューニングは難しいでしょうが、今回のような目的には非常に使い勝手が良かったと思います。

今後

現状はHartree FockのRHF計算のみ実装しましたので、今後はUHF計算をまずは実装したいと思います。 あとは、スパコンが使える今のうちに、少しくらいMPIチューニングなどで遊んでみるのも良いかもしれません。 その次は、DFTか微分計算ではありますが、おそらくUHFを書いてみるころには学生生活が終了している気がします...苦笑

2017.07.15 追記

UHFの計算はSzaboの本を見ながらすぐにかけてしまいました。 今はDFTを実装するべくいろいろ調べていますが、交換・相関ポテンシャルの空間積分のためのグリッド生成がなかなか難しいです。またここにメモをする予定です。

追記その2 Qiitaにメモしました。まだ積分は書いていなくて、ハードコーディングしてます。 qiita.com

Modern Quantum Chemistry: Introduction to Advanced Electronic Structure Theory (Dover Books on Chemistry)

Modern Quantum Chemistry: Introduction to Advanced Electronic Structure Theory (Dover Books on Chemistry)

新しい量子化学―電子構造の理論入門〈上〉

新しい量子化学―電子構造の理論入門〈上〉

密度汎関数法の基礎 (KS物理専門書)

密度汎関数法の基礎 (KS物理専門書)

*1:電子相関(≒電子同士の避け合いの効果)を考慮しない

*2:基底関数以外に経験的なパラメータを含めていない

*3:すいません、嘘言っているかも。全然この辺りの詳細は知りません

*4:どうしても2電子積分の値がズレてしまい放置...

VASPで計算したバンド構造の可視化のためのスクリプト

固体のバンド構造計算をこれまでは主にQuantum Espressoというフリーのコードを使って行っていたのですが、最近、少々vaspを利用した計算を行っていました。

vaspでバンド構造を計算すると、EIGENVALというファイルが作成されて、そこに各k点での固有値が記録されます。 このファイルを何らかのツールを使って可視化する必要があるわけですが、検索してみると主に出てくるのはp4vaspというプログラムを利用する、というもののようです。

しかし、p4vaspはなんだかバンド図の見た目がイケてない上に、磁性を考慮した場合にどうもupとdownのスピンのバンドを別々にカラーリングする方法がわからず、バンド図の可視化には使いづらいように思えました。 そこで、gnuplotによる可視化向けにEIGENVALSファイルを読み込んで整形するスクリプトRubyで作成しました。

github.com

本当はここに結果を載せたいのですが、一応vaspはライセンスの縛りがあるので、ここに結果のグラフを載せるのは避けておきます。

使い方は、EIGENVALSのあるディレクトリで、

ruby read_eigenvals.rb -n (各kpathの点の数)

とするだけです。オプションを適当に渡せば、スピンを考慮した計算結果のファイルや、gnuplot用のスクリプトも一緒に作成します。

理論化学の本を紹介していく(電子軌道の概念)

量子化学の一番初めは、とりあえず波動関数とは何ぞや?を抑えることだとすると、その後で来るのは、今度は軌道の概念を知ることだと思っています。 私は軌道のことは、友田修司先生の本で勉強をしました。

基礎量子化学―軌道概念で化学を考える

基礎量子化学―軌道概念で化学を考える

この本で使っている計算法は拡張ヒュッケル法などの現代の量子化学コミュニティではあまり使われなくなってしまった古い計算方法です。

しかし、軌道の概念をきちんと理解するためには、むしろこのような手法の方が変に原子軌道の混成を考えなくてよくなるため、定性的な理解には良いと思います。

この本はそれほど大きくない分子を用いて、

  • フロンティア軌道

  • 配位子場分裂・結晶場分裂(金属に配位したイオンなどにおいてみられる3d軌道の分裂)

などを議論しており、これらの概念を掴むには良いと思っています(ただし、すごく深いわけでもないと思います)。

また、溶解エネルギーなどの熱力学的な量に関する解説もある程度されており、数表としても私は使っています(ただし、そんなに量は無いので、本当に必要な時は化学便覧を見るのが良いです)。

私は指導教員によく言われるのは、「計算機で計算する数値データは、コンピュータとか計算手法が発展するので10年たったら役にたたないが、軌道を見た時の定性的な傾向(たとえば、フロンティア軌道は分子のどのあたりに局在しているか、など)は別の量子力学的効果を入れたりしない限り変わることは無いので、それをきちんと把握すること」ということです。

この本はその感覚をつかむために最適なのではないかと思っています。

ちなみに、軌道の概念をさらによく知るためには、洋書ですが、

Orbital Interactions in Chemistry

Orbital Interactions in Chemistry

かと思います。

化学屋から見たバンド理論 その2

今の時代、バンド計算も多くの場合、VASPやQuantum Espressoといった密度汎関数法による計算をすることになるかと思います。 バンド計算においては、必要なパラメータとして、

  • 汎関数(大体はPBE汎関数が多い模様。必要に応じて、PBE0やHSE06のようなハイブリッド汎関数を使うが計算が重くなる)

  • 擬ポテンシャル(PAWかUltrasoftを使えば大抵の場合いける、他にはノルム保存型などの選択肢がある)

  • カットオフエネルギー

  • k点メッシュ

  • スメアリング(ブロードニング)

を必ず指定する必要があります。

とりあえずこれらの意味が分からない時は、以下の本が良いと思います。

密度汎関数理論入門: 理論とその応用

密度汎関数理論入門: 理論とその応用

この本は、DFTの理論的な詳細には立ち入らず、あくまで使う側に必要な情報の提供に絞られています。 また、VASPやQuantum ESPRESSO, Ab Initなどのサンプルinputファイルが書かれているわけでもないので、 まずはネット上のチュートリアルなどからinputを探してきて、それの意味を一つ一つ理解していくのに使うべきと思います。

ちなみに、この本には、バンド構造(E-k分散図)は触れられていません。バンド構造の理解に関しては、前の記事で紹介した、ホフマンやCoxの本で理解のコツを覚えていくしかないと思っています。 また、バンド構造を計算するには、まずはブリルアンゾーンを指定する必要がありますが、私は

AFLOW.ORG: a distributed materials properties repository from high-throughput ab initio calculation

を使ってそれを指定しています(Kpath in the reciprocal spaceというのを指定して、出てきた結果をvaspの場合はKPOINTSファイルに指定)。

化学屋から見たバンド理論

化学系の学部・学科で理論計算を行う人は、gaussianやGAMESSといった量子化学計算プログラムを使って解析を行うケースが殆どだと思います。

しかし、シュレーディンガー方程式の解釈の仕方はそれだけではなくて、バンド理論というのも一つの見方です。

vaspやQuantum Espressoといったプログラムを用いて計算したバンド構造図は、初めての人には非常にわかりづらく、また、何度か見せられた場合でもよく言われる特徴(バンドギャップの有無、それが直接遷移か間接遷移か)くらいしかわからない人もいるかと思います。

化学分野の人がこれを理解しなければならない、というときは、以下の本が良かったので、ここに載せておきます。

固体と表面の理論化学

固体と表面の理論化学

固体の電子構造と化学

固体の電子構造と化学

一つ目は、福井謙一先生と一緒にノーベル賞を受賞されたホフマンの本で、バンドを直感的に理解するには非常にわかりやすい本であると言えます。バンドは実際には一つ一つが分子軌道と同じように考えることができます。

ただ、この本は絶版なので、図書館で探すしかないと思います。

二冊目は一冊目の本をもう少し現代的に書いた本で、4章でバンド理論のホフマンの本の要点は一通りまとめられてはいます(が、ホフマンの方が個人的にはわかりやすかった)。 実際にバンド理論を考えるには、3d軌道のように、電子が局在する最外殻軌道を考える必要があります。その時に出てくるのがHubbard-Uパラメータですが、このあたりについても説明がなされているので、必要に応じて読むのが良いと思います。