You're 
th visitor.
gnuplot 5.0以降で導入された、C言語で書かれた関数をgnuplotで取り込んで使うためのコマンド「import」の使用例を紹介します。ここでは、Bardeen-Cooper-Schrieffer（BCS）理論による超伝導ギャップと超伝導比熱の温度依存性を計算する関数をgnuplotで取り込めるようにします。これらの関数は初等関数で表すことはできず、関数の値を得るには、積分方程式を数値的に解いたり数値微分をしたりする必要があります。以下では、「C言語の関数を取り込む（import）その3 GSLを使う」のGSLの導入までは完了していると仮定します。
BCS理論によると、超伝導ギャップ\(\Delta(T)\)の温度依存性は、超伝導ギャップ方程式 \[ 1 = VN(0)\int^{E_{\mathrm{c}}}_0 d\xi\frac{1}{\sqrt{\xi^2+\Delta^2}} \tanh\left(\frac{\sqrt{\xi^2+\Delta^2}}{2k_{\mathrm{B}}T_{\mathrm{c}}} \right) \] の解として与えられます。ここで、\(V\)は対形成相互作用、\(N(0)\)はFermiエネルギー\(E_{\mathrm{F}}\)近傍での状態密度、\(\xi=\hbar^2 k^2/2m^\ast - E_{\mathrm{F}}\)は\(E_{\mathrm{F}}\)を基準とした常伝導電子の運動エネルギー、\(E_{\mathrm{c}}\)はカットオフエネルギー（ふつうはデバイエネルギーに取る）です。この方程式は、温度がゼロか\(T_{\mathrm{c}}\)に十分近い場合でないと解析的に解くことはできませんが、数値計算を使って比較的簡単に全温度にわたって解くことができます。
上述の式を、プログラムしやすいように無次元化しておきます。エネルギーの次元を持つ量は\(k_{\mathrm{B}}T_{\mathrm{c}}\)で規格化します。積分では\(\xi \to x = \xi/(k_{\mathrm{B}}T_{\mathrm{c}}) \)と変数変換し、超伝導ギャップも規格化して\(\delta = \Delta/(k_{\mathrm{B}}T_{\mathrm{c}}) \)とします。また、カットオフエネルギーも\( e_{\mathrm{c}} = E_{\mathrm{c}}/(k_{\mathrm{B}}T_{\mathrm{c}}) \)と規格化します（カットオフエネルギーは、\(k_{\mathrm{B}}T_{\mathrm{c}}\)の数百倍程度にしておけば正しい値が得られます）。さらに、\(\Delta = 0\)（すなわち\(T=T_{\mathrm{c}}\)）とした時のギャップ方程式から、\(1/VN(0)\)は、無次元量\(q\)を用いて \[ \frac{1}{VN(0)} = q \equiv \ln\left(\frac{2e^\gamma e_{\mathrm{c}}}{\pi}\right) \] と表せます。ただし、\(\gamma = 0.57721...\)はEular定数です。
上述の規格化と変形の下で、ギャップ方程式は以下のようになります: \[ \int^{e_{\mathrm{c}}}_0 dx\frac{1}{\sqrt{x^2+\delta^2}} \tanh\left(\frac{\sqrt{x^2+\delta^2}}{2} \right) - q =0 \]
同じくBCS理論によると、超伝導状態における比熱\(C\)の温度依存性は、 \[ C(T) = \frac{2}{T}\sum_{\boldsymbol{k}}\left(-\frac{\partial f_{\boldsymbol{k}} }{\partial E_{\boldsymbol{k}}}\right) \left(E_{\boldsymbol{k}}^2 - \frac{T}{2}\frac{\partial \Delta^2}{\partial T}\right) \] で与えられます。ここで\(E_{\boldsymbol{k}} = \sqrt{\xi_{\boldsymbol{k}}^2 + \Delta^2}\)は準粒子励起エネルギー、\(f_{\boldsymbol{k}} = 1/[\exp(-E_{\boldsymbol{k}}/k_{\mathrm{B}}T) + 1]\)は\(E_{\boldsymbol{k}}\)に関するFermi分布関数です。また、\(\boldsymbol{k}\)に関する和は、\(\xi_{\boldsymbol{k}} = \hbar^2\boldsymbol{k}^2/2m^\ast - \mu\)がカットオフ\(\pm \epsilon_{\mathrm{c}}\)の範囲に入る場合について取ります。
これもギャップ方程式と同様に規格化します。エネルギーを表す量は\(e_{\boldsymbol{k}} = E_{\boldsymbol{k}}/(k_{\mathrm{B}}T_{\mathrm{c}})\)、\(\delta_{\boldsymbol{k}} = \Delta_{\boldsymbol{k}}/(k_{\mathrm{B}}T_{\mathrm{c}})\)、\(e_{\mathrm{F}} = E_{\mathrm{F}}/(k_{\mathrm{B}}T_{\mathrm{c}})\)、\(e_{\mathrm{c}} = E_{\mathrm{c}}/(k_{\mathrm{B}}T_{\mathrm{c}})\)のように\(k_{\mathrm{B}}T_{\mathrm{c}}\)で規格化します。また、温度\(T\)も\(t = T/T_{\mathrm{c}}\)と規格化します。また、波数\(\boldsymbol{k}\)は、\(\hat{\boldsymbol{k}}^2 = e_{\mathrm{c}}-e_{\mathrm{F}} \)となるような\(\hat{\boldsymbol{k}}\)へと規格化します（この規格化の影響は、以下に述べるように\(C_{\mathrm{n}}\)を用いて比熱を規格化する際にキャンセルされます）。これらを用いると、上記の比熱の式は \[ {C} = \frac{2k_{\mathrm{B}}}{t}\sum_{\hat{\boldsymbol{k}}}\left(-\frac{\partial f_{\hat{\boldsymbol{k}}} }{\partial e_{\hat{\boldsymbol{k}}}}\right) \left(e_{\hat{\boldsymbol{k}}}^2 - \frac{t}{2}\frac{\partial \delta^2}{\partial t}\right) \] となります。さらに和を積分で書くと、 \[ {C} = \frac{2k_{\mathrm{B}}}{t}\int_{-k_{\mathrm{cutoff}}}^{k_{\mathrm{cutoff}}}4\pi k^2d k\left(-\frac{\partial f_{\hat{\boldsymbol{k}}} }{\partial e_{\hat{\boldsymbol{k}}}}\right) \left(e_{\hat{\boldsymbol{k}}}^2 - \frac{t}{2}\frac{\partial \delta^2}{\partial t}\right) \] また、比熱は\(T_{\mathrm{c}}\)における常伝導状態の電子比熱\(C_{\mathrm{n}} = \gamma_{\mathrm{e}}T_{\mathrm{c}}\)で規格化したいのですが、\(C_{\mathrm{n}}\)は\(T_{\mathrm{c}}\)での上記積分を実行することで得ることができます。
以下に、\(t = T/T_{\mathrm{c}}\)におけるBCS超伝導ギャップの値を返す関数BCS_gapと超伝導比熱の値を返す関数BCS_Cpを定義し、それらをgnuplotで読み込めるようにするためのC言語プログラムを示します。
#include "gnuplot_plugin.h"
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_roots.h>
#define EPSABS 1e-8
#define EPSREL 1e-10
#define MAXTIMES 400
#define INTEG_WORKSPACE_SIZE 200
#define INTEG_SUBINTERVAL_LIMIT 199
#define DIFF_EPS 1e-5 /*微分用の微小数。小さくしすぎるとこける*/
/* 関数に渡す引数を格納するための構造体 */
struct my_param { double a; double b;};
struct my_param3 { double a1; double a2; double a3;};
struct my_param5 { double a1; double a2; double a3; double a4; double a5;};
/* 超伝導ギャップ計算の際の被積分関数 xがエネルギー、tは規格化温度、dはギャップ */
double gap_integrand(double x, void *p){
struct my_param *params = (struct my_param *)p;
double d = (params -> a);
double t = (params -> b);
return tanh( sqrt(x*x + d*d)/(2.0*t) ) / sqrt(x*x+d*d) ;
}
/* 超伝導ギャップ計算のためのエネルギー積分を実行する関数 */
double gap_integral(double x, void *p){
/* xが規格化したギャップ、tが規格化した温度に対応、
cutoffはcutoffエネルギー（デバイ温度に取る場合が多い）に対応する。 */
struct my_param *Iparams = (struct my_param *)p;
double cutoff = (Iparams -> a);
double t = (Iparams -> b);
double q = log(2*exp(M_EULER)*cutoff/M_PI);
double uplim;
double result;
double abserr;
int neval;
gsl_function G;
struct my_param Gparams = {x, t};
gsl_integration_workspace *workspace = gsl_integration_workspace_alloc (INTEG_WORKSPACE_SIZE);
double gap_integrand(double, void *);
G.function = &gap_integrand;
G.params = &Gparams;
/* 積分の上限 */
uplim = cutoff;
gsl_integration_qag(&G, 0.0, uplim, EPSABS, EPSREL, INTEG_SUBINTERVAL_LIMIT, 1, workspace, &result, &abserr);
gsl_integration_workspace_free (workspace);
return result - q;
}
/* 超伝導ギャップの計算、tは規格化した温度T/Tc、cutoffはカットオフエネルギー。*/
double gap_main(double t, double cutoff){
double gap;
double gap_integral(double, void *);
int times;
int status;
double x_l, x_r;
gsl_function F;
struct my_param Fparams;
/* 解を数値的に求める際のWorkspaceの確保 */
const gsl_root_fsolver_type *T = gsl_root_fsolver_bisection;
gsl_root_fsolver *workspace_f = gsl_root_fsolver_alloc (T);
F.function = &gap_integral;
(&Fparams)->a = cutoff;
(&Fparams)->b = t;
F.params = &Fparams;
gsl_root_fsolver_set (workspace_f, &F, 0.001, 2.0);
/* 解を求めるループ */
for(times = 0; times < MAXTIMES; times++)
{
status = gsl_root_fsolver_iterate(workspace_f);
x_l = gsl_root_fsolver_x_lower(workspace_f);
x_r = gsl_root_fsolver_x_upper(workspace_f);
status = gsl_root_test_interval(x_l, x_r, EPSABS, EPSREL);
/* 集束したらループ停止 */
if(status != GSL_CONTINUE)
{
gap = (x_l+x_r)*0.5;
break;
}
}
/* free */
gsl_root_fsolver_free(workspace_f);
return gap;
}
/* 超伝導体での準粒子エネルギーE(k)を求める関数 kは波数、efはFermi Energy、dはGap */
double ek(double k, double ef, double d){
return sqrt((k*k-ef)*(k*k-ef) + d*d);
}
/* 比熱計算の際の被積分関数 */
double specific_heat_integrand(double x, void *p){
double ek(double, double, double);
struct my_param5 *Cparams = (struct my_param5 *)p;
double e;
double ef = (Cparams -> a1); /* Fermi energy */
double t = (Cparams -> a2); /* 規格化温度 */
double cutoff = (Cparams -> a3); /* Cutoff */
double d = (Cparams -> a4); /* gap */
double d2d = (Cparams -> a5); /* gapの2乗の温度微分 */
e = ek(x, ef, d); /* 波数xの準粒子エネルギー */
return x*x*( e*e - t*0.5*d2d ) * ( exp(-e/t)/t ) / ( (exp(-e/t)+1.0) * (exp(-e/t)+1.0) );
}
/* 超伝導状態の比熱を返す関数 tは規格化温度、cutoffはcutoffエネルギー、efはFermi Energy（10000程度） */
double specific_heat_main(double t, double cutoff, double ef){
double d, d1, dt, d2d;
double uplim, lowlim, result, abserr;
gsl_function C;
struct my_param5 Cparams;
double specific_heat_integrand(double, void *);
double gap_main(double, double);
/* 積分するためのWorkspaceの確保 */
gsl_integration_workspace *workspace = gsl_integration_workspace_alloc (INTEG_WORKSPACE_SIZE);
if (t<1.0){
d = gap_main(t, cutoff);
dt = DIFF_EPS;
d1 = gap_main(t+dt, cutoff);
// gapの2乗の温度微分
d2d = (d1*d1-d*d)/dt;
}
else {
// Normal state (t>1) ではGapはゼロ
d = 0.0;
d2d = 0.0;
}
C.function = &specific_heat_integrand;
(&Cparams)->a1 = ef;
(&Cparams)->a2 = t;
(&Cparams)->a3 = cutoff;
(&Cparams)->a4 = d;
(&Cparams)->a5 = d2d;
C.params = &Cparams;
/* 積分の上限と下限 */
uplim = sqrt(ef + cutoff );
lowlim = sqrt(ef - cutoff );
/* 積分の実行 */
gsl_integration_qag(&C, lowlim, uplim, EPSABS, EPSREL, INTEG_SUBINTERVAL_LIMIT, 1, workspace, &result, &abserr);
/* Workspaceを開放する */
gsl_integration_workspace_free (workspace);
/* 積分結果を温度で割ると比熱になる */
return result/t;
}
/* gnuplotでBCSギャップ関数を使うための関数
引数は2つで、1つ目が規格化温度、2つ目はカットオフエネルギー（普通の用途では500程度にしておく）*/
DLLEXPORT struct value BCS_gap(int nargs, struct value *arg, void *p){
double gap_main(double, double);
double t;
double cutoff;
struct value r;
/* 引数の数とタイプのチェック */
RETURN_ERROR_IF_WRONG_NARGS(r, nargs, 2);
RETURN_ERROR_IF_NONNUMERIC(r, arg[0]);
RETURN_ERROR_IF_NONNUMERIC(r, arg[1]);
/* 引数の代入 */
t = RVAL(arg[0]);
cutoff = fabs(RVAL(arg[1]));
/* 結果を返すための構造体のタイプを指定 */
r.type = CMPLX;
/* 虚部はゼロ */
r.v.cmplx_val.imag = 0.0;
/* 実部に結果を代入 */
if (t >= 1.0){
/* Normal状態ではギャップゼロ */
r.v.cmplx_val.real = 0.0;
}
else if (t < 0.001){
/* tがゼロに近いとエラーになるので、tが小さい場合はt=0.001での値で代用 */
r.v.cmplx_val.real = gap_main(0.001,cutoff);
}
else {
/* 上記以外ではギャップを計算 */
r.v.cmplx_val.real = gap_main(t,cutoff);
}
return r;
}
/* gnuplotでBCSの比熱(Normal状態の電子比熱係数で規格化)を使うための関数
引数は3つで、一つ目は規格化温度t、二つ目はcutoff（500程度にする）、三つめはFermi energy（10000程度にする）*/
DLLEXPORT struct value BCS_Cp(int nargs, struct value *arg, void *p){
double specific_heat_main(double, double, double);
struct value r;
double t;
double cutoff;
double ef;
/* 引数の数とタイプのチェック */
RETURN_ERROR_IF_WRONG_NARGS(r, nargs, 3);
RETURN_ERROR_IF_NONNUMERIC(r, arg[0]);
RETURN_ERROR_IF_NONNUMERIC(r, arg[1]);
RETURN_ERROR_IF_NONNUMERIC(r, arg[2]);
/* 引数の代入 */
t = RVAL(arg[0]);
cutoff = fabs(RVAL(arg[1]));
ef = fabs(RVAL(arg[2]));
/* 結果を返すための構造体のタイプを指定 */
r.type = CMPLX;
/* 虚部はゼロ */
r.v.cmplx_val.imag = 0.0;
/* 実部に結果を代入 */
if (t >= 1.0){
/* Normal状態の場合は温度に比例する比熱
（Normal状態の電子比熱係数で規格化するので、温度そのものを返せばいい*/
r.v.cmplx_val.real = t;
}
else if (t < 0.001){
/* 温度が十分低い場合は比熱ゼロ い*/
r.v.cmplx_val.real = 0.0;
}
else {
/* それ以外は比熱を計算しt=1の時の値で割る。*/
r.v.cmplx_val.real = specific_heat_main(t,cutoff,ef)/specific_heat_main(1.0,cutoff,ef);
}
return r;
}上記のプログラムをgap.cという名前で保存し、msys2/mingw-W64のMinGW64シェルまたはMinGW32シェルから以下のコマンドでコンパイルします。
gcc -static -shared -o gap.dll gap.c -lgsl -DHAVE_CONFIG_H -I.コンパイルがうまくいくと、gap.dllというdllファイルができます。
このdllを用いて、c言語で定義した関数をプロットしてみます。また、計算がうまくできていることを確かめるために、昔の論文：B. Mühlschlegel, Z. Phys. 155, 313 (1959).に掲載されている以下の数表と比較してみます（tが1以上のデータは筆者が追加）。このデータはBCS理論で得られる、超伝導ギャップ、エントロピー、自由エネルギー、臨界磁場、侵入長、比熱をt = T/Tcに対して計算したものです。
#Bernhard Muehlschlegel, Z. Phys. 155, 313 (1959).
#Cn=gamma*Tc
#t E_T/E_0 S/Cn -F/(Cn*Tc) Hc^2/(4*pi*Cn*Tc) 1-Lam_0/Lam_T C/Cn
1.4 0 1.0 0.5 0.0 1.0 1.4
1.35 0 1.0 0.5 0.0 1.0 1.35
1.3 0 1.0 0.5 0.0 1.0 1.3
1.25 0 1.0 0.5 0.0 1.0 1.25
1.2 0 1.0 0.5 0.0 1.0 1.2
1.15 0 1.0 0.5 0.0 1.0 1.15
1.1 0 1.0 0.5 0.0 1.0 1.1
1.05 0 1.0 0.5 0.0 1.0 1.05
1.001 0 1.0 0.5 0.0 1.0 1.001
1 0 1.0 0.5 0.0 1.0 2.4261
0.98 0.2436 0.9519 0.4805 0.0003 0.9601 2.3314
0.96 0.3416 0.9048 0.4619 0.0011 0.9206 2.2378
0.94 0.4148 0.8587 0.4443 0.0025 0.8814 2.1454
0.92 0.4749 0.8136 0.4276 0.0044 0.8425 2.0541
0.90 0.5263 0.7694 0.4117 0.0067 0.8041 1.9639
0.88 0.5715 0.7263 0.3968 0.0096 0.7660 1.8750
0.86 0.6117 0.6842 0.3827 0.0129 0.7283 1.7874
0.84 0.6480 0.6432 0.3694 0.0166 0.6911 1.7010
0.82 0.6810 0.6032 0.3569 0.0207 0.6544 1.6159
0.80 0.7110 0.5643 0.3453 0.0253 0.6182 1.5321
0.78 0.7386 0.5266 0.3344 0.0302 0.5826 1.4498
0.76 0.7640 0.4900 0.3242 0.0354 0.5475 1.3689
0.74 0.7874 0.4546 0.3148 0.0410 0.5131 1.2894
0.72 0.8089 0.4203 0.3060 0.0468 0.4793 1.2115
0.70 0.8288 0.3873 0.2979 0.0529 0.4463 1.1352
0.68 0.8471 0.3554 0.2905 0.0593 0.4140 1.0605
0.66 0.8640 0.3249 0.2837 0.0659 0.3825 0.9874
0.64 0.8796 0.2956 0.2775 0.0727 0.3518 0.9162
0.62 0.8939 0.2676 0.2719 0.0797 0.3221 0.8467
0.60 0.9070 0.2410 0.2668 0.0868 0.2933 0.7792
0.58 0.9190 0.2157 0.2622 0.0940 0.2656 0.7136
0.56 0.9299 0.1918 0.2582 0.1014 0.2389 0.6501
0.54 0.9399 0.1693 0.2545 0.1087 0.2133 0.5888
0.52 0.9488 0.1482 0.2514 0.1162 0.1890 0.5298
0.50 0.9569 0.1285 0.2486 0.1236 0.1660 0.4731
0.48 0.9641 0.1103 0.2462 0.1310 0.1442 0.4190
0.46 0.9704 0.0937 0.2442 0.1384 0.1239 0.3675
0.44 0.9760 0.0784 0.2425 0.1457 0.1055 0.3188
0.42 0.9809 0.0646 0.2410 0.1528 0.0878 0.2731
0.40 0.9850 0.0524 0.2399 0.1599 0.0721 0.2305
0.38 0.9885 0.0416 0.2389 0.1667 0.0580 0.1913
0.36 0.9915 0.0322 0.2382 0.1734 0.0456 0.1555
0.34 0.9938 0.0243 0.2376 0.1798 0.0348 0.1233
0.32 0.9957 0.0177 0.2372 0.1860 0.0257 0.0950
0.30 0.9971 0.0124 0.2369 0.1919 0.0182 0.0706
0.28 0.9982 0.0082 0.2367 0.1975 0.0123 0.0502
0.26 0.9989 0.0051 0.2366 0.2028 0.0078 0.0338
0.24 0.9994 0.0030 0.2365 0.2077 0.0046 0.0212
0.22 0.9997 0.0016 0.2365 0.2123 0.0024 0.0121
0.20 0.9999 0.0007 0.2364 0.2164 0.0011 0.0061
0.18 1.0000 0.0003 0.2364 0.2202 0.0005 0.0027
0.16 1.0000 0.0001 0.2364 0.2236 0.0001 0.0009
0.14 1.0000 0.0000 0.2364 0.2266 0.0000 0.0002このデータをMuehlschlegel.datとして保存しておき、以下のようなスクリプトを実行します。これは、カットオフエネルギーを500（×kBTc）としたときのBCSギャップ方程式の解を、先行研究の結果と比べることになります。
import BCS_gap(x,cutoff) from "gap:BCS_gap"
set term wxt font "Arial,12"
set samples 1000
set xrange [0:1]
set xlabel "T / Tc"
set ylabel "Delta"
plot BCS_gap(x,500), "Muehlschlegel.dat" u 1:($2*1.764) title "literature"すると、以下のような図が得られます。先行研究の計算結果を、我々のBCS_gap関数がしっかりと再現していることがわかります。
より詳細に誤差を検討するために、相対誤差：（文献値 − 計算値）/ 文献値を計算してみます。
import BCS_gap(x,cutoff) from "gap:BCS_gap"
set term wxt font "Arial,12"
set xrange [0:1]
set xlabel "T / Tc"
set ylabel "Relative Error (%)"
plot "BCS_thermodynamic_quantity.dat" u 1:($2*1.764-BCS_gap($1,500))/($2*1.764)*100 title "relative error"この実行結果は、以下のようになり、相対誤差は0.01%程度と、実用上問題ない程度になっています。
次に、以下のようなスクリプトで、BCS理論から得られる超伝導体の電子比熱（を温度で割ったもの）の温度依存性を先行研究と比較することにします。先述のようにカットオフエネルギーを500（×kBTc）とし、Fermiエネルギーは100000（×kBTc）に取ります。
import BCS_Cp(x,cutoff,eF) from "gap:BCS_Cp"
set key left Left
set xrange [0.0:1.2]
set xlabel "T / Tc"
set ylabel "C/gamma T"
plot BCS_Cp(x,500,100000)/x, \
"Muehlschlegel.dat" u 1:($7/$1) title "literature"
この実行結果は以下のようになり、またしても先行研究の結果をうまく再現していることがわかります。
ギャップ関数の時と同様に相対誤差を計算してみます。
import BCS_Cp(x,cutoff,eF) from "gap:BCS_Cp"
set xrange [0.0:1.2]
set xlabel "T / Tc"
set ylabel "Relative Error (%)"
plot "BCS_thermodynamic_quantity.dat" u 1:($7/$1-BCS_Cp($1,500,100000)/$1)/($7/$1)*100 notitleこの実行結果は以下（下の図はx軸近辺の拡大図）のようになり、特異点となるT=Tcちょうどと、比熱が小さくなる低温では誤差がやや大きくなっていますが、それ以外の温度域では0.5%以下の精度が得られています。
コンパイルしたdllを以下に置いておくので、（自己責任で）ご自由にお使いください。