#!/usr/bin/python import numpy as np # This allows to use "np" instead of "numpy" import glob from scipy import interpolate from matplotlib import pyplot as plt from matplotlib.colors import LinearSegmentedColormap from mpl_toolkits.mplot3d import Axes3D from matplotlib import cm from matplotlib import gridspec # resolution and light dev = 100 # number of devisions ligth_intensity = 1.5 # ligth_intensity E [0.5(light) 3(dark)] #colN = [ (50,205,50), (130,130,255), (255,69,0), (245,222,179), (160,82,45) ] # 'limegreen', 'deeppink', 'orangered' colN = [ (0,0,250)] col = [(240,245,247), (131,131,139), (73,73,74) ] # 'limegreen', 'deeppink', 'orangered' (210,180,190), (150,120,120), (150,70,70) RGB r_Ru = 0.24 # ========================================================================= ### function to constrcut the coordinates of sphere ---------------------- def mysphere(radius, pos): global x, y, z u, v = np.mgrid[0:2*np.pi:40j, 0:np.pi/2:10j] sphx = (max(x)-min(x))/(max(y)-min(y))*radius*np.cos(u)*np.sin(v) + pos[0] sphy = radius*np.sin(u)*np.sin(v) + pos[1] sphz = radius*np.cos(v) return sphx, sphy, sphz ### function to find the coordinates of sphere --------------------------- def mybond(rbond,vi,vf): global x, y, z normalvec = (vf-vi)/np.sqrt( sum((vf-vi)**2) ) # unit vector n' in the direction of vf-vi xyvec = np.array( [-normalvec[1], normalvec[0], 0] ) # e_z x n' xyvec = xyvec/np.sqrt(sum(xyvec**2)) if xyvec[0]==0: xyvec[0]=1e-5 ph = np.arctan( xyvec[1]/xyvec[0] ) th = np.arccos( normalvec[2] ) if xyvec[0] < 0.0 : ph = np.pi + ph a11 = np.cos(ph) a12 = -np.sin(ph)*np.cos(th) a13 = np.sin(ph)*np.sin(th) a21 = np.sin(ph) a22 = np.cos(th)*np.cos(ph) a23 = -np.sin(th)*np.cos(ph) a31 = 0.0 a32 = np.sin(th) a33 = np.cos(th) u, v = np.mgrid[0:2*np.pi:40j, 0:np.sqrt(sum((vf-vi)**2)):10j] cylx = rbond * np.cos(u) cyly = rbond * np.sin(u) cylz = v return a11*cylx+a12*cyly+a13*cylz+vi[0], a21*cylx+a22*cyly+a23*cylz+vi[1], a31*cylx+a32*cyly+a33*cylz+vi[2] ### function to create a colormap ----------------------------------------- def make_cmap(colors, position=None, bit=False): bit_rgb = np.linspace(0,1,256) if position == None: position = np.linspace(0,1,len(colors)) else: if len(position) != len(colors): sys.exit("position length must be the same as colors") elif position[0] != 0 or position[-1] != 1: sys.exit("position must start with 0 and end with 1") if bit: for i in range(len(colors)): colors[i] = (bit_rgb[colors[i][0]], bit_rgb[colors[i][1]], bit_rgb[colors[i][2]]) cdict = {'red':[], 'green':[], 'blue':[]} for pos, color in zip(position, colors): cdict['red'].append((pos, color[0], color[0])) cdict['green'].append((pos, color[1], color[1])) cdict['blue'].append((pos, color[2], color[2])) cmap = LinearSegmentedColormap('my_colormap',cdict,256) return cmap ### function to plot atoms ----------------------------------------------- def atomplot(Rrad, R, colorrgb, ligth_intensity, zorder): global dev colors = [(0.0,0.0,0.0), (int(colorrgb[0]/2),int(colorrgb[1]/2),int(colorrgb[2]/2)), colorrgb, (int((255+colorrgb[0])/2),int((255+colorrgb[1])/2),int((255+colorrgb[2])/2))] my_cmap = make_cmap(colors, bit=True) xplot, yplot, zplot = mysphere(Rrad, R) levels=[] for i in range(0,dev): levels.append( np.sqrt( 1.0 - (i*1.0/(dev-1))**2. ) ) levels=np.sort(levels) levels=Rrad*levels*ligth_intensity plt.contourf( xplot, yplot, zplot, levels, zdir='z', zorder=zorder, offset=0.1, cmap=my_cmap) plt.contour( xplot, yplot, zplot, levels, zdir='z', zorder=zorder+1, linewidth=1, offset=0.1, cmap=my_cmap) ### function to plot bonds ----------------------------------------------- def bondplot(rbond, Ri, Rf, colorrgb, ligth_intensity, zorder): global dev devp = dev/2 colors = [(0.0,0.0,0.0), colorrgb, (255,238,235)] my_cmap = make_cmap(colors, bit=True) xplot, yplot, zplot = mybond(rbond, Ri, Rf) levels=[] for i in range(0,devp): levels.append( np.sqrt( 1.0 - (i*1.0/(devp-1))**2) ) levels=np.sort(levels) levels=rbond*np.sort(levels)*ligth_intensity plt.contourf( xplot, yplot, zplot, levels, zdir='z', zorder=zorder, offset=0.1, cmap=my_cmap) plt.contour( xplot, yplot, zplot, levels, zdir='z', zorder=zorder+1, linewidth=1, offset=0.1, cmap=my_cmap) # =========================================================================================== ############################################################### Start PLOTS fig = plt.figure(figsize=(8.,8.)) gs = gridspec.GridSpec(2, 2) gs.update(wspace=0., hspace=0., left=0.08, right=0.82, bottom=0.17, top=0.9) # set the spacing between axes. ############################################################### def myplot(fname,sublabel,xlab,ylab,phi,origin): global x, y, z data=np.loadtxt(fname,dtype=float) x=data[:,0] y=data[:,1] z=data[:,2] zmin=np.min(z) z=z-zmin z[x<0.43]=max(z) # ***************************************************************** # interpolation # ***************************************************************** #M=100j #xnew, ynew= np.mgrid[min(x):max(x):M, min(y):max(y):M] #from scipy import interpolate as si #f=si.SmoothBivariateSpline(x, y, z, w=None, bbox=[None, None, None, None], kx=3, ky=3, s=0.1, eps=None) #znew=f(np.mgrid[min(x):max(x):M], np.mgrid[min(y):max(y):M]) ## define a treshold #zz=zz-np.min(znew) #znew=znew-np.min(znew) ## for atoms and geometry plot a = np.array([1.0,0.0]) b = np.array([np.cos(np.pi/3),(max(y)-min(y))/(max(x)-min(x))*np.sin(np.pi/3)]) shift = np.array([1.14, 3.]) # np.array([1.68, 3.8]) rs, rs2, rs3 = [], [], [] # rs2: position fcc & rs2: position hcp for i in 1./3., 2./3.: for j in 1./3., 2./3.: rs.append( (i*a + j*b) + shift ) if i!=2./3. and j<2./3.: rs2.append( (i*a + j*b)+1./9.*(a+b) + shift ) rs3.append( (i*a + j*b)+2./9.*(a+b) + shift ) rs=np.array(rs) rs2=np.array(rs2) rs3=np.array(rs3) #N1 = 0.5*(rs[0,:]+rs[3,:])+xp + 0.07*np.array([np.cos(phi*np.pi/180), (max(y)-min(y))/(max(x)-min(x))*np.sin(phi*np.pi/180)]) #N2 = 0.5*(rs[0,:]+rs[3,:])+yp - 0.07*np.array([np.cos(phi*np.pi/180), (max(y)-min(y))/(max(x)-min(x))*np.sin(phi*np.pi/180)]) N1 = origin[0]*a + origin[1]*b + shift + 0.07*np.array([np.cos(phi*np.pi/180), (max(y)-min(y))/(max(x)-min(x))*np.sin(phi*np.pi/180)]) N2 = origin[0]*a + origin[1]*b + shift - 0.07*np.array([np.cos(phi*np.pi/180), (max(y)-min(y))/(max(x)-min(x))*np.sin(phi*np.pi/180)]) xnew=np.zeros((101,101)) ynew=np.zeros((101,101)) znew=np.zeros((101,101)) count=0 for i in range(0,np.shape(xnew)[0]): for j in range(0,np.shape(xnew)[1]): xnew[i,j]=x[count] ynew[i,j]=y[count] znew[i,j]=z[count] count=count+1 levels = np.linspace(-0.5,2.5,25) cax = plt.contourf(xnew, ynew, znew, levels=levels, zdir='z', offset=0, cmap=cmap) CS = plt.contour(xnew, ynew, znew, levels, zdir='z', offset=0, linewidths=1.0, colors='k') manual_locations = [(1.15, 4.0),(0.88, 1.4),(0.66, 4.0), (0.6, 2.5),(0.77,1.8),(1.0, 3.0),(0.7, 2.0),(1.2, 3.0), (1.2, 2.2), (1.4, 1.8), (2.2,1.6)] plt.clabel(CS, inline=1, fmt='%1.1f', fontsize=12, manual=manual_locations) for k in range(0,len(rs)): atomplot(r_Ru, rs[k,:], col[0], ligth_intensity, 3) for k in range(0,len(rs2)): atomplot(r_Ru, rs2[k,:], col[1], ligth_intensity, 3) for k in range(0,len(rs3)): atomplot(r_Ru, rs3[k,:], col[2], ligth_intensity, 2) #vii, vff = np.array([ N1[0], N1[1], 0 ]), np.array([ N2[0], N2[1], 0 ]) #vi, vf = vii + 0.9*2*r_Ru/3*(vff-vii)/np.sqrt(sum((vff-vii)**2)), vff - 0.9*2*r_Ru/3*(vff-vii)/np.sqrt(sum((vff-vii)**2)) #bondplot(0.25*r_Ru/3, vi, vf, colN[0], 1.2, 4) atomplot(2*r_Ru/3, N1, colN[0], ligth_intensity, 4) atomplot(2*r_Ru/3, N2, colN[0], ligth_intensity, 4) # setting minor and major ticks ############################## minorLocator_x = np.mgrid[0:2.5:11j] minorLocator_y = np.mgrid[0:7.0:15j] ax.set_xticks(minorLocator_x, minor = True) ax.set_yticks(minorLocator_y, minor = True) plt.tick_params(axis = 'both', which = 'major', width=1, length=5, labelsize = 14) plt.tick_params(axis = 'both', which = 'minor', width=0.5, length=3.5) # ------------------------------------------------------------ plt.axis([0.4, 2.29, 0.2, 5.5]) if xlab=="No": plt.xticks([],[]) if ylab=="No": plt.yticks([],[]) if xlab=="Yes": plt.xlabel('r ($\AA$)', fontsize=15, color='black') if ylab=="Yes": plt.ylabel('Z ($\AA$)', fontsize=15, color='black') plt.text(0.45,0.7,sublabel,fontsize=15) return cax ####################################################################################### Subplots cmap = cm.jet # subplot 1 ax = fig.add_subplot(gs[0]) cax= myplot('RES','(a)',"No","Yes",0, [1./3.,1./3.]) #ax = point.plot(x='1.51', y='1.57', ax=ax, style='r-', label='point') plt.plot([0.769], [2.202], marker='o', markersize=7, color='white') plt.plot([1.096], [1.549], marker='o', markersize=7, color='white') # subplot 2 ax = fig.add_subplot(gs[1]) cax= myplot('RES2','(b)',"No","No",30, [1./3.+2./9.,1./3.+2./9.]) plt.plot([0.874], [1.586], marker='o', markersize=7, color='white') # subplot 1 ax = fig.add_subplot(gs[2]) cax= myplot('RES3','(c)',"Yes","Yes",0, [1./3.+1./6.,1./3.]) plt.plot([0.837], [1.777], marker='o', markersize=7, color='white') # subplot 2 ax = fig.add_subplot(gs[3]) cax= myplot('RES4','(d)',"Yes","No",120,[1./3.+5./18.,1./3.+5./18.]) plt.plot([0.837], [1.676], marker='o', markersize=7, color='white') fig.subplots_adjust(right=0.8) cbar_ax = fig.add_axes([0.85, 0.17, 0.03, 0.73]) #cbar_ax.set_ylabel('Potential (eV)') cbar = fig.colorbar(cax, cax=cbar_ax, ticks=[-0.5,0.0,0.5,1.0,1.5,2.0,2.5]) cbar.ax.tick_params(labelsize=19) cbar.ax.set_ylabel ('Potential (eV)', fontsize=19) plt.show()