ASEでQEのデータを作ってCar-ParrinelloMDをやってASEで可視化する

PythonのASE (Atomic Simulation Environment) でQuantum ESPRESSOのデータファイルを作って、cp.xによる第一原理MD計算をやってみます。

入力データ作成と計算実行

以下の流れで進みます。

1. ASEで入力構造を作成し、第一段階と第二段階の入力ファイルを作成。
2. 第一段階の計算（原子を動かさずに電子状態のみを計算）
3. 第二段階の計算（原子を動かすMD計算）

詳しい原理や計算キーワードについては他の詳しい記事などを参考にしてください。本記事では、とりあえず動けばいいというレベルで解説します。

ASEで入力構造を作成する部分については以下の記事がベースになっています。

ASEでQEのデータを作って第一原理MDをやる

PythonのASE (Atomic Simulation Environment) でQuantum ESPRESSOのデータファイルを作って、pw.xによる第一原理MD計算をやってみます。 cp.xのほうの計算はこちら モデル作成 メタ...

tech.gijukatsu.com

2023.12.29

そして第一段階の入力ファイルは以下です。

&CONTROL
   calculation='cp'
   restart_mode='from_scratch'
   dt=5.0
   nstep=50
   pseudo_dir='./'
/
&SYSTEM
   ntyp             = 2
   nat              = 5
   ibrav            = 0
   nosym=.true.
   nr1b=16
   nr2b=16
   nr3b=16
   ecutwfc=20.0
/
&ELECTRONS
   electron_dynamics='damp'
   electron_damping=0.2
/
&IONS
   ion_dynamics='none'
/
&CELL
/

ATOMIC_SPECIES
C 12.011 C.pbe-n-rrkjus_psl.1.0.0.UPF
H 1.008 H.pbe-rrkjus_psl.1.0.0.UPF

CELL_PARAMETERS angstrom
10.00000000000000 0.00000000000000 0.00000000000000
0.00000000000000 10.00000000000000 0.00000000000000
0.00000000000000 0.00000000000000 10.00000000000000

ATOMIC_POSITIONS angstrom
C 5.0000000000 5.0000000000 5.0000000000 
H 5.6291180000 5.6291180000 5.6291180000 
H 4.3708820000 4.3708820000 5.6291180000 
H 5.6291180000 4.3708820000 4.3708820000 
H 4.3708820000 5.6291180000 4.3708820000 

そして第二段階のファイルでは以下の部分を書き換えます。

• &CONTROL
  • restart_modeを’from_scratch’から’reset_counters’にする。（第一段階の計算ステップ数カウンターをリセットして、第一段階の計算の続きから計算する）
  • nstepにMD計算ステップ数を設定する。
• &ELECTRONS
  • electron_dynamicsを’verlet’にする。（ほかの設定もあるが詳しくは知らない）
• &IONS
  • ion_dynamicsを’none’から’verlet’にする。（第一段階では原子を動かさなかったが第二段階ではVerlet法で動かす）
  • ion_temperatureを’nose’にする。（温度制御を能勢の方法で行う）
  • tempwに任意の温度を設定する。（デフォルトでは300）

それ以外は第一段階と同じだと思いますが、エラーが出る場合やより計算内容を厳密に設定する場合があると思いますので、公式のドキュメントやexampleのファイルなどを参考にして書き換えてください。

第一段階の続きから計算するわけですが、計算フォルダ内にcp_50.saveのようや名前のフォルダができると思います。続きから計算するためのデータはここに入っています。

またcp.posなどのファイルは、第二段階の結果が第一段階の結果に追記される形になります。第一段階の結果が不要なのであれば、第二段階の計算前にcp.*を削除するといいと思います。

あとはこれらのファイルをcp.x -in input.in のような感じで実行します。

出力データの可視化

ASEで可視化をしたいと思いますが、ASEはQEのcp.xの出力データには対応していません。
ですがASEではQEのinputファイルには対応しているので、初期構造をinputファイルから読み込み、それ以降のトラジェクトリをcp.xのposファイルとcelファイルから作るコードを書きました。

from ase import Atoms
from ase.io import read, write
from ase.io.trajectory import Trajectory
from ase.units import Bohr

#input
qe_input  = "input_2.in"
cel_file  = "cp.cel"
pos_file  = "cp.pos"
#output
traj_file = "ase.traj"

input_atoms = read(qe_input, format="espresso-in")
symbols = input_atoms.get_chemical_symbols()
Natoms = len(symbols)

def read_file(filename, num_data):
    f = open(filename, "r")
    lines = f.readlines()
    f.close()
    data1    = []
    data2    = []
    data2tmp = []
    data2_count = 0
    for line in lines:
        line = line.replace("\n", "")
        values = line.split()
        if len(values) == 2:
            data1.append(values)
        if len(values) ==3:
            data2tmp.append(values)
            data2_count += 1
            if data2_count == num_data:
                data2.append(data2tmp)
                data2_count = 0
                data2tmp = []
    return data1, data2

step_time, cell_H = read_file(cel_file, 3)
cel = []
for H in cell_H:
    a = float(H[0][0]) * Bohr
    b = float(H[1][1]) * Bohr
    c = float(H[2][2]) * Bohr
    cel.append([a, b, c])

step_time, positions = read_file(pos_file, Natoms)

traj = Trajectory(traj_file, "w")
traj.write(atoms=input_atoms)

for i, pos in enumerate(positions):
    xyz = []
    for j in range(Natoms):
        x, y, z = pos[j]
        x = float(x) * Bohr
        y = float(y) * Bohr
        z = float(z) * Bohr
        xyz.append([x, y, z])
    snapshot = Atoms(symbols, xyz)
    snapshot.set_cell(cel[i])
    snapshot.set_pbc(True)
    traj.write(atoms=snapshot)

traj.close()

難しいことはやっていません。posやcelファイルはどちらも「ステップ数」と「データの中身」を繰り返し記録している構造なので、それを読み込む関数をdefしています。
posファイルとcelファイルには初期構造が記録されていないので、最初のトラジェクトリは入力ファイルのものを使用しています。
あとは最後のforループでASEのトラジェクトリを作って、ファイルをcloseします。
注意点として、長さの単位がposファイルとcelファイルではボーアになっているので、それをÅに直す定数を掛け算しています。

これを実行して最終的に出力されるase.trajファイルを以下のコマンドで読むことができます。

ase gui ase.traj

結果はこんな感じでした。

それっぽい動きをしていますね。またASEでの描画も問題なさそうです。