力技で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.
という結果が出ても偽の値が出てくる。

(構造最適化チェック|Supporting Information作成)プログラム

以前書いたAWKのプログラム
moki.hatenablog.com
の改良版@python。Gaussian09のアウトプットファイルを読み込んで、欲しい情報だけ抜き出す。
大学の共用サーバーを4種類使い分けてるのだが、全部に共通して入ってるpython2.5上で、numpy、scipyなしで動くように書いた。またpython3じゃ動かない(print文の文法に全くこだわっていないためなのだが)。また、他のプログラムを流用して作ったのでいろんな関数名やフラグ管理が適当である。
Logreader.pyと名前をつけて使っている。

"python Logreader.py 計算結果.log "

と実行すると構造最適化のチェック、

"python Logreader.py 計算結果.log SI "

と実行するとSupporting Information用の結果を出力する。

Logreader.py

#!/bin/python
import linecache
import sys

#read commandline
param = sys.argv

#setup parameters
number_atom  = {'1':'H','2':'He','3':'Li','4':'Be','5':'B','6':'C','7':'N','8':'O','9':'F','10':'Ne',
                '11':'Na','12':'Mg','13':'Al','14':'Si','15':'P','16':'S','17':'Cl','18':'Ar','19':'K','20':'Ca',
                '21':'Sc','22':'Ti','23':'V','24':'Cr','25':'Mn','26':'Fe','27':'Co','28':'Ni','29':'Cu','30':'Zn',
                '31':'Ga','32':'Ge','33':'As','34':'Se','35':'Br','36':'Kr','37':'Rb','38':'Sr','39':'Y','40':'Zr',
                '41':'Nb','42':'Mo','43':'Tc','44':'Ru','45':'Rh','46':'Pd','47':'Ag','48':'Cd','49':'In','50':'Sn',
                '51':'Sb','52':'Te','53':'I','54':'Xe','55':'Cs','56':'Ba'}

linenum  = sum(1 for line in open(param[1]))
#tempfile = open('templog','w')

Maximum_Force        = "Null"
RMS_Force            = "Null"
Maximum_Displacement = "Null"
RMS_Displacement     = "Null"
Firstfreq            = "Null"
Secondfreq           = "Null"
OptChk  = 1
FreqChk = 1
OptFlag      = 0
PrintFlag    = 0
PrintedFlag  = 0
ChargeFlag   = 1

try:
    if param[2] == "SI":
        print("-------------paste below in SI--------------")
        for line in range(linenum):
            Linedata = linecache.getline(param[1],line).split()
            try:
                if Linedata[0] == "Frequencies" and Linedata[1] == "--" and FreqChk == 1:
                    Firstfreq  = Linedata[2]
                    Secondfreq = Linedata[3]
                    if float(Firstfreq) > 0:
                        print " Number of imaginary frequency = 0"
                    elif float(Firstfreq) < 0 and float(Secondfreq) > 0:
                        print " Number of imaginary frequency = 1"
                    else:
                        print " Number of imaginary frequency >= 2"
                    FreqChk   = 0
                if Linedata[0] == "Zero-point" and Linedata[1] == "correction=":
                    PrintFlag = 1
                if PrintFlag == 1:
                    Line = linecache.getline(param[1],line)
                    print Line,
                if Linedata[4] == "thermal" and Linedata[5] == "Free":
                    PrintFlag = 0
            except:
                pass

        EE = 0
        for line in range(linenum):
            Linedata = linecache.getline(param[1],line).split()
            try:
                if Linedata[0] == "SCF" and Linedata[1] == "Done:":
                    EE = Linedata[4]
            except:
                pass           
            if line == linenum-1:
                print " Electric Energy =", EE
                print ""
        for line in range(linenum):
            Linedata = linecache.getline(param[1],line).split()
            try:
                if Linedata[1] == "Optimized" and Linedata[2] == "Parameters":
                    OptFlag = 1
                    if OptChk == 1:
                        OptChk = 0
                if Linedata[0] == "Standard" and Linedata[1] ==  "orientation:":
                    PrintFlag = OptFlag
                if Linedata[0] == "Rotational":
                    if PrintFlag == 1:
                        PrintedFlag = 1
                    PrintFlag = 0
                if PrintFlag == 1 and Linedata[0].isdigit() == True and PrintedFlag == 0:
                    print(" "+number_atom[Linedata[1]]+","+Linedata[2]+","+ Linedata[3]+","+ Linedata[4]+","+Linedata[5])
            except: 
                pass
                
