In [34]:

    

import matplotlib.pyplot as plt
%matplotlib inline
from pymatgen.util.plotting_utils import get_publication_quality_plot


    

In [35]:

    

import numpy as np
import pandas as pd


    

In [52]:

    

ls


    

    




10.26_Bohr/                      10.30_Bohr/                      Q5.ipynb                         results_1026_1030_combine.csv    results_1026_1030_rank_nkpt.csv

In [53]:

    

datafile2 = pd.read_csv('results_1026_1030_rank_nkpt.csv')


    

In [55]:

    

datafile2.head(10)


    

    
Out[55]:






  

    

      

      
filename(1026)
      
ecut(1026)
      
nkpts(1026)
      
alat(1026)
      
energy(1026)
      
filename(1030)
      
ecut(1030)
      
nkpts(1030)
      
alat(1030)
      
energy(1030)
      
energy_difference (meV/atom)
    
  
  

    

      
0
      
  Si_10_5_10.26.out
      
  10
      
 10
      
 10.26
      
-93.443092
      
  Si_10_5_10.3.out
      
  10
      
 10
      
 10.3
      
-93.437631
      
-37.148030
    
    

      
1
      
  Si_20_5_10.26.out
      
  20
      
 10
      
 10.26
      
-93.447850
      
  Si_20_5_10.3.out
      
  20
      
 10
      
 10.3
      
-93.448221
      
  2.525353
    
    

      
2
      
  Si_30_5_10.26.out
      
  30
      
 10
      
 10.26
      
-93.449064
      
  Si_30_5_10.3.out
      
  30
      
 10
      
 10.3
      
-93.449400
      
  2.283239
    
    

      
3
      
  Si_40_5_10.26.out
      
  40
      
 10
      
 10.26
      
-93.449238
      
  Si_40_5_10.3.out
      
  40
      
 10
      
 10.3
      
-93.449574
      
  2.283579
    
    

      
4
      
  Si_50_5_10.26.out
      
  50
      
 10
      
 10.26
      
-93.449229
      
  Si_50_5_10.3.out
      
  50
      
 10
      
 10.3
      
-93.449563
      
  2.275348
    
    

      
5
      
  Si_60_5_10.26.out
      
  60
      
 10
      
 10.26
      
-93.449239
      
  Si_60_5_10.3.out
      
  60
      
 10
      
 10.3
      
-93.449572
      
  2.270926
    
    

      
6
      
  Si_70_5_10.26.out
      
  70
      
 10
      
 10.26
      
-93.449241
      
  Si_70_5_10.3.out
      
  70
      
 10
      
 10.3
      
-93.449574
      
  2.266980
    
    

      
7
      
  Si_80_5_10.26.out
      
  80
      
 10
      
 10.26
      
-93.449245
      
  Si_80_5_10.3.out
      
  80
      
 10
      
 10.3
      
-93.449578
      
  2.268681
    
    

      
8
      
  Si_90_5_10.26.out
      
  90
      
 10
      
 10.26
      
-93.449249
      
  Si_90_5_10.3.out
      
  90
      
 10
      
 10.3
      
-93.449582
      
  2.268341
    
    

      
9
      
 Si_100_5_10.26.out
      
 100
      
 10
      
 10.26
      
-93.449251
      
 Si_100_5_10.3.out
      
 100
      
 10
      
 10.3
      
-93.449584
      
  2.268409
    
  

In [36]:

    

datafile = pd.read_csv('results_1026_1030_combine.csv')


    

In [21]:

    

datafile.head(10)


    

    
Out[21]:






  

    

      

      
filename(1026)
      
ecut(1026)
      
nkpts(1026)
      
alat(1026)
      
energy(1026)
      
filename(1030)
      
ecut(1030)
      
nkpts(1030)
      
alat(1030)
      
energy(1030)
      
energy_difference (meV/atom)
    
  
  

    

      
0
      
  Si_10_5_10.26.out
      
 10
      
 10
      
 10.26
      
-93.443092
      
  Si_10_5_10.3.out
      
 10
      
 10
      
 10.3
      
-93.437631
      
 37.148030
    
    

      
1
      
  Si_10_6_10.26.out
      
 10
      
 16
      
 10.26
      
-93.445716
      
  Si_10_6_10.3.out
      
 10
      
 16
      
 10.3
      
-93.440330
      
 36.638769
    
    

      
2
      
  Si_10_7_10.26.out
      
 10
      
 20
      
 10.26
      
-93.447019
      
  Si_10_7_10.3.out
      
 10
      
 20
      
 10.3
      
