Monologue-ひとりごと

日々の考えたこと

ラズパイPico + MicroPythonでのプログラミング-4

前回の記事から少し間が空きましたが、ラズパイPico+MicroPythonを用いた電子工作入門を再開していきます

tamurata435.hatenablog.com

今回からシリアル通信を試していきたいと思います。初めにUARTによる通信を試していきます。

UARTの概要

UART は I2C、SPI と並んで、マイコンで最も一般的に利用されるシリアル通信方式です。 シリアル通信は、1 本の信号線を用いて(双方向通信の場合は 2 本)、データを 1 ビットずつ送受信する方式です。送信側が 1 ビットを送り、受信側がそれを受け取るという単純な動作を繰り返すことでデータを伝送します。 このような通信を行うには、受信側・送信側ともに共通の信号(クロック)に従って動作する必要があります。 クロック信号を送受信側で共有する方式を同期通信、共有しない方式を非同期通信と呼び、UART は名称の “Universal Asynchronous Receiver Transmitter” の通り非同期通信に分類されます。

UARTのメリット

非同期通信である UART ではクロック信号を共有しないため、送受信間の配線は信号線と GND のみで済みます。 実装がシンプルで、クロック同期を気にする必要がないため長距離通信にも向いています。 I2C や SPI と比べると通信速度は遅いものの、扱いやすさから産業用の計測機器などで未だに使用する場面があります。

UARTの通信方式

クロックを共有しないため、送受信側であらかじめ1ビットの時間幅(ボーレート)やデータ長を一致させておく必要があります。 また、通信の開始と終了を示すために、データにスタートビットとストップビットを付与する必要があり、オーバーヘッドが大きく効率はあまり良くありません。 一般的には、スタートビット 1bit、データビット 8bit(パリティを含む場合は 7bit)、ストップビット 1〜2bit といった構成がよく使われます。

一般的に無信号区間は"1"の状態を維持するため、スタートビットは"0", ストップビットは"1"で固定となります。

ラズパイPicoによるUART

ラズパイPicoに搭載されたマイコンRP2040は2個のUARTコントローラーを持ち、以下のGPIOピンにマッピングすることが出来ます。

id TX RX
0 GPIO 0, 12, 16 GPIO 1, 13, 17
1 GPIO 4, 8 GPIO 5, 9

UARTのid番号、ボーレート、TX, RXとして使用するPinオブジェクトなどを引数にしてUARTオブジェクトを作成して通信を開始します。

uart1 = UART(0, baudrate=9600, tx=Pin(0), rx=Pin(1))

サンプルプログラム

2 台の Pico を使って送受信を確認するのが分かりやすいのですが、今回は都合により 1 台の Pico 内で通信を行います。uart0 と uart1 の 2 つの UART オブジェクトを作成し、uart0 から送信したデータを uart1 で受信します。 例として、Pico 内蔵の CPU 温度センサの値を読み取り、その値を UART で送受信してみます。 ジャンパピンの接続は、2 つの UART 間で TX と RX がクロスするように接続します(2 台の Pico を使う場合は GND を共通にする必要があります)。

# Raspberry Pi PicoのCPU温度を読み込み、UARTで送信/受信する

import time, utime
from machine import Pin, ADC, UART

temp_sensor = ADC(4)
led = Pin("LED", Pin.OUT, value=0)
uart0 = UART(0, baudrate=9600, tx=Pin(0), rx=Pin(1))
uart1 = UART(1, baudrate=9600, tx=Pin(4), rx=Pin(5))


def send_temp():
    read_value = temp_sensor.read_u16() * (3.3 / 65535)
    temp = 27 - (read_value - 0.706) / 0.001721
    message = f"temperature: {temp:.1f} degree\n"
    uart0.write(message)
    time.sleep(1)


def receive_temp():
    if uart1.any():
        led.value(1)
        print(uart1.readline())
        utime.sleep_ms(50)
        led.value(0)
    else:
        utime.sleep_ms(50)


while True:
    send_temp()
    receive_temp()

1秒ごとに温度が受信できればOKです!

