計算した結果のファイル出力が途中で終わる、または出力されない
解決したいこと
計算結果が正しくファイル出力されない。
計算結果はあっている。
発生している問題・エラー
↓のコードでいうとf_4のファイルへの出力が全く行われません。
またf_3のファイル入力も動いていない。
該当するソースコード
#include <stdio.h>
#include <math.h>
#include <complex.h>
#include <stdlib.h>
#define M 6 // 存在するサイト数
#define N 12 //サイト数×２ A+Bの個数
#define L 24 //対角化に必要 N*２
#define P 400
#define EPS 0.000001 //収束範囲
double b_pi;
double complex t_1, t_2, t_3, t_4, t_5, t_6, t_7, t_8, T_1, T_2, T_3, T_4,T_int; //,t_方向
double e[P+1],k_x[P+1];
double B_1[N][N], B_2[N][N],D[N][N];
double complex C[N][N];
double a[L][L],A[N];
double Energy[P+1][N];
double DOS[P+1][P+1][N];
double ddd[P+1][P+1],dd[P+1];
double fermi_1[P+1][N],fermi_2[P+1][N],fermisum_1[P+1],fermisum_2[P+1],f_2[P+1];
double complex u[P+1][N][N];
double u_2[P+1][N],u_3[P+1],u_4[P+1];
int main(int argc, char** argv){
b_pi=M_PI;
double k_xx = -b_pi;
double ene=-3.0;
double delta_E=2.75*pow(10,-2); //エネルギー分解能
double kt=273*1.38*pow(10,-23); //０℃＊ボルツマン定数
double meow=0.0; //科学ポテンシャル（フェルミエネルギー）
double U=3.0; //相互作用エネルギーの平均
double t=2.75; //飛び移り積分
double n_down[P+1][N]; //代入する値
double n_upper[P+1][N];
double n_updw[P+1][N];
FILE* fp;
FILE* fp_1;
FILE* fp_2;
FILE* fp_3;
FILE* fp_4;
FILE* fp_5;
if ((fp = fopen("zigintene.txt", "w")) == NULL)
{
printf("Cannot open the file\n");
exit(1);
}
if ((fp_1 = fopen("zigintdos.txt", "w")) == NULL)
{
printf("Cannot open the file\n");
exit(1);
}
if ((fp_2 = fopen("zigintupper.txt", "w")) == NULL)
{
printf("Cannot open the file\n");
exit(1);
}
if ((fp_4 = fopen("zigintdown.txt", "w")) == NULL)
{
printf("Cannot open the file\n");
exit(1);
}
/*fp_3 = fopen("zigintupper.txt", "r");
if(fp_3 == NULL) {
printf("ファイルを開くことが出来ませんでした．¥n");
}
for(int l=0; l<=P; l++){
for(int i=0;i<N;i++){
fscanf(fp_3, "%lf", &(n_upper[l][i]));
}
}*/
/*fp_5 = fopen("zigintdown.txt", "r");
if(fp_5 == NULL) {
printf("ファイルを開くことが出来ませんでした．¥n");
}
for(int l=0; l<=P; l++){
for(int i=0;i<N;i++){
fscanf(fp_5, "%lf", &(n_down[l][i]));
}
}*/
for(int i=0;i<=P;i++){
e[i]=ene+6.0/P*i;
k_x[i]=k_xx+2.0*b_pi/P*i;
}
for(int r=0;r<=P;r++){
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_7 = -1 * t_4;
t_8 = -1 * t_1;
T_1 = t_5 + t_6;
T_2 = t_2 + t_3;
T_3 = -1 * T_1;
T_4 = -1 * T_2;
for(int i=0;i<N;i++){
n_upper[r][i]=6.5;
n_down[r][i]=N-n_upper[r][i];
if(i%2==0){
n_updw[r][i]=0.5;
}
else{
n_updw[r][i]=n_upper[r][i]-n_down[r][i];
}
}
for (int k = 0;k < N;k++) {
for (int h = 0;h < N;h++) {
T_int=-U/t*n_updw[r][k];
for (int i = 0;i < M - 2;i++) {
if (h < M) {
if (k < M) {
C[h][k] = 0;
C[k][k]=-T_int;
}
else {
C[0][M] = T_1;
C[M - 1][N - 1] = T_1;
C[M - 1][N - 2] = t_4;
C[i + 1][M + i] = t_4;
C[i + 1][M + i + 1] = T_1;
}
}
else {
if (k < M) {
C[M][0] = T_2;
C[M][1] = t_1;
C[N - 1][M - 1] = T_2;
C[M + i + 1][i + 1] = T_2;
C[M + i + 1][i + 2] = t_1;
}
else {
C[h][k] = 0;
C[k][k]=-T_int;
}
}
}
}
}
int k, h;
for (k = 0;k < N;k++) {
for (h = 0;h < N;h++) {
D[h][k] = creal(C[h][k]);
B_1[h][k] = -1 * cimag(C[h][k]);
B_2[h][k] = cimag(C[h][k]);
}
}
for (k = 0;k < L;k++) {
for (h = 0;h < L;h++) {
if (h < N) {
if (k < N) {
a[h][k] = D[h % N][k % N];
}
else {
a[h][k] = B_1[h % N][k % N];
}
}
else {
if (k < N) {
a[h][k] = B_2[h % N][k % N];
}
else {
a[h][k] = D[h % N][k % N];
}
}
}
}
for(int i=0;i<N;i++){
for(int j=0;j<N;j++){
u[r][i][j]=0.0;
}
}
for(int i=0;i<N;i++){
u[r][i][i]=1.0;
}
while (1) {
double alpha, beta, gamma;
double s, c, w;
double wa, wb, wc;
double max;
int i, j, p, q, x, y;
//最大要素の行と列を検索
max = 0.0;
for (i = 0;i < L - 1;i++) {
for (j = i + 1;j < L;j++)
if (fabs(a[i][j]) > max) {
p = i;
q = j;
max = fabs(a[i][j]);
}
}
//収束したら解答打ち出し
if (max < EPS) break;
//sin, cos 計算
wa = a[p][p];
wb = a[p][q];
wc = a[q][q];
alpha = -wb;
beta = 0.5 * (wa - wc);
gamma = fabs(beta) / sqrt(alpha * alpha + beta * beta);
s = sqrt(0.5 * (1.0 - gamma));
if (alpha * beta < 0) s = -s;
c = sqrt(1.0 - s * s);
//直行変換
for (j = 0;j < L;j++) {
w = a[p][j] * c - a[q][j] * s;
a[q][j] = a[p][j] * s + a[q][j] * c;
a[p][j] = w;
}
for (j = 0;j < L;j++) {
a[j][p] = a[p][j];
a[j][q] = a[q][j];
}
w = 2.0 * wb * s * c;
a[p][p] = wa * c * c + wc * s * s - w;
a[q][q] = wa * s * s + wc * c * c + w;
a[p][q] = 0;
a[q][p] = 0;
for (int i = 0;i < N;i++)
{
A[i] = a[i][i];
}
double tmp;
for (int i=0; i<N; i++) {
for (int j=i+1; j<N ;j++) {
if (A[i] > A[j]) {
tmp = A[i];
A[i] = A[j];
A[j] = tmp;
}
}
}
for(i=0;i<N;i++){
w=u[r][i][p]*c-u[r][i][q]*s;
u[r][i][q]=u[r][i][p]*s+u[r][i][q]*c;
u[r][i][p]=w;
}
}
for(int i=0;i<N;i++){
Energy[r][i]=A[i];
}
}
//状態密度
for(int i=0;i<=P;i++){
for(int j=0;j<=P;j++){
for(int k=0;k<N;k++){
DOS[i][j][k]=1/pow(b_pi*2,0.5)/delta_E/b_pi*exp(-1*pow(e[i]-Energy[j][k],2)/2/delta_E/delta_E);
}
}
}
for(int i=0;i<=P;i++){
for(int j=0;j<=P;j++){
for(int k=0;k<N;k++){
ddd[i][j]+=DOS[i][j][k];
}
}
}
for(int i=0;i<=P;i++){
for(int j=0;j<=P;j++){
dd[i]+=ddd[i][j];
}
}
//波動関数の二乗＝確率
for(int l=0;l<=P;l++){
for(int i=0;i<N;i++){
for(int j=0;j<N;j++){
u_2[l][i]+=pow(cabs(u[l][i][j]),2);
}
}
}
for(int l=0;l<=P;l++){
for(int i=0;i<N;i++){
u_3[l]+=u_2[l][i];
}
}
for(int l=0;l<=P;l++){
u_4[l]=u_3[l]/N;
//printf( "%7.4lf\n", u_4[l]); //１
}
//個数をフェルミ関数から計算
for(int l=0;l<=P;l++){
for(int i=0;i<N;i++){
if(Energy[l][i]<=0){
fermi_1[l][i]=1/(1+exp((Energy[l][i]-meow)/kt)); //電子（↓側バンド）のフェルミ関数
}
if(Energy[l][i]>=0){
fermi_2[l][i]=1/(1+exp((-fabs(Energy[l][i])-meow)/kt)); //正孔（上側バンド）のフェルミ関数
}
}
}
for(int l=0;l<=P;l++){
for(int i=0;i<N;i++){
fermisum_1[l]+=fermi_1[l][i];
fermisum_2[l]+=fermi_2[l][i];
f_2[l]=fermisum_1[l]+fermisum_2[l];
}
}
//printf( "%2f\n",fermi_2[0]); //粒子数（サイト数）
for(int i=0;i<=P;i++){
for(int j=0;j<N;j++){
fprintf(fp, "%2f %lf\n",k_x[i],Energy[i][j]);
//fprintf(fp, "%2f %lf %lf\n",k_x[i],Energy[i][M-1],Energy[i][M]); //最低バンド
}
}
for(int i=0;i<=P;i++){
fprintf(fp_1, "%2f %lf\n",e[i],dd[i]); //状態密度
}
for(int l=0;l<=P;l++){
fprintf(fp_2,"%2f\n",fermisum_1[l]); //down
fprintf(fp_4,"%2f\n",fermisum_2[l]); //upper
}
fclose(fp);
fclose(fp_1);
fclose(fp_2);
fclose(fp_3);
fclose(fp_4);
fclose(fp_5);
return 0;
}
自分で試したこと
fprintf(fp_4,"%2f\n",fermisum_2[l])のかたちを変えたり、計算結果と見比べたりなど
0 likes
