まぃふぇいばりっと

機械学習やってます.Julia大好きです.勉強したことの殴り書きです.

Zigzag graphene nanoribbons のエネルギー分散関係

大昔に研究で使っていた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)を表しています f:id:physics303:20210215102702p:plain