コアダンプする箇所がわからない。
解決したいこと
コアダンプしていると出るが、その原因がわからない。
発生している問題・エラー
Segmentation fault (コアダンプ)
該当するソースコード
//グラフェンナノリボンエネルギーバンド 相互作用 反強磁性
#include <stdio.h>
#include <math.h>
#include <complex.h>
#include <stdlib.h>
#define M 10 // 存在するサイト数
#define N 20 //サイト数×２ A+Bの個数
#define P 100 //リボン長 全電子数 P*L
#define Q 1 //回数
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_2[2*N]; //代入する値
double meow_updw; //up&dwのフェルミ
double Energy[P+1][2*N];
double complex u[P+1][2*N][2*N]; //固有ベクトル
double complex fermi_1[P+1][2*N];
double gap[P+1]; //Eg,フェルミエネルギーバンド
double u_5_1[P+1][2*N][2*N],u_6_1[P+1][2*N],u_7_1[2*N],u_8_1[2*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("ziginteneantiferromagnetism M=4.txt", "w")) == NULL)
{
printf("Cannot open the file\n");
exit(1);
}
if ((fp_1 = fopen("zigintantiferromagnetism.txt", "w")) == NULL)
{
printf("Cannot open the file\n");
exit(1);
}
if ((fp_2 = fopen("zigintantiferromagnetism_1.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]; //全エネルギー
//刻み幅
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
//初期条件
for(int i=0;i<2*N;i++){
n_2[i]=0.4;
}
int z=0;
while(z<Q){
//行列計算
for(int r=0;r<=P;r++){
// 複素行列(エルミート)の対角化 (入力行列の配列は対角化後にユニタリ行列になる)
char jobz = 'V' ;// 固有ベクトルを計算する
char uplo = 'U' ;// Aに上三角行列を格納
int n = 2*N ;// 対角化する正方行列のサイズ
double _Complex A[2*N*2*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<2*N*2*N;i++){
A[i]=0;
}
for(int i=0;i<2*N;i++){
if(i<N){
A[i*2*N+N+i]=U/t*n_2[i];
}
else{
A[2*N*(i+N)+i]=U/t*n_2[i];
}
}
for(int j=1;j<M-1;j++){
A[2*N*j+M+j-1]=t_4; A[j*2*N+M+j]=T_1;
A[(j+M)*2*N+j]=T_2; A[(j+M)*2*N+j+1]=t_1;
A[2*N*(j+N)+N+M+j-1]=t_4; A[2*N*(j+N)+N+M+j]=T_1;
A[(j+N+M)*2*N+N+j]=T_2; A[(j+N+M)*2*N+N+j+1]=t_1;
}
A[M]=T_1;
A[2*N*(M-1)+N-2]=t_4; A[2*N*(M-1)+N-1]=T_1;
A[2*N*(N+1)-M]=T_1;
A[2*N*(N+M)-2]=t_4; A[2*N*(M+N)-1]=T_1;
A[2*N*M]=T_2; A[2*N*M+1]=t_1;
A[2*N*(N-1)+M-1]=T_2;
A[2*N*(N+M)+N]=T_2; A[2*N*(N+M)+N+1]=t_1;
A[2*N*2*N-M-1]=T_2;
/*for(int i=0;i<2*N*2*N;i++){
fprintf(fp_2,"%2f\n",A[i]);
}*/
double w[2*N] ;// 実固有値が小さい順に入る
int lda = 2*N ;// 対角化する正方行列のサイズ
double _Complex work[6*2*N] ;// 対角化する際に使用するメモリ
int lwork = 6*2*N ;// workの次元
double rwork[3*2*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<2*N;i++){
for(int j=0;j<2*N;j++){
u[r][i][j]=0.0;
u_5_1[r][i][j]=0.0;
}
u_6_1[r][i]=0.0;
u_7_1[i]=0.0;
u_8_1[i]=0.0;
}
for(int j=0;j<2*N;j++){
Energy[r][j]=w[j];
}
for(int i=0;i<2*N;i++){
for(int j=0;j<2*N;j++){
u[r][i][j]=A[i*2*N+j];
}
}
}
//フェルミエネルギー up&dw
for(int i=0;i<=P;i++){
for(int j=0;j<2*N;j++){
Energy_all[i*2*N+j]=Energy[i][j];
}
}
double tmp_ene;
for (int i=0; i<=P*2*N; i++) {
for (int j=i+1; j<=P*2*N; 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-2])/2;
/*for(int i=0;i<=P*N*2;i++){
printf("%f\n",Energy_all[i]);
}
printf("%f\n",meow_updw);*/
//Eg
double tmp;
for (int i=0; i<=P; i++) {
gap[i]=Energy[i][N]-Energy[i][N-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 l=0;l<=P;l++){
for(int i=0;i<2*N;i++){
fermi_1[l][i]=1/(1+exp((Energy[l][i]-meow_updw)/kt));
}
}
for(int i=0;i<2*N;i++){
for(int l=0;l<=P;l++){
for(int j=0;j<2*N;j++){
u_5_1[l][i][j]=pow(cabs(u[l][j][i]),2)*fermi_1[l][j];
u_6_1[l][i]+=u_5_1[l][i][j];
}
u_7_1[i]+=u_6_1[l][i];
}
u_8_1[i]=u_7_1[i]/(P+1);
}
for(int i=0;i<2*N;i++){
n_2[i]=u_8_1[i];
}
z++;
}
for(int i=0;i<=P;i++){
for(int j=0;j<2*N;j++){
fprintf(fp, "%2f %lf %lf %lf\n",k_x[i],Energy[i][j]-meow_updw,Energy[i][N]-meow_updw,Energy[i][N-1]-meow_updw);
// fprintf(fp, "%2f %lf %lf\n",k_x[i],Energy[i][N+1]-meow_updw,Energy[i][N-2]-meow_updw);
}
}
for(int i=0;i<2*N;i++){
fprintf(fp_2,"%2f\n",n_2[i]);
}
for(int i=0;i<M;i++){
fprintf(fp_1,"%2f\n",n_2[i]);
fprintf(fp_1,"%2f\n",n_2[i+M]);
}
printf("リボン長：%2f\n", ribbon_length);
printf("最小ギャップ：%f\n",gap[0]/(U/t));
fclose(fp);
fclose(fp_1);
fclose(fp_2);
return 0;
}
自分で試したこと
その前にn_2[2*N]をn_2[N]としたときはコンパイルが通ったのですが、意味として合うのは前者です。
n のところを変えると、別のエラーが出ました。
0
