Q&A
開発環境、コンパイラおよびそのバージョンを教えて下さい。
前提
状態密度を計算したい。
実現したいこと
そのために、eとk_xの変化ごとにsum_densityを計算して、最後にeごとにsum_dosの値をファイルに書きたい。
発生している問題・エラーメッセージ
一旦、k_x,eの変化ごとにsum_densityを表示しているが
-4.618141 -3.500000 1.636411 -4.618141 -3.500000 1.636411 -4.618141 -3.500000 1.636411 -4.618141 -3.500000 1.636411 -4.618141 -3.500000 1.636411 -4.618141 -3.500000 1.636411 -4.618141 -3.500000 1.636411 -4.523893 -3.500000 1.491209 -4.523893 -3.500000 1.491209 -4.523893 -3.500000 1.491209 -4.523893 -3.500000 1.491209 -4.523893 -3.500000 1.491209 -4.523893 -3.500000 1.491209 -4.523893 -3.500000 1.491209 -4.523893 -3.500000 1.491209
のように真ん中のeに変化がない。ここの原因が知りたい。
またこのように計算すればすべてのk_xについてsum_dosが最終的にさせているのかも知りたい。
該当のソースコード
C
1#include <stdio.h> 2#include <math.h> 3#include <complex.h> 4#include <stdlib.h> 5 6#define M 5 //単位法の数 7#define N 10 //次数設定 8#define L 20 9 10#define EPS 0.000001 //収束範囲 11 12int n = 100; //100000 //電子数＝2N 13int c = 1; 14double 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; //方向 15double b_pi; 16double a[L][L]; 17double D[N][N]; 18double complex C[N][N]; 19double A[L], B_1[N][N], B_2[N][N]; 20double E[L]; 21 22 23 24 25int main(int argc, char** argv) { 26 27 b_pi = M_PI; 28 double complex a_cube = 1;//= 0.246 * pow(10, -9); //格子定数[nm] 29 double k_x = -1.5*b_pi / a_cube; //-2*b_pi/a_cube/pow(3,0.5); //↓のk_xの値を変えるのも忘れずに！ K+,K-で値が変わるよ！ 30 double k_y = 0; 31 double e=-3.5; 32 double sum_dos[n]; 33 34 FILE* fp; 35 36 if ((fp = fopen("aaa N=5.txt", "w")) == NULL) 37 { 38 printf("Cannot open the file\n"); 39 exit(1); 40 } 41 42 for(int h=0;h<n;h++){ 43 double sum_density[n]; 44 e=e+7/n; 45 k_x = k_x + b_pi * 3/ a_cube / n; 46 47 for (int r = 0;r < n; r++) { 48 49 t_4 = cexp(-1 * I * (-1 * a_cube / pow(3, 0.5) * k_y)); 50 t_5 = cexp(-1 * I * (-a_cube / 2 * k_x - a_cube / 2 / pow(3, 0.5) * k_y)); //cos(a_cube / 2 / pow(3, 0.5) * k_x + a_cube / 2 * k_y) + I * sin((a_cube / 2 / pow(3, 0.5) * k_x + a_cube / 2 * k_y)); //cexp(I * (a_cube / 2 / pow(3, 0.5) * k_x + a_cube / 2 * k_y)); //expの中の符号を計算したあらかじめ 51 t_6 = cexp(-1 * I * (a_cube / 2 * k_x - a_cube / 2 / pow(3, 0.5) * k_y)); //cos(a_cube / 2 / pow(3, 0.5) * k_x - a_cube / 2 * k_y) + I * sin(a_cube / 2 / pow(3, 0.5) * k_x - a_cube / 2 * k_y); //cexp(I * (a_cube / 2 / pow(3, 0.5) * k_x - a_cube / 2 * k_y)); // 52 t_1 = cexp(I * a_cube / pow(3, 0.5) * k_y); 53 t_2 = cexp(I * (-a_cube / 2 * k_x - a_cube / 2 / pow(3, 0.5) * k_y)); 54 t_3 = cexp(I * (a_cube / 2 * k_x - a_cube / 2 / pow(3, 0.5) * k_y)); 55 t_7 = -1 * t_4; 56 t_8 = -1 * t_1; 57 T_1 = t_5 + t_6; 58 T_2 = t_2 + t_3; 59 T_3 = -1 * T_1; 60 T_4 = -1 * T_2; 61 62 for (int k = 0;k < N;k++) { 63 for (int h = 0;h < N;h++) { 64 for (int i = 0;i < M - 2;i++) { 65 if (h < M) { 66 if (k < M) { 67 C[h][k] = 0; 68 } 69 else { 70 C[0][M] = T_1; 71 C[M - 1][N - 1] = T_1; 72 C[M - 1][N - 2] = t_4; 73 C[i + 1][M + i] = t_4; 74 C[i + 1][M + i + 1] = T_1; 75 } 76 } 77 else { 78 if (k < M) { 79 C[M][0] = T_2; 80 C[M][1] = t_1; 81 C[N - 1][M - 1] = T_2; 82 C[M + i + 1][i + 1] = T_2; 83 C[M + i + 1][i + 2] = t_1; 84 85 } 86 else { 87 C[h][k] = 0; 88 } 89 } 90 } 91 } 92 } 93 94 int k, h; 95 for (k = 0;k < N;k++) { 96 for (h = 0;h < N;h++) { 97 D[h][k] = creal(C[h][k]); 98 B_1[h][k] = -1 * cimag(C[h][k]); 99 B_2[h][k] = cimag(C[h][k]); 100 } 101 } 102 103 104 for (k = 0;k < L;k++) { 105 for (h = 0;h < L;h++) { 106 if (h < N) { 107 if (k < N) { 108 a[h][k] = D[h % N][k % N]; 109 } 110 else { 111 a[h][k] = B_1[h % N][k % N]; 112 } 113 } 114 else { 115 if (k < N) { 116 a[h][k] = B_2[h % N][k % N]; 117 } 118 else { 119 a[h][k] = D[h % N][k % N]; 120 } 121 } 122 } 123 } 124 125 while (1) { 126 //double u[L][L]; //単位行列 127 double alpha, beta, gamma; 128 double s, c, w; 129 double wa, wb, wc; 130 double max; 131 int i, j, p, q, x, y; 132 133 //最大要素の行と列を検索 134 max = 0.0; 135 for (i = 0;i < L - 1;i++) { 136 for (j = i + 1;j < L;j++) 137 if (fabs(a[i][j]) > max) { 138 p = i; 139 q = j; 140 max = fabs(a[i][j]); 141 } 142 } 143 144 //収束したら解答打ち出し 145 if (max < EPS) break; 146 147 //sin, cos 計算 148 wa = a[p][p]; 149 wb = a[p][q]; 150 wc = a[q][q]; 151 alpha = -wb; 152 beta = 0.5 * (wa - wc); 153 gamma = fabs(beta) / sqrt(alpha * alpha + beta * beta); 154 s = sqrt(0.5 * (1.0 - gamma)); 155 if (alpha * beta < 0) s = -s; 156 c = sqrt(1.0 - s * s); 157 //直行変換 158 for (j = 0;j < L;j++) { 159 w = a[p][j] * c - a[q][j] * s; 160 a[q][j] = a[p][j] * s + a[q][j] * c; 161 a[p][j] = w; 162 } 163 164 for (j = 0;j < L;j++) { 165 a[j][p] = a[p][j]; 166 a[j][q] = a[q][j]; 167 } 168 169 170 171 w = 2.0 * wb * s * c; 172 a[p][p] = wa * c * c + wc * s * s - w; 173 a[q][q] = wa * s * s + wc * c * c + w; 174 a[p][q] = 0; 175 a[q][p] = 0; 176 177 for (int i = 0;i < L;i++) 178 { 179 A[i] = a[i][i]; 180 } 181 182 } 183 184 double tmp; 185 for (int i=0; i<L; i++) { 186 for (int j=i+1; j<L ;j++) { 187 if (A[i] > A[j]) { 188 tmp = A[i]; 189 A[i] = A[j]; 190 A[j] = tmp; 191 } 192 } 193 } 194 double ave, squ, v, s; 195 int data_num; 196 double sum = 0.0, sumv = 0.0; 197 198 ave = 0.0; 199 squ = 0.0; 200 v = 0.0; 201 s = 0.0; 202 data_num = L; 203 sum = 0.0; 204 sumv = 0.0; 205 for(int i = 0; i < data_num; i++){ 206 sum += A[i]; // データの総和 207 sumv += A[i] * A[i]; // 二乗和 208 } 209 210 // 平均・平方和・分散・標準偏差を求める 211 ave = sum / data_num; // 平均 212 squ = sumv - (sum * sum / data_num); // 平方和 213 v = squ / (data_num - 1); // 分散 214 s = sqrt(v); // 標準偏差 215 216 sum_density[r]=0.0; 217 for (int i=0;i<L;i++){ 218 E[i]=1/b_pi*exp(-pow((e-A[i]),2)/2/v); 219 sum_density[r] += E[i]; 220 } 221 222 //printf("%7.4lf\n",sum_density); 223 224 /*for (int i = 0;i < L;i++) { 225 fprintf(fp, "%2f %lf\n", k_x, A[i]); 226 }*/ 227 228 229 //fprintf(fp, "%2f %lf %lf\n",k_x,e,sum_density[r]); 230 231 232 } 233 234 for(int i=0;i<n;i++){ 235 sum_dos[h] +=sum_density[i]; 236 } 237 238 239 //printf("%lf\n",sum_dos); 240 fprintf(fp, "%2f %lf %lf\n",e,sum_density[h],sum_dos[h]); 241 } 242 243 fclose(fp); 244 return 0; 245} 246
試したこと
e+=7/n
e=e+7/n
と書いても変わりなかった。
補足情報（FW/ツールのバージョンなど）
ここにより詳細な情報を記載してください。
7/nがどういう値になるか理解できないってことですか。
回答1件
0
ベストアンサー
e=e+7/n;
この部分式 7 / n では '7' も 'n' も整数なので、整数除算が行われ結果は 0 です。
e = e + 7.0 / n; のように書いて、double の除算が行われるようにしましょう。
投稿2022/11/08 01:21
編集2022/11/08 01:22
総合スコア20845
ありがとうございます。
もう一つ質問なのですが、k_xは問題なく増えていくのですがそれはb_piがdouble型だからでしょうか。ここもnがintで3も整数値として入っていると思うのですが何が違うのでしょう。
> k_x = k_x + b_pi * 3/ a_cube / n;
b_pi * 3 は b_pi が double なので結果 double、
その / a_cube はもちろん double、
その / n は割られる数が double なので結果 double です。
ありがとうございます。変数の定義に立ち戻ります。
あなたの回答
tips
プレビュー
