import os
import re
import csv
import numpy as np
from math import log, exp
import scipy.special
from scipy import optimize
from scipy.optimize import minimize
from matplotlib import pyplot as plt
pi = 3.14159265358979323846
h = 6.6260755e-34 # Js";
hbar = 1.05459e-34 # "Js";
c = 2.99792458e8 # m/s";
e = 1.60218e-19 # C";
e0 = 8.854418782e-12; # C2N-1m-2";
kB = 1.380658e-23 # JK-1";
me = 9.1093897e-31 # kg";
nm2m = 1.0e-9
kMR = e * e / pi / h # convert delta_sigma to sheet conductance S
# parameters
# mode: 'sim', 'fit'
mode = 'sim'
# simulation B range
Bmin = -5.0
Bmax = 5.0
nB = 101
Bth = 1.0e-10
Bstep = (Bmax - Bmin) / (nB - 1)
# Files. Resistivity is in ohm*m
file_xx = 'MRxx.csv'
file_xy = 'MRxy.csv'
basename = os.path.basename(file_xx)
dirname = os.path.dirname(file_xx)
header, ext = os.path.splitext(file_xx)
filebody = os.path.basename(file_xx)
output_csv = os.path.join(dirname, filebody + '-fitting.csv')
output_txt = os.path.join(dirname, filebody + '-fitting.txt')
thickness = 5.0 # in nm
headerxx = []
Bxx = [] # magnetic field in T
MRxx = [] # resistivity in ohm*m
MCxx = [] # sheet conductivity in S
# fitting B range
Bxxmin = -5.0
Bxxmax = 5.0
# optimization parameters
sigma00 = 0.0 # sheet conductance at B = 0 in S
alpha0 = 1.0 # prefactor of weak (anti-)localization model. 1 for WL, -1 for AWL
mfp0 = 20.0 # Mean free path in nm
# initial parameter list
ai0 = [sigma00, alpha0, mfp0]
optid = [ 1, 0, 1]
#=============================================
# scipy.optimize.minimizeで使うアルゴリズム
#=============================================
#nelder-mead Downhill simplex
#powell Modified Powell
#cg conjugate gradient (Polak-Ribiere method)
#bfgs BFGS法
#newton-cg Newton-CG
#trust-ncg 信頼領域 Newton-CG 法
#dogleg 信頼領域 dog-leg 法
#L-BFGS-B’ (see here)
#TNC’ (see here)
#COBYLA’ (see here)
#SLSQP’ (see here)
#trust-constr’(see here)
#dogleg’ (see here)
#trust-exact’ (see here)
#trust-krylov’ (see here)
method = "nelder-mead"
#method = 'cg'
#method = 'powell'
#method = 'bfgs'
maxiter = 1000
tol = 1.0e-6
h_diff = 1.0e-3
# graph variabls
fig = None
# 実数値に変換できない文字列をfloat()で変換するとエラーになってプログラムが終了する
# この関数は、変換できなかったらNoneを返すが、プログラムは終了させない
def pfloat(str):
try:
return float(str)
except:
return None
# pfloat()のint版
def pint(str):
try:
return int(str)
except:
return None
# 起動時引数を取得するsys.argリスト変数は、範囲外のindexを渡すとエラーになってプログラムが終了する
# egtarg()では、範囲外のindexを渡したときは、defvalを返す
def getarg(position, defval = None):
try:
return sys.argv[position]
except:
return defval
# 起動時引数を実数に変換して返す
def getfloatarg(position, defval = None):
return pfloat(getarg(position, defval))
# 起動時引数を整数値に変換して返す
def getintarg(position, defval = None):
return pint(getarg(position, defval))
def usage(sigma0 = sigma00, alpha = alpha0, mfp = mfp0):
global Bmin, Bmax, nB
global file_xx, file_xy, thickness
global Bxxmin, Bxxmax
print("")
print("Usage: Variables in () are optional")
print(" MFP in nm B in T thickness in nm")
print(" (i) python {} sim (thickness sigma0 alpha mfp Bmin Bmax nB"
.format(sys.argv[0]))
print(" ex: python {} sim {} {} {} {} {} {} {}"
.format(sys.argv[0], thickness, sigma0, alpha, mfp, Bmin, Bmax, nB))
print(" (ii) python {} fit (file_xx thickness sigma0 id alpha id MFP id Bxxmin Bxxmax)"
.format(sys.argv[0]))
print(" ex: python {} fit {} {} {} {} {} {} {} {} {} {}"
.format(sys.argv[0], file_xx, thickness, sigma0, optid[0], alpha, optid[1], mfp, optid[2], Bxxmin, Bxxmax))
def terminate(message = None):
if message is not None:
print("")
print(message)
usage()
print("")
exit()
# 起動時引数を使ってglobal変数を更新する
def updatevars():
global mode
global Bmin, Bmax, Bstep, nB
global file_xx, thickness, Bxxmin, Bxxmax
global sigma00, alpha0, mfp0
global ai0, optid
argv = sys.argv
if len(argv) == 1:
terminate()
mode = getarg( 1, mode)
if mode == 'sim':
thickness = getfloatarg( 2, thickness)
sigma00 = getfloatarg( 3, sigma00)
alpha0 = getfloatarg( 4, alpha0)
mfp0 = getfloatarg( 5, mfp0)
Bmin = getfloatarg( 6, Bmin)
Bmax = getfloatarg( 7, Bmax)
nB = getintarg ( 8, nB)
Bstep = (Bmax - Bmin) / (nB - 1)
elif mode == 'fit':
file_xx = getarg ( 2, file_xx)
thickness = getfloatarg( 3, thickness)
sigma00 = getfloatarg( 4, sigma00)
optid[0] = getintarg ( 5, optid[0])
alpha0 = getfloatarg( 6, alpha0)
optid[1] = getintarg ( 7, optid[1])
mfp0 = getfloatarg( 8, mfp0)
optid[2] = getintarg ( 9, optid[2])
Bxxmin = getfloatarg(10, Bmin)
Bxxmax = getfloatarg(11, Bmax)
else:
terminate("Error: Invalid mode [{}]".format(mode))
ai0 = [sigma00, alpha0, mfp0]
#=============================
# other functions
#=============================
def savecsv(outfile, header, datalist):
try:
print("Write to [{}]".format(outfile))
f = open(outfile, 'w')
except:
# except IOError:
print("Error: Can not write to [{}]".format(outfile))
else:
fout = csv.writer(f, lineterminator='\n')
fout.writerow(header)
# fout.writerows(data)
for i in range(0, len(datalist[0])):
a = []
for j in range(len(datalist)):
a.append(datalist[j][i])
fout.writerow(a)
f.close()
def read_csv(infile, xmin = None, xmax = None, delimiter = ','):
data = []
try:
infp = open(infile, "r")
f = csv.reader(infp, delimiter = delimiter)
header = next(f)
print("header=", header)
for j in range(len(header)):
data.append([])
for row in f:
# print("row=", row)
x = pfloat(row[0])
if xmin is not None and xmin <= x <= xmax:
y = pfloat(row[1])
data[0].append(x)
data[1].append(y)
except:
terminate("Error: Can not read [{}]".format(infile))
return header, data[0], data[1]
def recover_parameters(ais, optid):
global ai0, sigma00, alpha0, mfp0
ai = []
c = 0
for i in range(len(optid)):
if optid[i] == 1:
ai.append(ais[c])
c += 1
else:
ai.append(ai0[i])
return ai
def choose_parameters(ais, optid):
global ai0, sigma00, alpha0, mfp0
ai = []
c = 0
for i in range(len(optid)):
if optid[i] == 1:
ai.append(ais[i])
return ai
# 弱局在による磁気伝導率変化 Δσ(B)の計算
def sigma_WL(B, l_in, alpha):
x = l_in * l_in * 4.0 * e * abs(B) / hbar
dsigma = scipy.special.digamma(0.5 + 1.0 / x) + log(x) #縦軸 e*e/pi/hbar単位
return alpha * kMR * dsigma
def sigma(B, l_in, alpha, sigma0):
return sigma0 + sigma_WL(B, l_in, alpha)
def calS2(xk):
global model
global Bxx, MRxx, MCxx
global alpha0
sigma0, alpha, mfp = recover_parameters(xk, optid)
sigma0 = max([0.0, sigma0])
mfp = max([0.0, mfp])
l_in = mfp * nm2m
S2 = 0.0
for i in range(len(Bxx)):
B = Bxx[i]
if abs(B) < Bth:
continue
sigmacal = sigma(B, l_in, alpha0, sigma0)
S2 += pow(sigmacal - MCxx[i], 2)
return S2 / len(Bxx)
# First derivatives to be used e.g. for cg method
# Approximate by forward difference method with the delta h = h_diff
def diff1(ai):
n = len(ai)
f0 = calS2(ai)
df = np.empty(n)
for i in range(n):
aii = ai
aii[i] = ai[i] + h_diff
df[i] = (calS2(aii) - f0) / h_diff
return df
# Callback function for scipy.optimize.minimize()
# Print variables every iteration, and update graph for every graphupdateinterval iterationsおｔ
iter = 0
def callback(fig, xk):
global iter
S2 = calS2(xk)
sigma0, alpha, mfp = recover_parameters(xk, optid)
# パラメータと残差二乗和を表示
print("callback {}: sigma0={:14.4g} alpha={:10.4g} mfp={:10.4g} S2={:12.6g}".format(iter, sigma0, alpha, mfp, S2))
iter += 1
def fitting():
global file_xx
global mfp0, alpha0, sigma00
global Bxxmin, Bxxmax
global headerxx, Bxx, MRxx, MCxx
global fig
print("")
print("MR_xx file : ", file_xx)
print(" Read [{}]".format(file_xx))
headerxx, Bxx, MRxx = read_csv(file_xx, Bxxmin, Bxxmax)
print(" Bxx range: {} - {}".format(Bxxmin, Bxxmax))
print(" thickness: {} nm".format(thickness))
print("")
print("Parameters")
print(" sigma0: {} S".format(sigma00))
print(" alpha: {}".format(alpha0))
print(" mfp: {} nm".format(mfp0))
print("")
print("Optimization configuration")
print(" method: ", method)
print(" maxiter: ", maxiter)
print(" tol : ", tol)
print("")
MCxx = []
for i in range(len(Bxx)):
MCxx.append(1.0 / MRxx[i] * thickness * nm2m) # S/m * thickness[m]
print("")
print("Calculate initial data")
l_in = mfp0 * nm2m
yMCxxcal = []
print("{:4s}\t{:8s}\t{:10s}\t{:10s}".format('i', 'B(T)', 'Ssheet(exp)(S)', 'Ssheet(cal)(S)'))
for i in range(len(Bxx)):
B = Bxx[i]
if abs(B) < Bth:
continue
sigmacal = sigma(B, l_in, alpha0, sigma00)
print("{:4d}\t{:8.4f}\t{:10.4g}\t{:10.4g}".format(i, B, MCxx[i], sigmacal))
yMCxxcal.append(sigmacal)
#=============================
# Optimization
#=============================
print("")
print("Optimization by {}".format(method))
print(" tol=", tol)
print(" ai0 =", ai0)
print(" optid=", optid)
ai = choose_parameters(ai0, optid)
print(" ai =", ai)
ret = minimize(calS2, ai, method = method, jac = diff1, tol = tol,
callback = lambda xk: callback(fig, xk),
options = {'maxiter':maxiter, "disp":True})
print("ret=", ret)
if method == 'nelder-mead':
simplex = ret['final_simplex']
ai = simplex[0][0]
fmin = ret['fun']
elif method == 'cg':
ai = ret['x']
fmin = ret['fun']
elif method == 'powell':
ai = ret['x']
fmin = ret['fun']
elif method == 'bfgs':
ai = ret['x']
fmin = calS2(ai)
sigma0, alpha, mfp = recover_parameters(ai, optid)
print("Optimized at S2={:12.6g}:".format(fmin))
print(" sigma0={:12.6g} S".format(sigma0))
print(" alpha={:12.6g}".format(alpha))
print(" mfp={:12.6g} nm".format(mfp))
print("")
print("Calculate final data")
l_in = mfp * nm2m
yMCxxfit = []
print("sigma0: {} S".format(sigma0))
print("{:4s}\t{:8s}\t{:10s}\t{:10s}".format('i', 'B(T)', 'dSsheet(exp)(S)', 'dSsheet(cal)(S)'))
for i in range(len(Bxx)):
B = Bxx[i]
if abs(B) < Bth:
continue
sigmafit = sigma(B, l_in, alpha, sigma0)
print("{:4d}\t{:8.4f}\t{:10.4g}\t{:10.4g}".format(i, B, MCxx[i] - sigma0, sigmafit - sigma0))
yMCxxfit.append(sigmafit)
#=============================
# Save fitting results
#=============================
print("")
print("Save fitting results to [{}]".format(output_csv))
savecsv(output_csv, ["B(T)", "rhoxx(exp)(ohm*m)", "Ssheet(exp)(S)", "Ssheet(ini)(S)", "Ssheet(fin)(S)"],
[Bxx, MRxx, MCxx, yMCxxcal, yMCxxfit])
print("")
print("Save fitting results to [{}]".format(output_txt))
f = open(output_txt, 'w')
f.write("file_xx={}\n".format(file_xx))
f.write("file_xy={}\n".format(file_xy))
f.write("thickness={} nm\n".format(thickness))
f.write("Bxxmin={} T\n".format(Bxxmin))
f.write("Bxxmax={} T\n".format(Bxxmax))
f.write("sigma0={} S\n".format(sigma0))
f.write("alpha={}\n".format(alpha))
f.write("mfp={} nm\n".format(mfp))
f.close()
#=============================
# Plot graphs
#=============================
print("")
print("plot")
fig = plt.figure(figsize = (8, 8))
ax1 = fig.add_subplot(2, 2, 1)
ax2 = fig.add_subplot(2, 2, 3)
ax3 = fig.add_subplot(2, 2, 4)
ax1.plot(Bxx, MRxx, label = 'MRxx', linestyle = 'none', marker = 'o', markersize = 1.0)
# ax1.plot(Bxx, yds, label = 'MRxx(WL,cal)', color = 'red')
ax1.set_xlabel("B(T)")
ax1.set_ylabel("rhoxx(ohm*m)")
ax1.set_title(file_xx)
ax1.legend()
ax2.plot(Bxx, MCxx, label = 'exp', linestyle = 'none', marker = 'o', markersize = 1.0)
ax2.plot(Bxx, yMCxxcal, label = 'initial', linestyle = '-', color = 'red', linewidth = 1.0)
ax2.plot(Bxx, yMCxxfit, label = 'final', linestyle = '-', color = 'blue', linewidth = 1.0)
ax2.set_xlabel("B(T)")
ax2.set_ylabel("Ssheet,xx(S)")
ax2.set_title(file_xx)
ax2.legend()
ax3.plot(Bxx, yMCxxfit, label = 'Ssheet,xx(WL,cal)')
ax3.set_xlabel("B(T)")
ax3.set_ylabel("Ssheet,xx(S)")
ax3.set_title("l_in = {:12.6g} nm".format(mfp))
ax3.legend()
plt.tight_layout()
plt.pause(0.1)
print("")
usage()
print("")
print("Next suggestion:")
print("python {} fit {} {} {:12.6g} {} {:12.6g} {} {:12.6g} {} {} {}"
.format(sys.argv[0], file_xx, thickness, sigma0, optid[0], alpha, optid[1], mfp, optid[2], Bxxmin, Bxxmax))
print("")
print("Press ENTER to exit>>", end = '')
input()
exit()
def simulate():
global nm2m
global mfp0, algpha0, sigma00
print("B range: {} - {} at {} step, {} points".format(Bmin, Bmax, Bstep, nB))
l_in = mfp0 * nm2m
xB = []
xx = []
yMCxxcal = []
yMRxxcal = []
print("{:4s}\t{:8s}\t{:10s}\t{:10s}\t{:10s}".format('i', 'B(T)', 'x', 'sigma(cal)(S)', 'rho(cal)(ohm*m)'))
for i in range(nB):
B = Bmin + i * Bstep
if abs(B) < Bth:
continue
x = l_in * l_in * 4.0 * e * B / hbar # dsの計算のためには必要はないが、printとxxリストへの追加に必要
dsigma = sigma(B, l_in, alpha0, sigma00)
rho = 1.0 / dsigma * (thickness * nm2m) # to ohm*m
print("{:4d}\t{:8.4f}\t{:10.4g}\t{:10.4g}\t{:10.4g}".format(i, B, x, dsigma, rho))
xB.append(B)
xx.append(x)
yMCxxcal.append(dsigma)
yMRxxcal.append(rho)
print("")
print("plot")
fig = plt.figure(figsize = (8, 8))
ax1 = fig.add_subplot(2, 2, 1)
ax2 = fig.add_subplot(2, 2, 2)
ax3 = fig.add_subplot(2, 2, 3)
ax4 = fig.add_subplot(2, 2, 4)
plt.title("mfp = {} nm".format(mfp0))
ax1.plot(xB, yMRxxcal, label = 'rho-B')
ax1.set_xlabel("B(T)")
ax1.set_ylabel("Ssheet(S)")
ax1.legend()
ax2.plot(xx, yMRxxcal, label = 'rho-x')
# ax2.set_ylim([-10, 10])
ax2.set_xlabel("x = l_in^2*4*e*B/hbar")
ax2.set_ylabel("Ssheet(S)")
ax2.legend()
ax3.plot(xB, yMCxxcal, label = 'Ssheet-B')
ax3.set_xlabel("B(T)")
ax3.set_ylabel("Ssheet(S)")
ax3.legend()
ax4.plot(xx, yMCxxcal, label = 'Ssheet-x')
# ax4.set_ylim([-10, 10])
ax4.set_xlabel("x = l_in^2*4*e*B/hbar")
ax4.set_ylabel("Ssheet(S)")
ax4.legend()
plt.tight_layout()
plt.pause(0.1)
print("")
usage()
print("")
print("Press ENTER to exit>>", end = '')
input()
def main():
global mode
global mfp
global Bmin, Bmax, Bstep, nB
updatevars()
print("")
print("mode: ", mode)
print("sigma0 : {} S".format(sigma00))
print("alpha : {} nm".format(alpha0))
print("Mean free path: {} nm".format(mfp0))
if mode == 'sim':
simulate()
elif mode == 'fit':
fitting()
else:
terminate("Error: Invalide mode [{}]".format(mode))
terminate()
if __name__ == '__main__':
main()
