-
Notifications
You must be signed in to change notification settings - Fork 5
/
draw_plane3b_unit_average.py
106 lines (83 loc) · 3.08 KB
/
draw_plane3b_unit_average.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
import numpy as np
import sys
import matplotlib.pyplot as plt
import argparse
SMALL_SIZE = 12
MEDIUM_SIZE = 14
BIGGER_SIZE = 16
plt.rc('font', size=BIGGER_SIZE) # controls default text sizes
plt.rc('axes', titlesize=BIGGER_SIZE) # fontsize of the axes title
plt.rc('axes', labelsize=MEDIUM_SIZE) # fontsize of the x and y labels
plt.rc('xtick', labelsize=SMALL_SIZE) # fontsize of the tick labels
plt.rc('ytick', labelsize=SMALL_SIZE) # fontsize of the tick labels
plt.rc('legend', fontsize=SMALL_SIZE) # legend fontsize
plt.rc('figure', titlesize=BIGGER_SIZE) # fontsize of the figure title
#plt.figure(figsize=(4,6))
parser = argparse.ArgumentParser(description='Plot planer average potential')
parser.add_argument('--filename', metavar='f',type=str,dest='filename',default='PLANAR_AVERAGE.dat')
parser.add_argument('--percent1', metavar='p1',type=float,dest='percent1',default=0.)
parser.add_argument('--percent2', metavar='p2',type=float,dest='percent2',default=1.)
args = parser.parse_args()
#filename= 'PLANAR_AVERAGE.dat'
filename = args.filename
percent1 = args.percent1
percent2 = args.percent2
## load data
data=np.loadtxt(filename) # this one works
x=data[:,0]
y=data[:,1]
#####################################################################
labeltype=1 # for 'PLANAR_AVERAGE.1e.a.dat'
if labeltype == 1:
# set up labels for plotting plane potential 'PLANAR_AVERAGE.1e.a.dat'
legend=''
#plt.legend()
title='Planar average potential'
xlabel='Distance'
xlabelunit=' ($\AA$)'
ylabel='Potential'
ylabelunit=' (eV $\AA$/e)'
color='tab:blue'
unittype=1 # use the default angstrom unit
#unittype=2 # use the bohr unit
if unittype==2:
x=x*1.88973 # 1A = 1.88973 Bohr
xlabelunit=' (bohr)'
#####################################################################
### vacuum level
arg=np.argmax(y)
### electro-static average in bulk
length = len(x)
len1 = int(percent1*length)
len2 = int(percent2*length)
## here len2 can == length
selected = y[len1:len2]
### use average
integrated = np.sum(selected)
print('Integration by summing all=%s' % integrated )
selected_ymin = selected.min()
selected_ymax = selected.max()
average = (selected_ymin + selected_ymax) / 2
ax = plt.gca()
## here len2 must be < length
if len2 >=length :
## to avoid the x[len2] outside range
len2 = length - 1
ax.fill_betweenx(np.linspace(selected_ymin,selected_ymax,50), x[len1], x[len2], color='green', alpha=0.3 )
#plt.vlines(x[len1],selected_ymin,selected_ymax,color='purple')
#plt.vlines(x[len2],selected_ymin,selected_ymax,color='purple')
title='average=%.5f' % (average)
print('The average electro-static potential is %s eV' % (average))
plt.hlines(average, x[0],x[-1],color='purple')
# black, tab:blue, tab:green, tab:orange
#plt.plot(x, y, color=color, label=legend)
plt.scatter(x, y, color=color, label=legend)
plt.title(title) # add title
plt.xlabel(xlabel+xlabelunit) # x label
plt.ylabel(ylabel+ylabelunit) # y label
plt.tight_layout()
figname= filename[:-3] + 'unit.average.pdf'
print('\tfigure saved as \n%s' % (figname))
plt.savefig(figname,dpi=600,transparent=True)
plt.show()
plt.close()