| Revision | b6fff3bfa7f0384dc31131ff086df2ffdfcd9ed4 (tree) |
|---|---|
| Zeit | 2018-02-15 13:59:43 |
| Autor |  mitsuaki1987 <kawamitsuaki@gmai...> |
| Commiter | mitsuaki1987 |
Backup
src/calc_nodeline.c (diff) src/fermi_patch.c (diff) src/free_patch.c (diff) src/initialize.c (diff) src/kumo.c (diff) src/read_file.c (diff)
src/variable.h (diff)| @@ -66,15 +66,16 @@ void calc_nodeline() { | ||
| 66 | 66 | #pragma omp for |
| 67 | 67 | for (itri = 0; itri < ntri[ib]; ++itri) { |
| 68 | 68 | |
| 69 | - eigsort(3, matp[ib][itri], sw); | |
| 69 | + eigsort(3, matp[ib][itri][0], sw); | |
| 70 | 70 | for (i = 0; i < 3; ++i) { |
| 71 | 71 | for (j = 0; j < 3; ++j) { |
| 72 | - a[i][j] = (0.f - matp[ib][itri][sw[j]]) / (matp[ib][itri][sw[i]] - matp[ib][itri][sw[j]]); | |
| 72 | + a[i][j] = (0.f - matp[ib][itri][sw[j]][0]) | |
| 73 | + / (matp[ib][itri][sw[i]][0] - matp[ib][itri][sw[j]][0]); | |
| 73 | 74 | }/*for (j = 0; j < 3; ++j)*/ |
| 74 | 75 | }/*for (i = 0; i < 3; ++i)*/ |
| 75 | 76 | |
| 76 | - if ((matp[ib][itri][sw[0]] < 0.0 && 0.0 <= matp[ib][itri][sw[1]]) | |
| 77 | - || (matp[ib][itri][sw[0]] <= 0.0 && 0.0 < matp[ib][itri][sw[1]])) { | |
| 77 | + if (( matp[ib][itri][sw[0]][0] < 0.0 && 0.0 <= matp[ib][itri][sw[1]][0]) | |
| 78 | + || (matp[ib][itri][sw[0]][0] <= 0.0 && 0.0 < matp[ib][itri][sw[1]][0])) { | |
| 78 | 79 | if (query == 0) { |
| 79 | 80 | for (i = 0; i < 3; ++i) { |
| 80 | 81 | kvnl[ib][nnl0][0][i] = kvp[ib][itri][sw[0]][i] * a[0][1] + kvp[ib][itri][sw[1]][i] * a[1][0]; |
| @@ -83,8 +84,8 @@ void calc_nodeline() { | ||
| 83 | 84 | }/*if (query == 0)*/ |
| 84 | 85 | nnl0 += 1; |
| 85 | 86 | }/*else if (matp[ib][itri][sw[0]] < 0.0 && 0.0 <= matp[ib][itri][sw[1]])*/ |
| 86 | - else if ((matp[ib][itri][sw[1]] < 0.0 && 0.0 <= matp[ib][itri][sw[2]]) | |
| 87 | - || (matp[ib][itri][sw[1]] <= 0.0 && 0.0 < matp[ib][itri][sw[2]])) { | |
| 87 | + else if ((matp[ib][itri][sw[1]][0] < 0.0 && 0.0 <= matp[ib][itri][sw[2]][0]) | |
| 88 | + || (matp[ib][itri][sw[1]][0] <= 0.0 && 0.0 < matp[ib][itri][sw[2]][0])) { | |
| 88 | 89 | if (query == 0) { |
| 89 | 90 | for (i = 0; i < 3; ++i) { |
| 90 | 91 | kvnl[ib][nnl0][0][i] = kvp[ib][itri][sw[0]][i] * a[0][2] + kvp[ib][itri][sw[2]][i] * a[2][0]; |
| @@ -81,13 +81,13 @@ static void triangle( | ||
| 81 | 81 | int ib, //!<[in] The band index |
| 82 | 82 | int *ntri0, //!<[inout] Index of triangle patch |
| 83 | 83 | int nbr, //!<[in] Bragg plane |
| 84 | - GLfloat mat1[3], //!<[in] The matrix element | |
| 84 | + GLfloat mat1[3][3], //!<[in] The matrix element | |
| 85 | 85 | GLfloat kvec1[3][3], //!<[in] @f$k@f$-vector of corners |
| 86 | 86 | GLfloat vf1[3][3] //!<[in] @f$v_f@f$-vector of corners |
| 87 | 87 | ) |
| 88 | 88 | { |
| 89 | 89 | int ibr, i, j, sw[3]; |
| 90 | - GLfloat prod[3], thr, thr2 = 0.001f, mat2[3], kvec2[3][3], | |
| 90 | + GLfloat prod[3], thr, thr2 = 0.001f, mat2[3][3], kvec2[3][3], | |
| 91 | 91 | vf2[3][3], a[3][3], bshift, vfave[3], norm[3]; |
| 92 | 92 | /* |
| 93 | 93 | If the area is nearly 0, it is ignored. |
| @@ -130,9 +130,9 @@ static void triangle( | ||
| 130 | 130 | All corners are outside of the Bragg plane |
| 131 | 131 | */ |
| 132 | 132 | for (i = 0; i < 3; ++i) { |
| 133 | - mat2[i] = mat1[sw[i]]; | |
| 134 | 133 | for (j = 0; j < 3; ++j) { |
| 135 | 134 | kvec2[i][j] = kvec1[sw[i]][j] + bshift * bragg[ibr][j]; |
| 135 | + mat2[i][j] = mat1[sw[i]][j]; | |
| 136 | 136 | vf2[i][j] = vf1[sw[i]][j]; |
| 137 | 137 | } |
| 138 | 138 | } |
| @@ -143,28 +143,30 @@ static void triangle( | ||
| 143 | 143 | /* |
| 144 | 144 | Single corner (#0) is inside of the Bragg plane |
| 145 | 145 | */ |
| 146 | - mat2[0] = mat1[sw[0]]; | |
| 147 | - mat2[1] = mat1[sw[0]] * a[0][1] + mat1[sw[1]] * a[1][0]; | |
| 148 | - mat2[2] = mat1[sw[0]] * a[0][2] + mat1[sw[2]] * a[2][0]; | |
| 149 | - for (i = 0; i < 3; ++i) { | |
| 146 | + for (i = 0; i < 3; ++i) { | |
| 150 | 147 | kvec2[0][i] = kvec1[sw[0]][i]; |
| 151 | 148 | kvec2[1][i] = kvec1[sw[0]][i] * a[0][1] + kvec1[sw[1]][i] * a[1][0]; |
| 152 | 149 | kvec2[2][i] = kvec1[sw[0]][i] * a[0][2] + kvec1[sw[2]][i] * a[2][0]; |
| 153 | 150 | |
| 151 | + mat2[0][i] = mat1[sw[0]][i]; | |
| 152 | + mat2[1][i] = mat1[sw[0]][i] * a[0][1] + mat1[sw[1]][i] * a[1][0]; | |
| 153 | + mat2[2][i] = mat1[sw[0]][i] * a[0][2] + mat1[sw[2]][i] * a[2][0]; | |
| 154 | + | |
| 154 | 155 | vf2[0][i] = vf1[sw[0]][i]; |
| 155 | 156 | vf2[1][i] = vf1[sw[0]][i] * a[0][1] + vf1[sw[1]][i] * a[1][0]; |
| 156 | 157 | vf2[2][i] = vf1[sw[0]][i] * a[0][2] + vf1[sw[2]][i] * a[2][0]; |
| 157 | 158 | }/*for (i = 0; i < 3; ++i)*/ |
| 158 | 159 | triangle(ib, ntri0, ibr + 1, mat2, kvec2, vf2); |
| 159 | 160 | |
| 160 | - mat2[0] = mat1[sw[0]] * a[0][1] + mat1[sw[1]] * a[1][0]; | |
| 161 | - mat2[1] = mat1[sw[1]]; | |
| 162 | - mat2[2] = mat1[sw[0]] * a[0][2] + mat1[sw[2]] * a[2][0]; | |
| 163 | 161 | for (i = 0; i < 3; ++i) { |
| 164 | 162 | kvec2[0][i] = kvec1[sw[0]][i] * a[0][1] + kvec1[sw[1]][i] * a[1][0]; |
| 165 | 163 | kvec2[1][i] = kvec1[sw[1]][i]; |
| 166 | 164 | kvec2[2][i] = kvec1[sw[0]][i] * a[0][2] + kvec1[sw[2]][i] * a[2][0]; |
| 167 | 165 | |
| 166 | + mat2[0][i] = mat1[sw[0]][i] * a[0][1] + mat1[sw[1]][i] * a[1][0]; | |
| 167 | + mat2[1][i] = mat1[sw[1]][i]; | |
| 168 | + mat2[2][i] = mat1[sw[0]][i] * a[0][2] + mat1[sw[2]][i] * a[2][0]; | |
| 169 | + | |
| 168 | 170 | vf2[0][i] = vf1[sw[0]][i] * a[0][1] + vf1[sw[1]][i] * a[1][0]; |
| 169 | 171 | vf2[1][i] = vf1[sw[1]][i]; |
| 170 | 172 | vf2[2][i] = vf1[sw[0]][i] * a[0][2] + vf1[sw[2]][i] * a[2][0]; |
| @@ -173,14 +175,15 @@ static void triangle( | ||
| 173 | 175 | kvec2[i][j] += bshift * bragg[ibr][j]; |
| 174 | 176 | triangle(ib, ntri0, ibr + 1, mat2, kvec2, vf2); |
| 175 | 177 | |
| 176 | - mat2[0] = mat1[sw[2]]; | |
| 177 | - mat2[1] = mat1[sw[1]]; | |
| 178 | - mat2[2] = mat1[sw[0]] * a[0][2] + mat1[sw[2]] * a[2][0]; | |
| 179 | 178 | for (i = 0; i < 3; ++i) { |
| 180 | 179 | kvec2[0][i] = kvec1[sw[2]][i]; |
| 181 | 180 | kvec2[1][i] = kvec1[sw[1]][i]; |
| 182 | 181 | kvec2[2][i] = kvec1[sw[0]][i] * a[0][2] + kvec1[sw[2]][i] * a[2][0]; |
| 183 | 182 | |
| 183 | + mat2[0][i] = mat1[sw[2]][i]; | |
| 184 | + mat2[1][i] = mat1[sw[1]][i]; | |
| 185 | + mat2[2][i] = mat1[sw[0]][i] * a[0][2] + mat1[sw[2]][i] * a[2][0]; | |
| 186 | + | |
| 184 | 187 | vf2[0][i] = vf1[sw[2]][i]; |
| 185 | 188 | vf2[1][i] = vf1[sw[1]][i]; |
| 186 | 189 | vf2[2][i] = vf1[sw[0]][i] * a[0][2] + vf1[sw[2]][i] * a[2][0]; |
| @@ -194,14 +197,15 @@ static void triangle( | ||
| 194 | 197 | /* |
| 195 | 198 | Two corners (#0, #1) are inside of the Bragg plane |
| 196 | 199 | */ |
| 197 | - mat2[0] = mat1[sw[0]] * a[0][2] + mat1[sw[2]] * a[2][0]; | |
| 198 | - mat2[1] = mat1[sw[1]] * a[1][2] + mat1[sw[2]] * a[2][1]; | |
| 199 | - mat2[2] = mat1[sw[2]]; | |
| 200 | 200 | for (i = 0; i < 3; ++i) { |
| 201 | 201 | kvec2[0][i] = kvec1[sw[0]][i] * a[0][2] + kvec1[sw[2]][i] * a[2][0]; |
| 202 | 202 | kvec2[1][i] = kvec1[sw[1]][i] * a[1][2] + kvec1[sw[2]][i] * a[2][1]; |
| 203 | 203 | kvec2[2][i] = kvec1[sw[2]][i]; |
| 204 | 204 | |
| 205 | + mat2[0][i] = mat1[sw[0]][i] * a[0][2] + mat1[sw[2]][i] * a[2][0]; | |
| 206 | + mat2[1][i] = mat1[sw[1]][i] * a[1][2] + mat1[sw[2]][i] * a[2][1]; | |
| 207 | + mat2[2][i] = mat1[sw[2]][i]; | |
| 208 | + | |
| 205 | 209 | vf2[0][i] = vf1[sw[0]][i] * a[0][2] + vf1[sw[2]][i] * a[2][0]; |
| 206 | 210 | vf2[1][i] = vf1[sw[1]][i] * a[1][2] + vf1[sw[2]][i] * a[2][1]; |
| 207 | 211 | vf2[2][i] = vf1[sw[2]][i]; |
| @@ -210,28 +214,30 @@ static void triangle( | ||
| 210 | 214 | kvec2[i][j] += bshift * bragg[ibr][j]; |
| 211 | 215 | triangle(ib, ntri0, ibr + 1, mat2, kvec2, vf2); |
| 212 | 216 | |
| 213 | - mat2[0] = mat1[sw[0]]; | |
| 214 | - mat2[1] = mat1[sw[1]]; | |
| 215 | - mat2[2] = mat1[sw[0]] * a[0][2] + mat1[sw[2]] * a[2][0]; | |
| 216 | 217 | for (i = 0; i < 3; ++i) { |
| 217 | 218 | kvec2[0][i] = kvec1[sw[0]][i]; |
| 218 | 219 | kvec2[1][i] = kvec1[sw[1]][i]; |
| 219 | 220 | kvec2[2][i] = kvec1[sw[0]][i] * a[0][2] + kvec1[sw[2]][i] * a[2][0]; |
| 220 | 221 | |
| 222 | + mat2[0][i] = mat1[sw[0]][i]; | |
| 223 | + mat2[1][i] = mat1[sw[1]][i]; | |
| 224 | + mat2[2][i] = mat1[sw[0]][i] * a[0][2] + mat1[sw[2]][i] * a[2][0]; | |
| 225 | + | |
| 221 | 226 | vf2[0][i] = vf1[sw[0]][i]; |
| 222 | 227 | vf2[1][i] = vf1[sw[1]][i]; |
| 223 | 228 | vf2[2][i] = vf1[sw[0]][i] * a[0][2] + vf1[sw[2]][i] * a[2][0]; |
| 224 | 229 | }/*for (i = 0; i < 3; ++i)*/ |
| 225 | 230 | triangle(ib, ntri0, ibr + 1, mat2, kvec2, vf2); |
| 226 | 231 | |
| 227 | - mat2[0] = mat1[sw[1]] * a[1][2] + mat1[sw[2]] * a[2][1]; | |
| 228 | - mat2[1] = mat1[sw[1]]; | |
| 229 | - mat2[2] = mat1[sw[0]] * a[0][2] + mat1[sw[2]] * a[2][0]; | |
| 230 | 232 | for (i = 0; i < 3; ++i) { |
| 231 | 233 | kvec2[0][i] = kvec1[sw[1]][i] * a[1][2] + kvec1[sw[2]][i] * a[2][1]; |
| 232 | 234 | kvec2[1][i] = kvec1[sw[1]][i]; |
| 233 | 235 | kvec2[2][i] = kvec1[sw[0]][i] * a[0][2] + kvec1[sw[2]][i] * a[2][0]; |
| 234 | 236 | |
| 237 | + mat2[0][i] = mat1[sw[1]][i] * a[1][2] + mat1[sw[2]][i] * a[2][1]; | |
| 238 | + mat2[1][i] = mat1[sw[1]][i]; | |
| 239 | + mat2[2][i] = mat1[sw[0]][i] * a[0][2] + mat1[sw[2]][i] * a[2][0]; | |
| 240 | + | |
| 235 | 241 | vf2[0][i] = vf1[sw[1]][i] * a[1][2] + vf1[sw[2]][i] * a[2][1]; |
| 236 | 242 | vf2[1][i] = vf1[sw[1]][i]; |
| 237 | 243 | vf2[2][i] = vf1[sw[0]][i] * a[0][2] + vf1[sw[2]][i] * a[2][0]; |
| @@ -261,18 +267,18 @@ static void triangle( | ||
| 261 | 267 | |
| 262 | 268 | if (prod[0] < 0.0f) { |
| 263 | 269 | for (i = 0; i < 3; ++i) { |
| 264 | - matp[ib][*ntri0][i] = mat1[2-i]; | |
| 265 | 270 | for (j = 0; j < 3; ++j) { |
| 266 | 271 | kvp[ib][*ntri0][i][j] = kvec1[2-i][j]; |
| 272 | + matp[ib][*ntri0][i][j] = mat1[2 - i][j]; | |
| 267 | 273 | nmlp[ib][*ntri0][i][j] = vf1[2-i][j]; |
| 268 | 274 | } |
| 269 | 275 | }/*for (i = 0; i < 3; ++i)*/ |
| 270 | 276 | } |
| 271 | 277 | else { |
| 272 | 278 | for (i = 0; i < 3; ++i) { |
| 273 | - matp[ib][*ntri0][i] = mat1[i]; | |
| 274 | 279 | for (j = 0; j < 3; ++j) { |
| 275 | 280 | kvp[ib][*ntri0][i][j] = kvec1[i][j]; |
| 281 | + matp[ib][*ntri0][i][j] = mat1[i][j]; | |
| 276 | 282 | nmlp[ib][*ntri0][i][j] = vf1[i][j]; |
| 277 | 283 | } |
| 278 | 284 | }/*for (i = 0; i < 3; ++i)*/ |
| @@ -308,14 +314,14 @@ static void tetrahedron( | ||
| 308 | 314 | int ib, //!< [in] The band index |
| 309 | 315 | int *ntri0, //!< [in,out] Counter of trixngle |
| 310 | 316 | GLfloat eig1[8], //!< [in] Orbital energies @f$\varepsilon_{n k}@f$ |
| 311 | - GLfloat mat1[8], //!< [in] Matrix elements @f$\Delta_{n k}@f$ | |
| 317 | + GLfloat mat1[8][3], //!< [in] Matrix elements @f$\Delta_{n k}@f$ | |
| 312 | 318 | GLfloat kvec1[8][3], //!< [in] @f$k@f$-vectors |
| 313 | 319 | GLfloat vf1[8][3] //!< [in] @f$v_f@f$-vectors |
| 314 | 320 | ) |
| 315 | 321 | { |
| 316 | 322 | int it, i, j, sw[4]; |
| 317 | - GLfloat eig2[4], mat2[4], kvec2[4][3], vf2[4][3], a[4][4], | |
| 318 | - kvec3[3][3], mat3[3], vf3[3][3], vol, thr = 0.000f; | |
| 323 | + GLfloat eig2[4], mat2[4][3], kvec2[4][3], vf2[4][3], a[4][4], | |
| 324 | + kvec3[3][3], mat3[3][3], vf3[3][3], vol, thr = 0.000f; | |
| 319 | 325 | |
| 320 | 326 | for (it = 0; it < 6; ++it) { |
| 321 | 327 | /* |
| @@ -323,8 +329,10 @@ static void tetrahedron( | ||
| 323 | 329 | */ |
| 324 | 330 | for (i = 0; i < 4; ++i) { |
| 325 | 331 | eig2[i] = eig1[corner[it][i]]; |
| 326 | - mat2[i] = mat1[corner[it][i]]; | |
| 327 | - for (j = 0; j < 3; ++j) vf2[i][j] = vf1[corner[it][i]][j]; | |
| 332 | + for (j = 0; j < 3; ++j) { | |
| 333 | + mat2[i][j] = mat1[corner[it][i]][j]; | |
| 334 | + vf2[i][j] = vf1[corner[it][i]][j]; | |
| 335 | + } | |
| 328 | 336 | /* |
| 329 | 337 | Fractional -> Cartecian |
| 330 | 338 | */ |
| @@ -349,14 +357,15 @@ static void tetrahedron( | ||
| 349 | 357 | kvec3[1][i] = kvec2[sw[0]][i] * a[0][2] + kvec2[sw[2]][i] * a[2][0]; |
| 350 | 358 | kvec3[2][i] = kvec2[sw[0]][i] * a[0][3] + kvec2[sw[3]][i] * a[3][0]; |
| 351 | 359 | |
| 360 | + mat3[0][i] = mat2[sw[0]][i] * a[0][1] + mat2[sw[1]][i] * a[1][0]; | |
| 361 | + mat3[1][i] = mat2[sw[0]][i] * a[0][2] + mat2[sw[2]][i] * a[2][0]; | |
| 362 | + mat3[2][i] = mat2[sw[0]][i] * a[0][3] + mat2[sw[3]][i] * a[3][0]; | |
| 363 | + | |
| 352 | 364 | vf3[0][i] = vf2[sw[0]][i] * a[0][1] + vf2[sw[1]][i] * a[1][0]; |
| 353 | 365 | vf3[1][i] = vf2[sw[0]][i] * a[0][2] + vf2[sw[2]][i] * a[2][0]; |
| 354 | 366 | vf3[2][i] = vf2[sw[0]][i] * a[0][3] + vf2[sw[3]][i] * a[3][0]; |
| 355 | 367 | } |
| 356 | - mat3[0] = mat2[sw[0]] * a[0][1] + mat2[sw[1]] * a[1][0]; | |
| 357 | - mat3[1] = mat2[sw[0]] * a[0][2] + mat2[sw[2]] * a[2][0]; | |
| 358 | - mat3[2] = mat2[sw[0]] * a[0][3] + mat2[sw[3]] * a[3][0]; | |
| 359 | - | |
| 368 | + | |
| 360 | 369 | vol = a[1][0] * a[2][0] * a[3][0]; |
| 361 | 370 | triangle(ib, ntri0, 0, mat3, kvec3, vf3); |
| 362 | 371 | } |
| @@ -366,14 +375,15 @@ static void tetrahedron( | ||
| 366 | 375 | kvec3[1][i] = kvec2[sw[0]][i] * a[0][3] + kvec2[sw[3]][i] * a[3][0]; |
| 367 | 376 | kvec3[2][i] = kvec2[sw[1]][i] * a[1][2] + kvec2[sw[2]][i] * a[2][1]; |
| 368 | 377 | |
| 378 | + mat3[0][i] = mat2[sw[0]][i] * a[0][2] + mat2[sw[2]][i] * a[2][0]; | |
| 379 | + mat3[1][i] = mat2[sw[0]][i] * a[0][3] + mat2[sw[3]][i] * a[3][0]; | |
| 380 | + mat3[2][i] = mat2[sw[1]][i] * a[1][2] + mat2[sw[2]][i] * a[2][1]; | |
| 381 | + | |
| 369 | 382 | vf3[0][i] = vf2[sw[0]][i] * a[0][2] + vf2[sw[2]][i] * a[2][0]; |
| 370 | 383 | vf3[1][i] = vf2[sw[0]][i] * a[0][3] + vf2[sw[3]][i] * a[3][0]; |
| 371 | 384 | vf3[2][i] = vf2[sw[1]][i] * a[1][2] + vf2[sw[2]][i] * a[2][1]; |
| 372 | 385 | } |
| 373 | - mat3[0] = mat2[sw[0]] * a[0][2] + mat2[sw[2]] * a[2][0]; | |
| 374 | - mat3[1] = mat2[sw[0]] * a[0][3] + mat2[sw[3]] * a[3][0]; | |
| 375 | - mat3[2] = mat2[sw[1]] * a[1][2] + mat2[sw[2]] * a[2][1]; | |
| 376 | - | |
| 386 | + | |
| 377 | 387 | vol = a[1][2] * a[2][0] * a[3][0]; |
| 378 | 388 | triangle(ib, ntri0, 0, mat3, kvec3, vf3); |
| 379 | 389 | /**/ |
| @@ -382,13 +392,14 @@ static void tetrahedron( | ||
| 382 | 392 | kvec3[1][i] = kvec2[sw[0]][i] * a[0][3] + kvec2[sw[3]][i] * a[3][0]; |
| 383 | 393 | kvec3[2][i] = kvec2[sw[1]][i] * a[1][2] + kvec2[sw[2]][i] * a[2][1]; |
| 384 | 394 | |
| 395 | + mat3[0][i] = mat2[sw[1]][i] * a[1][3] + mat2[sw[3]][i] * a[3][1]; | |
| 396 | + mat3[1][i] = mat2[sw[0]][i] * a[0][3] + mat2[sw[3]][i] * a[3][0]; | |
| 397 | + mat3[2][i] = mat2[sw[1]][i] * a[1][2] + mat2[sw[2]][i] * a[2][1]; | |
| 398 | + | |
| 385 | 399 | vf3[0][i] = vf2[sw[1]][i] * a[1][3] + vf2[sw[3]][i] * a[3][1]; |
| 386 | 400 | vf3[1][i] = vf2[sw[0]][i] * a[0][3] + vf2[sw[3]][i] * a[3][0]; |
| 387 | 401 | vf3[2][i] = vf2[sw[1]][i] * a[1][2] + vf2[sw[2]][i] * a[2][1]; |
| 388 | 402 | } |
| 389 | - mat3[0] = mat2[sw[1]] * a[1][3] + mat2[sw[3]] * a[3][1]; | |
| 390 | - mat3[1] = mat2[sw[0]] * a[0][3] + mat2[sw[3]] * a[3][0]; | |
| 391 | - mat3[2] = mat2[sw[1]] * a[1][2] + mat2[sw[2]] * a[2][1]; | |
| 392 | 403 | |
| 393 | 404 | vol = a[1][3] * a[3][0] * a[2][1]; |
| 394 | 405 | triangle(ib, ntri0, 0, mat3, kvec3, vf3); |
| @@ -399,13 +410,14 @@ static void tetrahedron( | ||
| 399 | 410 | kvec3[1][i] = kvec2[sw[3]][i] * a[3][1] + kvec2[sw[1]][i] * a[1][3]; |
| 400 | 411 | kvec3[2][i] = kvec2[sw[3]][i] * a[3][2] + kvec2[sw[2]][i] * a[2][3]; |
| 401 | 412 | |
| 413 | + mat3[0][i] = mat2[sw[3]][i] * a[3][0] + mat2[sw[0]][i] * a[0][3]; | |
| 414 | + mat3[1][i] = mat2[sw[3]][i] * a[3][1] + mat2[sw[1]][i] * a[1][3]; | |
| 415 | + mat3[2][i] = mat2[sw[3]][i] * a[3][2] + mat2[sw[2]][i] * a[2][3]; | |
| 416 | + | |
| 402 | 417 | vf3[0][i] = vf2[sw[3]][i] * a[3][0] + vf2[sw[0]][i] * a[0][3]; |
| 403 | 418 | vf3[1][i] = vf2[sw[3]][i] * a[3][1] + vf2[sw[1]][i] * a[1][3]; |
| 404 | 419 | vf3[2][i] = vf2[sw[3]][i] * a[3][2] + vf2[sw[2]][i] * a[2][3]; |
| 405 | 420 | } |
| 406 | - mat3[0] = mat2[sw[3]] * a[3][0] + mat2[sw[0]] * a[0][3]; | |
| 407 | - mat3[1] = mat2[sw[3]] * a[3][1] + mat2[sw[1]] * a[1][3]; | |
| 408 | - mat3[2] = mat2[sw[3]] * a[3][2] + mat2[sw[2]] * a[2][3]; | |
| 409 | 421 | |
| 410 | 422 | vol = a[0][3] * a[1][3] * a[2][3]; |
| 411 | 423 | triangle(ib, ntri0, 0, mat3, kvec3, vf3); |
| @@ -452,7 +464,7 @@ void fermi_patch() | ||
| 452 | 464 | private(ib,j0,i0,i1,ithread) |
| 453 | 465 | { |
| 454 | 466 | int ntri0, i, j, i2, j1, j2, ii0, ii1, ii2; |
| 455 | - GLfloat kvec1[8][3], mat1[8], eig1[8], vf1[8][3]; | |
| 467 | + GLfloat kvec1[8][3], mat1[8][3], eig1[8], vf1[8][3]; | |
| 456 | 468 | |
| 457 | 469 | ithread = get_thread(); |
| 458 | 470 | for (ib = 0; ib < nb; ++ib) { |
| @@ -519,16 +531,16 @@ void fermi_patch() | ||
| 519 | 531 | eig1[6] = eig[ib][ii0][ii1][i2] - EF; |
| 520 | 532 | eig1[7] = eig[ib][ii0][ii1][ii2] - EF; |
| 521 | 533 | /**/ |
| 522 | - mat1[0] = mat[ib][i0][i1][i2]; | |
| 523 | - mat1[1] = mat[ib][i0][i1][ii2]; | |
| 524 | - mat1[2] = mat[ib][i0][ii1][i2]; | |
| 525 | - mat1[3] = mat[ib][i0][ii1][ii2]; | |
| 526 | - mat1[4] = mat[ib][ii0][i1][i2]; | |
| 527 | - mat1[5] = mat[ib][ii0][i1][ii2]; | |
| 528 | - mat1[6] = mat[ib][ii0][ii1][i2]; | |
| 529 | - mat1[7] = mat[ib][ii0][ii1][ii2]; | |
| 530 | - /**/ | |
| 531 | 534 | for (j = 0; j < 3; j++) { |
| 535 | + mat1[0][j] = mat[ib][i0][i1][i2][j]; | |
| 536 | + mat1[1][j] = mat[ib][i0][i1][ii2][j]; | |
| 537 | + mat1[2][j] = mat[ib][i0][ii1][i2][j]; | |
| 538 | + mat1[3][j] = mat[ib][i0][ii1][ii2][j]; | |
| 539 | + mat1[4][j] = mat[ib][ii0][i1][i2][j]; | |
| 540 | + mat1[5][j] = mat[ib][ii0][i1][ii2][j]; | |
| 541 | + mat1[6][j] = mat[ib][ii0][ii1][i2][j]; | |
| 542 | + mat1[7][j] = mat[ib][ii0][ii1][ii2][j]; | |
| 543 | + /**/ | |
| 532 | 544 | vf1[0][j] = vf[ib][i0][i1][i2][j]; |
| 533 | 545 | vf1[1][j] = vf[ib][i0][i1][ii2][j]; |
| 534 | 546 | vf1[2][j] = vf[ib][i0][ii1][i2][j]; |
| @@ -571,24 +583,25 @@ void fermi_patch() | ||
| 571 | 583 | /* |
| 572 | 584 | Allocation of triangler patches |
| 573 | 585 | */ |
| 574 | - matp = (GLfloat***)malloc(nb * sizeof(GLfloat**)); | |
| 586 | + matp = (GLfloat****)malloc(nb * sizeof(GLfloat***)); | |
| 575 | 587 | clr = (GLfloat**)malloc(nb * sizeof(GLfloat*)); |
| 576 | 588 | kvp = (GLfloat****)malloc(nb * sizeof(GLfloat***)); |
| 577 | 589 | nmlp = (GLfloat****)malloc(nb * sizeof(GLfloat***)); |
| 578 | 590 | kvp_rot = (GLfloat**)malloc(nb * sizeof(GLfloat*)); |
| 579 | 591 | nmlp_rot = (GLfloat**)malloc(nb * sizeof(GLfloat*)); |
| 580 | 592 | for (ib = 0; ib < nb; ++ib) { |
| 581 | - matp[ib] = (GLfloat**)malloc(ntri[ib] * sizeof(GLfloat*)); | |
| 593 | + matp[ib] = (GLfloat***)malloc(ntri[ib] * sizeof(GLfloat**)); | |
| 582 | 594 | clr[ib] = (GLfloat*)malloc(12 * ntri[ib] * sizeof(GLfloat)); |
| 583 | 595 | kvp[ib] = (GLfloat***)malloc(ntri[ib] * sizeof(GLfloat**)); |
| 584 | 596 | nmlp[ib] = (GLfloat***)malloc(ntri[ib] * sizeof(GLfloat**)); |
| 585 | 597 | kvp_rot[ib] = (GLfloat*)malloc(9 * ntri[ib] * sizeof(GLfloat)); |
| 586 | 598 | nmlp_rot[ib] = (GLfloat*)malloc(9 * ntri[ib] * sizeof(GLfloat)); |
| 587 | 599 | for (i0 = 0; i0 < ntri[ib]; ++i0) { |
| 588 | - matp[ib][i0] = (GLfloat*)malloc(3 * sizeof(GLfloat)); | |
| 600 | + matp[ib][i0] = (GLfloat**)malloc(3 * sizeof(GLfloat*)); | |
| 589 | 601 | kvp[ib][i0] = (GLfloat**)malloc(3 * sizeof(GLfloat*)); |
| 590 | 602 | nmlp[ib][i0] = (GLfloat**)malloc(3 * sizeof(GLfloat*)); |
| 591 | 603 | for (i1 = 0; i1 < 3; ++i1) { |
| 604 | + matp[ib][i0][i1] = (GLfloat*)malloc(3 * sizeof(GLfloat)); | |
| 592 | 605 | kvp[ib][i0][i1] = (GLfloat*)malloc(3 * sizeof(GLfloat)); |
| 593 | 606 | nmlp[ib][i0][i1] = (GLfloat*)malloc(3 * sizeof(GLfloat)); |
| 594 | 607 | }/*for (i1 = 0; i1 < 3; ++i1)*/ |
| @@ -45,6 +45,7 @@ void free_patch() { | ||
| 45 | 45 | for (i0 = 0; i0 < ntri[ib]; ++i0) { |
| 46 | 46 | for (i1 = 0; i1 < 3; ++i1) { |
| 47 | 47 | free(nmlp[ib][i0][i1]); |
| 48 | + free(matp[ib][i0][i1]); | |
| 48 | 49 | free(kvp[ib][i0][i1]); |
| 49 | 50 | } |
| 50 | 51 | free(nmlp[ib][i0]); |
| @@ -123,10 +124,11 @@ void free_patch() { | ||
| 123 | 124 | */ |
| 124 | 125 | void max_and_min() { |
| 125 | 126 | int itri, ierr, ithread; |
| 126 | - GLfloat matmax, matmin, *max_th, *min_th; | |
| 127 | + GLfloat matmax, matmin, absmax, *max_th, *min_th, *abs_th; | |
| 127 | 128 | |
| 128 | 129 | max_th = (GLfloat*)malloc(nthreads * sizeof(GLfloat)); |
| 129 | 130 | min_th = (GLfloat*)malloc(nthreads * sizeof(GLfloat)); |
| 131 | + abs_th = (GLfloat*)malloc(nthreads * sizeof(GLfloat)); | |
| 130 | 132 | |
| 131 | 133 | printf("\n"); |
| 132 | 134 | if (fcscl == 1) printf(" ## Full color scale mode #############\n"); |
| @@ -139,20 +141,25 @@ void max_and_min() { | ||
| 139 | 141 | |
| 140 | 142 | if (fcscl == 1 || fcscl == 2 || fcscl == 7 || fcscl == 8) { |
| 141 | 143 | #pragma omp parallel default(none) \ |
| 142 | -shared(nb,ntri,matp,max_th,min_th) private(itri,ithread) | |
| 144 | +shared(nb,ntri,matp,max_th,min_th,abs_th) private(itri,ithread) | |
| 143 | 145 | { |
| 144 | 146 | int i, ib; |
| 147 | + GLfloat abs; | |
| 145 | 148 | |
| 146 | 149 | ithread = get_thread(); |
| 147 | 150 | max_th[ithread] = -1.0e10f; |
| 148 | 151 | min_th[ithread] = 1.0e10f; |
| 152 | + abs_th[ithread] = -1.0e10f; | |
| 149 | 153 | |
| 150 | 154 | for (ib = 0; ib < nb; ib++) { |
| 151 | 155 | #pragma omp for |
| 152 | 156 | for (itri = 0; itri < ntri[ib]; ++itri) { |
| 153 | 157 | for (i = 0; i < 3; ++i) { |
| 154 | - if (matp[ib][itri][i] > max_th[ithread]) max_th[ithread] = matp[ib][itri][i]; | |
| 155 | - if (matp[ib][itri][i] < min_th[ithread]) min_th[ithread] = matp[ib][itri][i]; | |
| 158 | + abs = sqrtf(matp[ib][itri][i][0] * matp[ib][itri][i][0] | |
| 159 | + + matp[ib][itri][i][1] * matp[ib][itri][i][1]); | |
| 160 | + if (matp[ib][itri][i][0] > max_th[ithread]) max_th[ithread] = matp[ib][itri][i][0]; | |
| 161 | + if (matp[ib][itri][i][0] < min_th[ithread]) min_th[ithread] = matp[ib][itri][i][0]; | |
| 162 | + if (abs > abs_th[ithread]) abs_th[ithread] = abs; | |
| 156 | 163 | } |
| 157 | 164 | }/*for (itri = 0; itri < ntri[ib]; ++itri)*/ |
| 158 | 165 | }/*for (ib = 0; ib < nb; ib++)*/ |
| @@ -160,12 +167,15 @@ shared(nb,ntri,matp,max_th,min_th) private(itri,ithread) | ||
| 160 | 167 | /**/ |
| 161 | 168 | matmax = max_th[0]; |
| 162 | 169 | matmin = min_th[0]; |
| 170 | + absmax = abs_th[0]; | |
| 163 | 171 | for (ithread = 1; ithread < nthreads; ithread++) { |
| 164 | 172 | if (max_th[ithread] > matmax) matmax = max_th[ithread]; |
| 165 | - if (min_th[ithread] > matmin) matmin = min_th[ithread]; | |
| 173 | + if (min_th[ithread] < matmin) matmin = min_th[ithread]; | |
| 174 | + if (abs_th[ithread] > absmax) absmax = abs_th[ithread]; | |
| 166 | 175 | } |
| 167 | 176 | printf(" Max. value : %f\n", matmax); |
| 168 | 177 | printf(" Min. value : %f\n\n", matmin); |
| 178 | + printf(" Max. absolute value : %f\n", absmax); | |
| 169 | 179 | /**/ |
| 170 | 180 | if (fcscl == 2 || fcscl == 8) { |
| 171 | 181 | printf(" Set min. value : "); |
| @@ -234,7 +244,7 @@ private(itri) | ||
| 234 | 244 | for (itri = 0; itri < ntri[ib]; ++itri) { |
| 235 | 245 | for (i = 0; i < 3; ++i) { |
| 236 | 246 | /**/ |
| 237 | - mat2 = (matp[ib][itri][i] - matmin) / (matmax - matmin); | |
| 247 | + mat2 = (matp[ib][itri][i][0] - matmin) / (matmax - matmin); | |
| 238 | 248 | mat2 = mat2 * 4.0f; |
| 239 | 249 | /**/ |
| 240 | 250 | if (mat2 <= 1.0) { |
| @@ -321,7 +331,7 @@ private(itri) | ||
| 321 | 331 | for (itri = 0; itri < ntri[ib]; ++itri) { |
| 322 | 332 | for (i = 0; i < 3; ++i) { |
| 323 | 333 | /**/ |
| 324 | - mat2 = matp[ib][itri][i] / 6.283185307f; | |
| 334 | + mat2 = matp[ib][itri][i][0] / 6.283185307f; | |
| 325 | 335 | mat2 = mat2 - floorf(mat2); |
| 326 | 336 | mat2 = mat2 * 6.0f; |
| 327 | 337 | /**/ |
| @@ -405,7 +415,7 @@ private(itri) | ||
| 405 | 415 | for (itri = 0; itri < ntri[ib]; ++itri) { |
| 406 | 416 | for (i = 0; i < 3; ++i) { |
| 407 | 417 | /**/ |
| 408 | - mat2 = (matp[ib][itri][i] - matmin) / (matmax - matmin); | |
| 418 | + mat2 = (matp[ib][itri][i][0] - matmin) / (matmax - matmin); | |
| 409 | 419 | /**/ |
| 410 | 420 | for (j = 0; j < 4; ++j) clr[ib][j + 4 * i + 12 * itri] = wgray[j] * mat2 + bgray[j] * (1.0f - mat2); |
| 411 | 421 | }/*for (i = 0; i < 3; ++i)*/ |
| @@ -246,15 +246,15 @@ void max_and_min_bz() { | ||
| 246 | 246 | for (ib = 0; ib < nb; ib++) { |
| 247 | 247 | eigmax = eig0[0][0][0][0]; |
| 248 | 248 | eigmin = eig0[0][0][0][0]; |
| 249 | - matmax = mat0[0][0][0][0]; | |
| 250 | - matmin = mat0[0][0][0][0]; | |
| 249 | + matmax = mat0[0][0][0][0][0]; | |
| 250 | + matmin = mat0[0][0][0][0][0]; | |
| 251 | 251 | for (i0 = 0; i0 < ng0[0]; ++i0) { |
| 252 | 252 | for (i1 = 0; i1 < ng0[1]; ++i1) { |
| 253 | 253 | for (i2 = 0; i2 < ng0[2]; ++i2) { |
| 254 | 254 | if (eig0[ib][i0][i1][i2] > eigmax) eigmax = eig0[ib][i0][i1][i2]; |
| 255 | 255 | if (eig0[ib][i0][i1][i2] < eigmin) eigmin = eig0[ib][i0][i1][i2]; |
| 256 | - if (mat0[ib][i0][i1][i2] > matmax) matmax = mat0[ib][i0][i1][i2]; | |
| 257 | - if (mat0[ib][i0][i1][i2] < matmin) matmin = mat0[ib][i0][i1][i2]; | |
| 256 | + if (mat0[ib][i0][i1][i2][0] > matmax) matmax = mat0[ib][i0][i1][i2][0]; | |
| 257 | + if (mat0[ib][i0][i1][i2][0] < matmin) matmin = mat0[ib][i0][i1][i2][0]; | |
| 258 | 258 | }/*for (i2 = 0; i2 < ng0[2]; ++i2)*/ |
| 259 | 259 | }/*for (i1 = 0; i1 < ng0[1]; ++i1)*/ |
| 260 | 260 | }/*for (i0 = 0; i0 < ng0[0]; ++i0)*/ |
| @@ -77,6 +77,7 @@ void interpol_energy() { | ||
| 77 | 77 | for (i1 = 0; i1 < ng[1]; i1++) { |
| 78 | 78 | for (i2 = 0; i2 < ng[2]; i2++) { |
| 79 | 79 | free(vf[ib][i0][i1][i2]); |
| 80 | + free(mat[ib][i0][i1][i2]); | |
| 80 | 81 | } |
| 81 | 82 | free(eig[ib][i0][i1]); |
| 82 | 83 | free(mat[ib][i0][i1]); |
| @@ -94,17 +95,18 @@ void interpol_energy() { | ||
| 94 | 95 | /**/ |
| 95 | 96 | for (ib = 0; ib < nb; ib++) { |
| 96 | 97 | eig[ib] = (GLfloat***)malloc(ng[0] * sizeof(GLfloat**)); |
| 97 | - mat[ib] = (GLfloat***)malloc(ng[0] * sizeof(GLfloat**)); | |
| 98 | + mat[ib] = (GLfloat****)malloc(ng[0] * sizeof(GLfloat***)); | |
| 98 | 99 | vf[ib] = (GLfloat****)malloc(ng[0] * sizeof(GLfloat***)); |
| 99 | 100 | for (i0 = 0; i0 < ng[0]; i0++) { |
| 100 | 101 | eig[ib][i0] = (GLfloat**)malloc(ng[1] * sizeof(GLfloat*)); |
| 101 | - mat[ib][i0] = (GLfloat**)malloc(ng[1] * sizeof(GLfloat*)); | |
| 102 | + mat[ib][i0] = (GLfloat***)malloc(ng[1] * sizeof(GLfloat**)); | |
| 102 | 103 | vf[ib][i0] = (GLfloat***)malloc(ng[1] * sizeof(GLfloat**)); |
| 103 | 104 | for (i1 = 0; i1 < ng[1]; i1++) { |
| 104 | 105 | eig[ib][i0][i1] = (GLfloat*)malloc(ng[2] * sizeof(GLfloat)); |
| 105 | - mat[ib][i0][i1] = (GLfloat*)malloc(ng[2] * sizeof(GLfloat)); | |
| 106 | + mat[ib][i0][i1] = (GLfloat**)malloc(ng[2] * sizeof(GLfloat*)); | |
| 106 | 107 | vf[ib][i0][i1] = (GLfloat**)malloc(ng[2] * sizeof(GLfloat*)); |
| 107 | 108 | for (i2 = 0; i2 < ng[2]; i2++) { |
| 109 | + mat[ib][i0][i1][i2] = (GLfloat*)malloc(3 * sizeof(GLfloat)); | |
| 108 | 110 | vf[ib][i0][i1][i2] = (GLfloat*)malloc(3 * sizeof(GLfloat)); |
| 109 | 111 | } |
| 110 | 112 | }/*for (i1 = 0; i1 < ng[1]; i1++)*/ |
| @@ -117,8 +119,10 @@ void interpol_energy() { | ||
| 117 | 119 | shared(nb,ng0,ng,eig,eig0,mat,mat0,interpol) \ |
| 118 | 120 | private (ib,i0,i1,i2,ii) |
| 119 | 121 | { |
| 120 | - int j0, j1, j2; | |
| 121 | - GLfloat coef[4], mat1[4][4][4], eig1[4][4][4], mat2[4][4], eig2[4][4], mat3[4], eig3[4]; | |
| 122 | + int j0, j1, j2, jj; | |
| 123 | + GLfloat coef[4], | |
| 124 | + mat1[4][4][4][2], mat2[4][4][2], mat3[4][2], | |
| 125 | + eig1[4][4][4], eig2[4][4], eig3[4]; | |
| 122 | 126 | |
| 123 | 127 | for (ib = 0; ib < nb; ib++) { |
| 124 | 128 | # pragma omp for nowait |
| @@ -132,9 +136,11 @@ void interpol_energy() { | ||
| 132 | 136 | eig1[j0][j1][j2] = eig0[ib][modulo(i0 + j0 - 1, ng0[0])] |
| 133 | 137 | [modulo(i1 + j1 - 1, ng0[1])] |
| 134 | 138 | [modulo(i2 + j2 - 1, ng0[2])]; |
| 135 | - mat1[j0][j1][j2] = mat0[ib][modulo(i0 + j0 - 1, ng0[0])] | |
| 136 | - [modulo(i1 + j1 - 1, ng0[1])] | |
| 137 | - [modulo(i2 + j2 - 1, ng0[2])]; | |
| 139 | + for (jj = 0; jj < 2; jj++) { | |
| 140 | + mat1[j0][j1][j2][jj] = mat0[ib][modulo(i0 + j0 - 1, ng0[0])] | |
| 141 | + [modulo(i1 + j1 - 1, ng0[1])] | |
| 142 | + [modulo(i2 + j2 - 1, ng0[2])][jj]; | |
| 143 | + } | |
| 138 | 144 | }/*for (j2 = 0; j2 < 4; j2++)*/ |
| 139 | 145 | }/*for (j1 = 0; j1 < 4; j1++)*/ |
| 140 | 146 | }/*for (i2 = 0; i2 < ng0[2]; i2++)*/ |
| @@ -143,10 +149,11 @@ void interpol_energy() { | ||
| 143 | 149 | for (j1 = 0; j1 < 4; j1++) { |
| 144 | 150 | for (j2 = 0; j2 < 4; j2++) { |
| 145 | 151 | eig2[j1][j2] = 0.0; |
| 146 | - mat2[j1][j2] = 0.0; | |
| 152 | + for (jj = 0; jj < 2; jj++) mat2[j1][j2][jj] = 0.0; | |
| 147 | 153 | for (ii = 0; ii < 4; ii++) { |
| 148 | 154 | eig2[j1][j2] += coef[ii] * eig1[ii][j1][j2]; |
| 149 | - mat2[j1][j2] += coef[ii] * mat1[ii][j1][j2]; | |
| 155 | + for (jj = 0; jj < 2; jj++) | |
| 156 | + mat2[j1][j2][jj] += coef[ii] * mat1[ii][j1][j2][jj]; | |
| 150 | 157 | }/*for (ii = 0; ii < 4; ii++)*/ |
| 151 | 158 | }/*for (j2 = 0; j2 < 4; j2++)*/ |
| 152 | 159 | }/*for (j1 = 0; j1 < 4; j1++)*/ |
| @@ -154,10 +161,11 @@ void interpol_energy() { | ||
| 154 | 161 | kumo_coef(j1, &coef[0]); |
| 155 | 162 | for (j2 = 0; j2 < 4; j2++) { |
| 156 | 163 | eig3[j2] = 0.0; |
| 157 | - mat3[j2] = 0.0; | |
| 164 | + for (jj = 0; jj < 2; jj++) mat3[j2][jj] = 0.0; | |
| 158 | 165 | for (ii = 0; ii < 4; ii++) { |
| 159 | 166 | eig3[j2] += coef[ii] * eig2[ii][j2]; |
| 160 | - mat3[j2] += coef[ii] * mat2[ii][j2]; | |
| 167 | + for (jj = 0; jj < 2; jj++) | |
| 168 | + mat3[j2][jj] += coef[ii] * mat2[ii][j2][jj]; | |
| 161 | 169 | }/*for (ii = 0; ii < 4; ii++)*/ |
| 162 | 170 | }/*for (j2 = 0; j2 < 4; j2++)*/ |
| 163 | 171 | for (j2 = 0; j2 < interpol; j2++) { |
| @@ -165,16 +173,18 @@ void interpol_energy() { | ||
| 165 | 173 | eig[ib][i0*interpol + j0] |
| 166 | 174 | [i1*interpol + j1] |
| 167 | 175 | [i2*interpol + j2] = 0.0; |
| 168 | - mat[ib][i0*interpol + j0] | |
| 169 | - [i1*interpol + j1] | |
| 170 | - [i2*interpol + j2] = 0.0; | |
| 176 | + for (jj = 0; jj < 2; jj++) | |
| 177 | + mat[ib][i0*interpol + j0] | |
| 178 | + [i1*interpol + j1] | |
| 179 | + [i2*interpol + j2][jj] = 0.0; | |
| 171 | 180 | for (ii = 0; ii < 4; ii++) { |
| 172 | 181 | eig[ib][i0*interpol + j0] |
| 173 | 182 | [i1*interpol + j1] |
| 174 | 183 | [i2*interpol + j2] += coef[ii] * eig3[ii]; |
| 175 | - mat[ib][i0*interpol + j0] | |
| 176 | - [i1*interpol + j1] | |
| 177 | - [i2*interpol + j2] += coef[ii] * mat3[ii]; | |
| 184 | + for (jj = 0; jj < 2; jj++) | |
| 185 | + mat[ib][i0*interpol + j0] | |
| 186 | + [i1*interpol + j1] | |
| 187 | + [i2*interpol + j2][jj] += coef[ii] * mat3[ii][jj]; | |
| 178 | 188 | }/*for (ii = 0; ii < 4; ii++)*/ |
| 179 | 189 | }/*for (j2 = 0; j2 < interpol; j2++)*/ |
| 180 | 190 | }/*for (j1 = 0; j1 < interpol; j1++)*/ |
| @@ -27,6 +27,7 @@ THE SOFTWARE. | ||
| 27 | 27 | #include <stdlib.h> |
| 28 | 28 | #include <stdio.h> |
| 29 | 29 | #include <math.h> |
| 30 | +#include <ctype.h> | |
| 30 | 31 | #include "variable.h" |
| 31 | 32 | #include "basic_math.h" |
| 32 | 33 | #if defined(HAVE_CONFIG_H) |
| @@ -46,6 +47,8 @@ void read_file( | ||
| 46 | 47 | { |
| 47 | 48 | int ib, i, j, i0, i1, i2, ii0, ii1, ii2, ierr; |
| 48 | 49 | FILE *fp; |
| 50 | + char* ctemp1; | |
| 51 | + char ctemp2[256], rc; | |
| 49 | 52 | /* |
| 50 | 53 | Open input file. |
| 51 | 54 | */ |
| @@ -62,7 +65,19 @@ void read_file( | ||
| 62 | 65 | /* |
| 63 | 66 | k-point grid |
| 64 | 67 | */ |
| 65 | - ierr = fscanf(fp, "%d%d%d", &ng0[0], &ng0[1], &ng0[2]); | |
| 68 | + ctemp1 = fgets(ctemp2, 256, fp); | |
| 69 | + ierr = sscanf(ctemp2, "%d%d%d%s", &ng0[0], &ng0[1], &ng0[2], &rc); | |
| 70 | + rc = tolower(rc); | |
| 71 | + if (ierr <= 3) rc = 'r'; | |
| 72 | + if (rc == 'r') | |
| 73 | + printf(" Real k-dependent quantity.\n"); | |
| 74 | + else if (rc == 'c') | |
| 75 | + printf(" Complex k-dependent quantity.\n"); | |
| 76 | + else { | |
| 77 | + printf(" Error! r or c is allowed. Input %c\n", rc); | |
| 78 | + exit(-1); | |
| 79 | + } | |
| 80 | + | |
| 66 | 81 | if (ierr == 0) printf("error ! reading ng"); |
| 67 | 82 | printf(" k point grid : %d %d %d \n", ng0[0], ng0[1], ng0[2]); |
| 68 | 83 | for (i = 0; i < 3; i++) ng[i] = ng0[i]; |
| @@ -127,30 +142,33 @@ void read_file( | ||
| 127 | 142 | Allocation of Kohn-Sham energies $ matrix elements |
| 128 | 143 | */ |
| 129 | 144 | eig0 = (GLfloat****)malloc(nb * sizeof(GLfloat***)); |
| 130 | - mat0 = (GLfloat****)malloc(nb * sizeof(GLfloat***)); | |
| 131 | 145 | eig = (GLfloat****)malloc(nb * sizeof(GLfloat***)); |
| 132 | - mat = (GLfloat****)malloc(nb * sizeof(GLfloat***)); | |
| 146 | + mat0 = (GLfloat*****)malloc(nb * sizeof(GLfloat****)); | |
| 147 | + mat = (GLfloat*****)malloc(nb * sizeof(GLfloat****)); | |
| 133 | 148 | vf = (GLfloat*****)malloc(nb * sizeof(GLfloat****)); |
| 134 | 149 | for (ib = 0; ib < nb; ib++) { |
| 135 | 150 | eig0[ib] = (GLfloat***)malloc(ng0[0] * sizeof(GLfloat**)); |
| 136 | - mat0[ib] = (GLfloat***)malloc(ng0[0] * sizeof(GLfloat**)); | |
| 137 | 151 | eig[ib] = (GLfloat***)malloc(ng0[0] * sizeof(GLfloat**)); |
| 138 | - mat[ib] = (GLfloat***)malloc(ng0[0] * sizeof(GLfloat**)); | |
| 152 | + mat0[ib] = (GLfloat****)malloc(ng0[0] * sizeof(GLfloat***)); | |
| 153 | + mat[ib] = (GLfloat****)malloc(ng0[0] * sizeof(GLfloat***)); | |
| 139 | 154 | vf[ib] = (GLfloat****)malloc(ng0[0] * sizeof(GLfloat***)); |
| 140 | 155 | for (i0 = 0; i0 < ng0[0]; i0++) { |
| 141 | 156 | eig0[ib][i0] = (GLfloat**)malloc(ng0[1] * sizeof(GLfloat*)); |
| 142 | - mat0[ib][i0] = (GLfloat**)malloc(ng0[1] * sizeof(GLfloat*)); | |
| 143 | 157 | eig[ib][i0] = (GLfloat**)malloc(ng0[1] * sizeof(GLfloat*)); |
| 144 | - mat[ib][i0] = (GLfloat**)malloc(ng0[1] * sizeof(GLfloat*)); | |
| 158 | + mat0[ib][i0] = (GLfloat***)malloc(ng0[1] * sizeof(GLfloat**)); | |
| 159 | + mat[ib][i0] = (GLfloat***)malloc(ng0[1] * sizeof(GLfloat**)); | |
| 145 | 160 | vf[ib][i0] = (GLfloat***)malloc(ng0[1] * sizeof(GLfloat**)); |
| 146 | 161 | for (i1 = 0; i1 < ng0[1]; i1++) { |
| 147 | 162 | eig0[ib][i0][i1] = (GLfloat*)malloc(ng0[2] * sizeof(GLfloat)); |
| 148 | - mat0[ib][i0][i1] = (GLfloat*)malloc(ng0[2] * sizeof(GLfloat)); | |
| 149 | 163 | eig[ib][i0][i1] = (GLfloat*)malloc(ng0[2] * sizeof(GLfloat)); |
| 150 | - mat[ib][i0][i1] = (GLfloat*)malloc(ng0[2] * sizeof(GLfloat)); | |
| 164 | + mat0[ib][i0][i1] = (GLfloat**)malloc(ng0[2] * sizeof(GLfloat*)); | |
| 165 | + mat[ib][i0][i1] = (GLfloat**)malloc(ng0[2] * sizeof(GLfloat*)); | |
| 151 | 166 | vf[ib][i0][i1] = (GLfloat**)malloc(ng0[2] * sizeof(GLfloat*)); |
| 152 | - for (i2 = 0; i2 < ng0[2]; ++i2) | |
| 167 | + for (i2 = 0; i2 < ng0[2]; ++i2) { | |
| 168 | + mat0[ib][i0][i1][i2] = (GLfloat*)malloc(3 * sizeof(GLfloat)); | |
| 169 | + mat[ib][i0][i1][i2] = (GLfloat*)malloc(3 * sizeof(GLfloat)); | |
| 153 | 170 | vf[ib][i0][i1][i2] = (GLfloat*)malloc(3 * sizeof(GLfloat)); |
| 171 | + } | |
| 154 | 172 | } |
| 155 | 173 | } |
| 156 | 174 | } |
| @@ -185,11 +203,17 @@ void read_file( | ||
| 185 | 203 | for (i2 = 0; i2 < ng0[2]; ++i2) { |
| 186 | 204 | if (lshift != 0) ii2 = i2; |
| 187 | 205 | else ii2 = modulo(i2 + (ng0[2] + 1) / 2, ng0[2]); |
| 188 | - ierr = fscanf(fp, "%e", &mat0[ib][ii0][ii1][ii2]); | |
| 189 | - } | |
| 190 | - } | |
| 191 | - } | |
| 192 | - } | |
| 206 | + if (rc == 'r') { | |
| 207 | + ierr = fscanf(fp, "%e", &mat0[ib][ii0][ii1][ii2][0]); | |
| 208 | + mat0[ib][ii0][ii1][ii2][1] = 0.0; | |
| 209 | + } | |
| 210 | + else | |
| 211 | + ierr = fscanf(fp, "%e%e", &mat0[ib][ii0][ii1][ii2][0], &mat0[ib][ii0][ii1][ii2][2]); | |
| 212 | + mat0[ib][ii0][ii1][ii2][2] = 0.0; | |
| 213 | + }/*for (i2 = 0; i2 < ng0[2]; ++i2)*/ | |
| 214 | + }/*for (i1 = 0; i1 < ng0[1]; ++i1)*/ | |
| 215 | + }/*for (i0 = 0; i0 < ng0[0]; ++i0)*/ | |
| 216 | + }/*for (ib = 0; ib < nb; ++ib)*/ | |
| 193 | 217 | fclose(fp); |
| 194 | 218 | } /* read_file */ |
| 195 | 219 |
| @@ -49,13 +49,13 @@ int nb; //!< The number of Bands | ||
| 49 | 49 | GLfloat avec[3][3]; //!< Direct lattice vector |
| 50 | 50 | GLfloat bvec[3][3]; //!< Reciprocal lattice vector |
| 51 | 51 | GLfloat ****eig0; //!< Eigenvalues @f$\varepsilon_{n k}@f$[::nb][::ng0[0]][::ng0[1]][::ng0[2]] |
| 52 | -GLfloat ****mat0; //!< Matrix element [::nb][::ng0[0]][::ng0[1]][::ng0[2]] | |
| 52 | +GLfloat *****mat0; //!< Matrix element [::nb][::ng0[0]][::ng0[1]][::ng0[2]][3] | |
| 53 | 53 | /* |
| 54 | 54 | Interpolation |
| 55 | 55 | */ |
| 56 | 56 | int ng[3]; //!< @b Interpolated @f$k@f$-grids |
| 57 | 57 | GLfloat ****eig; //!< Eigenvalues @f$\varepsilon_{n k}@f$ [::nb][::ng[0]][::ng[1]][::ng[2]] |
| 58 | -GLfloat ****mat; //!< Matrix element @f$\delta_{n k}@f$ [::nb][::ng[0]][::ng[1]][::ng[2]] | |
| 58 | +GLfloat *****mat; //!< Matrix element @f$\delta_{n k}@f$ [::nb][::ng[0]][::ng[1]][::ng[2]][3] | |
| 59 | 59 | GLfloat *****vf; //!< Matrix element @f$\{\bf v}_{{\rm f} n k}@f$ [::nb][::ng[0]][::ng[1]][::ng[2]][3] |
| 60 | 60 | int interpol; //!< Ratio of interpolation |
| 61 | 61 | /* |
| @@ -96,7 +96,7 @@ GLfloat ****nmlp; //!< Normal vector of patchs [::nb][::ntri][3][3] | ||
| 96 | 96 | GLfloat ****kvp; //!< @f$k@f$-vectors of points [::nb][::ntri][3][3] |
| 97 | 97 | GLfloat **nmlp_rot; //!< Normal vector of patchs [::nb][::ntri*3*3] |
| 98 | 98 | GLfloat **kvp_rot; //!< @f$k@f$-vectors of points [::nb][::ntri*3*3] |
| 99 | -GLfloat ***matp; //!< Matrix elements of points [::nb][::ntri][3] | |
| 99 | +GLfloat ****matp; //!< Matrix elements of points [::nb][::ntri][3][3] | |
| 100 | 100 | GLfloat **clr; //!< Colors of points [::nb][::ntri*3*4] |
| 101 | 101 | int itet; //!< Counter for tetrahedron |
| 102 | 102 | GLfloat side; //!< Which side is lighted |
fermisurfer
Fork