« 少し量子化学に近い話:maximaとgeomviewで遊ぶ(2) | トップページ | 再びmaxima »

2009年7月22日 (水)

少し量子化学に近い話:さらにHaskell

さて、昨日まで書いたような電子密度のグラフィックを描くためには、水素電子の波動関数が必要です。いろいろな参考書に出ているとは思うのですが、主量子数の高いところまではあまり出ていません。

とりあえず動径方向の関数の係数を計算したいと思って、Haskellでプログラムを書きました(マイナー過ぎますか)。

こんな感じです。

import System

main = do args <- getArgs
print $ reverse $ gcd1 $ ratio (read (head args)) $ reduce (2 * (read (head (tail args))) + 1) $ reduceExp $ orb ((read (head args)) + (read (head (tail args))))

orb :: Int -> [Integer]
orb n = (replicate n 0) ++ [1]

reduceExp :: [Integer] -> [Integer]
reduceExp xs = reduceExp1 ((length xs) - 1) xs

reduceExp1 :: Int -> [Integer] -> [Integer]
reduceExp1 n xs = if (n > 0) then (reduceExp1 (n-1) (reduceExp2 xs)) else (xs)

reduceExp2 :: [Integer] -> [Integer]
reduceExp2 xs = zipWith (-) (zipWith (*) [1..] ((tail xs) ++ [0])) xs

reduce :: Integer -> [Integer] -> [Integer]
reduce n xs = if (n > 0) then (reduce (n-1) (reduce1 xs)) else (xs)

reduce1 :: [Integer] -> [Integer]
reduce1 xs = zipWith (*) [1..] (tail xs)

ratio :: Integer -> [Integer] -> [Integer]
ratio n xs = zipWith (*) (ratio1 (length xs) n 2) xs

ratio1 :: Int -> Integer -> Integer -> [Integer]
ratio1 l x y = if (l == 1) then [1] else [x * (head (ratio1 (l-1) x y))] ++ (map (\n -> n * y) (ratio1 (l-1) x y))

gcd1 :: [Integer] -> [Integer]
gcd1 xs= gcd2 xs (gcd3 xs)

gcd2 :: [Integer] -> Integer -> [Integer]
gcd2 xs m = map (\n -> div n m) xs

gcd3 :: [Integer] -> Integer
gcd3 (x:xs) = if (xs == []) then x else (gcd x (gcd3 xs))

これをorb.hs(Orbitの略のつもり)として保存します(しばらく前に書いたものなので、今読みなおしても、何やってたんだかあまりよく理解できません)。

もちろん、先にHaskellシステムのインストールが必要でした。Glasgow Haskell Compilation Systemを入れることにして、

#apt-get install ghc6

としてインストールします(もちろんDebianとかの場合)。

コンパイルは以下の通り。

$ ghc orb.hs -o orb

でできたorbプログラムに、主量子数と軌道量子数を引数として与えて計算させます。

$ ./orb 3 0
[-2,18,-27]

こんな感じです。これは3Sの動径方向の関数の係数に対応します。

後、全体の係数とプロット範囲とかも求めなければならないのですが、それはまた後で。

こうやって作った6S軌道を載せておきます。

#6S

plot3d(
[p, q,
(exp(-sqrt(p*p+q*q)/6)*(((sqrt(p*p+q*q) -90) * (p*p+q*q) + 2700 * sqrt(p*p+q*q) - 32400) * (p*p+q*q) + 145800 * sqrt(p*p+q*q) - 174960))^2*(p*p+q*q)*0.000000000744904*0.5],
[p,-83.4725,83.4725],[q,-83.4725,83.4725],['grid,160,160]);

さっきのorbの出力はこうでした。

$ ./orb 6 0
[1,-90,2700,-32400,145800,-174960]

では今日はここまで。

|

« 少し量子化学に近い話:maximaとgeomviewで遊ぶ(2) | トップページ | 再びmaxima »

コメント

この記事へのコメントは終了しました。

トラックバック


この記事へのトラックバック一覧です: 少し量子化学に近い話:さらにHaskell:

« 少し量子化学に近い話:maximaとgeomviewで遊ぶ(2) | トップページ | 再びmaxima »