Quantum Espressoで第一原理計算#4 スピン軌道相互作用

今回は、スピン軌道相互作用(Spin Orbit Coupling: SOC)を考慮した計算を行ってみます。

スピン軌道相互作用(SOC)について

原子核の周囲を回る電子について、電子が静止した座標系を考えてみます。 この座標系では、電子の周囲を原子核が回転しているように見なすことができます。 正に帯電した原子核の回転による環状電流により、原子核周囲には磁場が生じます。

この磁場と電子スピンが相互作用する現象がスピン軌道相互作用です。 一般に重元素ほど効果が大きく、軽元素での効果は無視できる場合が多いとされています。

SOCを考慮した計算

SOCを考慮した計算を行うには、相対論的な効果が含まれた(full relativisticな)擬ポテンシャルを使用する必要があります。 適切な擬ポテンシャルを用意した上で、&SYSTEMブロックに以下のパラメータを追記します。

&SYSTEM
    noncolin = .TRUE.
    lspinorb = .TRUE.

"lspinorb"を.TRUE.にすることで、スピン軌道相互作用を考慮した計算が行えます。

"noncolin"を.TRUE.にすることで、任意の方向のスピンを取ることができるようになります。 ただ、このパラメータの物理的な意味についてはまだあまり理解できていません……。 非磁性な物質ではスピンの向きはランダムで全体として打ち消し合った状態となるため、そのような状況を扱うための設定なのでしょうか。

バンド分散図の計算

今回はSi, GaAsについてSOCを考慮したバンド分散図の計算を行い、SOCなしの場合と比較しました。 SOCなしの結果は、過去の計算結果を流用しています。

tamurata435.hatenablog.com

まず、SiのSOC有無での分散図を比較してみます。 大きな変化は見られず、SOCの効果が小さいことが分かります。 実際にはVBM付近のバンドに数十meVの分裂が生じるようですが、実用上の影響はほぼ無視できると言えそうです。

Siのバンド分散図

次に、GaAsについて見ていきましょう。 こちらは明確に変化が見られます。 VBMに注目すると、SOC無しでは縮退していたバンドが分裂しています。 また伝導帯についても同様に分裂が生じています。 以上より、重元素におけるSOCの効果が確認できました。

GaAsのバンド分散図

計算メモ

GaAsのバンド分散図計算において、nscf計算がなかなか収束せず苦労しました。 conv_thrの値を緩和することで(1.0d-8→1.0d-6)、なんとか計算を完了させることができました。

Quantum Espressoで第一原理計算#3 構造最適化

今回は、構造最適化についての理解を深めます。

構造最適化を行うには、入力ファイルの&CONTROLにおいて、calculationを'relax'もしくは'vc-relax'と指定します。 relaxは格子定数は固定した状態で単位胞内の原子位置のみ最適化、vc-relaxは原子位置に加えて格子定数も最適化する際に用います(vcはvariable-cellの略らしい)。 具体的な使い分けが理解できていないのですが、前者は回折実験で正確な格子定数が得られている場合に使うイメージでしょうか。

今回は、Siを計算対象としてvc-relaxによる構造最適化を行いました。 格子定数の初期値を変化させ、構造最適化を行った結果を見ていきます。

入力ファイルの作成

以下の入力ファイルを作成し、pw.xで計算を実行します。 今回は格子定数の初期値を5.2, 5.3, 5.4, 5.5, 5.6[Å]と変化させて計算しました。シェルスクリプトでLATTICE_CONSTANTの部分を置換し、処理を実行させていきます。

&CONTROL
    calculation = 'vc-relax',
    prefix = 'si',
    outdir = './tmp',
    pseudo_dir = './pseudo'
    etot_conv_thr = 1.0d-5,
    forc_conv_thr = 1.0d-4,
/

&SYSTEM
    ibrav = 2,
    a = LATTICE_CONSTANT,
    nat = 2,
    ntyp = 1,
    ecutwfc = 50,
    ecutrho = 200,
    occupations = 'fixed'
/

&ELECTRONS
    conv_thr = 1.0d-8
