久々に少々長いプログラムを書いた

論文を投稿したのだが、１ヶ月前に４人のレフェリー全員から好意的な結果が返ってきた。誰もrejectを示唆しなかったのはこれまでの経験でも殆ど無いし嬉しい。
とはいうものの、いずれのコメントでも修正点や質問が沢山あり、答えにくい質問も含まれていた。ちゃんとしたデータを出して認めてもらうことが必要であり、ないがしろにすると蹴られる。

大体の問題はボスの知識と１ヶ月に渡る文献調査とGaussian09を使った簡単な計算で解決されたのだが、質問のうちの１つがどうしても既存の文献だとサポートしきれず、新しくプログラムを書き下ろして解析することとなった。結果的に全部で1000行弱のものになった。

今回あったことと雑感

・今まではずっとPython 2系を使って書いてきたのだが、流石にサポートの問題とかもあるため、前々からPython 3に切り替えたいと思っていた。いい機会なので練習がてら使ってみることにした。print文が関数になっており戸惑ったり、rangeをリストとして利用できないなどの点に困ったが、慣れると案外使いやすい。

・初めて最適化アルゴリズムを採用した。この記事（最適化アルゴリズム - sonoshouのまじめなブログ）に従った。webにこういう記事があるありがたさを感じる。最初は焼きなまし法を採用したのだが、確率的に他の不安定状態にうつる箇所があるとどうしても上手く行かなかったため、ステップ数の増加に従ってパラメータの変位が小さくなる点を残してあとは自己流で適当に作った。Local minimumに落ちるのを回避するため、同じ計算を数度繰り返し、直近数回の誤差が１００万分の１とかになった時のみ採用するという条件を課した。39次元の変数の最適化を行ったが、この方法だと案外うまくいった。

・自分の不注意によるミスに多く出くわした。特に三角関数の微分で時間の係数を三角関数にかけたとき、ラジアンの変換に使う2πをかけ忘れててとんでもない値が出てた時は間違いに気づくまでに一晩かかった（寝て起きたら気づいた）。

割と楽しかったのでこれからもプログラミングを多用して研究していきたい。

ミネソタ汎関数に分散力補正を入れること等についての雑感

最近になって、ミネソタ汎関数に分散力補正を入れた論文を何報か立て続けに見た。違和感を覚えたのでメモしておきたい。
そもそも分散力補正は開発当時の既存の汎関数(B3LYPなど)が分子間力をうまく取り入れることができないという欠点を補うために開発された方法であり、すでに分子間力をうまく表現できるミネソタ汎関数に追加する必要はないはずである。
（ミネソタ汎関数の分子間力の再現性については下記の論文が分かりやすい
http://pubs.acs.org/doi/abs/10.1021/ar700111a）

おそらく、無理に分散力補正をつけた方法では、分子間力を大きく見積もりすぎて分子を無理やり近づける結果になっていることであろう。構造、エネルギーともにその信頼性に疑問が持たれる。
可能ならば、ベンチマークテストを行って結果を比較したいと思っている。

力技でSCF計算

SCF計算がまともに収束しないので力技で網羅的に計算を行った。
計算対象は伏せておくが、使ったオプションと結果は以下のとおり。
論文発表前だし、何かあると面倒なのでここに挙げたデータはちょっといじってある。
2つ計算値がある場合は下のほうが最終的な計算値。上はアルゴリズム切り替え前の計算値。

オプション：
scf=(xqc, maxconv=1~7,vtl)
scf=(xqc, maxconv=1~7)
scf=(xqc,maxconv=1,novaracc)
scf=yqc
scf=ssd
scf=(xqc,maxconv=1,vtl,vshift=50 or 1000 or 10000)


結果：
# UM06/6-31+G** scf=(xqc,maxconv=1,vtl)
SCF Done: E(UM06) = -555.736771470 A.U. after 2 cycles
SCF Done: E(UM06) = -555.994301760 a.u. after 19 cycles
Job cpu time: 0 days 0 hours 11 minutes 10.9 seconds.

# UM06/6-31+G** scf=(xqc,maxconv=2,vtl)
SCF Done: E(UM06) = -555.815032856 A.U. after 3 cycles
SCF Done: E(UM06) = -555.992542926 a.u. after 21 cycles
Job cpu time: 0 days 0 hours 7 minutes 59.8 seconds.

