wvogel日記

自分用の技術備忘録が多めです.

Haskellで書く3D結晶単位胞ビューア

暑すぎて熱中症になりそうです.
私は自宅ですが,三密が満たされている場所ではもっと暑そうです.
密といえば,立方最密構造ですよね!!
三密が禁忌とされる今,Haskellの3DCGライブラリnot-glossを使って,金属単結晶とイオン結晶体のビューアを作成します.

ソースコードは下記にあります.
GitHub - wvogel00/minerals

結晶構造

自然界での事象はその多くが安定な状態に収束する傾向にあり,
身近なところでは金属結晶などの集合体もその良い例です.
自然界の金を例にとると,自由電子を介する金属結合によって最も安定な状態で結晶として産出され,結晶構造は最も充填率が高い立方最密構造(または面心立方構造)をしています.
金属は単結晶では主に,立方最密構造,六方最密構造,体心立方構造で安定します.
10剤の結晶は無数の原子が連なっていますが,その繰り返しの最小単位を単位胞といい,単位胞が決まればその物質の結晶構造は説明されたことになります.
そして,ブラべ格子によって,三次元での結晶がとりうる結晶構造の種類は14種類であることが証明されています.
つまり,最大14種類に対応すれば良いわけですが,step by stepということで,
今回は金属結晶とNaCl型イオン結晶のみ実装しました.

XMLパーサ

金属によって安定な結晶状態は変わってきますので,XML形式で保存した元素データベースを作成し,ParsecをもとにしたパーサライブラリTrifectaでパースしていきます.
try関数に頼りまくっています.空白処理もありあまり美しくありませんが,Trifectaは割と使いやすいのでもし使うときには参考にしてください.

xmlの書式は下記の通りですが,まだ全てのデータを取り込んだ訳ではないので,
タグの読み込みに失敗することも想定してパース関数を作りました.
Teは六方最密構造をとるので格子定数にcがありますが,立方最密の場合にはcフィールドはなく,aのみが記載されています.

<elements>
    <element>
        <name>Te</name>
        <structure>hcp</structure>
        <a>0.446</a>
        <c>0.593</c>
        <density>6.236</density>
        <ion>
            <e>-2</e>
            <radius>0.221</radius>
        </ion>
    </element>
</elements>

このようにフィールドがない時にはMaybeが便利なのでとりあえず包みます.
また,元素によっては複数の価数イオンをもつので,XMLを格納する代数的データ型は以下のように定義しています.

data XML = XML{
    name :: Element, 
    structure :: Maybe Structure,
    _a :: Maybe Float,
    _c :: Maybe Float,
    _density :: Maybe Float,
    _ions :: [Ion]
    }
    deriving (Eq,Show)

化学式パーサ

XMLとは別に,実行ファイルに結晶化学式を引数として渡して,その結晶構造を表示させる必要がありますので,ここも簡単にパーサを書きました.あまりにも簡単すぎて書く必要もなかったかもしれません.まあ練習と思って.とても短いのでここにも掲載します.

parseCF :: Parser [(Maybe Element, Int)]
parseCF = (:) <$> element <*> many element

element :: Parser (Maybe Element, Int)
element = (tuple <$> upper <*> try (many lower) <*> try (many digit)) where
    tuple u l n = (findElem $ u:l, safeRead n)
    safeRead "" = 1
    safeRead x = read x

findElem :: String -> Maybe Element
findElem elem = find((==elem).show) [H ..]

結晶中の元素とその個数を返し,

"NaCl" -> [(Na,1), (Cl,1)]
"FeS2" -> [(Fe,1),(S,2)]

のように動作します.

ポーリングの法則とマーデルングエネルギー

まだ実装中ですが,イオン結晶を形成するには,各原子同士が単位胞内に収まる必要があり,幾何学的観点からはポーリングの法則に,静電エネルギー観点からはマーデルングエネルギーに従います.ポーリングは陽イオン・陰イオンの原子半径や結合の優先順位について説明しており,マーデルングエネルギーはイオン間距離と,結晶構造から決まるマーデルング定数によって決定されいます.
一応ポーリングの法則の判定関数は実装したのですが,マーデルングは結晶構造と関わってくるので先延ばしにしています.

例えば,ポーリングの法則に従うと,半径r1の陽イオンから半径r2の陰イオンへの配位数は下記関数に従います

coordinateN' r1 r2
    | r1/r2 < 0.155 = 2
    | r1/r2 < 0.225 = 3
    | r1/r2 < 0.414 = 4
    | r1/r2 < 0.732 = 6 -- NaCl型
    | otherwise = 8     -- CsCl型

陽イオンから陰イオンへの配位数がわかれば,陰イオンから陽イオンへの配位数は,ポーリングの第二法則によって決めることができます.

3Dライブラリ not-gloss

以上の知識を踏まえて,早速not-glossを使ってビューアを作成していきます.
not-glossについては割と詳しく下記で説明されています.
not-glossでお手軽3Dグラフィック描画 - Qiita
ただ,角度回転を使う際にはSpatialMathライブラリが必要ですので,そこだけご注意ください.(コンパイルすればすぐわかりますが)
また,text表示は日本語には対応していません.画像にして表示するのが良いと思います.
本当はテキスト入力とかパラメータ変更をできるUIにしたいのですが,その場合はwxHaskellとか使った方がいい気がしている.
今回はdisplay関数を使ってぐるぐる回せるものを作成しただけなので扱いで特にハマったところはありません.
アニメーションもいつか挑戦してもいいかもしれない.

金属の単結晶の場合は割と簡単で,xmlに元素ごとの安定な構造が記載されていますからそれにしたがって場合分けします

monoCrystal :: XML -> [VisObject Float]
monoCrystal xml = case structure xml of
    Just HCP -> hcpMonoStructure (fromJust $ _a xml) (fromJust $ _c xml)
    Just CCP -> ccpMonoStructure (fromJust $ _a xml)
    Just FCC -> ccpMonoStructure (fromJust $ _a xml)
    Nothing -> [drawText (show (name xml) ++" have no mono-crystal")]

今回,trace関数をはじめてまともに使いました.便利ですね
場合分けができたら結晶構造ごとにフレームを作成して元素を配置してokです.
グラフ理論では少し記述が難しいように感じたので今回全部決め打ちでコードにしています.ただ,この部分も外部ファイルに追い出しても良いかも.どなたかグラフでの表現含め,アイデアあればご教示ください.

いかに,六方最密構造の例としてMg,また立方最密構造の例としてAuを掲載します.

f:id:wvogel00:20200511183416p:plain
六方最密構造 Mg
f:id:wvogel00:20200511183507p:plain
立方最密構造 Au

イオン結晶

現段階ではまだ価数やポーリングの法則を使った場合分けができていないので,(陽->陰イオン配位数,陰->陽イオン配位数)= (6,4)の
NaCl型のものしか表示できませんが,その結果も掲載します.

f:id:wvogel00:20200511183725p:plain
NaCl型イオン結晶 NaCl

なんでこんなものを作り始めたのか謎ですが,結構はまってしまいました.