except:
        for line in range(linenum):
            Linedata = linecache.getline(param[1],line).split()
            try:
                if Linedata[0] == "Maximum" and Linedata[1] == "Force":
                    if Linedata[4] == 'YES':
                        Maximum_Force        = "YES"

                if Linedata[0] == "RMS" and Linedata[1] == "Force":
                    if Linedata[4] == 'YES':
                        RMS_Force            = "YES"

                if Linedata[0] == "Maximum" and Linedata[1] == "Displacement":
                    if Linedata[4] == 'YES':
                        Maximum_Displacement = "YES"
                if Linedata[0] == "RMS" and Linedata[1] == "Displacement":
                    if Linedata[4] == 'YES':
                        RMS_Displacement     = "YES"
                if Linedata[0] == "Input" and Linedata[1] ==  "orientation:":
                    Maximum_Force        = "NO"
                    RMS_Force            = "NO"
                    Maximum_Displacement = "NO"
                    RMS_Displacement     = "NO"
                if Linedata[0] == "Frequencies" and Linedata[1] == "--" and FreqChk == 1:
                    Firstfreq  = Linedata[2]
                    Secondfreq = Linedata[3]
                    print ""
                    if float(Firstfreq) > 0:
                        print "Number of imerginary frequency = 0"
                    elif float(Firstfreq) < 0 and float(Secondfreq) > 0:
                        print "Number of imerginary frequency = 1"
                    else:
                        print "Number of imerginary frequency >= 2"         
                    FreqChk   = 0
                if Linedata[1] == "Optimized" and Linedata[2] == "Parameters":
                    OptFlag = 1
                    if OptChk == 1:
                        print "Maximum_Force       ", Maximum_Force
                        print "RMS_Force           ", RMS_Force
                        print "Maximum_Displacement", Maximum_Displacement
                        print "RMS_Displacement    ", RMS_Displacement 
                        print ""
                        print "---Optimized Geometry---"
                        OptChk = 0
                if Linedata[0] == "Standard" and Linedata[1] ==  "orientation:":
                    PrintFlag = OptFlag
                if Linedata[0] == "Rotational":
                    if PrintFlag == 1:
                        PrintedFlag = 1
                    PrintFlag = 0
                if PrintFlag == 1 and Linedata[0].isdigit() == True and PrintedFlag == 0:
                    print number_atom[Linedata[1]]," ", Linedata[3]," ", Linedata[4]," ", Linedata[5]
            except:
                pass

何かの参考にどうぞ

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

<概略>
ボスの英語で聞き取れない点がある

概略は分かるが、その場でひとつひとつの単語を何度も聞き返すわけにはいかず、
わかったふりをして要点だけメモする

後でキーワードを手がかりにして成書なりインターネットなりを使って言われたことを調べる

会話の文脈を補完するために関連する基礎分野も調べる

いつの間にか説明された内容以上に体系的な知識がついてる

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

筋トレには学びがある

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

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

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

自転車を漕いでいると追い風には気付き難い。

漕ぐのをやめた時、風が吹いていた事に気づくんだなと思った。

雑感

こっちに来てから割と生活が楽になった。
学生の頃よりは実験時間は減ったし、それでなおインパクトのある成果が出やすい。

ボスが優秀すぎて、自分がやったことが成果に結びつきやすいというだけなのだが、
ここで出した成果を自分の能力だと思ってると後で痛い目に遭うだろう。

また、半年前に自分が「この分野の将来の課題」だと思っていたものがあったが、
今現在やってる研究とちょうど結びついてしまったので取り組むことになった。

これまでも研究の真似事をしてはいたが、最先端に来るのは初めてかもしれない。
あと残りの時間でどこまでできるだろうか。

技術を獲得すると色眼鏡がかかる

随分前に読んだmedtoolz氏のエントリーが、自分の中で年々大きくなっていくのを感じる。

偏見を獲得すること - レジデント初期研修用資料

