29062009 music
3曲目がやばい
HARD ISLANDSが、the sky was pinkとかDarkest starでのざらついた高揚感は感じられなかったな。
Darkest Star [12 inch Analog]
pythonでdeconvolution
演習本にデコンボリューション法が載っていたのだけど、実装までは載ってなかったので、調べたら、論文が見つかった。
ファーマコキネティクス―演習による理解GOTOとIFの組み合わせがwhileなのかforなのか悩むとこが幾つかあったが、pythonで書きなおした。結構、fortranっぽさが残っている気もするが。
#!/usr/bin/env python
def calc_eps():
eps = 1.0
while eps + 1.0 > 1.0: eps = eps/2.0
return eps
def sign(a1,a2):
if a2 >= 0: return abs(a1)
else: return -abs(a1)
def deconv(a,div,b,dose):
n = len(a); m = len(b); l = n + m -1
eps = calc_eps()
u = [0.0 for i in range(l)]
v = [0.0 for i in range(l)]
gamma = [0.0 for i in range(20)]
g = [0.0 for i in range(20)]
a.sort(lambda x,y:cmp(x[1],y[1]))
for II in range(n-1):
x1 = a[II][1] * -1
x2 = a[II+1][1] * -1
qx1 = 0.0
qx2 = 0.0
for i in range(n):
s1 = 1.0
s2 = 1.0
for j in range(n):
if j != i:
s1 = s1*(a[j][1] + x1)
s2 = s2*(a[j][1] + x2)
qx1 = qx1 + a[i][0] * s1
qx2 = qx2 + a[i][0] * s2
x3 = x1
qx3 = qx1
d = x2 -x1
e = d
while(1):
if abs(qx3) < abs(qx2):
x1 = x2; x2 = x3; x3 = x1
qx1 = qx2; qx2 = qx3; qx3 = qx1
tol = 4.0*eps*abs(x2)
xm = 0.5*(x3-x2)
if abs(xm) <= tol or qx2 == 0.0: break
if abs(e) >= tol and abs(qx1) > abs(qx2):
if x1 == x3:
s = qx2/qx1; p = 2.0*xm*s; q = 1.0-s
else:
q = qx1/qx3; r = qx2/qx3; s = qx2/qx1
p = s*(2.0*xm*q*(q-r) - (x2-x1)*(r-1.0))
q = (q-1.0)*(r-1.0)*(s-1.0)
if p > 0.0: q = -q
p = abs(p)
if 2.0*p < 3.0*xm*q - abs(tol*q) and p < abs(0.5*e*q):
e = d
d = p/q
else:
d = xm
e = d
else:
d = xm
e = d
x1 = x2
qx1 =qx2
x2 = x2 + d
if abs(d) <= tol: x2 = x2 - d + sign(tol,xm)
qx2 = 0.0
for i in range(n):
s2 = 1.0
for j in range(n):
if j != i: s2 = s2*(a[j][1] +x2)
qx2 = qx2 + a[i][0] * s2
if qx2*(qx3/abs(qx3)) > 0.0:
x3 = x1
qx3 = qx1
d = x2 - x1
e = d
gamma[II] = x2
s2 = 0.0
for j in range(n):
s1 = 0.0
for k in range(n):
if k != j:
s1 = s1 + 1.0/(gamma[II] + a[k][1])
s2 = s2 + a[j][0]*s1/(gamma[II] + a[j][1])
g[II] = 1.0/s2
ak1 = 0.0
ak2 = 0.0
ak3 = 0.0
for i in range(n):
ak1 = ak1 + a[i][0];
ak2 = ak2 + a[i][0]*a[i][1]
if n > 1 and i != n-1: ak3 = ak3 + g[i]/gamma[i]
ak1 = 1.0/ak1
ak2 = ak1*ak1*ak2
ak3 = ak2 - ak3
ak4 = 100.0*div/dose
uz = 0.0
for i in range(l):
s = 0.0
if i > m-1:
v[i] = -gamma[i-m]
for j in range(m):
s = s + b[j][0]/(gamma[i-m] + b[j][1])
u[i] = ak4*g[i-m]*s/gamma[i-m]
else:
v[i] = b[i][1]
u[i] = ak4*b[i][0]*(ak1 - ak3/b[i][1])
if n != 1:
for j in range(n-1):
s = s + g[j]/(gamma[j]*(gamma[j] + b[i][1]))
u[i] = u[i] - ak4*b[i][0]*s
uz = uz - u[i]
return uz,u,v
if __name__ =='__main__':
print "=== Data set 1 ==="
a = [[1.0427,5.0526],[0.97617,0.97567]]
div = 1.0
b = [[-0.80547,3.6373],[0.80547,0.83033]]
dose = 0.6
print deconv(a,div,b,dose)
print "=== Data set 2 ==="
a = [[1.3377,3.1851],[0.52381,0.61065]]
div = 1.0
b = [[-2.2637,2.0800],[2.2637,1.2513]]
dose = 0.6
print deconv(a,div,b,dose)
print "=== Data set 3 ==="
a = [[1.0308,5.9131],[1.0487,1.0262]]
div = 1.0
b = [[-1.6103,3.3563],[1.6103,1.3917]]
dose = 0.6
print deconv(a,div,b,dose)
print "=== Data set 4 ==="
a = [[1.3377,3.1851],[0.52381,0.61065]]
div = 1.0
b = [[-5.5440,2.2488],[5.5440,1.7521]]
dose = 0.6
print deconv(a,div,b,dose)
Table.1 のテスト計算を実行。論文のデータセット4のa2はタイポでしょう。
=== Data set 1 ===
(103.38116460095375,
[99.756988719460651, -23.218031509886707, -179.92012181052769],
[3.6373000000000002, 0.83033000000000001, 2.9469592648362699])
=== Data set 2 ===
(94.012343334076377,
[212.39909284817011, 2395.320119577972, -2701.7315557602187],
[2.0800000000000001, 1.2513000000000001, 1.3350740721242431])
=== Data set 3 ===
(94.362883301636657,
[-1704.7606134762168, 73.013373899587776, 1537.3843562749923],
[3.3563000000000001, 1.3916999999999999, 3.4906828227939406])
=== Data set 4 ===
(91.15932635164495,
[370.51731307178676, -1111.1884882257646, 649.51184880233291],
[2.2488000000000001, 1.7521, 1.3350740721242431])
26062009 q-chem
The Fragment Molecular Orbital Method
こんなの出てるのか
The Fragment Molecular Orbital Method: Practical Applications to Large Mokecular SystemsCrc Pr I Llc / ￥ 12,324 ()
通常9～13日以内に発送
買わねば
あと、インプットはここにあったので、論文読みながら実行すると吉。
25062009 q-chem
「基礎量子化学—軌道概念で化学を考える」を読んだ
良書だと思う。多少量子化学の知識は必要だけど。
分子軌道法（ＭＯ法）とは，原子に対する原子軌道の考え方を分子に対して適用し，分子の電子状態や反応機構を理解しようとするものである．この量子化学的視点を用いることによって，あらゆる化学現象を統一的に捉えて説明する．化学を学ぶもの必携の一冊．
基礎量子化学―軌道概念で化学を考える複合体結晶構造とかの結合の要因を探る時には、どうしても有機化学とか量子化学の知識が必要になってくる。且つドッキングというのは化学反応の特殊な一形態であり、それゆえ、弱い力に関する理論的な理解が必須となる（であろうと考えている）。
結局ドッキングというのは、反応により生成物が変わらない化学反応と捉えることができる。(結合を形成するようなものは自殺基質とか呼んで特殊な扱いをするね)
まぁそう考えれば、量子化学的な化学反応に対する理解は重要だろうと自然に導かれるけど。
- 第７章 分子構造の支配因子
- 第８章 分子間相互作用
- 第９章 分子構造と立体化学
- 第10章 官能基の性質と酸・塩基の強度
あたりは、コンフォメーションとか掘っていく際に役に立つはず。ただ、弱い相互作用に関する話は水素結合くらいしかなかったけど、結構詳しく載っていたので役立った。
その先のもっと弱い結合をどう評価していくかがテーマなんだけど、やっぱいまのところは論文を地道に追いかけていくしかないかな。
有機化学のための分子間力入門これも良い本なんだけど実例が基礎っぽすぎて、プラクティカルさが少し足りない。
24062009 life
プライムナンバーズ
「魅惑的で楽しい素数の事典」と副題にあるように、素数に関するアレコレが簡潔に書いてある。
プライムナンバーズ -魅惑的で楽しい素数の事典読んで理解するというよりは、ふーんそうなんだ的な。
まだ、Cまでしか読んでないけど、Project Eulerの問題にあった友愛数(amicable number)とか循環素数(circular prime)とか出てきて楽しい。
Project Eulerは50問解いたとこで休憩したままだけどまた再開しようかな。
24062009 q-chem
Energy Decomposition Analysis
複合体結晶構造から重要な残基とのインタラクションを見積もるのに、リガンドと特定の残基のシンプルなモデルを切り出してきて、きちんと計算すればよかろうと。
とか書いておいて、FMOの結果から相互作用の強いアミノ酸残基とリガンドを切り出してきて、morokuma-analysisしてみたけど、ちっとも収束しなかった。
うーん、これは困ったと、切り出したモデルでオプティマイズかけて再度計算してみたけどやっぱり収束しなかった。
結晶構造の座標みたいにあまり安定でないデータからエナジー出して、その成分をELECTROSTATIC,EXCHANGE REPULSION,POLARIZATION,CHARGE TRANSFERに分解する技術はきちんと押さえておく必要があるんだけど、どうやればいいんだろうか？
やっぱPIEDAの論文読んでおく必要ある？
43 Pair interaction energy decomposition analysis, D. G. Fedorov, K. Kitaura, J. Comp. Chem. 28 (2007) 222-237.
23062009 life
「DS幼児の脳トレ」を買った
娘がDSで遊びたがるので買ってみた。
考える力をぐんぐん伸ばす! DS幼児の脳トレ本当は中古で1000円くらいになっているのを期待していたのだけど、値崩れしないのね。娘は楽しそうに遊んでいるので良しとしたい。
あと、ちょっと触ってみた感じでは、ひらがな書き取り判定は結構シビアな気がする。僕が書き取りしても、書き順があってても、気を抜いて線からはみ出て「ぶーーーー」とか言われてちょっとむかついた。
【タイトルデザイン】修悦文字:東京新宿駅や日暮里駅などでガムテープ案内人として話題の佐藤修悦氏による制作。今回、はじめてこども向けの商品タイトルを手がける。
正直ガムテフォントが幼児の脳トレにとってどういう意味があるのかわからん。
GAMESSでMOROKUMA ANALYSIS
MOROKUMA ANALYSISのサンプル(EXAM28)を。
水とアンモニアのダイマー計算
$contrl scftyp=rhf runtyp=morokuma coord=zmt $end
$basis gbasis=n31 ngauss=4 $end
$guess guess=huckel $end
$morokm iatm(1)=3 $end
$data
water-ammonia dimer
Cs
H
O 1 rOH
H 2 rOH 1 aHOH
N 2 R 1 aHOH 3 0.0
H 4 rNH 3 aHNaxis 1 180.0
H 4 rNH 3 aHNaxis 5 +120.0
H 4 rNH 3 aHNaxis 5 -120.0
rOH=0.956
aHOH=105.2
rNH=1.0124
aHNaxis=112.1451 ! makes HNH=106.67
R=2.93
$end
結果
HARTREE KCAL/MOLE
ELECTROSTATIC ENERGY ES= -0.022336 -14.02
EXCHANGE REPULSION ENERGY EX= 0.014308 8.98
POLARIZATION ENERGY PL= -0.001786 -1.12
CHARGE TRANSFER ENERGY CT= -0.003780 -2.37
HIGH ORDER COUPLING ENERGY MIX= -0.000687 -0.43
TOTAL INTERACTION ENERGY, DELTA-E= -0.014282 -8.96
DECOMPOSITION OF CT
CHARGE TRANSFER ENERGY, MON= 1 CT= -0.000494 -0.31
CHARGE TRANSFER ENERGY, MON= 2 CT= -0.003286 -2.06
DECOMPOSITION OF PL
EPL, MON= 1 PL= -0.001068 -0.67
EPL, MON= 2 PL= -0.000584 -0.37
HIGH ORDER COUPLING FOR PL, PMIX= -0.000133 -0.08
複合体結晶構造から重要な残基とのインタラクションを見積もるのに、リガンドと特定の 残基のシンプルなモデルを切り出してきて、きちんと計算すればよかろうと。
すぐできる 量子化学計算ビギナーズマニュアル (KS化学専門書)