liboctave Wiki

量子系のダイナミクス

応用編

内容

シュレディンガー方程式を解いて、量子系の時間発展の様子を計算します。 量子系のハミルトニアンをHとすると、系の波動関数\psi(t)の時間発展はシュレディンガー方程式に従います。

i\frac{d}{dt}\psi(t)~=~H~\psi(t)

このシュレディンガー方程式を時間領域で差分化して、数値的に波動関数の時間発展を求めます。 時刻t+\Delta~tでの波動関数と時刻tでの波動関数との間には

\psi(t+\Delta~t)~=~e^{-iH\Delta~t}\psi(t)

という関係が成り立ちます。 ハミルトニアンの指数関数を求めてもいいのですが(行列の指数関数を求めたい)、 ここではこの関係式を、

(1+iH\Delta~t/2)\psi(t+\Delta~t)~=~(1-iH\Delta~t/2)\psi(t)

と近似し、\psi(t)から\psi(\Delta~t~+t)を 逐次的に求めてゆきます。ComplexMatrix::solve(ComplexColumnVector& b) メソッドを用いてこの方程式を解いています。 この近似法では確率保存則が保たれることが知られています。 http://www-cms.phys.s.u-tokyo.ac.jp/~naoki/CIPINTRO/CIP/wave.html のページの解説を参考にしました。詳しくはそちらを参照してください。

ヘッダファイル

fileQuantumDynamics.h

#ifndef _QuantumDynamics_h
#define _QuantumDynamics_h 

#include <iostream>
#include <vector>
#include <octave/config.h>
#include <octave/CColVector.h>
#include <octave/CMatrix.h>
#include <CDiagMatrix.h>

using namespace std;

//==================================================
struct Record {
  double time;
  ComplexColumnVector phi;
  Record(double itime, ComplexColumnVector iphi) :
    time(itime), phi(iphi) {};
};
typedef vector<Record> Data; 

//==================================================
class QuantumDynamics {
  ComplexMatrix H; // Hamiltonian
  ComplexColumnVector init_phi; // initival state
  double dt;
  double t_max;
  Data data;
  int skip;
public:
  QuantumDynamics(const ComplexMatrix& iH,  
		  const ComplexColumnVector& iinit_phi,  
		  double idt, double it_max, int iskip = 1) : 
    H(iH), init_phi(iinit_phi), dt(idt), t_max(it_max), skip(iskip) {};
  void calc();
  Data getData();
};

void QuantumDynamics::calc() {
  const Complex IU(0.0, 1.0);
  ComplexMatrix Unit(ComplexDiagMatrix(H.rows(),H.cols(),1.0));
  ComplexMatrix lhMatrix = Unit + 0.5 * dt * IU * H;
  ComplexMatrix rhMatrix = Unit - 0.5 * dt * IU * H;
  ComplexColumnVector phi = init_phi;
  ComplexColumnVector rhs;
 
  int count = 0;
  double t = 0.0;
  while (t < t_max) {
    if (count % skip == 0) {
      data.push_back(Record(t, phi));
    }
    rhs = rhMatrix * phi;
    phi = lhMatrix.solve(rhs);
    t += dt;
    count++;
  }
}

Data QuantumDynamics::getData() {
  return data;
}
 
#endif

サンプルプログラム

filemain2.cpp

#include "QuantumDynamics.h"

using namespace std; 

void print_probabilities(const Data& data) {
  for (int i = 0; i < data.size(); i++) {
    cout << i << " " <<data[i].time << " ";
    for (int j = 0; j < data[i].phi.rows(); j++) {
      cout << norm(data[i].phi(j)) << " ";
    }
    cout << endl;
  }
} 

int main() 
{
  // Hamiltonian
  double del = 0.0;
  double eps = 1.0;
  ComplexMatrix H(2,2);
  H(0,0) = 0.0; H(0,1) = eps;
  H(1,0) = eps; H(1,1) = del;

  // initial state
  ComplexColumnVector init_phi(2);
  init_phi(0) = 1.0;
  init_phi(1) = 0.0; 

  double dt = 1.e-2;
  // calculation of dynamics
  QuantumDynamics qDyn(H, init_phi, dt, 20.0, 1);
  qDyn.calc();
  print_probabilities(qDyn.getData());
}

解説

サンプルプログラムでは、エネルギー準位が縮退した2つの状態間にトンネルがあるような系の時間発展を求めています。 下に各状態の占有確率の時間発展を示します。 自由粒子の運動を求めるような、系の状態数が大きく、ハミルトニアンがスパースになる場合には、ここで用いたComplexMatrix::solve()を使った方法は効率的でないと思われます。 反復解法などを使うべきでしょう。 prob.png

コメント

&ref(): File not found: "wikilogo1.png" at page "量子系のダイナミクス";

[edit]