# UM06/6-31+G** scf=(xqc,maxconv=3,vtl)
SCF Done: E(UM06) = -555.957790071 A.U. after 4 cycles
SCF Done: E(UM06) = -555.992542926 a.u. after 20 cycles
Job cpu time: 0 days 0 hours 7 minutes 48.2 seconds.

# UM06/6-31+G** scf=(xqc,maxconv=4,vtl)
SCF Done: E(UM06) = -555.981373939 A.U. after 5 cycles
SCF Done: E(UM06) = -555.992542926 a.u. after 15 cycles
Job cpu time: 0 days 0 hours 7 minutes 14.7 seconds.

# UM06/6-31+G** scf=(xqc,maxconv=5,vtl)
SCF Done: E(UM06) = -555.983189973 A.U. after 6 cycles
SCF Done: E(UM06) = -555.992542926 a.u. after 18 cycles
Job cpu time: 0 days 0 hours 8 minutes 18.0 seconds.

# UM06/6-31+G** scf=(xqc,maxconv=6,vtl)
SCF Done: E(UM06) = -555.984044437 A.U. after 7 cycles
SCF Done: E(UM06) = -555.992542926 a.u. after 8 cycles
Job cpu time: 0 days 0 hours 7 minutes 13.7 seconds.

# UM06/6-31+G** scf=(xqc,maxconv=7,vtl)
SCF Done: E(UM06) = -555.984835916 A.U. after 8 cycles
SCF Done: E(UM06) = -555.992542926 a.u. after 8 cycles
Job cpu time: 0 days 0 hours 6 minutes 8.0 seconds.

# UM06/6-31+G** scf=(xqc,maxconv=1)
SCF Done: E(UM06) = -555.736771470 A.U. after 2 cycles
SCF Done: E(UM06) = -555.994301760 a.u. after 18 cycles
Job cpu time: 0 days 0 hours 7 minutes 33.3 seconds.

# UM06/6-31+G** scf=(xqc,maxconv=2)
SCF Done: E(UM06) = -555.815032856 A.U. after 3 cycles
SCF Done: E(UM06) = -555.992542926 a.u. after 21 cycles
Job cpu time: 0 days 0 hours 7 minutes 54.9 seconds.

# UM06/6-31+G** scf=(xqc,maxconv=3)
SCF Done: E(UM06) = -555.957790071 A.U. after 4 cycles
SCF Done: E(UM06) = -555.992542926 a.u. after 20 cycles
Job cpu time: 0 days 0 hours 8 minutes 4.7 seconds.

# UM06/6-31+G** scf=(xqc,maxconv=4)
SCF Done: E(UM06) = -555.981373939 A.U. after 5 cycles
SCF Done: E(UM06) = -555.992542926 a.u. after 15 cycles
Job cpu time: 0 days 0 hours 7 minutes 48.0 seconds.

# UM06/6-31+G** scf=(xqc,maxconv=5)
SCF Done: E(UM06) = -555.983189973 A.U. after 6 cycles
SCF Done: E(UM06) = -555.992542926 a.u. after 18 cycles
Job cpu time: 0 days 0 hours 7 minutes 49.5 seconds.

# UM06/6-31+G** scf=(xqc,maxconv=6)
SCF Done: E(UM06) = -555.984044437 A.U. after 7 cycles
SCF Done: E(UM06) = -555.992542926 a.u. after 8 cycles
Job cpu time: 0 days 0 hours 6 minutes 21.0 seconds.

# UM06/6-31+G** scf=(xqc,maxconv=7)
SCF Done: E(UM06) = -555.984835916 A.U. after 8 cycles
SCF Done: E(UM06) = -555.992542926 a.u. after 8 cycles
Job cpu time: 0 days 0 hours 6 minutes 28.1 seconds.

# UM06/6-31+G** scf=(xqc,maxconv=1,novaracc)
SCF Done: E(UM06) = -555.736771470 A.U. after 2 cycles
SCF Done: E(UM06) = -555.994301760 a.u. after 18 cycles
Job cpu time: 0 days 0 hours 9 minutes 44.7 seconds.

# UM06/6-31+G** scf=(yqc)
SCF Done: E(UM06) = -555.991797197 a.u. after 7 cycles
SCF Done: E(UM06) = -555.992542926 A.U. after 14 cycles
Job cpu time: 0 days 0 hours 3 minutes 8.9 seconds.

# UM06/6-31+G** scf=ssd
SCF Done: E(UM06) = -555.992542926 a.u. after 37 cycles
Job cpu time: 0 days 0 hours 14 minutes 13.7 seconds.

