kivantium活動日記

プログラムを使っていろいろやります

シリコンバレー観光記

留学に行くついでにシリコンバレーの観光をしたのでその記録を残しておきます。

青が1日目、緑が2日目、黄色が3日目の行き先です。

1日目(サンフランシスコ)

シリコンバレー観光といいながら1日目はサンフランシスコに行きました。この記事によると近年はシリコンバレーを離れてサンフランシスコに拠点を置くスタートアップが増えているそうです。

f:id:kivantium:20190806173857j:plain:w400
Twitter本社
f:id:kivantium:20190807111924j:plain:w400
Googleサンフランシスコ

BARTの券売機

空港からサンフランシスコへの移動はBARTという鉄道を利用したのですが、この券売機の仕様が分かりにくかったので記録しておきます。

f:id:kivantium:20190807132018j:plain:w600
インターフェース

切符を買うためにはまず券売機に貼られているシールを見て目的地までの切符料金を確認します。
次にその金額の切符を買うのですが、金額指定の方法が分かりにくかったです。例えば10ドル紙幣を持っているときに5.25ドルの切符を買いたいとすると、10ドル入れた後にBボタン (Subtract $1) を5回押したあと、Cボタン (Add 5¢) を5回押してCurrent valueを$5.25にしてからEボタン (Print $x.xx Ticket) を押して発券するという複雑な手順を要求されます。
クレジットカードを使う場合は手前のカードスロットでスキャンすると20ドル投入した扱いになるので同じ要領で値段を調整して購入します。
空港で最初にこのインターフェースを見たときはAdd $1が何を意味するのか分からなくて切符を買うことができませんでした。
空港のWi-Fiが通じるエリアまで戻ったあと「BART チケット 買い方」で検索してたどり着いたサンフランシスコのBARTチケットの購入方法|たびめもを見てようやく理解しました。
分かってしまえば簡単ですが、初見では分からないものです……

2日目(コンピュータ歴史博物館周辺)

シリコンバレーにはCaltrainを使って移動し、Airbnbで取った宿に宿泊しました。シリコンバレーの地価は異常に高いので写真のような4人部屋で1泊およそ5000円でした。

f:id:kivantium:20190809190840j:plain:w400
シリコンバレーの宿

宿から最初の目的地のコンピュータ歴史博物館にはMountain View Community Shuttleを利用して行きました。
地図を見た感じ徒歩で行けるような気がしてしまうのですが、アメリカはとにかく大きいので日本の感覚で地図を読むと読み誤ります。地図上の感覚に頼らず、徒歩で行くと何分かかるのか経路探索を使ってきちんと調べることが大事です。
Googleマップの導きに従ってバスに乗ったものの、降り方が分からず目的の停留所で降りることができませんでした。幸い次の停留所で降りる人がいたのでそこで一緒に降りて停留所一つぶん歩いて戻りました。アメリカのバスでは降りるバス停を紐を引っ張ることで伝えるのですが、そんなことは知らなかったのでずっとボタンを探していました。

f:id:kivantium:20190809135810j:plain:w400
バスの車内。分かりにくいが窓枠に沿って紐があり、この紐を引っ張ることで次の停留所で降りたい意思を運転手に伝えることができる。

コンピュータ歴史博物館は月曜日(と時期によっては火曜日)が休館日なので注意してください。(Hours & Admission)

f:id:kivantium:20190808103828j:plain:w400
コンピュータ歴史博物館の看板
f:id:kivantium:20190808105456j:plain:w400
最初の展示はそろばん
f:id:kivantium:20190808121511j:plain:w400
IBMの元となった会社のパンチカード処理機
f:id:kivantium:20190808125813j:plain:w400
ムーアの法則のグラフ
f:id:kivantium:20190808113352j:plain:w400
有名なティーカップ
f:id:kivantium:20190808133003j:plain:w400
座れるCLAY-1

毎日12時と14時から行われるガイド付きツアーに参加するのもおすすめです。展示に書いてあることよりも詳しく説明してくれます。

f:id:kivantium:20190808121715j:plain:w400
ガイドのおじさん。元CLAY従業員とのこと。

コンピュータ歴史博物館の後は徒歩で移動してMicrosoftGoogleを見学した後、Uberで宿に戻りました。

f:id:kivantium:20190808143552j:plain:w400
チューリング賞受賞者専用駐車スペースがあったことで有名なMicrosoft。その表示は現在存在しないようだった。

f:id:kivantium:20190808145722j:plain:w400
GoogleplexことGoogle本社。建物の中には入れないが全体が公園のようになっていて、観光客がたくさん来ていた。

