3. 研究手法
3.1 分子モデル
有機半導体分子は7ビーズの粗視化モデルとして構築した。ビーズ配列はT-T-C-C-C-T-Tであり、中央の3つのCビーズ（Core、タイプ1）は共役骨格のπ電子系を、両端の計4つのTビーズ（Tail、タイプ2）は側鎖アルキル基を表現する。
表1: 粗視化ビーズの特性
| ビーズタイプ | 物理的意味 | 個数/分子 | 質量 (m*) |
|---|---|---|---|
| C (Core, タイプ1) | 共役骨格 (π電子系) | 3 | 1.0 |
| T (Tail, タイプ2) | 側鎖 (アルキル基) | 4 | 1.0 |
3.2 力場パラメータ
全ての相互作用パラメータはLennard-Jones（LJ）還元単位系で定義した。
3.2.1 非結合相互作用
ビーズ間の非結合相互作用には12-6 LJ ポテンシャルを使用した：
$$U_{LJ}(r) = 4\varepsilon \left[\left(\frac{\sigma}{r}\right)^{12} - \left(\frac{\sigma}{r}\right)^{6}\right], \quad r < r_c$$
カットオフ距離 $r_c = 2.5\sigma$ とした。
表2: LJペアポテンシャルパラメータ
| ペア | $\varepsilon$ | $\sigma$ | 物理的意味 |
|---|---|---|---|
| C-C | 2.0 | 1.0 | 強いπ-πスタッキング |
| C-T | 1.2 | 1.0 | 骨格-側鎖間の中間的相互作用 |
| T-T | 1.0 | 1.0 | 弱い側鎖間ファンデルワールス力 |
C-Cペアのε値を2.0と大きく設定することで、有機半導体に特徴的なπ-πスタッキング相互作用を再現した。
3.2.2 結合相互作用
分子内結合にはFENE（Finitely Extensible Nonlinear Elastic）ポテンシャルを採用した：
$$U_{FENE}(r) = -\frac{1}{2}KR_0^2 \ln\left[1 - \left(\frac{r}{R_0}\right)^2\right] + 4\varepsilon\left[\left(\frac{\sigma}{r}\right)^{12} - \left(\frac{\sigma}{r}\right)^{6}\right] + \varepsilon$$
パラメータ: K = 30.0, $R_0$ = 1.5, $\varepsilon$ = 1.0, $\sigma$ = 1.0
3.2.3 角度相互作用
隣接3ビーズ間の角度ポテンシャルには調和型を使用し、分子の異なる部位に応じて剛性を変化させた：
$$U_{angle}(\theta) = K(\theta - \theta_0)^2$$
表3: 角度ポテンシャルパラメータ
| 角度タイプ | K (剛性) | $\theta_0$ (度) | 物理的意味 |
|---|---|---|---|
| T-T-C | 15.0 | 180 | 柔軟な側鎖結合 |
| T-C-C | 25.0 | 180 | 側鎖-骨格接合部 |
| C-C-C | 50.0 | 180 | 剛直な共役骨格 |
C-C-C角度の高いK値（50.0）により、共役骨格の面内剛性を表現した。
3.2.4 二面角相互作用
OPLS（Optimized Potentials for Liquid Simulations）型二面角ポテンシャルを適用した：
$$U_{dihedral}(\phi) = \frac{1}{2}\left[K_1(1+\cos\phi) + K_2(1-\cos 2\phi) + K_3(1+\cos 3\phi) + K_4(1-\cos 4\phi)\right]$$
表4: 二面角ポテンシャルパラメータ
| 二面角タイプ | $K_1$ | $K_2$ | $K_3$ | $K_4$ |
|---|---|---|---|---|
| T-T-C-C | 0.0 | 0.5 | 0.2 | 0.0 |
| T-C-C-C | 0.0 | 2.0 | 0.0 | 0.0 |
3.3 蒸着シミュレーション（PVDモデル）
3.3.1 シミュレーションボックスと境界条件
蒸着シミュレーションには、x方向・y方向に周期境界、z方向に固定境界（ff: fixed-fixed）を持つ直方体ボックス（$12 \times 12 \times 40\ \sigma^3$）を使用した。z方向にfixed境界を用いる理由は、shrink-wrap境界ではボックス上端が動的に変化し、上方から分子を挿入するcreate_atomsコマンドが正常に機能しないためである。
3.3.2 基板モデル
基板はLJ 9-3壁面ポテンシャル（wall/lj93）としてz = 0面に配置した：
$$U_{wall}(z) = \varepsilon_{wall}\left[\frac{2}{15}\left(\frac{\sigma}{z}\right)^{9} - \left(\frac{\sigma}{z}\right)^{3}\right], \quad z < z_c$$
デフォルトパラメータ: $\varepsilon_{wall}$ = 2.5, $\sigma$ = 1.0, cutoff = 2.5。基板を原子の集合体ではなく連続壁面ポテンシャルとして扱うことで、計算コストを大幅に削減しつつ、基板-分子間相互作用の本質を捉えた。ボックス上端にはreflect壁（wall/reflect）を設置し、蒸着中の分子の散逸を防止した。
3.3.3 蒸着プロトコル
蒸着は以下の逐次プロセスで実行した：
- 分子をランダムなxy位置、z軸周りにランダム回転した状態で、現在の膜表面より$3.0\sigma$上方に挿入
- 基板方向（-z方向）に$V_{dep} = 3.0$の初速度を付与
- NVE/limit積分器 + Langevin熱浴（T = 0.7, 減衰定数 = 1.0）で20,000ステップの緩和を実施
- 上記1-3を100分子について繰り返し
- 最終平衡化として50,000ステップのNVT緩和を実施
合計蒸着ステップ数: 100分子 × 20,000ステップ + 50,000ステップ = 2,050,000ステップ。時間刻み dt = 0.001で実施した。
3.3.4 周期系への変換
蒸着膜は壁面の影響を含むため、凝集法との公平な比較のために完全周期系への変換を行った。蒸着膜の原子座標を抽出し、全方向周期境界のボックスに配置した後、蒸着温度（0.7）でNPT平衡化し、その後AGG_FINAL_TEMP（0.4）まで冷却した。
3.4 凝集シミュレーション（バルク急冷モデル）
3.4.1 初期配置
100分子を$30 \times 30 \times 30\ \sigma^3$の完全周期ボックス内にランダムに分散配置した。
3.4.2 凝集プロトコル
凝集シミュレーションは4段階で実施した：
Phase 1 — ソフトプッシュオフ（10,000ステップ）: 分子の重なりを解消するため、反発力のみのソフトポテンシャルを使用。
Phase 2 — NVT平衡化（50,000ステップ）: 実際のLJ力場に切り替え、Langevin熱浴（T = 1.0）でNVT平衡化を実施。
Phase 3 — NPT圧縮（400,000ステップ）: 3段階のNPT圧縮で密な凝集体を形成。 - Stage 1: P = 0.1 → 0.5（150,000ステップ） - Stage 2: P = 0.5 → 1.0（150,000ステップ） - Stage 3: P = 1.0で平衡化（100,000ステップ）
Phase 4 — 急冷 + 最終平衡化（30,000ステップ）: T = 1.0 → 0.4の急冷（20,000ステップ）後、T = 0.4で10,000ステップの平衡化。
3.4.3 凝集法の特徴
凝集法ではバルク全体が同時にガラス転移を経験するため、全分子が同時に運動性を失う。これは蒸着法における逐次的表面最適化とは根本的に異なるメカニズムである。
3.5 熱スイープ解析
両手法で得られた薄膜構造に対して、同一のNPT熱スイーププロトコルを適用した。
- 温度範囲: $T^* = 0.3 \sim 2.0$（刻み幅 0.1、計18点）
- 各温度で10,000ステップの平衡化 + 20,000ステップのプロダクションラン
- 各温度で密度（$\rho$）とポテンシャルエネルギー（PE/atom）を計測
- ガラス転移温度（$T_g$）は密度の温度微分 $d(\rho)/dT$ の極小値から推定
3.6 パラメトリック研究
3.6.1 蒸着パラメトリック研究（50分子系）
計算コスト削減のため、パラメトリック研究では50分子系を使用した。以下の4パラメータについて系統的に検討した。
表5: 蒸着パラメトリック研究の設定
| パラメータ | 変数名 | 検討値 | デフォルト値 |
|---|---|---|---|
| 蒸着速度 | $V_{dep}$ | 0.5, 1.0, 3.0, 5.0, 10.0 | 3.0 |
| 基板温度 | $T_{dep}$ | 0.3, 0.5, 0.7, 0.9, 1.1 | 0.7 |
| 基板相互作用 | $\varepsilon_{wall}$ | 0.5, 1.5, 2.5, 4.0 | 2.5 |
| 成膜段階 | $N_{mol}$ | 5, 10, 20, 50 | 50 |
3.6.2 凝集パラメトリック研究（50分子系）
以下の5パラメータについて系統的に検討した。
表6: 凝集パラメトリック研究の設定
| パラメータ | 変数名 | 検討値 | デフォルト値 |
|---|---|---|---|
| 冷却速度 | quench_steps | 5,000 / 10,000 / 20,000 / 50,000 / 100,000 | 20,000 |
| 圧縮圧力 | $P_{final}$ | 0.1, 0.5, 1.0, 2.0, 5.0 | 1.0 |
| C-C相互作用 | $\varepsilon_{CC}$ | 1.0, 2.0, 3.0, 4.0 | 2.0 |
| 主鎖剛性 | $K_{CCC}$ | 10, 30, 50, 100 | 50 |
| アニーリング | 条件 | 急冷(10k) / 急冷(100k) / アニール(T=1.5) / アニール(T=2.0) | 急冷(20k) |
3.7 解析手法
- 密度: NPTシミュレーション中のボックス体積と総質量から算出
- ポテンシャルエネルギー: 系全体のPEを原子数で正規化（PE/atom）
- ガラス転移温度: 密度-温度曲線の温度微分の極小値として推定
- 動径分布関数（RDF）: dump fileからビーズ間距離の分布を解析
- 膜厚プロファイル: z方向の密度分布（成膜過程の解析用）