-93.441816
      
 35.389563
    
    

      
3
      
  Si_10_8_10.26.out
      
 10
      
 29
      
 10.26
      
-93.447441
      
  Si_10_8_10.3.out
      
 10
      
 29
      
 10.3
      
-93.442182
      
 35.770454
    
    

      
4
      
  Si_10_9_10.26.out
      
 10
      
 35
      
 10.26
      
-93.447557
      
  Si_10_9_10.3.out
      
 10
      
 35
      
 10.3
      
-93.442162
      
 36.703396
    
    

      
5
      
 Si_10_10_10.26.out
      
 10
      
 47
      
 10.26
      
-93.447692
      
 Si_10_10_10.3.out
      
 10
      
 47
      
 10.3
      
-93.442302
      
 36.669042
    
    

      
6
      
  Si_20_5_10.26.out
      
 20
      
 10
      
 10.26
      
-93.447850
      
  Si_20_5_10.3.out
      
 20
      
 10
      
 10.3
      
-93.448221
      
 -2.525353
    
    

      
7
      
  Si_20_6_10.26.out
      
 20
      
 16
      
 10.26
      
-93.450785
      
  Si_20_6_10.3.out
      
 20
      
 16
      
 10.3
      
-93.451192
      
 -2.763860
    
    

      
8
      
  Si_20_7_10.26.out
      
 20
      
 20
      
 10.26
      
-93.451882
      
  Si_20_7_10.3.out
      
 20
      
 20
      
 10.3
      
-93.452218
      
 -2.288001
    
    

      
9
      
  Si_20_8_10.26.out
      
 20
      
 29
      
 10.26
      
-93.452291
      
  Si_20_8_10.3.out
      
 20
      
 29
      
 10.3
      
-93.452645
      
 -2.413310
    
  

In [22]:

    

x = list(set(datafile['nkpts(1026)'].tolist()))
x


    

    
Out[22]:





[35, 10, 47, 16, 20, 29]

In [37]:

    

datafile10 =datafile[(datafile['ecut(1026)']==10)]
datafile20 =datafile[(datafile['ecut(1026)']==20)]
datafile30 =datafile[(datafile['ecut(1026)']==30)]
datafile40 =datafile[(datafile['ecut(1026)']==40)]
datafile50 =datafile[(datafile['ecut(1026)']==50)]
datafile60 =datafile[(datafile['ecut(1026)']==60)]
datafile70 =datafile[(datafile['ecut(1026)']==70)]
datafile80 =datafile[(datafile['ecut(1026)']==80)]
datafile90 =datafile[(datafile['ecut(1026)']==90)]
datafile100 =datafile[(datafile['ecut(1026)']==100)]


    

In [56]:

    

datakpt10 = datafile2[(datafile2['nkpts(1026)']==10)]
datakpt16 = datafile2[(datafile2['nkpts(1026)']==16)]
datakpt20 = datafile2[(datafile2['nkpts(1026)']==20)]
datakpt29 = datafile2[(datafile2['nkpts(1026)']==29)]
datakpt35 = datafile2[(datafile2['nkpts(1026)']==35)]
datakpt47 = datafile2[(datafile2['nkpts(1026)']==47)]


    

In [59]:

    

x_cut_10 =datakpt10['ecut(1026)'].tolist()
y_cut_conver_10 = datakpt10['energy_difference (meV/atom)'].tolist()
x_cut_16 =datakpt16['ecut(1026)'].tolist()
y_cut_conver_16 = datakpt16['energy_difference (meV/atom)'].tolist()
x_cut_20 =datakpt20['ecut(1026)'].tolist()
y_cut_conver_20 = datakpt20['energy_difference (meV/atom)'].tolist()
x_cut_29 =datakpt29['ecut(1026)'].tolist()
y_cut_conver_29 = datakpt29['energy_difference (meV/atom)'].tolist()
x_cut_35 =datakpt35['ecut(1026)'].tolist()
y_cut_conver_35 = datakpt35['energy_difference (meV/atom)'].tolist()
x_cut_47 =datakpt47['ecut(1026)'].tolist()
y_cut_conver_47 = datakpt47['energy_difference (meV/atom)'].tolist()


    

In [24]:

    

datafile10


    

    
Out[24]:






  

    

      

      
filename(1026)
      
ecut(1026)
      
nkpts(1026)
      
alat(1026)
      
energy(1026)
      
filename(1030)
      
ecut(1030)
      
nkpts(1030)
      
alat(1030)
      
energy(1030)
      
energy_difference (meV/atom)
    
  
  

    

      
0
      
  Si_10_5_10.26.out
      
 10
      
 10
      
 10.26
      
