[Python×Ansys Lumerical]Python APIでLumericalを制御してLiNbO3 光変調器の解析を実施してみよう

はじめに

本記事では、Ansys Lumerical DEVICE を Python API（lumapi）から制御し、薄膜 LiNbO₃（LN）位相変調器の電気光学特性を一括して評価するスクリプトについて解説します。

今回のスクリプトは、Ansys Optics の Application Gallery に公開されている「Thin Film Lithium Niobate Electro-Optic Phase Modulator(https://optics.ansys.com/hc/en-us/articles/19435937674387-Thin-Film-Lithium-Niobate-Electro-Optic-Phase-Modulator)」
の解析例をベースに作成したものです。元のモデル構成や解析フローは、公式ドキュメントに準拠しています。

今回のスクリプトの目的は、同一のLumericalのセッション内で CHARGE 解析と FEEM（電気光学モード解析）を連続実行し、複数波長における有効屈折率および、中心波長 1.55 µm における VπL、光損失、群屈折率などをpython APIを用いて一括して取得することです。
特に、同じ解析を波長ごとに何度も起動する無駄を避け、「一度の実行ですべての物理量を取得する」ことを意識してスクリプトを組んでみました。

スクリプト全体(LN_modulator.py)

こちらがスクリプトになります。
application galleryにあるLumericalのプロジェクトファイルとlsfファイルをpython側で呼び出し、
位相シフタの特性値の取得を行ってます。
群屈折率は複数(1.545 umと1.555 um)の波長の屈折率の差分をとることで、中心波長(1.55 um)の屈折率を求めています。

冒頭では、lumapi を利用するために Lumerical の Python API へのパスを追加し、数値処理・可視化・CSV 出力のために NumPy、Pandas、Matplotlib を読み込んでいます。その後、解析対象となる DEVICE プロジェクトファイル（.ldev）と、CHARGE および FEEM を実行するための LSF スクリプトのパスを定義しています。これにより、Python 側から GUI 操作と同等の解析を再現できるようになっています。

# -*- coding: utf-8 -*-
"""
Created on Sun Feb  2 21:45:30 2025

@author: babatake
"""

import sys

sys.path.append(r"C:\Program Files\Lumerical\v252\api\python")
import lumapi
import os
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd

# 📂 パス設定
project_dir   = r"D:/LN modulator"
project_file  = os.path.join(project_dir, "LN_phase_modulator.ldev").replace("\\", "/")
charge_script = os.path.join(project_dir, "LN_phase_modulator_CHARGE.lsf").replace("\\", "/")
feem_script   = os.path.join(project_dir, "LN_phase_modulator_FEEM.lsf").replace("\\", "/")

# ---------- 1D float 正規化ユーティリティ ----------
def _to_1d_float(name, arr):
    """任意の配列を 1D float64 に正規化（(N,1) 等や複素も吸収）"""
    if arr is None:
        raise ValueError(f"{name}: None")
    a = np.asarray(arr)
    a = np.real_if_close(a, tol=1e5)
    a = np.ravel(a)
    try:
        a = a.astype(np.float64, copy=False)
    except Exception as e:
        raise ValueError(f"{name}: cannot cast to float64; shape={a.shape}, dtype={a.dtype}") from e
    return a

def _safe_getv(device, key, dtype=float):
    """device.getv(key) の安全取得（失敗時 None）。"""
    try:
        val = device.getv(key)
        if val is None:
            return None
        return np.array(val, dtype=dtype)
    except Exception:
        return None

# ---------- 統合：3波長＋1.55µmの L_pi・loss・VπL を1回で取得 ----------
def run_all_once(params, wavelengths=(1.545e-6, 1.55e-6, 1.555e-6)):
    """
    同一 DEVICE セッションで、指定波長群を順に CHARGE→FEEM 実行。
    - 各波長の neff_TE を取得
    - 中心波長(1.55 µm)では L_pi, alpha_dB, Volt, VpiL_step も取得
    戻り値:
      neff_dict: {λ: neff_TE(np.ndarray)}
      ng_TE: 1.55 µm の群屈折率（Volt配列と同長）または None
      volt_1550: 電圧配列 (np.ndarray) または None
      Lpi_1550: 1.55 µm の L_pi [m] (np.ndarray) または None
      loss_1550: 1.55 µm の loss [dB/cm] (np.ndarray) ※alpha_dB（無ければ推定 or None）
      VpiL_step_1550: LSF中心差分の VπL [V·cm] または None
      VpiL_prod_1550: 積定義 V×Lπ の VπL [V·cm] または None
    """
    lam_center = 1.55e-6
    neff_dict = {}
    volt_1550 = None
    Lpi_1550  = None
    loss_1550 = None
    VpiL_step_1550 = None

    with lumapi.DEVICE() as device:
        device.load(project_file)
        device.switchtolayout()

        for wl in wavelengths:
            # layout モードで設定してから実行
            device.switchtolayout()
            device.setnamed("FEEM", "wavelength", wl)

            # CHARGE → FEEM
            with open(charge_script, "r", encoding="utf-8") as f:
                device.eval(f.read())
            with open(feem_script,  "r", encoding="utf-8") as f:
                device.eval(f.read())

            # neff 取得
            neff_TE = _safe_getv(device, "neff_TE", dtype=np.complex128)
            if neff_TE is None:
                raise RuntimeError(f"neff_TE not produced at wavelength={wl}")
            neff_dict[wl] = neff_TE

            # 中心波長なら各種も取得
            if abs(wl - lam_center) < 1e-12:
                Lpi_1550  = _safe_getv(device, "L_pi", dtype=float)
                loss_1550 = _safe_getv(device, "alpha_dB", dtype=float)  # dB/cm（無いこともある）
                # LSF 側に Volt / VpiL_step を保存してある前提（無ければフォールバック）
                volt_1550     = _safe_getv(device, "Volt", dtype=float)
                VpiL_step_1550 = _safe_getv(device, "VpiL_step", dtype=float)

                # Volt が取れない場合は等間隔で補う（neff 長に合わせる）
                if volt_1550 is None:
                    volt_1550 = np.linspace(0.0, 1.0, len(neff_TE))

    # 群屈折率 ng（両端差分：ng = neff - λ dneff/dλ）
    if (1.545e-6 in neff_dict) and (lam_center in neff_dict) and (1.555e-6 in neff_dict):
        neff_1545 = np.asarray(neff_dict[1.545e-6], dtype=np.complex128)
        neff_1550 = np.asarray(neff_dict[lam_center], dtype=np.complex128)
        neff_1555 = np.asarray(neff_dict[1.555e-6], dtype=np.complex128)
        dlam = 1.555e-6 - 1.545e-6
        dneff_dlam = (neff_1555 - neff_1545) / dlam
        ng_TE = np.real(neff_1550 - lam_center * dneff_dlam)
    else:
        print("⚠️ ng を計算するための 3 波長結果が不足")
        ng_TE = None

    # loss フォールバック（alpha_dB が無い場合、Im(neff) から推定）
    if (loss_1550 is None or loss_1550.size == 0) and (lam_center in neff_dict):
        k0  = 2 * np.pi / lam_center
        n_im = np.imag(np.asarray(neff_dict[lam_center], dtype=np.complex128))
        if n_im.size > 0:
            # power-loss[dB/m] = 8.686 * 2 * k0 * Im(n)
            loss_dB_per_m  = 8.686 * 2.0 * k0 * n_im
            loss_1550 = loss_dB_per_m * 0.01  # dB/cm
        else:
            loss_1550 = None

    # VπL（積の定義: V×Lπ）[V·cm]
    VpiL_prod_1550 = None
    if (volt_1550 is not None) and (Lpi_1550 is not None):
        VpiL_prod_1550 = np.asarray(volt_1550, dtype=float) * (np.asarray(Lpi_1550, dtype=float) * 100.0)

    return neff_dict, ng_TE, volt_1550, Lpi_1550, loss_1550, VpiL_step_1550, VpiL_prod_1550

# ==============================
#  🚀 メイン（スタンドアロン実行）
# ==============================
if __name__ == "__main__":
    # 幾何パラメータ（必要なら使用：今回は LSF 側の変数をそのまま利用）
    #params = {
    #    "wg_width_top": 0.4e-6,
    #    "wg_height": 0.5e-6,
    #    "Ground_Electrode_Left_xmin": -15e-6,
    #    "Ground_Electrode_Left_xmax": -11e-6,
    #    "Signal_Electrode_xmin": -2e-6,
    #    "Signal_Electrode_xmax": 2e-6,
    #    "Ground_Electrode_Right_xmin": 11e-6,
    #    "Ground_Electrode_Right_xmax": 15e-6
    #}

    # ★ これ1回だけで全取得（重複実行なし）
    (neff_dict, ng_TE, Volt, Lpi_m,
     loss_dB_per_cm, VpiL_step, VpiL_prod) = run_all_once(params)

    # --- 参考表示（neff と ng）---
    if (1.55e-6 in neff_dict) and (ng_TE is not None):
        print("\n✅ 解析結果（1.55 µm）:")
        print(f"neff_TE at 1.55 µm: {neff_dict[1.55e-6]}")
        print(f"Group Index ng_TE: {ng_TE}")
    else:
        print("⚠️ 1.55 µm の neff/ng が取得できず")

    # --- VπL & loss → CSV/プロット ---
    if (Volt is not None) and (Lpi_m is not None):
        # 1D 正規化
        Volt1 = _to_1d_float("Volt", Volt)
        Lpi1  = _to_1d_float("L_pi", Lpi_m)
        loss1 = None if loss_dB_per_cm is None else _to_1d_float("loss", loss_dB_per_cm)
        VpiL_step1 = None if VpiL_step is None else _to_1d_float("VpiL_step", VpiL_step)
        VpiL_prod1 = None if VpiL_prod is None else _to_1d_float("VpiL_prod", VpiL_prod)

        # 長さ合わせ（存在する配列の最小長）
        seqs = [Volt1, Lpi1]
        if loss1 is not None:      seqs.append(loss1)
        if VpiL_step1 is not None: seqs.append(VpiL_step1)
        if VpiL_prod1 is not None: seqs.append(VpiL_prod1)
        n = min(map(len, seqs))
        Volt1, Lpi1 = Volt1[:n], Lpi1[:n]
        if loss1 is not None:      loss1 = loss1[:n]
        if VpiL_step1 is not None: VpiL_step1 = VpiL_step1[:n]
        if VpiL_prod1 is not None: VpiL_prod1 = VpiL_prod1[:n]

        # 0V と非有限を除外（0V の L_pi / VπL は無くても良い前提）
        mask = (Volt1 > 0) & np.isfinite(Lpi1)
        if loss1 is not None:      mask &= np.isfinite(loss1)
        if VpiL_step1 is not None: mask &= np.isfinite(VpiL_step1)
        if VpiL_prod1 is not None: mask &= np.isfinite(VpiL_prod1)

        Volt1, Lpi1 = Volt1[mask], Lpi1[mask]
        if loss1 is not None:      loss1 = loss1[mask]
        if VpiL_step1 is not None: VpiL_step1 = VpiL_step1[mask]
        if VpiL_prod1 is not None: VpiL_prod1 = VpiL_prod1[mask]

        # CSV 保存（両定義の VπL を併記）
        out_csv = os.path.join(project_dir, "lumerical_loss_vpil.csv")
        data = {"Volt[V]": Volt1, "L_pi[m]": Lpi1}
        if VpiL_step1 is not None: data["VpiL_step[V*cm]"] = VpiL_step1
        if VpiL_prod1 is not None: data["VpiL_prod[V*cm]"] = VpiL_prod1
        if loss1 is not None:      data["loss[dB/cm]"] = loss1
        pd.DataFrame(data).to_csv(out_csv, index=False, float_format="%.12e")
        print(f"✅ Saved CSV: {out_csv}")
        print("[shapes]", {k: np.asarray(v).shape for k, v in data.items()})

        # プロット：VπL
        if (VpiL_step1 is not None) or (VpiL_prod1 is not None):
            plt.figure(figsize=(8,5))
            if VpiL_step1 is not None:
                plt.plot(Volt1, VpiL_step1, marker="o", label="VπL (center-diff step)")
            #if VpiL_prod1 is not None:
            #    plt.plot(Volt1, VpiL_prod1, marker="s", label="VπL (V × Lπ)")
            plt.xlabel("Voltage (V)"); plt.ylabel("VπL (V·cm)")
            plt.title("VπL vs Voltage @ 1.55 μm")
            plt.grid(True, which="both", linestyle="--"); plt.legend(); plt.show()

        # プロット：loss（あれば）
        if loss1 is not None:
            plt.figure(figsize=(8,5))
            plt.plot(Volt1, loss1, marker="o")
            plt.xlabel("Voltage (V)"); plt.ylabel("Loss (dB/cm)")
            plt.title("Optical Loss vs Voltage @ 1.55 μm")
            plt.grid(True, which="both", linestyle="--"); plt.show()
        else:
            print("⚠️ loss は取得できなかった（alpha_dB 未出力 or Im(neff) 未取得）")
    else:
        print("⚠️ VπL/loss の出力条件を満たしていない（1.55 µm データ不足）")

最後に、Matplotlib を用いて VπL–電圧特性および光損失–電圧特性をプロットしています。特に VπL は、電気光学変調器の性能指標として最も重要な量の一つであり、電圧依存性を可視化することで設計パラメータの妥当性を確認できます。

出力結果

こちらが出力結果になります。
求めた位相屈折率と群屈折率はpythonプロンプト上で表示させています。

こちらがVπLと光損失の計算結果の電圧依存性になります。

なおスクリプトの実行終了時にはVπLと光損失ともにCSVで出力させるようにしています。

※本記事は筆者個人の見解であり、所属組織の公式見解を示すものではありません。

問い合わせ

光学シミュレーションソフトの導入や技術相談、
設計解析委託をお考えの方はサイバネットシステムにお問合せください。

光学ソリューションサイトについては以下の公式サイトを参照：
👉 光学ソリューションサイト(サイバネット)

光学分野のエンジニアリングサービスについては以下の公式サイトを参照：
👉 光学エンジニアリングサービス(サイバネット)