物質科学シミュレーションのポータルサイト

次に擬ポテンシャルをダウンロードしたのち、pw.xを実行する。

$ wget https://www.physics.rutgers.edu/gbrv/ni_pbesol_v1.4.uspp.F.UPF
$ wget https://pseudopotentials.quantum-espresso.org/upf_files/O.pbesol-n-kjpaw_psl.0.1.UPF
$ pw.x < NiO.in > NiO.out

# plot_spin.py
import matplotlib
matplotlib.use('Agg')

import re
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# --- 設定項目 ---
# Quantum ESPRESSOの出力ファイルへのパスを指定してください
QE_OUTPUT_FILE = 'NiO.out'
OUTPUT_IMAGE_FILE = 'spin_structure.png' # 保存する画像ファイル名

# 描画設定
SPIN_SCALE_FACTOR = 1.8  # スピンの矢印の長さを調整する係数
ATOM_SIZE = 200          # 原子を表す点の大きさ
SPIN_THRESHOLD = 0.1     # この値(Bohr mag)以下の磁化を持つ原子のスピンは表示しない

# 元素ごとの色を定義 (計算に合わせて追加・変更してください)
ELEMENT_COLORS = {
    'Ni': 'darkcyan',
    'O': 'red',
}

# --- Quantum ESPRESSO出力解析 ---

def parse_qe_output(filepath):
    """
    Quantum ESPRESSO (pw.x) のSCF出力ファイルをパースする。
    最後の構造と磁化の情報を抽出する。
    """
    BOHR_TO_ANGSTROM = 0.529177210903
    alat, lattice_vectors_alat, atoms, magnetization = None, [], [], []

    try:
        with open(filepath, 'r') as f:
            lines = f.readlines()
    except FileNotFoundError:
        raise FileNotFoundError(f"エラー: ファイル '{filepath}' が見つかりません。")

    # 構造情報を抽出
    alat_line = next((line for line in lines if 'lattice parameter (alat)' in line), None)
    if alat_line:
        alat_bohr = float(alat_line.split('=')[1].split()[0])
        alat = alat_bohr * BOHR_TO_ANGSTROM

    axes_start_index = next((i for i, line in enumerate(lines) if 'crystal axes: (cart. coord. in units of alat)' in line), None)
    if axes_start_index is not None:
        for i in range(1, 4):
            line = lines[axes_start_index + i]
            vec_part = line.split('=')[1]
            vec_str = vec_part.replace('(', '').replace(')', '').strip()
            vec = [float(x) for x in vec_str.split()]
            lattice_vectors_alat.append(vec)

    pos_start_index = next((i for i, line in enumerate(lines) if 'site n.     atom                  positions (alat units)' in line), None)
    if pos_start_index is not None:
        for i in range(pos_start_index + 1, len(lines)):
            line = lines[i]
            if 'tau(' in line:
                parts = line.split()
                element = ''.join(filter(str.isalpha, parts[1]))
                coord_part = line.split('=')[1]
                coords_str = coord_part.replace('(', '').replace(')', '').strip()
                coords_alat = np.array([float(c) for c in coords_str.split()])
                atoms.append({'element': element, 'coords_alat': coords_alat})
            elif not line.strip():
                break

    mag_start_index = -1
    for i in range(len(lines) - 1, -1, -1): # 後ろからインデックスを探索
        if 'Magnetic moment per site' in lines[i]:
            mag_start_index = i
            break # 最初に見つかった(＝ファイルの一番後ろの)セクションで確定
            
    # 見つかった場合、そのセクションを読み込む
    if mag_start_index != -1:
        for i in range(mag_start_index + 1, len(lines)):
            line = lines[i]
            if 'atom' in line and 'magn=' in line:
                magn_z = float(line.split('magn=')[1].split()[0])
                magnetization.append(magn_z)
            # 磁化セクションの終わりを判定 (空行または別のセクションの開始)
            elif not line.strip() or 'total cpu time' in line:
                break # 磁化セクション後の空行などでループを抜ける

    # --- データ検証と整形 ---
    if not all([alat, lattice_vectors_alat, atoms, magnetization]):
        # デバッグ用にどの情報が欠けているかを表示
        print("--- デバッグ情報 ---")
        print(f"格子定数 (alat) 抽出成功: {'はい' if alat else 'いいえ'}")
        print(f"格子ベクトル (lattice_vectors_alat) 抽出成功: {'はい' if lattice_vectors_alat else 'いいえ'}")
        print(f"原子座標 (atoms) 抽出成功: {'はい' if atoms else 'いいえ'}")
        print(f"磁化 (magnetization) 抽出成功: {'はい' if magnetization else 'いいえ'}")
        print("--------------------")
        raise ValueError("エラー: 出力ファイルから必要な情報（格子、原子、磁化）をすべて抽出できませんでした。")

    if len(atoms) != len(magnetization):
        raise ValueError(f"エラー: 原子数({len(atoms)})と磁化の数({len(magnetization)})が一致しません。")

    lattice_vectors = np.array(lattice_vectors_alat) * alat
    final_atoms = [{'element': atom['element'], 'coords': atom['coords_alat'] * alat} for atom in atoms]
    spin_vectors = [[0.0, 0.0, mz] for mz in magnetization]

    return {
        'lattice_vectors': lattice_vectors,
        'atoms': final_atoms,
        'spin_vectors': spin_vectors,
    }

