U[i]の値が変化しない
Discussion
Closed
解決したいこと
数値計算上でU[i]の値をいろいろ設定しているのだが、全く反映されない
発生している問題・エラー
０の値しか出てこない
該当するソースコード
//グラフェンナノリボンエネルギーバンド 相互作用 強磁性
#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 5 //回数
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 e[P+1],k_x[P+1],U[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 Energy_sum,Energy_sum_1;
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];
double magnetic_0[P+1],magnetic_1[P+1],magnetic_5[P+1];
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;
// FILE* fp_3;
// FILE* fp_4;
if ((fp = fopen("ferrosite1M=10.txt", "w")) == NULL)
{
printf("Cannot open the file\n");
exit(1);
}
if ((fp_1 = fopen("ferrosite2M=10.txt", "w")) == NULL)
{
printf("Cannot open the file\n");
exit(1);
}
if ((fp_2 = fopen("ferrosite5M=10.txt", "w")) == NULL)
{
printf("Cannot open the file\n");
exit(1);
}
/*if ((fp_3 = fopen("zigintdosferromagnetism.txt", "w")) == NULL)
{
printf("Cannot open the file\n");
exit(1);
}
if ((fp_4 = fopen("zigintspin.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 kt=300*8.6171*pow(10,-5); //27℃＊ボルツマン定数(eV/K)
//相互作用エネルギーの平均(eV) u/t
double t=1.0; //2.75; //飛び移り積分(eV)
double ribbon_length; //リボン幅
double Energy_all[(P+1)*2*N*2]; //全エネルギー
double meu_B= 5.788*pow(10,-5); //ボーア磁子 [/eV/T]
//刻み幅
for(int i=0;i<=P;i++){
k_x[i]=k_xx+2.0*b_pi/P*i;
U[i]=10/(P+1)*i;
printf("%f\n",U[i]);
}
//リボン幅
ribbon_length=0.246/pow(3,0.5)/2*(3*M-2); //0.0710140
for(int h=0;h<=P;h++){
//初期条件
//端が違う
n_up[0]=0.01;
n_up[1]=0.01;
n_up[M]=0.6;
n_up[M+1]=0.6;
n_up[N-2]=0.99;
n_up[N-1]=0.99;
n_dw[0]=0.99;
n_dw[1]=0.99;
n_dw[M]=0.4;
n_dw[M+1]=0.4;
n_dw[N-2]=0.01;
n_dw[N-1]=0.01;
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[h]/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++){
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[h]/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++){
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+1)*N*2; i++) {
for (int j=i+1; j<(P+1)*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+1)*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][j][i]),2)*fermi_dw_1[l][j];
u_6_dw_1[l][i]+=u_5_dw_1[l][i][j];
u_5_up_1[l][i][j]=pow(fabs(u_up[l][j][i]),2)*fermi_up_1[l][j];
u_6_up_1[l][i]+=u_5_up_1[l][i][j];
}
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_up_1[i];
n_dw[i]=u_8_dw_1[i];
}
z++;
}
//magnetic
magnetic_0[h]=(n_up[0]-n_dw[0])*meu_B;
magnetic_1[h]=(n_up[1]-n_dw[1])*meu_B;
magnetic_5[h]=(n_up[4]-n_dw[4])*meu_B;
}
for(int i=0;i<=P;i++){
fprintf(fp, "%lf %lf\n",U[i], magnetic_0[i]);
fprintf(fp_1, "%lf %lf\n",U[i], magnetic_1[i]);
fprintf(fp_2, "%lf %lf\n",U[i], magnetic_5[i]);
}
printf("リボン長：%lf\n", ribbon_length);
fclose(fp);
fclose(fp_1);
fclose(fp_2);
// fclose(fp_3);
// fclose(fp_4);
return 0;
}
自分で試したこと
非常に長く専門性が高い内容になりますが、正直ミス自体は基本的な気がします。
環境はlinuxのコマンドでやっています。
0