diff --git a/examples/Cu/AMR-xy/AMR_rhotheta.py b/examples/Cu/AMR-xy/AMR_rhotheta.py index 86a16209..834f3537 100644 --- a/examples/Cu/AMR-xy/AMR_rhotheta.py +++ b/examples/Cu/AMR-xy/AMR_rhotheta.py @@ -28,7 +28,7 @@ def extract_rho_theta(f,dtheta,num_show,interval_line,choose,output_file, output if __name__ == '__main__': - theta_interval=5 + theta_interval=15 Btau_show=6 Btau_num=101 Btau_max=10 diff --git a/examples/Cu/AMR-xy/AMR_rhoxx.gnu b/examples/Cu/AMR-xy/AMR_rhoxx.gnu index 260ba0bd..717d6555 100644 --- a/examples/Cu/AMR-xy/AMR_rhoxx.gnu +++ b/examples/Cu/AMR-xy/AMR_rhoxx.gnu @@ -7,7 +7,7 @@ set ylabel '{/Symbol r}_{xx}*{/Symbol t} ({/Symbol W}*m*s)' set format y "%1.1e" set xlabel '{/Symbol \161}' set xrange [0:180] -set yrange [1.2e-21:2.5e-21] +#set yrange [1.2e-21:2.5e-21] set xtics 30 set key outside #unset key diff --git a/examples/Cu/AMR-xy/AMR_rhozz.gnu b/examples/Cu/AMR-xy/AMR_rhozz.gnu index 9e6f558c..6a5c4131 100644 --- a/examples/Cu/AMR-xy/AMR_rhozz.gnu +++ b/examples/Cu/AMR-xy/AMR_rhozz.gnu @@ -7,7 +7,7 @@ set ylabel '{/Symbol r}_{zz}*{/Symbol t} ({/Symbol W}*m*s)' set format y "%1.1e" set xlabel '{/Symbol \161}' set xrange [0:180] -set yrange [1.2e-21:2.7e-21] +#set yrange [1.2e-21:2.7e-21] set xtics 30 set key outside #unset key diff --git a/examples/Cu/AMR-xy/xyplane.sh b/examples/Cu/AMR-xy/xyplane.sh index 9d58a254..a6156404 100755 --- a/examples/Cu/AMR-xy/xyplane.sh +++ b/examples/Cu/AMR-xy/xyplane.sh @@ -3,11 +3,11 @@ # alpha is the angle between the magnetic field and the z' axis in the z'-b plane, where z' axis is # perpendicular to a and b axis. -for ((iphi=0; iphi<=36; iphi++)) +for ((iphi=0; iphi<=12; iphi++)) do theta=90 -phi=`echo "$iphi*5"|bc` +phi=`echo "$iphi*15"|bc` dir='Btheta'$theta'Bphi'$phi echo $theta $phi $dir mkdir $dir @@ -35,9 +35,9 @@ OmegaNum = 1 ! omega number OmegaMin = 0 ! energy interval OmegaMax = 0 ! energy interval E_i= OmegaMin+ (OmegaMax-OmegaMin)/(OmegaNum-1)*(i-1) EF_integral_range = 0.05 ! in eV, a broadening factor to choose the k points for integration -Nk1 =81 ! Kmesh(1) for KCUBE_BULK -Nk2 =81 ! Kmesh(2) for KCUBE_BULK -Nk3 =81 ! Kmesh(3) for KCUBE_BULK +Nk1 =41 ! Kmesh(1) for KCUBE_BULK +Nk2 =41 ! Kmesh(2) for KCUBE_BULK +Nk3 =41 ! Kmesh(3) for KCUBE_BULK BTauNum= 101 ! Number of B*tau we calculate BTauMax = 10.0 ! The maximum B*tau, starting from Btau=0. Tmin = 30 ! Temperature in Kelvin @@ -78,39 +78,22 @@ KCUBE_BULK EOF cat>$dir/wt-theta.sh<15s}{:>16s}{:>16s}{:>16s}{:>16s}{:>16s}{:>16s}{:>16s}{:>16s}{:>16s}\n'.format('BTau','xx','xy','xz','yx','yy','yz','zx','zy','zz')) + # for ibtau, btau in enumerate(new_btaulist): + # row = [btau] + [sigma_bands[imu, iT, ibtau, ib, i] for i in range(9)] + # sbandfile.write('{:>15.6f}{:>16.6e}{:>16.6e}{:>16.6e}{:>16.6e}{:>16.6e}{:>16.6e}{:>16.6e}{:>16.6e}{:>16.6e}\n'.format(*row)) + + # store Tlist, Omegalist, btaulist and sigma_bands in sigma_bands.npy in the form of dictionary + np.save('sigma_bands.npy',{'Tlist':Tlist, 'Omegalist':Omegalist, 'btaulist':new_btaulist, 'sigma_bands':sigma_bands}, allow_pickle=True) + + + # write sigma and rhotau file + rho_all = np.zeros((OmegaNum, NumT, BTauNum, 9)) + + for imu, mu in enumerate(Omegalist): + # sum the conductivity and resistivity of all bands and write them in + # sigma_total_mu_*.dat and rhotau_total_mu_*.dat + + sigmafile = open(f'sigma_total_mu_{mu:.3f}eV.dat', 'w') + rhofile = open(f'rhotau_total_mu_{mu:.3f}eV.dat', 'w') + for iT, T in enumerate(Tlist): + sigmafile.write(f'# T = {T}K\n') + sigmafile.write('#{:>15s}{:>16s}{:>16s}{:>16s}{:>16s}{:>16s}{:>16s}{:>16s}{:>16s}{:>16s}\n'.format('BTau','xx','xy','xz','yx','yy','yz','zx','zy','zz')) + rhofile.write(f'# T = {T}K\n') + rhofile.write('#{:>15s}{:>16s}{:>16s}{:>16s}{:>16s}{:>16s}{:>16s}{:>16s}{:>16s}{:>16s}\n'.format('BTau','xx','xy','xz','yx','yy','yz','zx','zy','zz')) + for ibtau, btau in enumerate(new_btaulist): + sigma = sum(sigma_bands[imu, iT, ibtau]) + row = [btau] + [sigma[i] for i in range(9)] + sigmafile.write('{:>15.6f}{:>16.6e}{:>16.6e}{:>16.6e}{:>16.6e}{:>16.6e}{:>16.6e}{:>16.6e}{:>16.6e}{:>16.6e}\n'.format(*row)) + if abs(np.linalg.det(sigma.reshape(3,3)))>1e-9: + rho = np.linalg.inv(sigma.reshape(3,3)) + rho_all[imu, iT, ibtau] = rho.reshape(1,9)[0] + row = [btau] + [rho.reshape(1,9)[0][i] for i in range(9)] + rhofile.write('{:>15.6f}{:>16.6e}{:>16.6e}{:>16.6e}{:>16.6e}{:>16.6e}{:>16.6e}{:>16.6e}{:>16.6e}{:>16.6e}\n'.format(*row)) + else: + rhofile.write('# error: sigma is zero since no k points contribute to the calculations of MR\n') + rho_all[imu, iT, ibtau] = np.NaN + sigmafile.write('\n') + rhofile.write('\n') + # store Tlist, Omegalist, btaulist and rho_all in rho_all.npy in the form of dictionary + np.save('rho_all.npy', {'Tlist':Tlist, 'Omegalist':Omegalist, 'btaulist':new_btaulist, 'rho_all':rho_all}, allow_pickle=True) + +def plot_cbar(component=['xx']): + + Ncomp = len(component) + + # read Tlist, Omegalist, btaulist and rho_all stored as dictionary from rho_all.npy + # Note that the true chemical potential is mu = Omega + 0.04 eV + Omegalist = np.load('rho_all.npy', allow_pickle=True).item()['Omegalist'] + Tlist = np.load('rho_all.npy', allow_pickle=True).item()['Tlist'] + btaulist = np.load('rho_all.npy', allow_pickle=True).item()['btaulist'] + rho_all = np.load('rho_all.npy', allow_pickle=True).item()['rho_all'] + + OmegaNum = len(Omegalist) + NumT = len(Tlist) + BTauNum = len(btaulist) + + print('Omegalist', Omegalist) + print('here debug', OmegaNum, NumT, BTauNum) + rho = np.zeros((OmegaNum, Ncomp, NumT, BTauNum)) + for imu, mu in enumerate(Omegalist): + for iT, T in enumerate(Tlist): + icomp = 0 + if 'xx' in component: + rho[imu, icomp, iT, :] = rho_all[imu, iT, :, 0] + icomp += 1 + if 'xy' in component: + rho[imu, icomp, iT, :] = rho_all[imu, iT, :, 1] + icomp += 1 + if 'xz' in component: + rho[imu, icomp, iT, :] = rho_all[imu, iT, :, 2] + icomp += 1 + if 'yx' in component: + rho[imu, icomp, iT, :] = rho_all[imu, iT, :, 3] + icomp += 1 + if 'yy' in component: + rho[imu, icomp, iT, :] = rho_all[imu, iT, :, 4] + icomp += 1 + if 'yz' in component: + rho[imu, icomp, iT, :] = rho_all[imu, iT, :, 5] + icomp += 1 + if 'zx' in component: + rho[imu, icomp, iT, :] = rho_all[imu, iT, :, 6] + icomp += 1 + if 'zy' in component: + rho[imu, icomp, iT, :] = rho_all[imu, iT, :, 7] + icomp += 1 + if 'zz' in component: + rho[imu, icomp, iT, :] = rho_all[imu, iT, :, 8] + icomp += 1 + assert icomp == Ncomp + + # The subplots will have the same width + # and there will be a horizontal space of 0.4 between the subplots + fig, axs = plt.subplots(figsize=(12,5), nrows=1, ncols=Ncomp,gridspec_kw={'width_ratios': [1 for i in range(Ncomp)], 'wspace': 0.4}) + print('plot data at mu = {:.3f} eV'.format(mu)) + + cmap = plt.get_cmap('jet',len(Tlist)) + for icomp, comp in enumerate(component): + for iT, T in enumerate(Tlist): + if T: + if Ncomp >1: + if comp == 'xx' or comp == 'yy' or comp =='zz': + + # uncomment this to plot resistivity + axs[icomp].plot(btaulist, rho[imu, icomp, iT, :], linewidth=2.5, color=cmap(iT), label='T=%.2f K'%(T)) + # plot MR + #axs[icomp].plot(btaulist, (rho[imu, icomp, iT, :]-rho[imu, icomp, iT, 0])/rho[imu, icomp, iT, 0]*100, linewidth=2.5, color=cmap(iT), label='T=%.2f K'%(T)) + else: + + # uncomment this to plot resistivity + axs[icomp].plot(btaulist, rho[imu, icomp, iT, :], linewidth=2.5, color=cmap(iT), label='T=%.2f K'%(T)) + # plot rho_xy-rho_xy(0) + # axs[icomp].plot(btaulist, rho[imu, icomp, iT, :]-rho[imu, icomp, iT, 0], linewidth=2.5, color=cmap(iT), label='T=%.2f K'%(T)) + else: + if comp == 'xx' or comp == 'yy' or comp =='zz': + axs.plot(btaulist, (rho[imu, icomp, iT, :]-rho[imu, icomp, iT, 0])/rho[imu, icomp, iT, 0]*100, linewidth=2.5, color=cmap(iT), label='T=%.2f K'%(T)) + else: + axs.plot(btaulist, rho[imu, icomp, iT, :]-rho[imu, icomp, iT, 0], linewidth=2.5, color=cmap(iT), label='T=%.2f K'%(T)) + + if Ncomp >1: + axs[icomp].set_xlabel(r'$B\tau\ (T\cdot ps)$') + if comp == 'xx' or comp == 'yy' or comp =='zz': + axs[icomp].set_ylabel(r'$ \rho_{%s}\tau\ (\Omega\cdot m\cdot s) $'%(comp)) + #axs[icomp].set_ylabel(r'$ MR_{%s} $'%(comp)) + else: + axs[icomp].set_ylabel(r'$ \rho_{%s}\tau\ (\Omega\cdot m\cdot s) $'%(comp)) + #axs[icomp].set_ylabel(r'$ (\rho_{%s}-\rho_0)\ (\Omega\cdot m) $'%(comp)) + else: + axs.set_xlabel(r'$B\tau\ (T\cdot ps)$') + if comp == 'xx' or comp == 'yy' or comp =='zz': + axs.set_ylabel(r'$ MR_{%s} $'%(comp)) + else: + axs.set_ylabel(r'$ (\rho_{%s}-\rho_0)\tau\ (\Omega\cdot m\cdot s) $'%(comp)) + + # Add title + fig.suptitle(r'$\mu=%.0f$ meV'%(mu*1000), ha='center', y=0.95) + + # Adjust the subplot margins + plt.subplots_adjust(top=0.8) + plt.subplots_adjust(bottom=0.15) + + norm = mpl.colors.Normalize(vmin=Tlist[0],vmax=Tlist[-1]) + sm = plt.cm.ScalarMappable(cmap=cmap, norm=norm) + sm.set_array([]) + clb = plt.colorbar(sm, ticks=np.linspace(Tlist[0], Tlist[-1], 5)) + clb.ax.set_title('T(K)') + plt.savefig('%.0fmeV.pdf'%((mu+0.000001)*1000), bbox_inches='tight') + #plt.show() + plt.close() + +def interpolate(component=['xx']): + # component is a list of components to be plotted + Ncomp = len(component) + + # read Tlist, Omegalist, btaulist, rho_all from rho_all.npy + # rho_all is a 4D array with shape (len(Omegalist), len(Tlist), len(btaulist), 9) + Omegalist = np.load('rho_all.npy', allow_pickle=True).item()['Omegalist'] + Tlist_interpolated = np.load('rho_all.npy', allow_pickle=True).item()['Tlist'] + btaulist = np.load('rho_all.npy', allow_pickle=True).item()['btaulist'] + rho_all = np.load('rho_all.npy', allow_pickle=True).item()['rho_all'] + NBtau = len(btaulist) + # print(NumT, NOmega, NBtau) + + # read concentration, Tlist, mu from fixdoping_mu.npy + # mulist is a 2D array with shape (len(concenlist), len(Tlist)) + concenlist = np.load('fixdoping_mu.npy', allow_pickle=True).item()['concentration'] + Tlist = np.load('fixdoping_mu.npy', allow_pickle=True).item()['Tlist'] + mulist = np.load('fixdoping_mu.npy', allow_pickle=True).item()['mu'] + # choose T in Tlist which is in the range of Tlist_interpolated + Tlist = [T for T in Tlist if (T > min(Tlist_interpolated) and T < max(Tlist_interpolated))] + Nconcen, NumT = len(concenlist), len(Tlist) + + + for comp in component: + # make sure the component is valid or print error message + if comp not in ['xx', 'xy', 'xz', 'yx', 'yy', 'yz', 'zx', 'zy', 'zz']: + print('Invalid component: %s'%(comp)) + return + + # generate the array for the component and use 'comp' as the prefix of the array name + if comp == 'xx': + rho_comp = rho_all[:, :, :, 0] + elif comp == 'xy': + rho_comp = rho_all[:, :, :, 1] + elif comp == 'xz': + rho_comp = rho_all[:, :, :, 2] + elif comp == 'yx': + rho_comp = rho_all[:, :, :, 3] + elif comp == 'yy': + rho_comp = rho_all[:, :, :, 4] + elif comp == 'yz': + rho_comp = rho_all[:, :, :, 5] + elif comp == 'zx': + rho_comp = rho_all[:, :, :, 6] + elif comp == 'zy': + rho_comp = rho_all[:, :, :, 7] + elif comp == 'zz': + rho_comp = rho_all[:, :, :, 8] + + # # replace the NaN in rho_comp with interpolated values + # if np.isnan(rho_comp).any(): + # print('NaN detected in {} component'.format(comp)) + # nan_mask = np.isnan(rho_comp) + # # nan_idx = np.where(nan_mask) + # # for idx in zip(nan_idx[0],nan_idx[1],nan_idx[2]): + # # print ('Omega, T, btau = ', '{:.3f}'.format(Omegalist[idx[0]]), '{:.3f}'.format(Tlist[idx[1]]), '{:.3f}'.format(btaulist[idx[2]])) + + # rho_comp[nan_mask] = np.interp(np.flatnonzero(nan_mask), np.flatnonzero(~nan_mask), rho_comp[~nan_mask]) + # # print(rho_comp[nan_mask]) + + + ################ define the interpolation function############# + # print(np.shape(rho_comp)) + interpolator = RegularGridInterpolator((Omegalist, Tlist_interpolated, btaulist), rho_comp) + ################ interpolate the resistivity at every concentration and temperature + ################ store the data in file 'rho{}_fixdoping.dat'.format{comp} + with open('rho{}_fixdoping.dat'.format(comp), 'w') as f: + print('write data of rho{}'.format(comp)) + for icon in range(Nconcen): + concentration = concenlist[icon] + # write a comment line to + print('write rho{} at n = {:.3e}'.format(comp, concentration)) + f.write('# n = {:.3e}\n'.format(concentration)) + for it in range(NumT): + T = Tlist[it] + mu = mulist[icon, it] + f.write('# T = {:.3f} K, mu = {:.3f} eV\n'.format(T,mu)) + f.write('#{:>15s}{:>16s}\n'.format('BTau',comp)) + for btau in btaulist: + rho_interp= interpolator((mu, T, btau)) + f.write('{:>15.6f}{:>16.6e}\n'.format(btau, rho_interp)) + + +def plot_interpolate(component=['xx']): + # component is a list of components to be plotted + Ncomp = len(component) + + # read Tlist, Omegalist, btaulist, rho_all from rho_all.npy + # rho_all is a 4D array with shape (len(Omegalist), len(Tlist), len(btaulist), 9) + Omegalist = np.load('rho_all.npy', allow_pickle=True).item()['Omegalist'] + Tlist_interpolated = np.load('rho_all.npy', allow_pickle=True).item()['Tlist'] + btaulist = np.load('rho_all.npy', allow_pickle=True).item()['btaulist'] + rho_all = np.load('rho_all.npy', allow_pickle=True).item()['rho_all'] + NBtau = len(btaulist) + # print(NumT, NOmega, NBtau) + + # read concentration, Tlist, mu from fixdoping_mu.npy + # mulist is a 2D array with shape (len(concenlist), len(Tlist)) + concenlist = np.load('fixdoping_mu.npy', allow_pickle=True).item()['concentration'] + Tlist = np.load('fixdoping_mu.npy', allow_pickle=True).item()['Tlist'] + mulist = np.load('fixdoping_mu.npy', allow_pickle=True).item()['mu'] + # choose T in Tlist which is in the range of Tlist_interpolated + Tlist = [T for T in Tlist if (T > min(Tlist_interpolated) and T < max(Tlist_interpolated))] + Nconcen, NumT = len(concenlist), len(Tlist) + + + ################# plot the resistivity at every concentration and temperature + ################# store the data in file 'rho{}_fixdoping.pdf'.format{comp} + cmap = plt.get_cmap('jet',len(Tlist)) + for icon in range(Nconcen): + print('plot data at n = {:.3e}'.format(concenlist[icon])) + concentration = concenlist[icon] + # plt.figure() + fig, axs = plt.subplots(figsize=(13,4), nrows=1, ncols=Ncomp+1 , + gridspec_kw={'width_ratios': [1 for i in range(Ncomp+1)], 'wspace': 0.5}) + fig.suptitle(r'$n = {:.1e}$'.format(concentration)+' cm$^{-3}$') + for icomp, comp in enumerate(component): + muvsT = [] + for it in range(NumT): + T = Tlist[it] + muvsT.append(mulist[icon, it]) + # print(1+icon*(NumT*(NBtau+2)+1)+it*(NBtau+2)) + rho = np.loadtxt('rho{}_fixdoping.dat'.format(comp), comments='#', skiprows=1+icon*(NumT*(NBtau+2)+1)+it*(NBtau+2), max_rows=len(btaulist)) + ############### the first subplot is the resistivity vs BTau + # axs[icomp].plot(rho[:,0], rho[:,1], color=cmap(it),label='T = {:.3f} K'.format(T)) + ############### the first subplot is the MR vs BTau + if 0:# comp == 'xx' or comp =='yy' or comp =='zz': + axs[icomp].plot(np.array(rho[:,0])*T, (rho[:,1]-rho[0,1])/rho[0,1], color=cmap(it),label='T = {:.3f} K'.format(T)) + else: + axs[icomp].plot(np.array(rho[:,0])*T, (rho[:,1]-rho[0,1]), color=cmap(it),label='T = {:.3f} K'.format(T)) + + if comp == 'xx' or comp =='yy' or comp =='zz': + # the xlabel of the first subplot is BTau, the xlabel of the second subplot is T(K) + axs[icomp].set_xlabel('BTau') + # axs[icomp].set_ylabel('MR{}'.format(comp)) + axs[icomp].set_ylabel('MR{}'.format(comp)) + # set xlim to [0,1] + axs[icomp].set_xlim([0,20]) + axs[icomp].set_ylim([0,1e-17]) + else: + # the xlabel of the first subplot is BTau, the xlabel of the second subplot is T(K) + axs[icomp].set_xlabel('BTau') + # axs[icomp].set_ylabel('MR{}'.format(comp)) + axs[icomp].set_ylabel(r'$\rho_{}$'.format(comp)) + axs[icomp].set_ylim([0,1e-17]) + axs[icomp].set_xlim([0,10]) + + + # add colorbar to the first subplot + norm = mpl.colors.Normalize(vmin=Tlist[0],vmax=Tlist[-1]) + sm = plt.cm.ScalarMappable(cmap=cmap, norm=norm) + sm.set_array([]) + clb = plt.colorbar(sm, ticks=np.linspace(Tlist[0], Tlist[-1], 5)) + clb.ax.set_title('T(K)') + + # the second subplot is the chemical potential vs T + axs[Ncomp].plot(Tlist, muvsT, color='k') + axs[Ncomp].set_ylabel(r'$\mu$ (eV)') + axs[Ncomp].set_xlabel('T(K)') + + plt.savefig('rho_{:.3e}.pdf'.format( concentration), bbox_inches='tight') + # plt.savefig('MR_{:.3e}.pdf'.format( concentration), bbox_inches='tight') + #plt.show() + plt.close() + + +def read_BoltzTrap(): + + # open file 'case.trace_fixdoping' and read the chemical potential at each temperature at each concentration + # the file contains the data of ten concentrations and every concentration are temperature dependent + # the temperature varies from 5k to 300k with step 5k + Tlist = np.arange(5, 305, 5) + concentration = np.array([-10e18, -7e18, -5E18, -3e18, -1e18, 1e18, 3e18, 5e18, 7e18, 10e18]) + NumT = len(Tlist) + Ncon = len(concentration) + mu = np.zeros((Ncon, NumT)) + + + with open('case.trace_fixdoping', 'r') as f: + # the first line is the comment line, skip it + f = f.readlines()[1:] + + # generate the loop to read data of every concentration + for icon in range(Ncon): + + # the temperature-dependent icon-th concentration is stored in the [icon*(NumT+1):icon*(NumT+1)+NumT] rows + data = f[icon*(NumT+1):icon*(NumT+1)+NumT] + + # generate the loop to read data of every temperature + # the first column is the temperature, the tenth column is the chemical potential + # Note that the chemical potential is in unit of Rydberg, so we need to multiply it by 13.6057 to convert it to eV + for it in range(NumT): + mu[icon, it] = float(data[it].split()[9])*13.6057 + + # store the concentration, Tlist, and mu as a dictionary in the file 'fixdoping_mu.npy' + np.save('fixdoping_mu.npy', {'concentration':concentration, 'Tlist':Tlist, 'mu':mu}, allow_pickle=True) + + +# main program +if __name__ == '__main__': + # define the band list + band_list=[6] + tau_list= [1] + + # band_list=['band50'] + component=['xx','yx','yy'] + # component=['xx'] + # component=['xz'] + + # remove all npy files in the folder + os.system('rm -f *.npy') + + readsigma_writerho(band_list, tau_list) + plot_cbar(component=component) + # read_BoltzTrap() + # interpolate(component=component) + # plot_interpolate(component=component) diff --git a/examples/Cu/Readme.txt b/examples/Cu/Readme.txt index 7bae15a4..a5e88b18 100644 --- a/examples/Cu/Readme.txt +++ b/examples/Cu/Readme.txt @@ -1,32 +1,65 @@ -# Added on Sep.05.2019 By QuanSheng Wu -This is a example for calculating ordinary magnetoresistance with given magnetic field direction. -0. unzip the hr.dat.tar.gz - $ tar xzvf wannier90_hr.dat_nsymm48.tar.gz -1. Calculate band structure and Fermi surface - $ cp wt.in-bands wt.in - $ mpiexec -np 4 wt.x& - $ gnuplot bulkek.gnu - $ xcrysden --bxsf FS3D.bxsf & - -2. Calculate MR for different magnetic fields along Theta=0, 18, 30, 45 degree and Phi=90. - $ cp wt.in-OHE-theta0 wt.in - $ mpiexec -np 4 wt.x & - $ cp sigma_band_6_mu_0.00eV_T_30.00K.dat sigma_band_6_mu_0.00eV_T_30.00K.dat-theta0 - $ cp wt.in-OHE-theta18 wt.in - $ mpiexec -np 4 wt.x & - $ cp sigma_band_6_mu_0.00eV_T_30.00K.dat sigma_band_6_mu_0.00eV_T_30.00K.dat-theta18 - $ cp wt.in-OHE-theta30 wt.in - $ mpiexec -np 4 wt.x & - $ cp sigma_band_6_mu_0.00eV_T_30.00K.dat sigma_band_6_mu_0.00eV_T_30.00K.dat-theta30 - $ cp wt.in-OHE-theta45 wt.in - $ mpiexec -np 4 wt.x & - $ cp sigma_band_6_mu_0.00eV_T_30.00K.dat sigma_band_6_mu_0.00eV_T_30.00K.dat-theta45 - $ gnuplot rho.gnu - $ evince rho.pdf - - -3. Get the evolution of momentum k under magnetic field - $ cp wt.in-OHE-evolve wt.in - $ mpiexec -np 4 wt.x & - Then open evolve_band_6.dat file and check whether there are k points evolution. - +0. unzip the hr.dat.tar.gz +$ tar xzvf wannier90_hr.dat_nsymm48.tar.gz + +1. Calculate band structure and Fermi surface +$ mkdir band +$ cp wt.in-bands band/wt.in +$ cp wannier90_hr.dat_nsymm48 band/ +$ cd band +$ mpirun -np 4 wt.x & +$ gnuplot bulkek.gnu +$ xcrysden --bxsf FS3D.bxsf + +2. Calculate MR for different magnetic fields along Theta=0, 18, 30, 45 degree and Phi=90. +$ mkdir OHE-theta0 +$ cp wannier90_hr.dat_nsymm48 OHE-theta0/ +$ cp wt.in-OHE-theta0 OHE-theta0/wt.in +$ cd OHE-theta0 +$ mpirun -np 4 wt.x & +$ cp rho_total_mu_0.00eV.dat ../rho_total_mu_0.00eV.dat-theta0 +$ cd .. + +repeat...... + +$ mkdir OHE-theta45 +$ cp wannier90_hr.dat_nsymm48 OHE-theta45/ +$ cp wt.in-OHE-theta0 OHE-theta45/wt.in +$ cd OHE-theta45 +$ mpirun -np 4 wt.x & +$ cp rho_total_mu_0.00eV.dat ../rho_total_mu_0.00eV.dat-theta45 +$ cd .. +$ gnuplot rhotheta.gnu + +3. Get the bulk Fermi surface in a fixed k plane +$ mkdir FS-contour +$ cp wannier90_hr.dat_nsymm48 FS-contour/ +$ cp wt.in-FS-contour FS-contour/wt.in +$ cd FS-contour +$ mpirun -np 4 wt.x & +$ gnuplot fs_kplane.gnu + +4. Get the evolution of momentum k under magnetic field +$ mkdir OHE-evolve +$ cp wannier90_hr.dat_nsymm48 OHE-evolve/ +$ cp wt.in-OHE-evolve OHE-evolve/wt.in +$ cd OHE-evolve +$ mpirun -np 4 wt.x & +Then open evolve_band_6.dat file and check whether there are k points evolution. + +5. Calculate AMR +$ cd AMR-xy +$ cp ../wannier90_hr.dat_nsymm48 ./ +$ chmod 777 xyplane.sh +replace the content between "cat>\$dir/wt-theta.sh<