大昔に研究で使っていたCのコードが発掘された.
グラフェンのジグザグナノリボンのバンド図を求めます。Nはユニットセル内の原子の数です。最近接のとびうつりのみを考えたタイトバインディングで計算しています。計算結果は/dataにdatファイルとして出力します。ハミルトニアンの対角化は,Lapackでやってます.
#include <stdio.h> #include <math.h> #include <complex.h> /**************************************************************** zigzagグラフェンナノリボンの計算。波数はky方向にのみ定義できる。 *****************************************************************/ int main(void) { int N =70; //ユニットセル内の原子数 FILE *saveband[N]; char bandname[50]; //エネルギーバンドを保存するファイル名 // Lapackに必要な変数 char jobz = 'V' ;// 固有ベクトルを計算する char uplo = 'U' ;// Aに上三角行列を格納 int n = N ;// 対角化する正方行列のサイズ double _Complex H[N*N] ;// 対角化する行列。対角化後は固有ベクトルが並ぶ double w[N] ;// 実固有値が小さい順に入る int lda = N ;// 対角化する正方行列のサイズ double _Complex work[6*N];// 対角化する際に使用するメモリ int lwork = 6*N ;// workの次元 double rwork[3*N-2] ;// 固定 int info ;// 成功すれば0、失敗すれば0以外を返す int i,j,p; double t=1; //transfer integral double a=1; //lattice constant double ky,kx; //wave vector kx=0 double fin=100; //計算の細かさ double V=0.0; //ユニットセル内の2つの原子のオンサイトポテンシャルの差 //保存するファイル名を決定する。 for(j=0;j<N;j++){ sprintf(bandname,"./data/band%d.dat",j); saveband[j] = fopen(bandname,"w"); } //各kyについてハミルトニアンを構成する for( p=-fin ; p<fin ; p++){ ky = 3.14*p/fin; for(i=0;i<N*N;i++){ H[i]=0.0; } for(j=0;j<N;j++){ if(2*j*(N+1) < N*N) H[2*j*(N+1)] = V ; if(2*j*(N+1)-1 > 0 && 2*j*(N+1) < N*N ) H[2*j*(N+1)-1] = t ; if(2*j*(N+1)+1 < N*N) H[2*j*(N+1)+1] = 2*t*cos(a*ky/2.0) ; if((2*j+1)*(N+1) < N*N) H[(2*j+1)*(N+1)] = -V; if((2*j+1)*(N+1)-1 > 0 && (2*j+1)*(N+1)-1 < N*N ) H[(2*j+1)*(N+1)-1] = 2*t*cos(a*ky/2.0) ; if((2*j+1)*(N+1)+1 < N*N) H[(2*j+1)*(N+1)+1] = t; } zheev_( &jobz, &uplo, &n, H, &lda, w, work, &lwork, rwork, &info ); for(j=0;j<N;j++){ fprintf(saveband[j],"%f\t%f\n",ky,w[j]); } } for(j=0;j<N;j++){ fclose(saveband[j]); } return 0; }
N=70の計算結果をgnuplotっでプロットさせたものです。横軸は波数、縦軸はエネルギー(eV)を表しています