-93.443092
      
  Si_10_5_10.3.out
      
 10
      
 10
      
 10.3
      
-93.437631
      
 37.148030
    
    

      
1
      
  Si_10_6_10.26.out
      
 10
      
 16
      
 10.26
      
-93.445716
      
  Si_10_6_10.3.out
      
 10
      
 16
      
 10.3
      
-93.440330
      
 36.638769
    
    

      
2
      
  Si_10_7_10.26.out
      
 10
      
 20
      
 10.26
      
-93.447019
      
  Si_10_7_10.3.out
      
 10
      
 20
      
 10.3
      
-93.441816
      
 35.389563
    
    

      
3
      
  Si_10_8_10.26.out
      
 10
      
 29
      
 10.26
      
-93.447441
      
  Si_10_8_10.3.out
      
 10
      
 29
      
 10.3
      
-93.442182
      
 35.770454
    
    

      
4
      
  Si_10_9_10.26.out
      
 10
      
 35
      
 10.26
      
-93.447557
      
  Si_10_9_10.3.out
      
 10
      
 35
      
 10.3
      
-93.442162
      
 36.703396
    
    

      
5
      
 Si_10_10_10.26.out
      
 10
      
 47
      
 10.26
      
-93.447692
      
 Si_10_10_10.3.out
      
 10
      
 47
      
 10.3
      
-93.442302
      
 36.669042
    
  

In [38]:

    

x_nkpts_10 = datafile10['nkpts(1026)'].tolist()
y_conver_10 = datafile10['energy_difference (meV/atom)'].tolist()
x_nkpts_20 = datafile20['nkpts(1026)'].tolist()
y_conver_20 = datafile20['energy_difference (meV/atom)'].tolist()
x_nkpts_30 = datafile30['nkpts(1026)'].tolist()
y_conver_30 = datafile30['energy_difference (meV/atom)'].tolist()
x_nkpts_40 = datafile40['nkpts(1026)'].tolist()
y_conver_40 = datafile40['energy_difference (meV/atom)'].tolist()
x_nkpts_50 = datafile50['nkpts(1026)'].tolist()
y_conver_50 = datafile50['energy_difference (meV/atom)'].tolist()
x_nkpts_60 = datafile60['nkpts(1026)'].tolist()
y_conver_60 = datafile60['energy_difference (meV/atom)'].tolist()
x_nkpts_70 = datafile70['nkpts(1026)'].tolist()
y_conver_70 = datafile70['energy_difference (meV/atom)'].tolist()
x_nkpts_80 = datafile80['nkpts(1026)'].tolist()
y_conver_80 = datafile80['energy_difference (meV/atom)'].tolist()
x_nkpts_90 = datafile90['nkpts(1026)'].tolist()
y_conver_90 = datafile90['energy_difference (meV/atom)'].tolist()
x_nkpts_100 = datafile100['nkpts(1026)'].tolist()
y_conver_100 = datafile100['energy_difference (meV/atom)'].tolist()


    

In [68]:

    

fig, (ax0,ax1)=plt.subplots(nrows=2)
fig.set_size_inches(10,10)

ax0.set_title('Energy Convergence VS K-point',fontsize=24, y=1.05)
ax0.set_ylabel('Energy Difference (meV/atom)',fontsize=18)
ax0.tick_params(labelsize = 18)
ax0.set_xlabel('# of K-points',fontsize=24)
# ax0.set_yscale('log')


#plot ecut 10 data
# ax0.scatter(x_nkpts_10,y_conver_10,marker='o',alpha=0.75,s=60,c='sandybrown',label='Encut: 10 (Ry)')
# ax0.plot(x_nkpts_10,y_conver_10,c='black')

#plot ecut 20 data
ax0.scatter(x_nkpts_20,y_conver_20,marker='o',alpha=0.75,s=60,c='lime',label='Encut: 20 (Ry)')
# ax0.plot(x_nkpts_20,y_conver_20,c='black')

#plot ecut 30 data
ax0.scatter(x_nkpts_30,y_conver_30,marker='o',alpha=0.75,s=40,c='dodgerblue',label='Encut: 30 (Ry)')
# ax0.plot(x_nkpts_30,y_conver_30,c='black')

#plot ecut 40 data
ax0.scatter(x_nkpts_40,y_conver_40,marker='o',alpha=0.75,s=40,c='tomato',label='Encut: 40 (Ry)')
# ax0.plot(x_nkpts_40,y_conver_40,c='black')