"学習とは、自らの偏見に基づいた、断片的な知識の再配列に他ならない。教科書の配列をそのまま受容することは「暗記」であって、学習とは違う。偏見を持たない人は、裏を返せば学べない。

漠然と読んだ知識は、居場所がないから残らない。あらかじめ知識の居場所を作っておいて、「この知識は、この知識と関連づけて、脳内地図のこの番地に置いておこう」という態度で教科書や論文を読むと、読んだ分だけ知識が残るし、漠然と読むよりもスピードが上がる。

「偏った」人と、「勉強が得意な」人との差異は、ある問題を解決するにあたって、それぞれが獲得してきた偏見が役に立つのかどうかで決定される。偏見に優劣は存在しないし、問題が変わればもちろん、「使える人」の居場所に立つ顔ぶれも変わってくる。"

 厳密に言うと引用箇所の主題となっている「学習」とは違うのだけれども、研究の進め方、問題の捉え方という点が類似しているように感じている。最近、ラボのメンバーと関わって感じることがあったので備忘録としてまとめておきたい。


 新しく入ってきた院生達は、計算・反応・合成と、それぞれ別のバックグラウンドを持っている。彼らの課題に対するアプローチはそれぞれ個性的で、バックグラウンドに大きく依存している。皆、実験も計算もやらなければいけないテーマを行っているんだけど、例えば、計算出身ならば手を動かす前に量子化学計算で解析するし、反応・合成出身ならとりあえず反応をかけている。
 彼らは自分の慣れ親しんだ技術からアプローチする。効率的だし理にかなっていると思う。


 課題へのアプローチの仕方の違いは、物事の捉え方の違いに由来する。
 自分は化学反応を見るとき、「既知の理論で説明可能な現象か否か」「反応素過程の速度はどのくらいか、現実的に測定可能か」「合成における有用性はあるか」などの観点から理解していく。いずれも自分のバックグラウンドと習得した技術に密接に関係している。また、化学物質それ自体に対しては、数学における変数のようなものだと考えていて、それ自体に対してあまり興味関心を持っていない傾向にある。これも経験に拠って作られた価値観である。
 一方で合成を専門にしている研究者は、物質を合成ルート解析の観点から見るだろう。彼らの頭にストックされた反応の知識と物質の構造が照らし合わされ、「原理的に合成が可能か」「困難か、簡易か」「ルートが既知と類似するか、完全に新規か」などを判断基準にして物事を捉えるかもしれない。専門として向き合ってきた対象と、自分の持つ技術の解決可能な範囲が、物事の解像度に影響を与える。
 もちろん、技術や専門以外にも周囲の常識とか社会的要請、現在の環境などの要素が複雑に絡み合って物事は捉えられ、アプローチの仕方が決定されるわけだけど、その中でも研究という個人に属する活動においては技術という切り口は大きいように思っている。


 この分野での経験が少しずつ溜まってきて、使える特殊技術が増えてきた。それと同時に、上記の化学反応の見方のような、独自の「偏見」が強くなってきたように感じる。偏見は、認識した対象を自分の得意なフィールドの中で理解させる。「この技術を持っている自分なら」どう解析するか、「この分野を概観する知識を(とりあえず)持っている自分なら」この結果に価値を見いだせるか、などなど。偏見という色眼鏡をかけることによって自分独自の価値観に照らし合わせて物事を見ることが可能になる。このことで、最初に挙げたエントリー中の「知識の居場所」のように、自分が目にしたデータや読んだ論文は、成形されてから自分の中に居座るようになる。


 色眼鏡を通して見た範囲が広がっていき、様々なものが「知識の居場所」に収まっていくと、それらが自分の中で有機的に結びついていく。この、「偏見を通して見た範囲」が十分に広がり、研究対象を丸ごと自分の価値観から把握できるに至ったとき、面白く・深みを帯びて見えるし、それに対する新しいアプローチを思いつくようになるんだろうなと思っている。現ボスを含むトップクラスの先生と話させていただいたときなど、多岐にわたる事象について独自の観点からの深い見解などが伺えることがあったが、それは専門分野の色眼鏡で世界を見渡した結果かもしれない。いずれにせよ、その境地は遠い。