# UM06/6-31+G** scf=(xqc,maxconv=1,vtl,vshift=1000)
SCF Done: E(UM06) = -555.736771470 A.U. after 2 cycles
SCF Done: E(UM06) = -555.992542926 a.u. after 37 cycles
Job cpu time: 0 days 0 hours 9 minutes 2.5 seconds.

# UM06/6-31+G** scf=(xqc,maxconv=1,vtl,vshift=10000)
SCF Done: E(UM06) = -555.736771470 A.U. after 2 cycles
SCF Done: E(UM06) = -555.992542927 a.u. after 381 cycles
Job cpu time: 0 days 1 hours 28 minutes 47.2 seconds.

# UM06/6-31+G** scf=(xqc,maxconv=1,vtl,vshift=50)
SCF Done: E(UM06) = -555.736771470 A.U. after 2 cycles
SCF Done: E(UM06) = -555.994301761 a.u. after 18 cycles
Job cpu time: 0 days 0 hours 11 minutes 28.8 seconds.



まず、最終的に得られた計算値がそれぞれの手法で異なっているのがわかる。これは極小値にハマって正確なSCF計算ができずに偽のエネルギー値を出しているからである。ここには載せてないがデフォルトのアルゴリズムだとSCF計算初期に現れた極小値に見事にハマってしまい、とんでもなく大きなエネルギーを出す。

scf=(xqc,maxconv=1,vtl)は今回の場合は信用できる手段のようだ（データを載せてないが、xqcの代わりにqc使っても高エネルギーで収束してしまう。）。というか、scf=(xqc,maxconv=1)以外の方法だと事実上このシステムは信頼できるエネルギーを計算できない。
エネルギーの信頼性は、構造を少し変化させつつポテンシャルエネルギー曲面全体をプロットすることで判断することができる。（間違った値が得られるとそこだけ曲線から外れる）

新アルゴリズムのyqcには期待してたがあんまり有効でないようだ。
また、vshiftの値を大きくすると折角maxconv=1を使っていても間違った値が出る。
デフォルト(vshift=100)か、あるいは小さめに設定するといい。



（追記）
Stableキーワードも有効っぽい
www.hpc.co.jp
しかしながら、現在自分が扱っているシステムでは
The wavefunction is stable under the perturbations considered.
The wavefunction is already stable.
という結果が出ても偽の値が出てくる。

英語ができない結果として専門知識が体系的に蓄積する

＜概略＞
ボスの英語で聞き取れない点がある
↓
概略は分かるが、その場でひとつひとつの単語を何度も聞き返すわけにはいかず、
わかったふりをして要点だけメモする
↓
後でキーワードを手がかりにして成書なりインターネットなりを使って言われたことを調べる
↓
会話の文脈を補完するために関連する基礎分野も調べる
↓
いつの間にか説明された内容以上に体系的な知識がついてる

こんな感じになってる。まあ聞き取れるに越したことはないんだけど、
ネイティブってスラング混ぜて手加減なしに喋るから多少は仕方ないよねと言った感じ。
まだ書籍化されていないような内容についてはお手上げで、後で再度聞きに行くハメになる。

筋トレには学びがある

筋肉を肥大させる際、闇雲にトレーニングすることは推奨されていない。
適切な方法で筋肉に負荷をかけ、十分な食事と休息をとらないと効果は現れない。
無理をすれば怪我をしてしまい、これまでの努力が無駄になることだってある。
また、努力しても遺伝的に素質のある人間にはかなわない。
周囲にいる骨格からして異なるアメリカ人と自分を比較しても意味はなく、比較すべきは過去の自分である。

自分はこの考え方が好きであるし、仕事なり何なりもこのような思想で運用するほうがいいと思っている。長時間何かに打ち込むことはスキルアップの上で必要だが、負担に感じた時は十分な食事や睡眠をとらないと精神や肉体を壊してしまい、これまで積み上げてきたものが崩れてしまう。とてつもない成果を出している超人の真似をする必要はないし、当然それを他人に強いてはいけない。彼らは誰も真似出来ないからこその超人なのだから。

一昨年、調子に乗って長距離を走っていたら軽いランナー膝になってしまったのだが、今でも膝に負荷がかかるとそれが再発することがある。
ひょっとすると、これから一生この痛みの再発への不安を抱えながら走ることになるかもしれない。
こういうことのないように適度にやって行きたい。