/

&IONS
/

&CELL
/

ATOMIC_SPECIES
    Si 28.0855 Si.pbe-n-kjpaw_psl.1.0.0.UPF
ATOMIC_POSITIONS {crystal}
    Si 0.00 0.00 0.00
    Si 0.25 0.25 0.25
K_POINTS {automatic}
    6 6 6 0 0 0

etot_conv_thr, forc_conv_thrは、それぞれ全エネルギーの変化量の閾値、原子に働く力の閾値を指定します。 &IONS, &CELLのブロックでは最適化に関する各種オプションを指定できます。オプションを指定しない場合もブロックを入力する必要があるので注意です。

結果の解析

出力を見ると、"Final enthalpy"の出力前後に以下の内容が記述されています。

bfgs converged in   5 scf cycles and   4 bfgs steps
     (criteria: energy <  1.0E-05 Ry, force <  1.0E-04 Ry/Bohr, cell <  5.0E-01 kbar)

     End of BFGS Geometry Optimization

     Final enthalpy           =     -93.4518849618 Ry

     File ./tmp/si.bfgs deleted, as requested
Begin final coordinates
     new unit-cell volume =    276.12823 a.u.^3 (    40.91798 Ang^3 )
     density =      2.27954 g/cm^3

CELL_PARAMETERS (alat= 10.01554846)
  -0.516042576  -0.000000000   0.516042576
  -0.000000000   0.516042576   0.516042576
  -0.516042576   0.516042576  -0.000000000

ATOMIC_POSITIONS (crystal)
Si              -0.0000000000       -0.0000000000        0.0000000000
Si               0.2500000000        0.2500000000        0.2500000000
End final coordinates

デフォルトでは、50回計算のループが回るようですが、今回は数回で収束したことが分かります。

計算結果を見ていきましょう。 "new unit-cell volume"から基本単位胞の体積、CELL_PARAMETERSから基本単位格子の単位ベクトル、ATOMIC_POSITIONSから最適化された原子位置が分かります。 今回の計算では、原子位置の変化はないようです。

格子定数の変化を見ていきましょう。構造最適化後の基本単位胞の体積は40.91798 Ang3となりました。 Siの場合、以下の式で慣用単位格子における格子定数が得られます(基本単位胞は慣用単位胞の体積の1/4となっているため)。

 (40.91798\times4)^\frac{1}{3} \simeq 5.47

もしくは、CELL_PARAMETERSの単位ベクトルからも格子定数が求められます。 alatに記載された値が基準の長さとなるので、各ベクトルにalatの値を掛けます。 なお、単位はBohrとなるので注意が必要です。

 \mathbf{a_{\mathnormal{1}}}=(-0.516042576,\  0,\  0.516042576)\times 10.01554846 \simeq (-5.168,\ 0,\ 5.168)

Åに変換するため、各成分に0.529を掛けます。

  \mathbf{a_{\mathnormal{1}}} \simeq (−2.734, 0, 2.734)

Siの基本単位格子の単位ベクトルの大きさと慣用単位格子における格子定数は以下の関係があります。

 |\mathbf{a_{\mathnormal{1}}}|=\frac{a}{\sqrt{2}}

従って格子定数は以下のように求められます。

 a=\sqrt{2} \times \sqrt{(-2.734)^{2} + (2.734)^{2}} \simeq 5.47

入力ファイルで与えた初期値と、最適化後の格子定数の関係をプロットしてみました。 いずれも5.47[Å]付近の値となっており、初期値に依存せず一定の値に収束していることが分かります。 また、実際の値(5.43Å)に近い値となっています。擬ポテンシャルを選ぶことで、より実験値に近い値が得られるかもしれません。

Quantum Espressoで第一原理計算#2 収束判定

これむでは特に気にせずに計算を進めましたが、全エネルギーの収束を確認することで、正しいSCF計算を行なったかを確認することが重要です。 ここでは、単体Siを例に、主要パラメータ(ecutwfc、k点分割数、ecutrho)を変化させた際のエネルギー収束を検証しました。