3日目(スタンフォード大学

3日目はスタンフォード大学を見学しました。

f:id:kivantium:20190809095342j:plain:w400
観光案内所のお姉さんが教えてくれた一番有名な建物Memorial Church。
f:id:kivantium:20190809103301j:plain:w400
Gates Computer Science。名前的にコンピュータサイエンス学科の建物だと思われる。館内にはコンピュータ歴史博物館と連携したComputer History Exhibitsという展示がある。
f:id:kivantium:20190809101405j:plain:w400
コンピュータサイエンス学科の建物の掲示板にはCofounderを求める張り紙がたくさん
f:id:kivantium:20190809120011j:plain:w400
キャンパスの雰囲気。奥に見えるのがHoover Tower。

f:id:kivantium:20190809112814j:plain:w400
1日2回キャンパスツアーが行われており、学生の話を聞きながらキャンパスを回ることができる。集合時間にVisitor Centerに集まるだけで申し込み不要・無料で参加できる。

その後バスで移動してFacebookの有名な看板を見てきました。

f:id:kivantium:20190809145518j:plain:w400
表側は観光客がひっきりなしに記念撮影を行っている
f:id:kivantium:20190809145629j:plain:w400
裏側はサンマイクロシステムズのロゴが書かれている。Facebookの現本社は以前サンマイクロシステムズの本社であり、没落した会社がどうなるかの戒めとしてわざとサンマイクロシステムズの看板を裏返して使っているとのこと。

その後はサンフランシスコ国際空港で飛行機に乗って留学先に向かいました。

その他

シリコンバレーは世界のテック企業が集まる技術の中心なので高層ビルが立ち並ぶ大都会を想像していましたが、実際に行ってみるとかなりの田舎であるという感想を持ちました。最寄りのコンビニに行くにも20分歩く必要があるなど、車社会であることを感じさせられました。日本と違って土地が広大だから建物を集積させる必要がないなどの事情があるのでしょうが、東京などに慣れてしまった日本人としてはかなり住みにくいと思いました。あと食べ物がおいしくないのはとてもつらいです。

f:id:kivantium:20190807162241j:plain:w400
f:id:kivantium:20190807175258j:plain:w400
このようなのどかな風景がずっと続いている

留学の資金源について

最初に述べた通り、この観光は留学のついでに行ったものなので留学の資金源についても述べておきます。
今回の留学は官民協働留学支援制度を主な資金源にする予定です。この奨学金は語学留学や単位取得ではなくインターンなどの実践活動を主な目的とする留学を対象としています。民間企業から集めた寄付を文部科学省が取りまとめて日本学生支援機構を経由して学生に給付するという形を取っているようです。留学計画に対して審査があり、僕が参加している11期では1939人中544人が採用されました現在12期の募集を行っているので、留学に興味がある人は応募してもいいかもしれません。

採用者から見たこの奨学金の良いところと悪いところについて書いておくので参考にしてください。

良いところ

給付型奨学金なので返済の必要がありません。

悪いところ

金額が少ない

留学先の地域によってもらえる奨学金の額が異なりますが、僕が行く北米では月額16万円が支給されます。この額はアメリカで生活する上では不十分です。今月は家賃として16万2000円支払ったのでそれだけで奨学金が全額溶けました。また、アメリカに長期間滞在する場合はビザが必要となりますが、労働に従事しない研究者用のビザを発行するためにはアメリカで労働する必要がないことを証明するために母国で十分な資金を受け取っていることを示す必要があります。この「十分な資金」のラインが年間2万5000ドルでした(これは学校によるらしいので僕の場合)。月16万円の奨学金を受け取っても年間1万9200ドル(1ドル100円の場合)にしかなりません。つまり、アメリカ合衆国からするとこの奨学金アメリカでの生活に十分な資金源だと認めることができないということです。
別に奨学金だけで留学資金をまかなう必要はありませんし、一人あたりの金額を増やすと採用人数が減ってしまうことを考えると一概に金額を増やせとは言えませんが、給付額が十分でないという点は事実です。

支給後に返還を求められる可能性がある

留学の途中で何らかの事情があって留学計画を変更することがありますが、計画変更が文科省によって承諾されなかった場合は既に支給された奨学金であっても返還を求められる可能性があります。留学計画の変更の審査には2ヶ月かかるらしく、2ヶ月以上前に変更申請するよう求められていますが、これは非常に難しいです。

  • 8/10以降に開始される留学を対象とした奨学金であるにも関わらず、変更申請を出せるのは7/19以降なので8/10の予定変更を申請することは不可能
  • 留学先の国からビザが許可されるかどうかによって留学計画は大きく変わるが、2ヶ月前にビザの可否が分かることは通常ない
  • 留学先機関が急に留学生を受け入れなくなって留学を中止することになるかもしれないが、2ヶ月前に分かることは通常ない

したがって、たとえば2つの留学先に行く留学計画を提出していたときに、1つめの留学先での滞在を終えたタイミングで2つめの留学先に行けなくなることが判明したが、申請が遅かったため計画変更が承認されず既に支給された奨学金を全額返還させられるというケースが起こりえます。実際自分も出発直前になってビザが間に合わないことが判明したので計画変更を行う必要があり、この留学の奨学金が支払われるのかまだ確定しない状態です。国家プロジェクトなので不正受給を防ぐ必要があるのはわかりますが、学生にとって非常に使いにくい奨学金になっていると思います。

結果通知が遅い

この奨学金は2019年8月10日以降に開始される留学を対象にしていますが、採用結果が通知されたのは6月18日でした。奨学金に採用されたら留学を行うことにしていたので結果が通知された日からビザなどの手続きを行ったのですが、受け入れ先の休暇などが重なり渡航開始には間に合いませんでした。そもそもアメリカ大使館によればビザの申請は渡航の3ヶ月前に行うべきだとのことで、結果通知が届いた時点で遅かったということです。奨学金が通るかどうかに関わらずビザの手続きを進めておくべきだったのでしょうが、留学が可能かどうか分からない時点で手続きを進めるのは現実的には難しいです。また、もし仮にこの奨学金が十分な額を支払っているなら本来その奨学金の額を大使館に提出して収入証明にするのが道理なので手続きとしても奨学金の結果通知を待つのが正当なはずですが、現状ではそれが不可能になっています。

研修がつらい

奨学金に採用された学生は文科省主催の2日間の研修を受ける必要があるのですが、非常につらい内容でした

  • 「自分の軸を見つめ直す」と称して、自分の人生を小学1年生から振り返るシートを作りグループディスカッションをさせられた
  • 留学先で日本文化を発信する方法を考えるグループディスカッションをさせられた(配布された模範的な発信方法の例: 外国人にけん玉を教える、折り紙を教える、流しそうめんパーティーを開く)
  • グローバルリーダーになるために必要なことについてグループディスカッションをさせられた
  • 奨学金の原資となる寄付金を集めるのかがいかに難しかったかなどの運営苦労話を聞かされた

などなど。少なくとも、海外の機関で研究を行うことを目的とする人に役立つ内容は皆無でした。(研究室でけん玉をする必要があるでしょうか?)全国の学生を交通費も支払わず東京に呼び集めて2日間拘束して行う教育の内容がこれでいいのか非常に疑問です。

事務手続きが煩雑

この奨学金を受け取るためには採用後も大量の書類を提出する必要があります。受け入れ証明書・受け入れ担当者の署名が入った在籍証明書(締め切りは毎月5日まで)・飛行機の半券・パスポートの入出国履歴のコピー・留学状況報告書……。そしてその書類は大学の事務を経由して文科省に送られるので、大学の事務担当者は学生と文科省の間の調整役としてさらに多くの事務作業を行っているはずで、事務手続きを行うたびに申し訳なくなっています。

この奨学金のほかに民間企業からの奨学金を受け取っているのですが、そちらは採用された後に銀行口座の情報を送るだけで振り込んでくれたので非常に使い勝手がよかったです。留学用奨学金は今のところ1円も振り込まれていないので飛行機代などは全てそこから支払いました。支援していただいたトヨタドワンゴには感謝しています。

本留学は官民協働留学支援制度の支援を受けている。この奨学金を受け取った学生は留学エヴァンジェリストとして日本の学生に留学の魅力を伝えることが求められており、この記事は筆者のエヴァンジェリスト活動の一環として書かれた。この記事が留学機運醸成の一助となることを願っている。

NWChemをコンパイルして動かす

なにも分からない。

NWChemは計算化学ソフトで、GSoCのOpenChemistryに参加していたので動かしてみたかった。

とりあえず動いた方法

環境はUbuntu 18.04。

ドキュメントのCompiling NWChemにあるNWChem 6.6 on Ubuntu 14.04 (Trusty Tahr)の項目を見ればいいと思っていたが微妙に動かなかった。そもそも現在の最新リリースバージョンは6.8.1なのでドキュメントが更新されてなさそう。以下のように改変したらとりあえず動くようになった。

必要なパッケージのインストール

python-dev gfortran libopenblas-dev libopenmpi-dev openmpi-bin tcsh make 

最新リリースをクローンする

git clone -b hotfix/release-6-8 https://github.com/nwchemgit/nwchem.git nwchem-6.8.1
cd nwchem-6.8.1/src

環境変数の設定

ドキュメントでは

export USE_MPI=y  
export NWCHEM_TARGET=LINUX64  
export USE_PYTHONCONFIG=y  
export PYTHONVERSION=2.7  
export PYTHONHOME=/usr 
export BLASOPT="-lopenblas -lpthread -lrt"  
export BLAS_SIZE=4  
export USE_64TO32=y

とするように指示されているが、BLASOPTを指定するとmake時にLAPACK_LIBを指定しろというエラーが発生する。しかし、その正しい値についてはドキュメントに言及がなく分からなかった。そこで

unset BLASOPT

を実行したところとりあえずエラーが消えたので良しとした。(絶対良くないんだよなぁ)

コンパイル

make nwchem_config NWCHEM_MODULES="all python"
make 64_to_32
make

定石どおり make install しようとしたら動かなかったのでコンパイル後にできたbin/以下をPATHに加えてインストール終了とした。

動作チェック

Getting Started にあるコードが動くことを確認した。(aptでインストールしたバージョンではこれが動かなかったのがそもそもコンパイルしないといけなくなった原因なのでNWChemのメンテナンスが全般的に間に合ってないっぽい)

ついでに閉殻Hartree-Fock法によるHeH+のエネルギー計算 - kivantium活動日記と同じ設定で走らせて結果がどうなるか見てみた。

title "HeH+ STO-3G SCF energy calculation"
charge 1
geometry  
  H 0 0 0
  He 0 0 1.4632
end
basis
  H  library sto-3g
  He library sto-3g
end
task scf energy

結果は

       Final RHF  results 
       ------------------ 

         Total SCF energy =     -2.825194209471
      One-electron energy =     -4.549711559836
      Two-electron energy =      1.001202357200
 Nuclear repulsion energy =      0.723314993165

        Time for solution =      0.0s

となり、前回の結果のTOTAL ENERGY = -2.860662163500と近い値が出ることが確認できた。

閉殻Hartree-Fock法によるHeH+のエネルギー計算

量子化学計算を勉強するために新しい量子化学―電子構造の理論入門〈上〉を読んでいます。
サンプルコードがFortranだったので勉強がてらC言語に移植しました。

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

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

理論

本で200ページくらいかけて説明されている内容をブログに書くのは大変なので省略します。

ちなみに、以下のツイートの「100ページかけて説明してる別の本」はこの本を指します。
難しい内容を易しく簡潔に説明しようとするのは大変で、難しい内容は難しいまま長い時間をかけて説明する必要があるのだなぁという学びがありました。


実装

本のFortran実装を素直にC言語で書き直しました。
本のコードを打ち込むのが面倒だったので、検索して見つけたページにあったコードと見比べながら書きました。

最終的に得られたエネルギーはFortran実装の結果と小数第5位まで一致しました。
この差がどこで生じているのか謎ですが、実装は正しいんじゃないでしょうか……

#include <math.h>
#include <stdio.h>

// プロトタイプ宣言
void HFCALC(int, int, double, double, double, double, double);
void INTGRL(int, int, double, double, double, double, double);
void COLECT(int, int, double, double, double, double, double);
void SCF(int, int, double, double, double, double, double);
double calcS(double, double, double);
double calcT(double, double, double);
double F0(double);
double calcV(double, double, double, double, double);
double TWOE(double, double, double, double, double, double, double);
void FORMG(void);
void DIAG(double [][2], double [][2], double [][2]);
void MULT(double [][2], double [][2], double [][2], int, int);
void MATOUT(double [][2], int, int, int, int, char*);

double S12, T11, T12, T22, V11A, V12A, V22A, V11B, V12B, V22B,
       V1111, V2111, V2121, V2211, V2221, V2222;

// Slater型関数に対する最小2乗近似の短縮係数と軌道指数
double COEF[3][3] = {{1.0, 0.678914, 0.444635}, 
                     {0.0, 0.430129, 0.535328},
                     {0.0, 0.0, 0.154329}};
double EXPON[3][3] = {{0.27095, 0.151623, 0.109818},
                      {0.0, 0.851819, 0.405771},
                      {0.0, 0.0, 2.22766}};
double D1[3], A1[3], D2[3], A2[3];
double PI = 3.1415926535898;

double S[2][2], X[2][2], XT[2][2], H[2][2], F[2][2], G[2][2], C[2][2],
       FPRIME[2][2], CPRIME[2][2], P[2][2], OLDP[2][2], TT[2][2][2][2], E[2][2];

int main(void) {
  int IOP = 2;
  // 基底関数STO-NG
  int N = 3;
  // 基底状態の平衡結合長(a.u.)
  double R = 1.4632;
  // Heに対するSTO-3G軌道指数値
  // He原子に対する最良の軌道指数は1.6875
  // HeH+では電子雲が縮んでいることを反映して1.24倍
  double ZETA1 = 2.0925;
  // Hに対するSTO-3G軌道指数値
  // H原子に対する最良の軌道指数は1.0
  // HeH+では電子雲が縮んでいることを反映して1.24倍
  double ZETA2 = 1.24;
  double ZA = 2.0;
  double ZB = 1.0;
  HFCALC(IOP, N, R, ZETA1, ZETA2, ZA, ZB);
}

/*
 * 1s基底関数に対応するSTO-NG最小基底を用いて2電子系2原子分子に対する
 * Hartree-Fock計算を行う関数
 *
 * IOP=0 何も表示しない
 * IOP=1 収束の結果のみ表示
 * IOP=2 全繰り返しの結果を表示
 * N STO-NG を使う (N=1, 2, 3)
 * R 結合長 (a.u.)
 * ZETA1 原子Aの基底関数のSlater軌道指数値
 * ZETA2 原子Bの基底関数のSlater軌道指数値
 * ZA 原子Aの原子番号
 * ZB 原子Bの原子番号
 */
void HFCALC(int IOP, int N, double R, double ZETA1, double ZETA2, double ZA,
            double ZB) {
  if (IOP != 0) {
    printf("   STO-%dG for atomic numbers %5.2lf and %5.2lf\n\n\n", N, ZA, ZB);
  }
  INTGRL(IOP, N, R, ZETA1, ZETA2, ZA, ZB);
  COLECT(IOP, N, R, ZETA1, ZETA2, ZA, ZB);
  SCF(IOP, N, R, ZETA1, ZETA2, ZA, ZB);
}

/* SCF計算に必要な積分を計算する関数*/
void INTGRL(int IOP, int N, double R, double ZETA1, double ZETA2, double ZA,
            double ZB) {
  double R2 = R * R;
  // 軌道指数のスケーリングと短縮係数の規格化
  for (int i = 0; i < N; ++i) {
    A1[i] = EXPON[i][N-1] * pow(ZETA1, 2.0);
    D1[i] = COEF[i][N-1] * pow(2.0 * A1[i] / PI, 0.75);
    A2[i] = EXPON[i][N-1] * pow(ZETA2, 2.0);
    D2[i] = COEF[i][N-1] * pow(2.0 * A2[i] / PI, 0.75);
  }
  S12 = 0.0;
  T11 = 0.0;
  T12 = 0.0;
  T22 = 0.0;
  V11A = 0.0;
  V12A = 0.0;
  V22A = 0.0;
  V11B = 0.0;
  V12B = 0.0;
  V22B = 0.0;
  V1111 = 0.0;
  V2111 = 0.0;
  V2121 = 0.0;
  V2211 = 0.0;
  V2221 = 0.0;
  V2222 = 0.0;

  // 1電子積分の計算
  // V12A = 原子核Aへの各引力の非対角項
  // RAP2 = 中心Aと中心Pの距離の2乗
  for (int i = 0; i < N; ++i) {
    for (int j = 0; j < N; ++j) {
      double RAP = A2[j] * R / (A1[i] + A2[j]);
      double RAP2 = pow(RAP, 2.0);
      double RBP2 = pow(R - RAP, 2.0);
      S12 += calcS(A1[i], A2[j], R2) * D1[i] * D2[j];
      T11 += calcT(A1[i], A1[j], 0.0) * D1[i] * D1[j];
      T12 += calcT(A1[i], A2[j], R2) * D1[i] * D2[j];
      T22 += calcT(A2[i], A2[j], 0.0) * D2[i] * D2[j];
      V11A += calcV(A1[i], A1[j], 0.0, 0.0, ZA) * D1[i] * D1[j];
      V12A += calcV(A1[i], A2[j], R2, RAP2, ZA) * D1[i] * D2[j];
      V22A += calcV(A2[i], A2[j], 0.0, R2, ZA) * D2[i] * D2[j];
      V11B += calcV(A1[i], A1[j], 0.0, R2, ZB) * D1[i] * D1[j];
      V12B += calcV(A1[i], A2[j], R2, RBP2, ZB) * D1[i] * D2[j];
      V22B += calcV(A2[i], A2[j], 0.0, 0.0, ZB) * D2[i] * D2[j];
    }
  }
  // 2電子積分の計算
  for (int i = 0; i < N; ++i) {
    for (int j = 0; j < N; ++j) {
      for (int k = 0; k < N; ++k) {
        for (int l = 0; l < N; ++l) {
          double RAP = A2[i] * R / (A2[i] + A1[j]);
          //double RBP = R - RAP;
          double RAQ = A2[k] * R / (A2[k] + A1[l]);
          double RBQ = R - RAQ;
          double RPQ = RAP - RAQ;
          double RAP2 = RAP * RAP;
          //double RBP2 = RBP * RBP;
          //double RAQ2 = RAQ * RAQ;
          double RBQ2 = RBQ * RBQ;
          double RPQ2 = RPQ * RPQ;
          V1111 += TWOE(A1[i], A1[j], A1[k], A1[l], 0.0, 0.0, 0.0) 
                   * D1[i] * D1[j] * D1[k] * D1[l];
          V2111 += TWOE(A2[i], A1[j], A1[k], A1[l], R2, 0.0, RAP2)
                   * D2[i] * D1[j] * D1[k] * D1[l];
          V2121 += TWOE(A2[i], A1[j], A2[k], A1[l], R2, R2, RPQ2)
                   * D2[i] * D1[j] * D2[k] * D1[l];
          V2211 += TWOE(A2[i], A2[j], A1[k], A1[l], 0.0, 0.0, R2)
                   * D2[i] * D2[j] * D1[k] * D1[l];
          V2221 += TWOE(A2[i], A2[j], A2[k], A1[l], 0.0, R2, RBQ2)
                   * D2[i] * D2[j] * D2[k] * D1[l];
          V2222 += TWOE(A2[i], A2[j], A2[k], A2[l], 0.0, 0.0, 0.0) * D2[i] *
                   D2[j] * D2[k] * D2[l];
        }
      }
    }
  }
  if (IOP != 0) {
    printf("   R          ZETA1      ZETA2      S12       T11\n\n");
    printf("%11.6lf%11.6lf%11.6lf%11.6lf%11.6lf\n\n\n", R, ZETA1, ZETA2, S12, T11);
    printf("   T12        T22        V11A       V12A        V22A\n\n");
    printf("%11.6lf%11.6lf%11.6lf%11.6lf%11.6lf\n\n\n", T12, T22, V11A, V12A, V22A);
    printf("   V11B       V12B       V22B       V1111        V2111\n\n");
    printf("%11.6lf%11.6lf%11.6lf%11.6lf%11.6lf\n\n\n", V11B, V12B, V22B, V1111,
           V2111);
    printf("   V2121      V2211       V2221     V2222\n\n");
    printf("%11.6lf%11.6lf%11.6lf%11.6lf\n", V2121, V2211, V2221, V2222);
  }
}

// 重なり積分
double calcS(double A, double B, double RAB2) {
  return pow(PI / (A + B), 1.5) * exp(-A * B * RAB2 / (A + B));
}

// 運動エネルギー積分
double calcT(double A, double B, double RAB2) {
  return A * B / (A + B) * (3.0 - 2.0 * A * B * RAB2 / (A + B)) *
         pow(PI / (A + B), 1.5) * exp(-A*B*RAB2/(A+B));
}

// F0関数
double F0(double ARG) {
  double ret;
  if (ARG < 1.0e-6) {
    ret = 1.0 - ARG / 3.0;
  } else {
    ret = sqrt(PI / ARG) * erf(sqrt(ARG)) / 2.0;
  }
  return ret;
}

// 核引力積分
double calcV(double A, double B, double RAB2, double RCP2, double ZC) {
  double v = 2.0 * PI / (A + B) * F0((A + B) * RCP2) * exp(-A * B * RAB2 / (A + B));
  return -v * ZC;
}

// 2電子積分
double TWOE(double A, double B, double C, double D, double RAB2, double RCD2,
            double RPQ2) {
  return 2.0 * pow(PI, 2.5) / ((A + B) * (C + D) * sqrt(A + B + C + D))
         * F0((A + B) * (C + D) * RPQ2 / (A + B + C + D))
         * exp(-A * B * RAB2 / (A + B) - C * D * RCD2 / (C + D));
}

// 積分の結果から行列を求める
void COLECT(int IOP, int N, double R, double ZETA1, double ZETA2, double ZA,
            double ZB) {
  // 核ハミルトニアン cf. 式(3.153)
  H[0][0] = T11 + V11A + V11B;
  H[0][1] = T12 + V12A + V12B;
  H[1][0] = H[0][1];
  H[1][1] = T22 + V22A + V22B;

  // 重なり行列 cf. 式(3.136)
  S[0][0] = 1.0;
  S[0][1] = S12;
  S[1][0] = S[0][1];
  S[1][1] = 1.0;

  // 正準直交化 cf. 式(3.169), p.191
  X[0][0] = 1.0 / sqrt(2.0 * (1.0 + S12));
  X[1][0] = X[0][0];
  X[0][1] = 1.0 / sqrt(2.0 * (1.0 - S12));
  X[1][1] = -X[0][1];

  // 変換行列の転置
  XT[0][0] = X[0][0];
  XT[0][1] = X[1][0];
  XT[1][0] = X[0][1];
  XT[1][1] = X[1][1];

  // 2電子積分の行列
  TT[0][0][0][0] = V1111;
  TT[1][0][0][0] = V2111;
  TT[0][1][0][0] = V2111;
  TT[0][0][1][0] = V2111;
  TT[0][0][0][1] = V2111;
  TT[1][0][1][0] = V2121;
  TT[0][1][1][0] = V2121;
  TT[1][0][0][1] = V2121;
  TT[0][1][0][1] = V2121;
  TT[1][1][0][0] = V2211;
  TT[0][0][1][1] = V2211;
  TT[1][1][1][0] = V2221;
  TT[1][1][0][1] = V2221;
  TT[1][0][1][1] = V2221;
  TT[0][1][1][1] = V2221;
  TT[1][1][1][1] = V2222;

  if (IOP != 0) {
    MATOUT(S, 2, 2, 2, 2, "S   ");
    MATOUT(X, 2, 2, 2, 2, "X   ");
    MATOUT(H, 2, 2, 2, 2, "H   ");
    printf("\n");
    for (int i = 0; i < 2; ++i) {
      for (int j = 0; j < 2; ++j) {
        for (int k = 0; k < 2; ++k) {
          for (int l = 0; l < 2; ++l) {
            printf("   (%2d%2d%2d%2d )%10.6lf\n", i+1, j+1, k+1, l+1, TT[i][j][k][l]);
          }
        }
      }
    }
  }
}

// SCFを実行する
void SCF(int IOP, int N, double R, double ZETA1, double ZETA2, double ZA,
         double ZB) {
  // 密度行列の収束規準
  double CRIT = 1.0e-4;
  // 繰り返しの最大回数
  int MAXIT = 25;
  // 核ハミルトニアンを最初の推定に用いる(つまりPを零行列で初期化)
  for (int i = 0; i < 2; ++i) {
    for (int j = 0; j < 2; ++j) {
      P[i][j] = 0.0;
    }
  }
  if (IOP == 2) {
    MATOUT(P, 2, 2, 2, 2, "P    ");
  }
  int ITER = 0;
  double EN;
  for (;;) {
    ITER++;
    if (IOP == 2) {
      printf("\n    START OF ITERATION NUMBER = %2d\n", ITER);
    }
    // Fock行列の2電子部分を計算
    FORMG();

    if (IOP == 2) {
      MATOUT(G, 2, 2, 2, 2, "G    ");
    }

    // 核ハミルトニアンを足してFock行列を得る
    for (int i = 0; i < 2; ++i) {
      for (int j = 0; j < 2; ++j) {
        F[i][j] = H[i][j] + G[i][j];
      }
    }

    // 電子エネルギーを計算する
    EN = 0.0;
    for (int i = 0; i < 2; ++i) {
      for (int j = 0; j < 2; ++j) {
        EN += 0.5 * P[i][j] * (H[i][j] + F[i][j]);
      }
    }

    if (IOP == 2) {
      MATOUT(F, 2, 2, 2, 2, "F   ");
      printf("\n\n\n    ELECTRONIC ENERGY = %20.12lf\n", EN);
    }

    // Gを用いてFock行列を変換
    MULT(F, X, G, 2, 2);
    MULT(XT, G, FPRIME, 2, 2);
    // Fock行列の対角化
    DIAG(FPRIME, CPRIME, E);
    // 固有ベクトルから行列Cを得る
    MULT(X, CPRIME, C, 2, 2);
    // 新しい密度行列を計算する
    for (int i = 0; i < 2; ++i) {
      for (int j = 0; j < 2; ++j) {
        OLDP[i][j] = P[i][j];
        P[i][j] = 0.0;
        for (int k = 0; k < 1; ++k) { // ここはk<1で本当にいいのか不明
          P[i][j] += 2.0 * C[i][k] * C[j][k];
        }
      }
    }
    if (IOP == 2) {
      MATOUT(FPRIME, 2, 2, 2, 2, "F'   ");
      MATOUT(CPRIME, 2, 2, 2, 2, "C'   ");
      MATOUT(E, 2, 2, 2, 2, "E   ");
      MATOUT(C, 2, 2, 2, 2, "C   ");
      MATOUT(P, 2, 2, 2, 2, "P   ");
    }
    // deltaを計算する
    double DELTA = 0.0;
    for (int i = 0; i < 2; ++i) {
      for (int j = 0; j < 2; ++j) {
        DELTA += pow(P[i][j] - OLDP[i][j], 2.0);
      }
    }
    DELTA = sqrt(DELTA / 4.0);
    if (IOP != 0) {
      printf("\n    DELTA(CONVERGENCE OF DENSITY MATRIX) = %10.6lf\n", DELTA);
    }

    // 収束判定
    if (DELTA < CRIT) break;
    if (ITER < MAXIT) continue;
    // 収束しなかった場合
    printf("    NO CONVERGENCE IN SCF");
    return;
  }

  double ENT = EN + ZA * ZB / R;
  if (IOP != 0) {
    printf("\n\n    CALCULATION CONVERGED"
           "\n\n    ELECTORONIC ENERGY = %20.12lf"
           "\n\n    TOTAL ENERGY =       %20.12f", EN, ENT);
  }
  if (IOP == 1) {
    MATOUT(G, 2, 2, 2, 2, "G   ");
    MATOUT(F, 2, 2, 2, 2, "F   ");
    MATOUT(E, 2, 2, 2, 2, "E   ");
    MATOUT(C, 2, 2, 2, 2, "C   ");
    MATOUT(P, 2, 2, 2, 2, "P   ");
  }
  MULT(P, S, OLDP, 2, 2);
  if (IOP != 0) {
    MATOUT(OLDP, 2, 2, 2, 2, "PS  ");
  }
}

void FORMG(void) {
  for (int i = 0; i < 2; ++i) {
    for (int j = 0; j < 2; ++j) {
      G[i][j] = 0.0;
      for (int k = 0; k < 2; ++k) {
        for (int l = 0; l < 2; ++l) {
          G[i][j] += P[k][l] * (TT[i][j][k][l] - 0.5 * TT[i][l][k][j]);
        }
      }
    }
  }
}

// Fを対角化し、Cに固有ベクトル、Eに固有値を入れる
void DIAG(double F[][2], double C[][2], double E[][2]) {
  double theta;
  if (fabs(F[0][0] - F[1][1]) > 1.0e-20) {
    theta = 0.5 * atan(2.0 * F[0][1] / (F[0][0] - F[1][1]));
  } else {
    theta = PI / 4.0;
  }
  C[0][0] = cos(theta);
  C[1][0] = sin(theta);
  C[0][1] = sin(theta);
  C[1][1] = -cos(theta);
  E[0][0] = F[0][0] * pow(cos(theta), 2.0) + F[1][1] * pow(sin(theta), 2.0)
            + F[0][1] * sin(2.0 * theta);
  E[1][1] = F[1][1] * pow(cos(theta), 2.0) + F[0][0] * pow(sin(theta), 2.0)
            - F[0][1] * sin(2.0 * theta);
  E[1][0] = 0.0;
  E[0][1] = 0.0;
  // 固有値と固有ベクトルを整列する
  if (E[1][1] > E[0][0]) return;
  double temp = E[1][1];
  E[1][1] = E[0][0];
  E[0][0] = temp;
  temp = C[0][1];
  C[0][1] = C[0][0];
  C[0][0] = temp;
  temp = C[1][1];
  C[1][1] = C[1][0];
  C[1][0] = temp;
}

// AとBの積をCに入れる
void MULT(double A[][2], double B[][2], double C[][2], int IM, int M) {
  for (int i = 0; i < M; ++i) {
    for (int j = 0; j < M; ++j) {
      C[i][j] = 0.0;
      for (int k = 0; k < M; ++k) {
        C[i][j] += A[i][k] * B[k][j];
      }
    }
  }
}

void MATOUT(double A[][2], int IM, int IN, int M, int N, char *LABEL) {
  int IHIGH = 0;
  for (;;) {
    int LOW = IHIGH + 1;
    IHIGH = IHIGH + 5;
    IHIGH = (IHIGH < N) ? IHIGH : N;
    printf("\n\n\n    THE %s ARRAY\n", LABEL);
    printf("               ");
    for (int i = LOW; i <= IHIGH; ++i) {
      printf("          %3d      ", i);
    }
    printf("\n");
    for (int i = 1; i <= M; ++i) {
      printf("%10d     ", i);
      for (int j = LOW; j <= IHIGH; ++j) {
        printf(" %18.10lf", A[i - 1][j - 1]);
      }
      printf("\n");
    }
    if (N <= IHIGH) break;
  }
}

最終的な出力は

    CALCULATION CONVERGED

    ELECTORONIC ENERGY =      -4.227529304014

    TOTAL ENERGY =            -2.860662163500


    THE PS   ARRAY
                           1                  2      
         1            1.5296370381       1.1199277981
         2            0.6424383868       0.4703629619

でした。参照したFortranコードの出力は

    CALCULATION CONVERGED

    ELECTRONIC ENERGY =  -0.422752913203D+01

    TOTAL ENERGY =       -0.286066199152D+01



    THE PS   ARRAY
                           1                  2
         1        0.1529637121D+01   0.1119927796D+01
         2        0.6424383099D+00   0.4703628793D+00

でした。

FortranからCに移植するときのtips

詰まったことを挙げておきます。

配列の添字

Fortranでは添字が1から始まるので配列にアクセスするときに添字をうまく調整する必要があります。

多次元配列の並び順

このページに詳しく書いてありますが、C言語の多次元配列はメモリ上でa[0][0], a[0][1], a[1][0], a[1][1], a[2][0], a[2][1]のような順番で並びますが、Fortranではa(1,1), a(2,1), a(3,1), a(1,2), a(2,2), a(3,2)のような順番に並びます。多次元配列の初期化を書き直すときに気をつける必要があります。

整合配列

Fortranではサブルーチンに配列を渡した配列のサイズをサブルーチン側で定義することができます。(整合配列という)

      dimension a(6)
      a(1) = 1
      a(2) = 2
      a(3) = 3
      a(4) = 4
      a(5) = 5
      a(6) = 6
      call sub(a,3,2)
      end

      subroutine sub(a,n,m)
      dimension a(n,m) <---- a(n,*) でも同じ結果を得る
      do i=1,n
         write(6,*) (a(i,j),j=1,m)
      end do
      return
      end

コードはFortran プログラミングの基礎知識から引用しました。

C99で導入された可変長配列の機能を使うと同じことができるそうです。cf. 3.3 整合配列、可変長配列
今回は配列のサイズが決め打ちできたのでvoid MULT(double A[][2], double B[][2], double C[][2], int IM, int M) { のように関数を定義しました。

C言語のabs

特にFortranとは関係がないのですが、C言語のabsはint型に対して使う関数で、double型に対してはfabsという別の関数が用意されています。(これはC言語で関数オーバーロードができないことが原因です。C11で導入されたGenericという機能を使うとそれっぽいことが実現できるらしいです。cf. Generic - C言語入門

2019年2月

2月もあっという間に終わってしまいました。

今月のトピック

今月読んだ本

今月は論文を書くのに忙しくてあまり本を読めませんでした(大嘘)

天冥の標X Part3

全ての力を尽くして天冥の標シリーズをオススメする - 基本読書を読んだのをきっかけに追い始めたシリーズの完結篇。
2013年1月の記事なので読み始めてから6年間経ったことになる。

2013年というと高校2年生で、当時の自分は医学部に行こうと思っていたはずのだが、人類の軌跡を宇宙規模の争いと絡めながら語っていくスケールの大きなお話を読んでしまった結果、自分がやりたいのは人類の科学を進めて歴史を動かすことであってどうせ死ぬ人間をほんの少し延命させることではないななどとうっかり思ってしまい(このシリーズの主人公は医者なのに!)、医学部志望をやめてよく分からない学科を受験してしまったので、私の人生は天冥の標に狂わされたと言っていいだろう。
そんな思い入れの深いシリーズが予定通り10巻で終わってしまうことになり完結篇を読むのが非常に怖かったのだが、楽しみにしていてよかったと思える結末に辿りついてくれていた。

1巻は導入としてだけではなく植民星での革命劇として純粋に面白かったし、2巻の感染症との戦いは素晴らしく、3巻の世界観や酸素要らずの価値観には人生を破壊され、4巻はおいといて、5巻のスケールの大きさには度肝を抜かれ、6巻での戦いに圧倒され、7巻での絶望に心を打たれ、8巻での伏線回収に感心し、9,10巻での物語の畳み方には舌を巻いた素晴らしいシリーズであった。最終巻のあとがきに「読者の何割かに多少のエネルギーと二、三個の発見を提供できたかもしれない」とあったが、少なくとも私には大変なエネルギーを与えてくれる作品だった。

かぐや様は告らせたい 1〜13巻

アニメ放映に合わせてBookwalkerで80%オフキャンペーンをやっていたので全巻購入した。「天才たちの恋愛頭脳戦」をやっていた最初の頃は正直しょうもないことやってるなと思っていまいちだったのだが、5巻収録の「花火の音は聞こえない」以降、かぐや様が天才設定を捨てたベタ惚れポンコツ乙女になってからは毎回素晴らしかった。恋愛ものはあまり好きじゃないと思っていたが、単に面白い作品に出会ってないだけだったのではと心変わりした。

五等分の花嫁 1〜8巻

かぐや様を読んでラブコメもいいなと思っていたところに80%オフキャンペーンをやっていたのを見たので購入。性格悪いヒロインが多くてあまり好きになれなかったのだが、7巻あたりで皆がデレ出してからは面白く、8巻での二乃の一転攻勢は素晴らしくてここまで読んでよかったと思える出来だった。(ここまで書いてきて思ったが単に自分はデレデレヒロインが好きなだけなのでは……)

山田エルフ大先生の恋する純真ごはん

kivantium花嫁修業の元ネタである山田エルフ先生が料理を作る漫画。連載されると聞いてからはマガジンWALKERに契約して毎月電撃大王を読んでいたので全て読んだ話ではあったがお布施。

みらいのふうふですけど?

みらいのふうふですけど?(2) (百合姫コミックス)

みらいのふうふですけど?(2) (百合姫コミックス)

中学生百合。二人のコミカルなやりとりが素晴らしい作品だったが2巻で完結してしまった。作者が途上のIPコンテンツと言っているので打ち切りっぽい。悲しい。

お兄ちゃんはおしまい!

お兄ちゃんはおしまい!   (2) (IDコミックス)

お兄ちゃんはおしまい! (2) (IDコミックス)

妹に女の子にされてしまったニートのお兄ちゃんの話。1巻を買ったので2巻も。

GUNSLINGER GIRL

セールをやっていたので1巻だけ購入。全体に漂う暗い雰囲気や諦念が好みじゃなかったので2巻以降はまたの機会に。

量子化学

量子化学 上巻

量子化学 上巻

ラボの人に勧められたので読んだ。仮定と近似がはっきり書かれていており、説明も省略が少なかったので分かりやすかった。

2019年1月

2018年12月分も含みます。

今月のトピック

行列の微分を求めるサイト

NIPS 2018で発表された成果らしい: Computing Higher Order Derivatives of Matrix and Tensor Expressions

Vimプラグイン ale

flaskなどのlinterをVim上で走らせることができる上に、保存したときにインデントを自動で修正することもできる。便利。

その他

https://twitter.com/tmaehara/status/937385784440582144



読んだ本

量子計算理論 量子コンピュータの原理

量子計算理論 量子コンピュータの原理

量子計算の理論に関する本。古典確率計算を量子計算のようなベクトル表示で表すことで対比する導入は面白かった。
物理的な実現方法やShorのアルゴリズムのような具体的な量子アルゴリズムは扱っておらず、「量子コンピュータを持っていると主張しているQ-Wave社が本当に量子計算を行っていると検証するにはどうすればよいか」のような理論について扱っていた。論文のダイジェストという感じですんなり読める感じではなかった。
以前読んだ量子コンピュータ入門(第2版)のほうが一冊目としては良いだろう。

ここまでの9巻で広げた風呂敷を回収していくのがすごく良い。

まほろまてぃっくの作者の新作。ママ(小学生)がかわいい。

三ツ星カラーズ6 (電撃コミックスNEXT)

三ツ星カラーズ6 (電撃コミックスNEXT)

ハミングガール (百合姫コミックス)

ハミングガール (百合姫コミックス)

2018年11月

今月のトピック

PRMLの無料公開

PRML以外で無料公開されている有名本・PDFを挙げておく

いろいろな法則とか

こういうのは名前を覚えておかないとドヤ顔で披露できないので記憶する必要がある。

Content-Disposition レスポンスヘッダー

ブラウザでPDFを開いたときに、ブラウザで開かれる場合とダウンロードされる場合があるが、この挙動はHTTPヘッダーで指定できる。
HTML5のdownload属性で指定される場合やJavaScriptでダウンロードされる場合もある)

その他



2018年10月

今月知った有益情報

Multi-service image search

複数サービスを横断して画像検索ができる

UCL Course on RL

http://www0.cs.ucl.ac.uk/staff/d.silver/web/Teaching.html
強化学習のオンライン講座。誰かが分かりやすいと言っていた。

The Architecture of Open Source Applications

著名なOSSプロジェクトの内部について簡単に説明していく本らしい。

古地図

いろんな古地図を見ることができる

Exploit-Exercises: Protostar (v2)

CTFのpwn問題の勉強におすすめの問題とのこと

電子サイコロの回路出典

4015を1つ使うだけでサイコロを作れる回路の出典は「初歩のラジオ1980年10月号に松本悟さんが発表された記事」とのこと

読んだ本

量子コンピュータ入門(第2版)

量子コンピュータ入門(第2版)

ゲート型量子コンピュータの解説本。今まで読んだ量子コンピュータ本の中で一番わかりやすかった。

青葉ちゃんと紅葉ちゃんのコンペ対決。人々が頑張るのはよい。

まちカドまぞく (4) (まんがタイムKRコミックス)

まちカドまぞく (4) (まんがタイムKRコミックス)

4コマオブザイヤーで3年連続1位になっている作品。不思議な魅力がある。4巻ではギャグ描写だと思っていたものが展開上重要なアイテムになる展開に感心した。

3巻までは能力を持たないナナが能力者たちと対決する話として非常によくできており、4巻では少し方向性を変えつつナナの過去や能力者との友情が描かれていた。面白かったが、展開や新キャラからそこはかとない車輪の国っぽさを感じてしまった。

零號琴

零號琴

飛 浩隆の16年ぶりの長編。巨大楽器の竣工を記念した仮面舞踏会という面白くなさそうなあらすじだったが、高い描写力や緻密な世界構築に圧倒された。さすがにこれはふざけすぎではとか展開に納得感がなかったなぁなどと思わないでもない。

RISC-V原典  オープンアーキテクチャのススメ

RISC-V原典 オープンアーキテクチャのススメ

RISC-Vの命令セットについての本。ちょっとだけお気持ちが書いてあるが基本的にデータシートと同じ内容が書いてあるので買うほどではないかもしれない。

特定商取引法に定められた事項は請求により遅滞なく提供する