Fitting -- Birch's formula

#  Young-Woo Son @ KIAS

#

#   input gnuplot-file for fitting vol(a.u.). vs. energy(a.u.) into

#   Birch's equation upto second order. Input & output B0 is given in Gpa

#

#  Inclusion of fourth order correction should improve B0,B0P.

#  I've found from the example of iron that the fourth order correction

#  gives almost equal result as Murnaghan's equation, but third order only

#  is slightly wrong. Note that the form of second and third order terms are

# invariant under inclusion of higher order corrections.

E(V) = E0 + 9.0/8.0*B0*V0*((V0/V)**(2.0/3.0)-1)**2 + 9.0/16.0*B0*(B0P-4)* \

       V0*((V0/V)**(2.0/3.0)-1.0)**3.0 + R4*((V0/V)**(2.0/3.0)-1.0)**4.0

B0P   = 2.60567 ; B0 = 0.111821

V0     = 167.834 ; E0 = 0.1            ; R4  = -73.033

set term postscript eps enhanced 10

set output 'Ca_bondinglength.eps'

set encoding iso_8859_1

set size 0.35,0.4

set xr[150:200]

set yr[0.0:0.24]

set xtics 150, 10, 200

set ytics 0.0,0.04,0.4

set format y "%.2f"

#set mytics 2

set xtics nomirror

set nokey

set title "Ca bulk (Ultra Soft PP, PW91)\n Birch eq. fit: a=5.516{\305}, B_0=17.9 GPa" font 'Helvetica'

set xlabel 'Volume({\305}^3)

set ylabel 'E_{total} (eV)'

fit E(x) 'Etot.dat'  u 1:($2*13.6058) via B0P,B0,V0,E0,R4

plot E(x),'Etot.dat' u 1:($2*13.6058) w p pt 7 ps 1

# converting into GPa (energy = eV, length = Ang)

# B0(GPa) = B0(eV/Ang^3) *1.6022*100.0

## In this example, cell volume in Ang^3 and energy is in Ryd (thus *13.6058 is necessary.)

## It is higly recommended to use V0 and E0 at the minimum position.

<fig.log> file

...

Final set of parameters            Asymptotic Standard Error

=======================            ==========================

B0P             = 2.60567          +/- 0.1231       (4.725%)

B0              = 0.111821         +/- 0.0005028    (0.4497%)

V0              = 167.834          +/- 0.02815      (0.01677%)

E0              = 0.0213131        +/- 0.0001264    (0.5932%)

R4              = -73.033          +/- 17.73        (24.28%)

correlation matrix of the fit parameters:

               B0P    B0     V0     E0     R4    

B0P             1.000

B0             -0.489  1.000

V0             -0.922  0.476  1.000

E0              0.108 -0.704 -0.137  1.000

R4              0.829 -0.863 -0.735  0.422  1.000