DFTでは電子密度を自己無撞着に求めますが、電子密度が正しい値に近づくほど全エネルギーは極小値に収束します。 そのため、パラメータを変化させたときの全エネルギーの変化を確認することで、適切なパラメータを判断できます。

使用した擬ポテンシャルはPAW型の"Si.pbe-n-kjpaw_psl.1.0.0.UPF"で、以下の入力ファイルをベースに計算を行いました。

&CONTROL
    calculation = 'scf',
    prefix = 'si',
    outdir = './tmp',
    pseudo_dir = './pseudo'
/

&SYSTEM
    ibrav = 2,
    celldm(1) = 10.26,
    nat = 2,
    ntyp = 1,
    ecutwfc = ECUTWFC,
    ecutrho = ECUTRHO,
    occupations = 'fixed'
/

&ELECTRONS
    conv_thr = 1.0d-8
/

ATOMIC_SPECIES
    Si 28.0855 Si.pbe-n-kjpaw_psl.1.0.0.UPF
ATOMIC_POSITIONS {crystal}
    Si 0.00 0.00 0.00
    Si 0.25 0.25 0.25
K_POINTS {automatic}
    KPOINTS KPOINTS KPOINTS 0 0 0

計算は、以下のようなシェルスクリプトを作成し、実行しました。

for ECUTWFC in 20 30 40 50 60 70; do
    sed "s/ECUTWFC/${ECUTWFC}/g" ./input/ecutwfc.in > ./input/ecutwfc_${ECUTWFC}.in
    pw.x < ./input/ecutwfc_${ECUTWFC}.in > ./output/ecutwfc_${ECUTWFC}.out
done

ecutwfcの決定

ecutwfcは、波動関数を平面波基底で展開する際に、どこまで短波長成分を取り込むかを決めるパラメータになります。 値を大きくすると波動関数の表現精度は向上しますが、その分計算コストも増大します。精度とコストのバランスを考慮して設定することが必要です。

下にecutwfcを20〜70 [Ry]の範囲で10刻みで変化させてSCF計算を行った際のエネルギーの変化を示します。 値の増加に伴い、エネルギーが一定値へと収束していくことが確認できました。 特に50 [Ry]付近でほぼ収束しており、擬ポテンシャルに記載された推奨値(44 [Ry])と整合性が取れた結果となりました。

k点分割数の決定

次に、k点分割数を2×2×2から10×10×10まで変化させて計算を行いました。 波数と波長が逆数関係にあることを踏まえると、k点分割数の増加はどこまで長波長成分を取り込むかを決めるパラメータとなります。

ecutwfcは50 [Ry]に固定したうえで計算した結果、k点分割数の増加に伴い全エネルギーが収束していく様子が確認できました。 6×6×6を超えたあたりでエネルギーがほぼ一定となり、この系に対して十分な分割数であると判断できます。

ecutrhoの決定

ecutrhoは、電荷密度およびポテンシャルの表現精度を決めるパラメータです。 計算コストへの影響は小さいものの、値が小さすぎるとSCF収束が不安定になるとのことです。

ecutwfc=50 [Ry]、k点分割数=6×6×6の条件下で、ecutrhoを100〜400 [Ry]の範囲で50刻みで変化させました。 ecutrhoが増加するにつれて全エネルギーが一定値へと収束し、200 [Ry]付近でほぼ安定していることが分かりました。 これは、擬ポテンシャルの推奨値(175 [Ry])とも一致した結果となっています。

エネルギー変化量に注目すると、ecutwfcやk点分割数を変化させた場合と比べて非常に小さく、ほぼ誤差範囲のようにも見えます。 一方で、構造最適化では電荷密度の精度が重要となるとのことで、やはり推奨値以上のecutrhoを設定することが望ましいでしょう。 なお、擬ポテンシャルの種類によって、適切な値は変化するそうです。一般的に、ノルム保存型擬ポテンシャルではecutwfcの約4倍、US/PAWでは約8倍以上が目安とのことです。

Quantum Espressoで第一原理計算#1 擬ポテンシャルについて

