In [14]:

    

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 [6]:

    

ls


    

    




Q3.ipynb                    Si_100_4_10.26.out          Si_20_4_10.26.out           Si_40_4_10.26.out           Si_60_4_10.26.out           Si_80_4_10.26.out           analyze.py
Si.pbe-n-kjpaw_psl.0.1.UPF  Si_100_4_10.26.pw.in        Si_20_4_10.26.pw.in         Si_40_4_10.26.pw.in         Si_60_4_10.26.pw.in         Si_80_4_10.26.pw.in         results.csv
Si.pw.in.template           Si_10_4_10.26.out           Si_30_4_10.26.out           Si_50_4_10.26.out           Si_70_4_10.26.out           Si_90_4_10.26.out           run_pw.py*
Si.pz-n-kjpaw_psl.0.1.UPF   Si_10_4_10.26.pw.in         Si_30_4_10.26.pw.in         Si_50_4_10.26.pw.in         Si_70_4_10.26.pw.in         Si_90_4_10.26.pw.in

In [20]:

    

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


    

In [21]:

    

datafile.head(0)


    

    
Out[21]:






  

    

      

      
filename
      
ecut
      
nkpts
      
alat
      
energy
      
total_force
      
Force (meV/angstrom)
      
Force Converge (meV/angstrom)
    
  
  

    

      
0
      
  Si_10_4_10.26.out
      
  10
      
 24
      
 10.26
      
-93.416654
      
 0.150311
      
 3864.652599
      
  0.000000
    
    

      
1
      
  Si_20_4_10.26.out
      
  20
      
 24
      
 10.26
      
-93.419639
      
 0.153839
      
 3955.361159
      
 90.708560
    
    

      
2
      
  Si_30_4_10.26.out
      
  30
      
 24
      
 10.26
      
-93.420874
      
 0.153770
      
 3953.587097
      
 -1.774062
    
    

      
3
      
  Si_40_4_10.26.out
      
  40
      
 24
      
 10.26
      
-93.421050
      
 0.153827
      
 3955.052626
      
  1.465529
    
    

      
4
      
  Si_50_4_10.26.out
      
  50
      
 24
      
 10.26
      
-93.421039
      
 0.153907
      
 3957.109510
      
  2.056883
    
    

      
5
      
  Si_60_4_10.26.out
      
  60
      
 24
      
 10.26
      
-93.421049
      
 0.153905
      
 3957.058088
      
 -0.051422
    
    

      
6
      
  Si_70_4_10.26.out
      
  70
      
 24
      
 10.26
      
-93.421052
      
 0.153907
      
 3957.109510
      
  0.051422
    
    

      
7
      
  Si_80_4_10.26.out
      
  80
      
 24
      
 10.26
      
-93.421057
      
 0.153895
      
 3956.800977
      
 -0.308533
    
    

      
8
      
  Si_90_4_10.26.out
      
  90
      
 24
      
 10.26
      
-93.421061
      
 0.153911
      
 3957.212354
      
  0.411377
    
    

      
9
      
 Si_100_4_10.26.out
      
 100
      
 24
      
 10.26
      
-93.421064
      
 0.153883
      
 3956.492445
      
 -0.719909
    
  

In [22]:

    

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


    

    
Out[22]:





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

In [23]:

    

y_energy = datafile['Force (meV/angstrom)'].tolist()
y_conver = datafile['Force Converge (meV/angstrom)'].tolist()


    

In [33]:

    

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(r'Force $meV/\AA$',fontsize=20)
ax0.set_ylim(3860,3970)
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('Force Converge $meV/\AA$',fontsize=18)
# ax1.set_yscale('log')
ax1.set_ylim(-10,100)
ax1.annotate(r'Convergence < 5 $meV/\AA$ Cutoff 30 Ry', xy=(32, 5),  xycoords='data',
                xytext=(0.95, 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()