同じような計算なのに結果が全く異なる。
解決したいこと
u_8_up_1[N]、u_8_dw_1[N]という数を計算しているのだが、ある程度のズレは理解できるが、大きく異なるポイントが出てしまう。
発生している問題・エラー
例えばN＝12ならu_8_up_1[N]、u_8_dw_1[N]をprintfすると
0.049526 0.049195
0.287440 0.291578
0.308701 0.312578
0.366130 0.385112
0.437588 0.422055
0.467727 0.517449
0.423547 0.423635
0.444844 0.444611
0.455366 0.454081
0.463847 0.464400
0.471478 0.462416
0.537887 0.353247
のように最後の数字u_8_up_1[12]、u_8_dw_1[12]で大きく異なった値となってしまう。
該当するソースコード
//グラフェンナノリボンエネルギーバンド 相互作用 強磁性
#include <stdio.h>
#include <math.h>
#include <complex.h>
#include <stdlib.h>
#define M 6 // 存在するサイト数 //M=24で5nmくらい
#define N 12 //サイト数×２ A+Bの個数
#define P 100 //リボン長 全電子数 P*L
#define Q 10 //回数
double b_pi,N_c;
double complex t_1, t_2, t_3, t_4, t_5, t_6, t_7, t_8, T_1, T_2; //t_方向
double k_x[P+1];
double n_up[N],n_dw[N]; //代入する値
double n_1_up[M],n_1_dw[M],n_updw;
double meow_updw; //up&dwのフェルミ
double Energy_dw[P+1][N],Energy_up[P+1][N];
double u_up[P+1][N][N],u_dw[P+1][N][N]; //固有ベクトル
double fermi_up_1[P+1][N],fermi_dw_1[P+1][N];
double gap[P+1],gap_up[P+1],gap_dw[P+1]; //Eg,フェルミエネルギーバンド
double u_5_up_1[P+1][N][N],u_6_up_1[P+1][N],u_7_up_1[N],u_8_up_1[N];
double u_5_dw_1[P+1][N][N],u_6_dw_1[P+1][N],u_7_dw_1[N],u_8_dw_1[N];
void zheev_( char*, char*, int*, double _Complex*, int*, double*, double _Complex*, int*, double*, int* );
int main(void)
{
FILE* fp;
FILE* fp_1;
FILE* fp_2;
if ((fp = fopen("ziginteneferromagnetism M=.txt", "w")) == NULL)
{
printf("Cannot open the file\n");
exit(1);
}
if ((fp_1 = fopen("zigintferromagnetism_spin.txt", "w")) == NULL)
{
printf("Cannot open the file\n");
exit(1);
}
if ((fp_2 = fopen("zigintferromagnetism_spindraw.txt", "w")) == NULL)
{
printf("Cannot open the file\n");
exit(1);
}
b_pi=M_PI;
N_c=(P+1)*N;
double k_xx = -b_pi;
double kt=0.00001*8.6171*pow(10,-5); //-273℃＊ボルツマン定数(eV/K)
double U=4.32; //相互作用エネルギーの平均(eV)
double t=2.75; //飛び移り積分(eV)
double ribbon_length; //リボン幅
double Energy_all[(P+1)*2*N*2]; //全エネルギー
//刻み幅
for(int i=0;i<=P;i++){
k_x[i]=k_xx+2.0*b_pi/P*i;
}
//リボン幅
ribbon_length=0.246/pow(3,0.5)/2*(3*M-2); //0.0710140
//初期条件
n_up[0]=0.01;
n_up[1]=0.99;
n_up[M]=0.4;
n_up[M-1]=0.6;
n_up[N-2]=0.01;
n_up[N-1]=0.99;
n_dw[0]=0.99;
n_up[1]=0.01;
n_dw[M]=0.6;
n_dw[M-1]=0.4;
n_up[N-2]=0.99;
n_dw[N-1]=0.01;
/*for(int i=0;i<N;i++){
n_up[i]=0.4;
n_dw[i]=0.4;
}*/
int z=0;
while(z<Q){
//行列計算
//dw
for(int r=0;r<=P;r++){
// 複素行列(エルミート)の対角化 (入力行列の配列は対角化後にユニタリ行列になる)
char jobz = 'V' ;// 固有ベクトルを計算する
char uplo = 'U' ;// Aに上三角行列を格納
int n = N ;// 対角化する正方行列のサイズ
double _Complex A[N*N] ;// 対角化する行列。対角化後は固有ベクトルが並ぶ
t_4 = 1;
t_5 = cexp(-1 * I * -1/ 2 * k_x[r]);
t_6 = cexp(-1 * I * 1/2 * k_x[r]);
t_1 =1;
t_2 = cexp(I * -1 / 2 * k_x[r]);
t_3 = cexp(I * 1 / 2 * k_x[r]);
T_1 = t_5 + t_6;
T_2 = t_2 + t_3;
for(int i=0;i<N*N;i++){
A[i]=0;
}
for(int i=0;i<N;i++){
A[i*N+i]=U/t*n_up[i];
}
for(int j=1;j<M-1;j++){
A[j*N+M+j]=t_4; A[j*N+M+j+1]=T_1;
A[(j+M)*N+j]=T_2; A[(j+M)*N+j+1]=t_1;
}
A[M]=T_1; A[M+1]=0; A[N*M-2]=t_4; A[N*M-1]=T_1;
A[N*M]=T_2; A[N*M+1]=t_1; A[N*N-1-M]=T_2; A[N*N-M]=0;
double w[N] ;// 実固有値が小さい順に入る
int lda = N ;// 対角化する正方行列のサイズ
double _Complex work[6*N] ;// 対角化する際に使用するメモリ
int lwork = 6*N ;// workの次元
double rwork[3*N-2] ;// 3*SIZE-2で固定
int info ;// 成功すれば0、失敗すれば0以外を返す
zheev_( &jobz, &uplo, &n, A, &lda, w, work, &lwork, rwork, &info );
for(int i=0;i<N;i++){
for(int j=0;j<N;j++){
u_dw[r][i][j]=0.0;
u_5_dw_1[r][i][j]=0.0;
}
u_6_dw_1[r][i]=0.0;
u_7_dw_1[i]=0.0;
u_8_dw_1[i]=0.0;
}
for(int j=0;j<N;j++){
Energy_dw[r][j]=w[j];
}
for(int i=0;i<N;i++){
for(int j=0;j<N;j++){
if(i*N<=i<(i+1)*N){
u_dw[r][i][j]=A[i*N+j];
}
}
}
}
//up
for(int r=0;r<=P;r++){
// 複素行列(エルミート)の対角化 (入力行列の配列は対角化後にユニタリ行列になる)
char jobz = 'V' ;// 固有ベクトルを計算する
char uplo = 'U' ;// Aに上三角行列を格納
int n = N ;// 対角化する正方行列のサイズ
double _Complex A[N*N] ;// 対角化する行列。対角化後は固有ベクトルが並ぶ
t_4 = 1;
t_5 = cexp(-1 * I * -1/ 2 * k_x[r]);
t_6 = cexp(-1 * I * 1/2 * k_x[r]);
t_1 =1;
t_2 = cexp(I * -1 / 2 * k_x[r]);
t_3 = cexp(I * 1 / 2 * k_x[r]);
T_1 = t_5 + t_6;
T_2 = t_2 + t_3;
for(int i=0;i<N*N;i++){
A[i]=0;
}
for(int i=0;i<N;i++){
A[i*N+i]=U/t*n_dw[i];
}
for(int j=1;j<M-1;j++){
A[j*N+M+j]=t_4; A[j*N+M+j+1]=T_1;
A[(j+M)*N+j]=T_2; A[(j+M)*N+j+1]=t_1;
}
A[M]=T_1; A[M+1]=0; A[N*M-2]=t_4; A[N*M-1]=T_1;
A[N*M]=T_2; A[N*M+1]=t_1; A[N*N-1-M]=T_2; A[N*N-M]=0;
double w[N] ;// 実固有値が小さい順に入る
int lda = N ;// 対角化する正方行列のサイズ
double _Complex work[6*N] ;// 対角化する際に使用するメモリ
int lwork = 6*N ;// workの次元
double rwork[3*N-2] ;// 3*SIZE-2で固定
int info ;// 成功すれば0、失敗すれば0以外を返す
zheev_( &jobz, &uplo, &n, A, &lda, w, work, &lwork, rwork, &info );
for(int i=0;i<N;i++){
for(int j=0;j<N;j++){
u_up[r][i][j]=0.0;
u_5_up_1[r][i][j]=0.0;
}
u_6_up_1[r][i]=0.0;
u_7_up_1[i]=0.0;
u_8_up_1[i]=0.0;
}
for(int j=0;j<N;j++){
Energy_up[r][j]=w[j];
}
for(int i=0;i<N;i++){
for(int j=0;j<N;j++){
if(i*N<=i<(i+1)*N){
u_up[r][i][j]=A[i*N+j];
}
}
}
}
//フェルミエネルギー up&dw
for(int i=0;i<=P;i++){
for(int j=0;j<N;j++){
Energy_all[i*N+j]=Energy_dw[i][j];
Energy_all[i*N+j+(P+1)*N]=Energy_up[i][j];
}
}
double tmp_ene;
for (int i=0; i<=P*N*2; i++) {
for (int j=i+1; j<=P*N*2; j++) {
if (Energy_all[i] > Energy_all[j]) {
tmp_ene = Energy_all[i];
Energy_all[i] = Energy_all[j];
Energy_all[j] = tmp_ene;
}
}
}
meow_updw=(Energy_all[(P+1)*N]+Energy_all[(P+1)*N+1])/2;
/*for(int i=0;i<=P*N*2;i++){
printf("%f\n",Energy_all[i]);
}
printf("%f\n",meow_updw);*/
for(int l=0;l<=P;l++){
for(int i=0;i<N;i++){
fermi_dw_1[l][i]=1/(1+exp((Energy_dw[l][i]-meow_updw)/kt));
fermi_up_1[l][i]=1/(1+exp((Energy_up[l][i]-meow_updw)/kt));
}
}
for(int i=0;i<N;i++){
for(int l=0;l<=P;l++){
for(int j=0;j<N;j++){
u_5_dw_1[l][i][j]=pow(fabs(u_dw[l][i][j]),2)*fermi_dw_1[l][i];
u_6_dw_1[l][i]+=u_5_dw_1[l][j][i];
u_5_up_1[l][i][j]=pow(fabs(u_up[l][i][j]),2)*fermi_up_1[l][i];
u_6_up_1[l][i]+=u_5_up_1[l][j][i];
}
u_7_dw_1[i]+=u_6_dw_1[l][i];
u_7_up_1[i]+=u_6_up_1[l][i];
}
u_8_dw_1[i]=u_7_dw_1[i]/(P+1);
u_8_up_1[i]=u_7_up_1[i]/(P+1);
}
for(int i=0;i<N;i++){
n_up[i]=u_8_dw_1[i];
n_dw[i]=u_8_up_1[i];
}
z++;
}
//Eg
double tmp,tmp_dw, tmp_up;
for (int i=0; i<=P; i++) {
gap_up[i]=Energy_up[i][M]-Energy_up[i][M-1];
gap_dw[i]=Energy_dw[i][M]-Energy_dw[i][M-1];
gap[i]=Energy_up[i][M]-Energy_dw[i][M-1];
// gap[i]=Energy_up[i][M]-Energy_up[i][M-1];
}
for (int i=0; i<=P; i++) {
for (int j=i+1; j<=P; j++) {
if (gap[i] > gap[j]) {
tmp = gap[i];
gap[i] = gap[j];
gap[j] = tmp;
}
}
}
for(int i=0;i<N;i++){
n_updw+=n_up[i]*n_dw[i];
}
for(int i=0;i<=P;i++){
for(int j=0;j<N;j++){
fprintf(fp, "%2f %lf %lf %lf %lf\n",k_x[i],Energy_up[i][j]-meow_updw,Energy_dw[i][j]-meow_updw,Energy_up[i][M]-meow_updw,Energy_up[i][M-1]-meow_updw);
// fprintf(fp, "%2f %lf %lf %lf %lf\n",k_x[i],Energy_dw[i][j]-meow_updw,Energy_up[i][j]-meow_updw,Energy_dw[i][M]-meow_updw,Energy_dw[i][M-1]-meow_updw);
// fprintf(fp, "%2f %lf %lf\n",k_x[i],Energy_up[i][j]-meow_updw-U/t*n_updw,Energy_dw[i][j]-meow_updw-U/t*n_updw);
// // fprintf(fp, "%2f %lf %lf\n",k_x[i],Energy_up[i][M]-meow_updw,Energy_dw[i][M-1]-meow_updw);
}
}
for(int i=0;i<N;i++){
fprintf(fp_1,"%lf %lf\n",n_up[i],n_dw[i]);
}
for(int i=0;i<M;i++){
fprintf(fp_2,"%lf %lf\n",n_up[i],n_dw[i]);
fprintf(fp_2,"%lf %lf\n",n_up[i+M],n_dw[i+M]);
}
printf("リボン長：%lf\n", ribbon_length);
printf("最小ギャップ：%lf\n",gap[0]/(U/t));
fclose(fp);
fclose(fp_1);
fclose(fp_2);
return 0;
}
自分で試したこと
非常に専門性が高く、長ったらしいプログラムですが、そこのバグが入る原因がわからないです。もしかすると仕様なのかもしれませんが、なぜここまで違う数字が入るのかの原因がわからないです。
環境はubuntuでコンパイルはicc -qmklを用いて計算しています。
0
