pythonでmdtrajからaseへ変換する方法

pythonで分子動力学法のトラジェクトリを扱うさい，よく知られたパッケージとしてmdtrajとaseがある．今回はmdtrajで得た構造をaseに変換する必要があったので，やり方について自分用にメモしておく．

大きな問題はmdtrajの原子は例えば同じ水素でもhやh1,h2というふうに結合によって名前が違うので，これを単にaseに入れてもエラーになってしまうということ．これを防ぐためには，以下のようにtopologyを取得する必要がある．

次に注意する点として単位の違いが挙げられる．mdtrajは分子動力学でよく使うスケールであるnm単位を採用しているが，aseはあくまで原子系の計算をメインターゲットにしておりAngstrom単位を採用している．そこで以下のようにmdtrajから変換する時には10をかけて補正する必要がある．

というわけでできたのがこちら．入力にmdtrajの1ステップを入れると，対応するaseを出力してくれる．

def convert_mdtraj_to_ase(mdtraj_snap):
    # idの取得
    table, bonds =mdtraj_snap.topology.to_dataframe()
    #
    from ase import Atoms 
    mols = Atoms(symbols=table['element'],
                 positions=snap.xyz.reshape([-1,3])*10,  # nm(gro) to ang (ase)
                 cell= snap.unitcell_vectors.reshape([3,3])*10,
                 pbc=[1, 1, 1]) #pbcは周期境界条件．これをxyz軸にかけることを表している．
    return mols   # ase object