# --- メインの描画処理 (変更なし) ---
def plot_spin_structure_from_qe(qe_outfile):
    try:
        data = parse_qe_output(qe_outfile)
    except (ValueError, FileNotFoundError) as e:
        print(e)
        return

    fig = plt.figure(figsize=(10, 10))
    ax = fig.add_subplot(111, projection='3d')

    for i, atom in enumerate(data['atoms']):
        x, y, z = atom['coords']
        element = atom['element']
        color = ELEMENT_COLORS.get(element, 'gray')
        ax.scatter([x], [y], [z], s=ATOM_SIZE, c=color, label=element, depthshade=False)
        
        mx, my, mz = data['spin_vectors'][i]
        if np.linalg.norm([mx, my, mz]) > SPIN_THRESHOLD:
            ax.quiver(x, y, z, mx, my, mz,
                      length=SPIN_SCALE_FACTOR,
                      normalize=True,
                      color='black',
                      arrow_length_ratio=0.3)

    vectors = data['lattice_vectors']
    origin = np.array([0,0,0])
    corners = [
        origin, origin + vectors[0], origin + vectors[1], origin + vectors[2],
        origin + vectors[0] + vectors[1], origin + vectors[0] + vectors[2],
        origin + vectors[1] + vectors[2], origin + vectors[0] + vectors[1] + vectors[2]
    ]
    edges = [
        (0, 1), (0, 2), (0, 3), (1, 4), (1, 5), (2, 4),
        (2, 6), (3, 5), (3, 6), (4, 7), (5, 7), (6, 7)
    ]
    for start_idx, end_idx in edges:
        p1, p2 = corners[start_idx], corners[end_idx]
        ax.plot3D(*zip(p1, p2), color="gray", linestyle='--')
    
    ax.set_xlabel('X (Å)')
    ax.set_ylabel('Y (Å)')
    ax.set_zlabel('Z (Å)')
    ax.set_title('Crystal Structure with Spin Vectors (from Quantum ESPRESSO)')
    
    handles, labels = ax.get_legend_handles_labels()
    by_label = dict(zip(labels, handles))
    ax.legend(by_label.values(), by_label.keys())

    x_lims, y_lims, z_lims = ax.get_xlim(), ax.get_ylim(), ax.get_zlim()
    ax.set_box_aspect([lim[1] - lim[0] for lim in [x_lims, y_lims, z_lims]])
    ax.set_anchor('C')

    plt.savefig(OUTPUT_IMAGE_FILE, dpi=300, bbox_inches='tight')
    print(f"画像 '{OUTPUT_IMAGE_FILE}' を保存しました。")
    plt.show()


if __name__ == '__main__':
    plot_spin_structure_from_qe(QE_OUTPUT_FILE)

pythoh3 plot_spin.py