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
"""
縦MR, 横MRから弱(反)局在、Drude型電子−正孔混合伝導のキャリア濃度、移動度をもとめる
最適化にはSIMPLEX法を使う
係数の比較のため、線形最小二乗法による多項式回帰も可能
ρxx(B)はシート抵抗Sxx(B)でフィッティング
ρxy(B)は横抵抗率でフィッティング
cf. Two carrier model for Hall resistance and magnetoresistance, APL 107, 182411 (2015)
"""
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 # convert nm to m
kMR = e * e / pi / h # convert delta_sigma to sheet conductance S
# parameters
# mode: 'sim', 'fit'
mode = 'sim'
# model: combination of 'poly', 'wl', 'wal', 'single', 'two', 'paraboric', 'same_mu', 'same_n'
model = 'wl'
# simulation B range
Bmin = -5.0
Bmax = 5.0
nB = 101
Bstep = (Bmax - Bmin) / (nB - 1)
xth = 1.0e-10 # minimum threshold to avoid x = zero calculation
# Files. The resistivity in the files are in ohm*m, converted to ohm*cm
file_xx = 'MRxx.csv'
file_xy = 'MRxy.csv'
#file_xy = '-'
thickness = 5.0 # in nm
# fitting B range
Bxxmin = -5.0
Bxxmax = 5.0
Bxymin = -5.0
Bxymax = 5.0
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_fittingxx_csv = os.path.join(dirname, header + '-xx-fitting.csv')
output_fittingxy_csv = os.path.join(dirname, header + '-xy-fitting.csv')
output_simulationxx_csv = os.path.join(dirname, header + '-xx-simulation.csv')
output_simulationxy_csv = os.path.join(dirname, header + '-xy-simulation.csv')
output_txt = os.path.join(dirname, header + '-fitting.txt')
output_bat = os.path.join(dirname, header + '.bat')
headerxx = []
headerxy = []
Bxx = [] # magnetic field in T
Bxy = [] # magnetic field in T
MRxx = [] # resistivity in ohm*m
MCxx = [] # sheet conductivity in S
MRxy = []
MCxy = []
yMRxxini = []
yMCxxini = []
yMRxyini = []
yMCxyini = []
# optimization parameters
# m-th order polynomial
norders = [0, 2, 4, 6]
coeffs = [0.0, 0.0, 1.0]
# Weak (anti)localization model
sigma00 = -1.0 # sheet conductance at B = 0 in S. For sigma00 < 0, initial sigma0 will be sought from MCxx
alpha0 = 1.0 # prefactor of weak (anti-)localization model. 1 for WL, -1 for AWL
mfp0 = 20.0 # Mean free path in nm
# Drude model
# hole density, hole mobility, electron density, electron mobility
nh0 = 3e23 # 2*1.35e24 # m^-3
muh0 = 0.006 # 2*0.200 # m^2/Vs
ne0 = 1.3e22 # 2*1.8e24 # m^-3
mue0 = 0.02 # 5*0.0300 # m^2/Vs
#ρxxとρxyの残差二乗和のweight。wxxは1.0に固定
wxx = 1.0
wxy = 1.0
# variable lists and initial parameter list
varname = ["sigma0", "alpha", "mfp", "nh", "muh", "ne", "mue"]
varunit = [ "S", "", "nm", "cm-3", "cm2/Vs", "cm-3", "cm2/Vs"]
ai0 = [ sigma00, alpha0, mfp0, nh0, muh0, ne0, mue0]
optid = [ 0, 0, 0, 0, 0, 0, 0]
#=============================================
# 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-4
h_diff = 1.0e-3
#=============================
# Graph variables
#=============================
fig = None
fontsize = 12
legend_fontsize = 16
graphupdateinterval = 10
# 実数値に変換できない文字列を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, nh = nh0, muh = muh0, ne = ne0, mue = mue0):
global Bmin, Bmax, nB
global file_xx, file_xy, thickness
global Bxxmin, Bxxmax, Bxymin, Bxymax
print("")
print("Usage: Variables in () are optional")
print(" MFP in nm B in T thickness in nm")
print(" model: combination of poly, wl, wal, two, same_mu, same_n, single")
print("")
print(" (i-a) python {} sim (model file_xx file_xy thickness sigma0 alpha mfp Bmin Bmax nB"
.format(sys.argv[0]))
print(" if file_xx and/or file_xy include valid B data, calculation will be performed for B[i]")
print(" Otherwise, Bmin, Bmax, and nB are used.")
print(" ex: python {} sim {} {} {} {} {} {} {} {} {} {}"
.format(sys.argv[0], model, '-', '-', thickness, sigma0, alpha, mfp, Bmin, Bmax, nB))
print(" ex: python {} sim {} {} {} {} {} {} {} {} {} {}"
.format(sys.argv[0], model, file_xx, file_xy, thickness, sigma0, alpha, mfp, Bmin, Bmax, nB))
print("")
print(" (i-b) python {} sim poly (file_xx file_xy thickness [a0 coeff0]*norders"
.format(sys.argv[0]))
print(" ex: python {} sim {} {} {} {} {} {} {} {}"
.format(sys.argv[0], 'poly', '-', '-', thickness, 0, 0.0, 2, 1.0))
print("")
print(" (ii-a) python {} fit poly (file_xx file_xy norder list)".format(sys.argv[0]))
print(" ex: python {} fit {} {} {} {} {} {} {}"
.format(sys.argv[0], 'poly', file_xx, file_xy, 0, 2, 4, 6))
print("")
print(" (ii-b) python {} fit (model file_xx file_xy thickness sigma0 id alpha id MFP id "
+ "nh id muh id ne id mue id wxx wxy Bxxmin Bxxmax maxiter tol)"
.format(sys.argv[0]))
print(" ex: python {} fit {} {} {} {}"
.format(sys.argv[0], model, file_xx, file_xy, thickness))
print(" ex: python {} fit {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {}"
.format(sys.argv[0], model, file_xx, file_xy, thickness,
sigma0, optid[0], alpha, optid[1], mfp, optid[2],
nh, optid[3], muh, optid[4], ne, optid[5], mue, optid[6],
wxx, wxy,
Bxxmin, Bxxmax, Bxymin, Bxymax, maxiter, tol))
print(" ex: python {} fit {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {}"
.format(sys.argv[0], 'wl+same_mu', file_xx, file_xy, thickness,
sigma0, optid[0], alpha, optid[1], mfp, optid[2],
nh, optid[3], muh, optid[4], ne, optid[5], mue, optid[6],
wxx, wxy,
Bxxmin, Bxxmax, Bxymin, Bxymax, maxiter, tol))
print(" ex: python {} fit {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {}"
.format(sys.argv[0], 'wl+same_n', file_xx, file_xy, thickness,
sigma0, optid[0], alpha, optid[1], mfp, optid[2],
nh, optid[3], muh, optid[4], ne, optid[5], mue, optid[6],
wxx, wxy,
Bxxmin, Bxxmax, Bxymin, Bxymax, maxiter, tol))
def terminate(message = None):
if message is not None:
print("")
print(message)
usage()
print("")
exit()
#=============================
# Treat argments
#=============================
# 起動時引数を使ってglobal変数を更新する
def updatevars():
global mode, model
global Bmin, Bmax, Bstep, nB
global file_xx, Bxxmin, Bxxmax
global file_xy, Bxymin, Bxymax
global thickness
global sigma00, alpha0, mfp0
global nh0, muh0, ne0, mue0
global wxx, wxy
global norders, coeffs
global maxiter, tol, graphupdateinterval
global ai0, optid
global output_fittingxx_csv, output_fittingxy_csv
global output_simulationxx_csv, output_simulationxy_csv
global output_txt
global output_bat
argv = sys.argv
if len(argv) == 1:
terminate()
mode = getarg(1, mode)
model = getarg(2, model)
if mode == 'sim':
file_xx = getarg ( 3, file_xx)
file_xy = getarg ( 4, file_xy)
thickness = getfloatarg( 5, thickness)
if model == 'poly':
n0 = getfloatarg( 6, None)
if n0 is not None:
norders = []
coeffs = []
i = 0
while 1:
n0 = getintarg( 6+2*i, None)
if n0 is None:
break
norders.append(n0)
coeffs.append(getfloatarg( 6+2*i+1, None))
i += 1
else:
sigma00 = getfloatarg( 6, sigma00)
alpha0 = getfloatarg( 7, alpha0)
mfp0 = getfloatarg( 8, mfp0)
Bmin = getfloatarg( 9, Bmin)
Bmax = getfloatarg(10, Bmax)
nB = getintarg (11, nB)
Bstep = (Bmax - Bmin) / (nB - 1)
elif mode == 'fit' and model == 'poly':
file_xx = getarg ( 3, file_xx)
file_xy = getarg ( 4, file_xy)
thickness = getfloatarg( 5, thickness)
n0 = getfloatarg( 6, None)
if n0 is not None:
norders = []
coeffs = []
i = 0
while 1:
n0 = getintarg( 6+i, None)
if n0 is None:
break
norders.append(n0)
coeffs.append(0.0)
i += 1
elif mode == 'fit':
file_xx = getarg ( 3, file_xx)
file_xy = getarg ( 4, file_xy)
thickness = getfloatarg( 5, thickness)
sigma00 = getfloatarg( 6, sigma00)
optid[0] = getintarg ( 7, optid[0])
alpha0 = getfloatarg( 8, alpha0)
optid[1] = getintarg ( 9, optid[1])
mfp0 = getfloatarg(10, mfp0)
optid[2] = getintarg (11, optid[2])
nh0 = getfloatarg(12, nh0)
optid[3] = getintarg (13, optid[3])
muh0 = getfloatarg(14, muh0)
optid[4] = getintarg (15, optid[4])
ne0 = getfloatarg(16, ne0)
optid[5] = getintarg (17, optid[5])
mue0 = getfloatarg(18, mue0)
optid[6] = getintarg (19, optid[6])
wxx = getfloatarg(20, wxx)
wxy = getfloatarg(21, wxy)
Bxxmin = getfloatarg(22, Bxxmin)
Bxxmax = getfloatarg(23, Bxxmax)
Bxymin = getfloatarg(24, Bxymin)
Bxymax = getfloatarg(25, Bxymax)
maxiter = getintarg (26, maxiter)
tol = getfloatarg(27, tol)
graphupdateinterval = getintarg(28, graphupdateinterval)
else:
terminate("Error: Invalid mode [{}]".format(mode))
if 'wl' in model:
alpha0 = abs(alpha0)
elif 'wal' in model:
alpha0 = -abs(alpha0)
if 'wl' in model or 'wal' in model:
if getarg( 7, None) is None:
optid[0] = 1
if getarg(11, None) is None:
optid[2] = 1
if 'two' in model:
sigma00 = 0.0
optid[0] = 0
if getarg(13, None) is None:
optid[3] = 1
if getarg(15, None) is None:
optid[4] = 1
if getarg(17, None) is None:
optid[5] = 1
if getarg(19, None) is None:
optid[6] = 1
if 'same_mu' in model:
sigma00 = 0.0
optid[0] = 0
if getarg(13, None) is None:
optid[3] = 1
if getarg(15, None) is None:
optid[4] = 1
if getarg(17, None) is None:
optid[5] = 1
optid[6] = 0
if 'same_n' in model:
print("same n")
sigma00 = 0.0
optid[0] = 0
if getarg(13, None) is None:
optid[3] = 1
if getarg(15, None) is None:
optid[4] = 1
if getarg(19, None) is None:
optid[6] = 1
optid[5] = 0
if 'single' in model:
sigma00 = 0.0
optid[0] = 0
if getarg(13, None) is None:
optid[3] = 1
if getarg(15, None) is None:
optid[4] = 1
optid[5] = 0
optid[6] = 0
ai0 = [ sigma00, alpha0, mfp0, nh0, muh0, ne0, mue0]
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_fittingxx_csv = os.path.join(dirname, header + '-xx-fitting.csv')
output_fittingxy_csv = os.path.join(dirname, header + '-xy-fitting.csv')
output_simulationxx_csv = os.path.join(dirname, header + '-xx-simulation.csv')
output_simulationxy_csv = os.path.join(dirname, header + '-xy-simulation.csv')
output_txt = os.path.join(dirname, header + '-fitting.txt')
output_bat = os.path.join(dirname, header + '.bat')
#=============================
# 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)
for j in range(len(header)):
data.append([])
for row in f:
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, aidef):
ai = []
c = 0
for i in range(len(optid)):
if optid[i] == 1:
ai.append(ais[c])
c += 1
else:
ai.append(aidef[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
flabel = ['1', 'sin(x)', 'cos(x)', 'sin(2x)', 'cos(2x)', 'sin(3x)', 'cos(3x)']
def lsqfunc(i, x):
global norders
try:
return pow(x, norders[i])
except:
print("Error in lsqfunc. Return 0.0")
return 0.0
def caly(x, norders, coeffs):
y = 0.0
for i in range(len(norders)):
y += coeffs[i] * lsqfunc(i, x)
return y
def mlsq(norders, x, y, IsPrint = 1):
m = len(norders)
n = len(x)
Si = np.empty([m, 1])
Sij = np.empty([m, m])
for l in range(0, m):
v = 0.0
for i in range(0, n):
v += y[i] * lsqfunc(l, x[i])
Si[l, 0] = v
for j in range(0, m):
for l in range(j, m):
v = 0.0
for i in range(0, n):
v += lsqfunc(j, x[i]) * lsqfunc(l, x[i])
Sij[j, l] = Sij[l, j] = v
ci = np.linalg.inv(Sij) @ Si
ci = ci.transpose().tolist()
if IsPrint:
print("Vector and Matrix:")
print("Si=\n", Si)
print("Sij=\n", Sij)
print("")
print("coeffs=", ci[0])
print("")
return ci[0]
def search_sigma0(B, sigma):
sigma0 = 0.0
minB = 1.0e10
for i in range(len(B)):
if abs(B[i]) < minB:
minB = abs(B[i])
sigma0 = sigma[i]
return sigma0
def plot(fig, Bxx, MRxx, MCxx, Bxy, MRxy, MRxxini, MCxxini, MRxyini, MRxxfit = None, MCxxfit = None, MRxyfit = None):
global plt
global graphupdateinterval
plt.clf()
ax11 = fig.add_subplot(2, 3, 1)
ax12 = fig.add_subplot(2, 3, 2)
ax13 = fig.add_subplot(2, 3, 3)
ax21 = fig.add_subplot(2, 3, 4)
ax22 = fig.add_subplot(2, 3, 5)
ax23 = fig.add_subplot(2, 3, 6)
if len(MRxx) > 0:
ax11.plot(Bxx, MRxx, label = 'raw', linestyle = 'none', marker = 'o', markersize = 1.0)
if len(MRxxini) > 0:
ax11.plot(Bxx, MRxxini, label = 'ini', linestyle = '-', color = 'red', linewidth = 1.0)
if MRxxfit:
ax11.plot(Bxx, MRxxfit, label = 'fit', linestyle = '-', color = 'blue', linewidth = 1.0)
if len(MRxxini) > 0:
ax12.plot(Bxx, MRxxini, label = 'ini', linestyle = '-', color = 'red', linewidth = 1.0)
if MCxx is not None and len(MCxx) > 0:
ax21.plot(Bxx, MCxx, label = 'raw', linestyle = 'none', marker = 'o', markersize = 1.0)
if MCxxini is not None and len(MCxxini) > 0:
ax21.plot(Bxx, MCxxini, label = 'ini', linestyle = '-', color = 'red', linewidth = 1.0)
if MCxxfit is not None and len(MCxxfit):
ax21.plot(Bxx, MCxxfit, label = 'fit', linestyle = '-', color = 'blue', linewidth = 1.0)
if len(MCxxini) > 0:
ax22.plot(Bxx, MCxxini, label = 'ini', linestyle = '-', color = 'red', linewidth = 1.0)
if len(MRxy) > 0:
ax13.plot(Bxy, MRxy, label = 'raw', linestyle = 'none', marker = 'o', markersize = 1.0)
if MRxyini is not None and len(MRxyini) > 0:
ax13.plot(Bxy, MRxyini, label = 'ini', linestyle = '-', color = 'red', linewidth = 1.0)
if MRxyfit is not None and len(MRxyfit) > 0:
ax13.plot(Bxy
, MRxyfit, label = 'fit', linestyle = '-', color = 'blue', linewidth = 1.0)
ax11.tick_params(labelsize = fontsize)
ax12.tick_params(labelsize = fontsize)
ax13.tick_params(labelsize = fontsize)
ax21.tick_params(labelsize = fontsize)
ax22.tick_params(labelsize = fontsize)
ax23.tick_params(labelsize = fontsize)
ax11.set_xlabel("B(T)")
ax11.set_ylabel("rhoxx(ohm*m)")
ax11.set_title(file_xx)
ax12.set_xlabel("B(T)")
ax12.set_ylabel("rhoxx(ohm*m)")
ax13.set_xlabel("B(T)")
ax13.set_ylabel("rhoxy(ohm*m)")
ax21.set_xlabel("B(T)")
ax21.set_ylabel("Ssheet,xx(S)")
ax21.set_title(file_xx)
ax22.set_xlabel("B(T)")
ax22.set_ylabel("Ssheet,xx(S)")
# ax3.set_title("l_in = {:12.6g} nm".format(mfp))
ax11.legend(fontsize = legend_fontsize)
ax12.legend(fontsize = legend_fontsize)
ax13.legend(fontsize = legend_fontsize)
ax21.legend(fontsize = legend_fontsize)
ax22.legend(fontsize = legend_fontsize)
# ax23.legend(fontsize = legend_fontsize)
plt.tight_layout()
plt.pause(0.1)
def rho_twocarrier(B, nh, muh, ne, mue):
global e
B2 = B * B
muh2mue2 = muh * muh * mue * mue
nhmuhnemue = nh * muh + ne * mue
nhmuenemuh = nh * mue + ne * muh
ndiff = nh - ne
nmu2diff = nh * muh * muh - ne * mue * mue
a1 = nhmuhnemue + nhmuenemuh * muh * mue * B2
a2 = nhmuhnemue * nhmuhnemue + ndiff * ndiff * muh2mue2 * B2
rhoxx = 1.0 / e * a1 / a2
rhoxx0 = 1.0 / e / nhmuhnemue
MRxx = rhoxx - rhoxx0
a3a = nmu2diff # nh*muh*muh - ne*mue*mue
a3b = ndiff * muh2mue2 * B2 # (nh - ne) * muh*muh*mue*mue * B2
a4a = nhmuhnemue * nhmuhnemue # pow(nh*muh + ne*mue, 2)
a4b = ndiff * ndiff * muh2mue2 * B2 #pow(nh - ne, 2) * muh*muh*mue*mue * B2
a3 = a3a + a3b
a4 = a4a + a4b
MRxy = B / e * a3 / a4
return rhoxx, rhoxx0, MRxx, MRxy
def rho_singlecarrier(B, n, mu):
global e
B2 = B * B
sigmaxx0 = e * n * mu
rhoxx0 = 1.0 / sigmaxx0
rhoxx = rhoxx0 * (1.0 + mu * mu * B2)
MRxx = rhoxx - rhoxx0
MRxy = B / e / n
return rhoxx, rhoxx0, MRxx, MRxy
def rho_samemu(B, nh, mu, ne):
return rho_twocarrier(B, nh, mu, ne, mu)
def rho_samen(B, n, muh, mue):
return rho_twocarrier(B, n, muh, n, mue)
# return in MKS
def rho_Drude(model, B, nh, muh, ne, mue):
if 'two' in model:
return rho_twocarrier(B,nh, muh, ne, mue)
elif 'single' in model:
return rho_singlecarrier(B, nh, muh)
elif 'same_mu' in model:
return rho_samemu(B, nh, muh, ne)
elif 'same_n' in model:
return rho_samen(B, nh, muh, mue)
print("")
print("Error in rho_Drude: Invalide model [{}]".format(model))
print("")
exit()
# 弱局在による磁気伝導率変化 Δσ(B)の計算。シート電導度
def sigma_WL(B, l_in, alpha):
x = l_in * l_in * 4.0 * e * abs(B) / hbar # MKS単位で無次元化
if abs(x) < xth:
x = xth
dsigma = scipy.special.digamma(0.5 + 1.0 / x) + log(x) #縦軸 e*e/pi/hbar単位
return alpha * kMR * dsigma
# シート電導度で返す
def cal_sigma(model, B, thickness, sigma0, alpha, mfp, nh, muh, ne, mue, norders, coeffs):
global nm2m
l_in = mfp * nm2m # m
sigma = 0.0
if 'wl' in model or 'wal' in model:
sigma += sigma0 + sigma_WL(B, l_in, alpha)
if 'two' in model or 'single' in model or 'same_mu' in model or 'same_n' in model:
rxx, rxx0, MR_xx, MR_xy = rho_Drude(model, B, nh, muh, ne, mue)
sigma += 1.0 / rxx * (thickness * nm2m)
if 'poly' in model:
sigma += caly(B, norders, coeffs)
return sigma
# in MKS
def cal_rhoxy_list(model, B, nh, muh, ne, mue, norders, coeffs, IsPrint = 1):
global MRxy
if IsPrint:
if len(MRxy) > 0:
print("{}\t{}\t{}".format("B (T)", "MRxy (ohm*m)(obs)", "MRxy (ohm*m)(cal)"))
else:
print("{}\t{}".format("B (T)", "MRxy (ohm*m)(cal)"))
MRxycal = []
for i in range(len(B)):
if 'two' in model or 'same_mu' in model or 'same_n' in model or 'single' in model:
rxx, rxx0, MR_xx, MR_xy = rho_Drude(model, B[i], nh, muh, ne, mue)
elif 'poly' in model:
MR_xy = caly(B[i], norders, coeffs)
# sigma_xy = caly(B[i], norders, coeffs)
# if sigma_xy == 0.0:
# MR_xy = 0.0
# else:
# MR_xy = 1.0 / sigma_xy * (thickness * nm2m)
# else:
# MR_xy = 0.0
MRxycal.append(MR_xy)
if IsPrint:
if len(MRxy) > 0:
print("{}\t{}\t{}".format(B[i], MRxy[i], MRxycal[i]))
else:
print("{}\t{}".format(B[i], MRxycal[i]))
if IsPrint:
print("")
return MRxycal
# シート電導度
def cal_sigmaxx_list(model, Bin, thickness, sigma0, alpha, mfp, nh, muh, ne, mue, norders, coeffs, IsPrint = 1):
MCxxcal = []
MRxxcal = []
if IsPrint:
if len(MRxx) == 0:
print("{:4s}\t{:8s}\t{:10s}\t{:10s}"
.format('i', 'B(T)', 'rho(cal)(ohm*m)', 'Ssheet(cal)(S)'))
else:
print("{:4s}\t{:8s}\t{:10s}\t{:10s}\t{:10s}\t{:10s}"
.format('i', 'B(T)', 'rho(exp)(ohm*m)', 'rho(cal)(ohm*m)', 'Ssheet(exp)(S)', 'Ssheet(cal)(S)'))
for i in range(len(Bin)):
B = Bin[i]
sigmacal = cal_sigma(model, B, thickness, sigma0, alpha, mfp, nh, muh, ne, mue, norders, coeffs)
if sigmacal == 0.0:
rhocal = 0.0
else:
rhocal = 1.0 / sigmacal * (thickness * nm2m)
if IsPrint:
if len(MRxx) == 0:
print("{:4d}\t{:8.4f}\t{:10.4g}\t{:10.4g}"
.format(i, B, rhocal, sigmacal))
else:
print("{:4d}\t{:8.4f}\t{:10.4g}\t{:10.4g}\t{:10.4g}\t{:10.4g}"
.format(i, B, MRxx[i], rhocal, MCxx[i], sigmacal))
MCxxcal.append(sigmacal)
MRxxcal.append(rhocal)
return MRxxcal, MCxxcal
def calS2(xk):
global model
global Bxx, MRxx, MCxx
global Bxy, MRxy
global alpha0
global ai0, optid
sigma0, alpha, mfp, nh, muh, ne, mue = recover_parameters(xk, optid, ai0)
sigma0 = max([0.0, sigma0])
mfp = max([0.0, mfp])
l_in = mfp * nm2m
S2 = 0.0
for i in range(len(Bxx)):
sigmacal = cal_sigma(model, Bxx[i], thickness, sigma0, alpha, mfp, nh, muh, ne, mue)
S2 += wxx * pow(sigmacal - MCxx[i], 2)
if 'two' in model or 'same_mu' in model or 'same_mu' in model or 'single' in model:
for i in range(len(Bxy)):
rxx, rxx0, MR_xx, MR_xy = rho_Drude(model, Bxy[i], nh, muh, ne, mue)
S2 += wxy * pow(MR_xy - MRxy[i], 2)
return S2 / (wxx * len(Bxx) + wxy * len(Bxy))
# 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
global model
global ai0, optid
global Bxx, MRxx, MCxx
global Bxy, MRxy
global yMRxxini, yMCxxini, yMRxyini
S2 = calS2(xk)
sigma0, alpha, mfp, nh, muh, ne, mue = recover_parameters(xk, optid, ai0)
if 'same_mu' in model:
mue = muh
if 'same_n' in model:
ne = nh
if 'single' in model:
ne = 0.0
mue = 0.0
# パラメータと残差二乗和を表示
print("callback {}: sigma0={:14.4g} alpha={:6.3f} mfp={:8.4g} nh={:10.4g} muh={:8.4g} ne={:10.4g} mue={:8.4g} S2={:12.6g}"
.format(iter, sigma0, alpha, mfp, nh, muh, ne, mue, S2))
iter += 1
# Update graphs
if iter % graphupdateinterval == 0:
yMRxxfit, yMCxxfit = cal_sigmaxx_list(model, Bxx, thickness, sigma0, alpha, mfp, nh, muh, ne, mue, 0)
yMRxyfit = cal_rhoxy_list(model, Bxy, nh, muh, ne, mue, 0)
plot(fig, Bxx, MRxx, MCxx, Bxy, MRxy, yMRxxini, yMCxxini, yMRxyini, yMRxxfit, yMCxxfit, yMRxyfit)
def fitting():
global file_xx
global Bxxmin, Bxxmax
global headerxx, Bxx, MRxx, MCxx
global yMRxxini, yMCxxini, yMRxyini
global sigma00, alpha0, mfp0
global nh0, muh0, ne0, mue0
global ai0, optid
global fig
fig = plt.figure(figsize = (12, 8))
if 'poly' in model:
print("Optimization configuration")
method = 'linear LSQ'
print(" method: ", method)
else:
print("")
print("wxx: ", wxx)
print("wxy: ", wxy)
print("")
print("Optimization configuration")
print(" method: ", method)
print(" maxiter: ", maxiter)
print(" tol : ", tol)
print("")
print("Calculate initial data")
yMRxxini, yMCxxini = cal_sigmaxx_list(model, Bxx, thickness, sigma00, alpha0, mfp0, nh0, muh0, ne0, mue0)
yMRxyini = cal_rhoxy_list(model, Bxy, nh0, muh0, ne0, mue0)
# Plot graphs
plot(fig, Bxx, MRxx, MCxx, Bxy, MRxy, yMRxxini, yMCxxini, yMRxyini)
#=============================
# Optimization
#=============================
coeffsxx = []
coeffsxy = []
print("")
print("Optimization by {}".format(method))
if 'poly' in model:
MCxx = []
MCxy = []
for i in range(len(MRxx)):
# MCxx.append(1.0 / MRxx[i])
MCxx.append(1.0 / MRxx[i] * (thickness * nm2m)) # S/cm * thickness[cm]
for i in range(len(MRxy)):
# MCxy.append(1.0 / MRxy[i])
MCxy.append(1.0 / MRxy[i] * (thickness * nm2m)) # S/cm * thickness[cm]
coeffsxx = mlsq(norders, Bxx, MCxx)
coeffsxy = mlsq(norders, Bxy, MRxy)
fmin = 0.0
print("")
print("Optimized at S2={:12.6g}:".format(fmin))
print(" norders = ", norders)
print(" coeffsxx for sheet conductance = ", coeffsxx)
print(" coeffsxy for resistivity = ", coeffsxy)
sigma0 = 0.0
alpha, mfp, nh, muh, ne, mue = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0
for i in range(len(norders)):
if norders[i] == 0:
sigma0 = coeffsxx[i]
break
else:
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, nh, muh, ne, mue = recover_parameters(ai, optid, ai0)
if 'same_mu' in model:
mue = muh
if 'same_n' in model:
ne = nh
if 'single' in model:
ne = 0.0
mue = 0.0
print("")
print("Optimized at S2={:12.6g}:".format(fmin))
# print("sigma0: {} S".format(sigma0))
for i in range(len(varname)):
print(" {:10s}: optid={}: {}"
.format(varname[i], optid[i], [sigma0, alpha, mfp, nh, muh, ne, mue][i], varunit[i]))
print("")
print("Calculate final data")
yMRxxfit, yMCxxfit = cal_sigmaxx_list(model, Bxx, thickness,
sigma0, alpha, mfp, nh, muh, ne, mue, norders, coeffsxx)
yMRxyfit = cal_rhoxy_list(model, Bxy, nh0, muh0, ne0, mue0, norders, coeffsxy)
#=============================
# Save fitting results
#=============================
print("")
print("Save fitting results to [{}]".format(output_fittingxx_csv))
if 'poly' in model:
savecsv(output_fittingxx_csv,
["B(T)",
"rhoxx(exp)(ohm.cm)", "rhoxx(fin)(ohm.cm)",
"Ssheet(exp)(S)", "Ssheet(fin)(S)"],
[Bxx, MRxx, yMRxxfit, MCxx, yMCxxfit])
else:
savecsv(output_fittingxx_csv,
["B(T)",
"rhoxx(exp)(ohm.cm)", "rhoxx(ini)(ohm.cm)", "rhoxx(fin)(ohm.cm)",
"Ssheet(exp)(S)", "Ssheet(ini)(S)", "Ssheet(fin)(S)"],
[Bxx, MRxx, yMRxxini, yMRxxfit, MCxx, yMCxxini, 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))
for i in range(len(varname)):
f.write("{}={} {}\n".format(varname[i], [sigma0, alpha, mfp, nh, muh, ne, mue][i], varunit[i]))
f.close()
#=============================
# Plot graphs
#=============================
print("")
print("plot")
plot(fig, Bxx, MRxx, MCxx, Bxy, MRxy, yMRxxini, yMCxxini, yMRxyini, yMRxxfit, yMCxxfit, yMRxyfit)
print("")
usage()
print("")
s = "python {} fit {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {}".format(
sys.argv[0], model, file_xx, file_xy, thickness,
"{:12.6g}".format(sigma0).strip(), optid[0],
"{:12.6g}".format(alpha).strip(), optid[1],
"{:12.6g}".format(mfp).strip(), optid[2],
"{:12.6g}".format(nh).strip(), optid[3],
"{:12.6g}".format(muh).strip(), optid[4],
"{:12.6g}".format(ne).strip(), optid[5],
"{:12.6g}".format(mue).strip(), optid[6],
wxx, wxy, Bxxmin, Bxxmax, Bxymin, Bxymax, maxiter, tol)
print(" Saved to [{}]".format(output_bat))
f = open(output_bat, 'w')
f.write(s + "\n")
f.close()
print("Next suggestion:")
print(s)
print("")
print("Press ENTER to exit>>", end = '')
input()
exit()
def simulate():
global nm2m
global mode, model
global Bmin, Bmax, Bstep, nB
global file_xx, headerxx, Bxx, MRxx, MCxy
global file_xy, headerxy, Bxy, MRxy, MCxy
global sigma00, alpha0, mfp0, nh0, muh0, ne0, mue0
global norders, coeffs
print("")
print("Calculate conductivity(B)")
print("B range: {} - {} at {} step, {} points".format(Bmin, Bmax, Bstep, nB))
if len(Bxx) > 0:
xBxx = Bxx
else:
xBxx = [Bmin + i * Bstep for i in range(nB)]
if len(Bxy) > 0:
xBxy = Bxy
else:
xBxy = [Bmin + i * Bstep for i in range(nB)]
l_in = mfp0 * nm2m
xx = [l_in * l_in * 4.0 * e * xBxx[i] / hbar for i in range(len(xBxx))]
yMRxxini, yMCxxini = cal_sigmaxx_list(model, xBxx, thickness,
sigma00, alpha0, mfp0, nh0, muh0, ne0, mue0,
norders, coeffs)
yMRxyini = cal_rhoxy_list(model, xBxy, nh0, muh0, ne0, mue0, norders, coeffs)
#=============================
# Save simulation results
#=============================
print("")
print("Save simulation results to [{}]".format(output_simulationxx_csv))
savecsv(output_simulationxx_csv,
["B(T)", "rhoxx(ini)(ohm.cm)", "Ssheet(ini)(S)"],
[xBxx, yMRxxini, yMCxxini])
if len(Bxy) > 0:
print("")
print("Save simulation results to [{}]".format(output_simulationxy_csv))
savecsv(output_simulationxy_csv,
["B(T)", "rhoxy(ini)(ohm.cm)"],
[xBxy, yMRxyini])
print("")
print("plot")
fig = plt.figure(figsize = (12, 8))
plot(fig, xBxx, MRxx, MCxx, xBxy, MRxy, yMRxxini, yMCxxini, yMRxyini)
print("")
usage()
print("")
print("Press ENTER to exit>>", end = '')
input()
def main():
global mode, model
global Bmin, Bmax, Bstep, nB
global file_xx, headerxx, Bxx, MRxx, MCxx
global file_xy, headerxy, Bxy, MRxy, MCxy
global sigma00, alpha0, mfp0, nh0, muh0, ne0, mue0
global norders, coeffs
updatevars()
print("")
print("mode : ", mode)
print("model: ", model)
print("")
print("Measured MR data (the resistivity in file are in ohm*m)")
print(" thickness: {} nm".format(thickness))
print(" Read [{}]".format(file_xx))
print(" Bxx range: {} - {}".format(Bxxmin, Bxxmax))
try:
headerxx, Bxx, MRxx = read_csv(file_xx, Bxxmin, Bxxmax)
print(" ndata: ", len(Bxx))
except:
print(" Warning: Could not read file [{}]".format(file_xx))
print(" Read [{}]".format(file_xy))
print(" Bxy range: {} - {}".format(Bxymin, Bxymax))
try:
headerxy, Bxy, MRxy = read_csv(file_xy, Bxymin, Bxymax)
print(" ndata: ", len(Bxy))
except:
print(" Warning: Could not read file [{}]".format(file_xy))
MCxx = []
for i in range(len(Bxx)):
MCxx.append(1.0 / MRxx[i] * (thickness * nm2m)) # S/cm * thickness[cm]
if sigma00 < 0.0:
sigma00 = search_sigma0(Bxx, MCxx)
print("")
print("sigma00 is updated from MRxx: {} S".format(sigma00))
ai0[0] = sigma00
print("")
print("parameters for model [{}]:".format(model))
if model == 'poly':
for i in range(len(norders)):
print(" {}: order = {} coeff = {}".format(i, norders[i], coeffs[i]))
else:
for i in range(len(varname)):
print(" {:10s}: optid={}: {}".format(varname[i], optid[i], ai0[i], varunit[i]))
if mode == 'sim':
simulate()
elif mode == 'fit':
fitting()
else:
terminate("Error in main: Invalide mode [{}]".format(mode))
terminate()
if __name__ == '__main__':
main()
