In [1]:

    

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


    

In [2]:

    

import numpy as np
import pandas as pd


    

In [3]:

    

ls


    

    




Converge_Image.png                                   Si_100_8_10.26.pw.in                                 Si_40_8_10.26.pw.in                                  Si_80_8_10.26.pw.in
NANO266 Lab 2 - Convergence of DFT calculations.pdf  Si_10_8_10.26.out                                    Si_50_8_10.26.out                                    Si_90_8_10.26.out
Q1.ipynb                                             Si_10_8_10.26.pw.in                                  Si_50_8_10.26.pw.in                                  Si_90_8_10.26.pw.in
README.md                                            Si_20_8_10.26.out                                    Si_60_8_10.26.out                                    analyze.py*
Si.pbe-n-kjpaw_psl.0.1.UPF                           Si_20_8_10.26.pw.in                                  Si_60_8_10.26.pw.in                                  results.csv
Si.pw.in.template                                    Si_30_8_10.26.out                                    Si_70_8_10.26.out                                    run_pw.py*
Si.pz-n-kjpaw_psl.0.1.UPF                            Si_30_8_10.26.pw.in                                  Si_70_8_10.26.pw.in
Si_100_8_10.26.out                                   Si_40_8_10.26.out                                    Si_80_8_10.26.out

In [4]:

    

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


    

In [5]:

    

datafile.head(0)


    

    
Out[5]:






  

    

      

      
filename
      
ecut (Ry)
      
nkpts
      
alat
      
energy
      
total_force
      
Total Energy (eV)
      
energy difference (meV/atom)
    
  
  

    

      
0
      
  Si_10_8_10.26.out
      
  10
      
 29
      
 10.26
      
-93.447441
      
 0
      
-1271.417145
      
  0.000000
    
    

      
1
      
  Si_20_8_10.26.out
      
  20
      
 29
      
 10.26
      
-93.452291
      
 0
      
-1271.483131
      
 32.993192
    
    

      
2
      
  Si_30_8_10.26.out
      
  30
      
 29
      
 10.26
      
-93.453525
      
 0
      
-1271.499924
      
  8.396549
    
    

      
3
      
  Si_40_8_10.26.out
      
  40
      
 29
      
 10.26
      
-93.453698
      
 0
      
-1271.502283
      
  1.179205
    
    

      
4
      
  Si_50_8_10.26.out
      
  50
      
 29
      
 10.26
      
-93.453688
      
 0
      
-1271.502151
      
 -0.065920
    
    

      
5
      
  Si_60_8_10.26.out
      
  60
      
 29
      
 10.26
      
-93.453698
      
 0
      
-1271.502282
      
  0.065647
    
    

      
6
      
  Si_70_8_10.26.out
      
  70
      
 29
      
 10.26
      
-93.453700
      
 0
      
-1271.502309
      
  0.013606
    
    

      
7
      
  Si_80_8_10.26.out
      
  80
      
 29
      
 10.26
      
-93.453705
      
 0
      
-1271.502373
      
  0.032041
    
    

      
8
      
  Si_90_8_10.26.out
      
  90
      
 29
      
 10.26
      
-93.453708
      
 0
      
-1271.502415
      
  0.020817
    
    

      
9
      
 Si_100_8_10.26.out
      
 100
      
 29
      
 10.26
      
-93.453710
      
 0
      
-1271.502449
      
  0.017075
    
  

In [6]:

    

x = datafile['ecut (Ry)'].tolist()
x


    

    
Out[6]:





[10, 20, 30, 40, 50, 60, 70, 80, 90, 100]

In [7]:

    

y_energy = datafile['Total Energy (eV)'].tolist()
y_conver = datafile['energy difference (meV/atom)'].tolist()


    

In [8]:

    

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('Kinetic Energy Cutoff (Ry)',fontsize=24)
ax0.set_ylabel('Total Energy (eV)',fontsize=20)
ax0.set_ylim( -1271.54,-1271.38)
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('Kinetic Energy Cutoff (Ry)',fontsize=24)
ax1.set_ylabel('Energy Converge (meV)',fontsize=18)
# ax1.set_yscale('log')
ax1.set_ylim(-10,80)
ax1.annotate('Convergence < 5 meV', xy=(41, 5),  xycoords='data',
                xytext=(0.8, 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()