#plot ecut 50 data
ax0.scatter(x_nkpts_50,y_conver_50,marker='o',alpha=0.75,s=40,c='sage',label='Encut: 50 (Ry)')
# ax0.plot(x_nkpts_50,y_conver_50,c='black')

#plot ecut 60 data
ax0.scatter(x_nkpts_60,y_conver_60,marker='o',alpha=0.75,s=40,c='blue',label='Encut: 60 (Ry)')
# ax0.plot(x_nkpts_60,y_conver_60,c='black')

#plot ecut 70 data
ax0.scatter(x_nkpts_70,y_conver_70,marker='o',alpha=0.75,s=40,c='darkviolet',label='Encut: 70 (Ry)')
# ax0.plot(x_nkpts_70,y_conver_70,c='black')

#plot ecut 80 data
ax0.scatter(x_nkpts_80,y_conver_80,marker='o',alpha=0.75,s=40,c='coral',label='Encut: 80 (Ry)')
# ax0.plot(x_nkpts_80,y_conver_80,c='black')

ax0.legend(loc='best')


ax1.set_title('Energy Convergence VS Ecut',fontsize=24, y=1.05)
ax1.set_ylabel('Energy Difference (meV/atom)',fontsize=18)
ax1.tick_params(labelsize = 18)
ax1.set_xlabel('Energy cutoff (Ry)',fontsize=24)

#plot kpts 10 data
ax1.scatter(x_cut_10,y_cut_conver_10,marker='o',alpha=0.75,s=60,c='sandybrown',label='Number of K-point: 10')
# ax1.plot(x_cut_10,y_cut_conver_10,c='black')

#plot kpts 16 data
ax1.scatter(x_cut_16,y_cut_conver_16,marker='o',alpha=0.75,s=60,c='lime',label='Number of K-point: 16')
# ax1.plot(x_cut_16,y_cut_conver_16,c='black')

#plot kpts 20 data
ax1.scatter(x_cut_20,y_cut_conver_20,marker='o',alpha=0.75,s=60,c='blue',label='Number of K-point: 20')
# ax1.plot(x_cut_20,y_cut_conver_20,c='black')

#plot kpts 29 data
ax1.scatter(x_cut_29,y_cut_conver_29,marker='o',alpha=0.75,s=60,c='maroon',label='Number of K-point: 29')
# ax1.plot(x_cut_29,y_cut_conver_29,c='black')

#plot kpts 35 data
ax1.scatter(x_cut_35,y_cut_conver_35,marker='o',alpha=0.75,s=60,c='darkviolet',label='Number of K-point: 35')
# ax1.plot(x_cut_35,y_cut_conver_35,c='black')

#plot kpts 47 data
ax1.scatter(x_cut_47,y_cut_conver_47,marker='o',alpha=0.75,s=60,c='red',label='Number of K-point: 47')
# ax1.plot(x_cut_10,y_cut_conver_47,c='black')

ax1.legend(loc='best')
ax1.set_ylim(1.8,3)

plt.subplots_adjust(hspace=0.5)
plt.show()


    

In [13]:

    

fig, (ax0,ax1)=plt.subplots(nrows=2)
fig.set_size_inches(10,10)

ax0.scatter(x, y_energy,marker='o',alpha=0.75,s=60,c='r')
ax0.plot(x,y_energy,c='black')
ax0.set_title('Energy',fontsize=30,y=1.05)
ax0.set_xlabel('# of K-points',fontsize=24)
ax0.set_ylabel('Total Energy (eV)',fontsize=20)
ax0.set_ylim( -1271.51,-1271.43)
ax0.tick_params(labelsize = 18)

ax1.scatter(x,y_conver,marker='o',alpha=0.75, s=60,c='b')
ax1.plot(x,y_conver,c='black')
ax1.set_title('Energy Difference',fontsize=24,y=1.05)
ax1.set_xlabel('# of K-points',fontsize=24)
ax1.set_ylabel('Energy Converge (meV/atom)',fontsize=18)
# ax1.set_yscale('log')
ax1.set_ylim(-5,30)
ax1.annotate('Convergence < 5 meV, # of K-point: 29 ', xy=(30, 5),  xycoords='data',
                xytext=(0.9, 0.95), textcoords='axes fraction',
                arrowprops=dict(facecolor='black', shrink=0.05),
                horizontalalignment='right', verticalalignment='top',fontsize=16
                )
ax1.tick_params(labelsize = 18)

plt.subplots_adjust(hspace=0.5)
plt.show()


    

In [17]:

    




    

    
Out[17]:





<matplotlib.text.Text at 0x10891e9d0>

In [18]:

    

plt.show()