これまで、半導体・金属・磁性体といった代表的な系について、一通りの基本的な計算を行ってきました。

今回からは、入門レベルから一歩進んで、各種パラメータの意味をより深く理解したり、応用的な計算にも挑戦していきます。 とはいえ、私自身まだ初心者なので、アドバイスやご指摘があればぜひ教えてください。

最初のテーマとして、これまで「なんとなく」で設定していたエネルギーカットオフや k 点数について、見直していきます。 これらの値は使用する擬ポテンシャルの種類によって適切な値が変わるため、導入として、擬ポテンシャルの概要を(自分用の整理として)まとめておきます。

擬ポテンシャルの概要

擬ポテンシャル法では、原子核および内殻電子の影響を有効ポテンシャルとして置き換え、そのポテンシャルの下で価電子の状態を計算します。 内殻電子が直接関与する現象(例えば XAS など)の計算には向きませんが、多くの物性では内殻電子の寄与は小さいため、非常に有効な手法となります。

なぜ擬ポテンシャルが必要となるのか。

第一原理計算では波動関数を求めますが、その表現方法として、原子軌道基底を使う方法や、平面波基底を使う方法などがあります。

平面波基底は固体の周期構造と相性が良く、数値的に扱いやすいという利点があります。 しかしながら、原子核付近の急峻な波動関数の変化を平面波で表現しようとすると、非常に多くの基底関数が必要となり、計算コストが大きくなってしまいます。 そこで、内殻電子を明示的に扱わず、ポテンシャルとして置き換えることで、計算コストを抑えつつ、物性に重要な価電子の挙動を記述するアプローチが擬ポテンシャル法です。 Quantum ESPRESSO は平面波基底を用いるコードであるため、計算には擬ポテンシャルの準備が必要です。

擬ポテンシャルの種類

Quantum Espressoで使用可能な擬ポテンシャルは、大別して3種類に分類されます。

  • ノルム保存型(Norm-Conserving: NC)
  • ウルトラソフト(Ultrasoft: US)
  • PAW(Projector Augmented-Wave)

いずれも、全電子計算で得られた波動関数をもとに、価電子・内殻電子を分離し、原子核付近の急峻な変化を取り除くことで構築されます。 具体的には、原子核の中心として球対称座標を考え、価電子と内殻電子の境界となる半径を r_{c}とし、 r \lt r_{c}の領域(内殻電子領域)で滑らかな擬波動関数を定義します。

ノルム保存型 (NC)

NCでは、 r \lt r_{c}の領域において、擬波動関数と全電子波動関数電荷量(ノルム)が一致するように、擬ポテンシャルを構築します。 もっとも基本的で分かりやすい手法ですが、ノルム保存という制約上、擬波動関数が急峻になりやすく、ecutwfcで指定する波動関数のカットオフは比較的高く設定する必要があります。

ウルトラソフト (US)

USでは、ノルム保存条件を緩和し、より滑らかな擬波動関数を許容します。 これにより、基底関数の数を減らすことができ、計算コストを下げることができます。 ただし、ノルム保存条件の緩和により、擬波動関数と全電子波動関数電荷量が一致しないため、補正電荷を導入して整合性を保ちます。 この補償電荷は空間的に局在化しているため、電荷密度の観点では変化が急峻となりやすく、ecutrhoで指定する電荷密度のカットオフは他の擬ポテンシャルと比較して、高い値に設定する必要があります。

PAW (Projector Augmented-Wave)

PAWにおいても、滑らかな擬波動関数を定義しますが、USと異なり補償電荷を導入するのではなく、 全電子波動関数に一致するように、 r \lt r_{c}の領域で全電子波動関数と擬波動関数の差分を足していきます。 価電子部分の記述は擬波動関数でも十分であるため、精度が不十分な原子核近傍においてのみ全電子波動関数を扱うというようなアプローチになります。 これにより、ecutwfcの値はUS並みに抑えつつも、USほど高いecutrhoの値も不要という、良い所取りができます。

精度と計算コストの兼ね合いから、PAWを使用するのが主流となっているようです。