TMM4100 Materialteknikk
Parameterstudie og kurvetilpasning for atombinding¶
Den generiske modellen for bindingsenergi fra Atombinding er bare en empirisk modell, og skal den ha noe nytte, må vi kunne finne/tilpasse de ulike parametrene slik at modellen passer med et virkelig system.
Oppsummering av relasjoner:
Energifunksjonen:
\begin{equation} E = -A r^{-m} + B r^{-n} \tag{1} \end{equation}Kraft er energi derivert med hensyn på $r$
\begin{equation} F = mA r^{-m-1} -nB r^{-n-1} \tag{2} \end{equation}Ved likevekt er $E = E_0$ når $r = r_0$ slik at (1) blir
\begin{equation} E_0 = -A r_0^{-m} + B r_0^{-n} \tag{3} \end{equation}Ved likevekt er $F=0$ og ligning (2) blir dermed
\begin{equation} 0 = mA r_0^{-m-1} -nB r_0^{-n-1} \tag{4} \end{equation}Dermed kan vi løse (3) og (4) med hensyn på $A$ og $B$:
Fra (4):
\begin{equation} B = \frac{m}{n} A r_0^{n-m} \tag{5} \end{equation}Substituerer $B$ fra (5) i (3):
\begin{equation} A = -\frac{E_0 r_0^m}{1-m/n} \tag{6} \end{equation}Bindingens stivhet $k$ ved likevektsavstanden $r_0$ er
\begin{equation} k = \frac{dF}{dr} \bigg |_{r=r_0} =-(m+1)mA r_0^{-m-2} + (n+1)nB r_0^{-n-2} \tag{7} \end{equation}Separasjonsavstanden $r_s$ er der kraften har et maksimum:
\begin{equation} \frac{dF}{dr}\bigg |_{r=r_s} = 0 \Rightarrow -(m+1)mA r_s^{-m-2} + (n+1)nB r_s^{-n-2} = 0 \Rightarrow r_s = \bigg( \frac{m(m+1)A}{n(n+1)B} \bigg)^\frac{1}{m-n} \tag{8} \end{equation}Maks kraft $F_s$ blir nå gitt av (2) innsatt for $r = r_s$ fra (8)
Numerisk eksempel:
import numpy as np
from matplotlib import pyplot as plt
E0 = -3
r0 = 0.25
m = 8
n = 4
A = -(E0*r0**m)/(1-m/n) # ligning (6)
B = (m/n)*A*r0**(n-m) # ligning (5)
r=np.linspace(0.2, 0.4, 1000)
E = -A*r**(-m) + B*r**(-n) # ligning (1)
F = m*A*r**(-m-1) -n*B*r**(-n-1) # ligning (2)
k = -(m+1)*m*A*r0**(-m-2) + (n+1)*n*B*r0**(-n-2) # ligning (7)
rs = ( (m*(m+1)*A)/(n*(n+1)*B) )**(1/(m-n)) # ligning (8)
Fs = m*A*rs**(-m-1) -n*B*rs**(-n-1) # ligning (2) for r=rs
fig,(ax1,ax2)=plt.subplots(nrows=1, ncols=2,figsize=(9,5))
ax1.plot(r,E ,color='blue',label='$E$')
ax1.set_xlabel('$r$ [nm]')
ax1.set_ylabel('$E$ [eV]')
ax1.legend(loc='best')
ax1.grid(True)
ax1.set_xlim(0.2,max(r))
ax1.set_ylim(E0*1.1,0)
ax1.set_title('$E_0$={}, $r_0$={}, $m$={:.1f}, $n$={:.1f}, $A$={:.2E}, $B$={:.2E}'.format(E0,r0,m,n,A,B), fontsize=10)
ax2.plot(r,F ,color='red',label='$F$')
ax2.set_xlabel('$r$ [nm]')
ax2.set_ylabel('$F$ [eV / nm]')
ax2.legend(loc='best')
ax2.grid(True)
ax2.set_xlim(0.2,max(r))
ax2.set_ylim(-Fs/3,Fs*1.1)
ax2.set_title('$F_s$={:.2f}, $r_s$={:.3f}, k={:.1f}'.format(Fs,rs,k), fontsize=10)
plt.tight_layout()
plt.draw()
Effekten av $m$ og $n$
Separasjonskraften $F_s$ og bindingens stivhet $k$ øker med økende verdier av $m$ og $n$ som vist i følgende eksempel:
E0 = -3 # holdes konstant
r0 = 0.25 # holdes konstant
fig,(ax1,ax2)=plt.subplots(nrows=1, ncols=2,figsize=(9,5))
r=np.linspace(0.2, 0.5, 1000)
for m,n in zip((2,2,2,3,3,3),(4,5,6,4,5,6)):
A = -(E0*r0**m)/(1-m/n) # ligning (6)
B = (m/n)*A*r0**(n-m) # ligning (5)
E = -A*r**(-m) + B*r**(-n) # ligning (1)
F = m*A*r**(-m-1) -n*B*r**(-n-1) # ligning (2)
k = -(m+1)*m*A*r0**(-m-2) + (n+1)*n*B*r0**(-n-2) # ligning (7)
ax1.plot(r,E ,label='m={}, n={}'.format(m,n))
ax1.set_xlabel('$r$ [nm]')
ax1.set_ylabel('$E$ [eV]')
ax1.legend(loc='best')
ax1.set_ylim(E0*1.1,0)
ax1.grid(True)
ax2.plot(r,F ,label='k={:.0f}'.format(k))
ax2.set_xlabel('$r$ [nm]')
ax2.set_ylabel('$F$ [eV]')
ax2.legend(loc='best')
ax2.set_ylim(-5,15)
ax2.grid(True)
plt.tight_layout()
plt.draw()
Parameter fra empiriske data
Anta at vi har, på en eller annen måte, funnet følgende data gjennom empiriske målinger:
$$E_0 = -4, \quad r_0 = 0.3, \quad k = 900, \quad F_s = 900$$Hva blir verdiene av $m$, $n$, $A$ og $B$?
Nå har vi tilstrekkelig med ligninger og grensebetingelser, men vi har et system av ligninger som ikke akkurat er lineære. Systemet kan løses numerisk med f.eks. scipy.optimize.fsolve, men denne metoden krever uansett god første gjetting på startverdier for iterasjonene, så vi kan likegodt gjøre en mindre avansert prosess:
E0 = -4 # holdes konstant
r0 = 0.3 # holdes konstant
k_empirisk = 900
Fs_empirisk = 16
for i in range(0,1000000): # kjører en million tester...
m = 1+10*np.random.random() # antar at n er mellom 1 og 10
n = 1+10*np.random.random() # antar at m er mellom 1 og 10
A = -(E0*r0**m)/(1-m/n) # ligning (6)
B = (m/n)*A*r0**(n-m) # ligning (5)
k = -(m+1)*m*A*r0**(-m-2) + (n+1)*n*B*r0**(-n-2) # ligning (7)
rs = ( (m*(m+1)*A)/(n*(n+1)*B) )**(1/(m-n)) # ligning (8)
Fs = m*A*rs**(-m-1) -n*B*rs**(-n-1) # ligning (2) for r=rs
if abs(k-k_empirisk) < 1 and abs(Fs-Fs_empirisk) < 0.1: # sjekker om resultater er tilstrekkelig nærme
print('m:{:.2f}, n:{:.2f}, k:{:.1f}, Fs:{:.2f}'.format(m,n,k,Fs))
Resultatet viser at det er to sett av parameter som tilfredstiller empiriske data:
$$m \approx 2.3, \quad n \approx 8.6$$og
$$m \approx 8.6, \quad n \approx 2.3$$... og de gir identiske kurver:
fig,(ax1,ax2)=plt.subplots(nrows=1, ncols=2,figsize=(9,5))
r=np.linspace(0.2, 0.5, 1000)
for m,n in zip((2.3, 8.6),(8.6,2.3)):
A = -(E0*r0**m)/(1-m/n) # ligning (6)
B = (m/n)*A*r0**(n-m) # ligning (5)
E = -A*r**(-m) + B*r**(-n) # ligning (1)
F = m*A*r**(-m-1) -n*B*r**(-n-1) # ligning (2)
k = -(m+1)*m*A*r0**(-m-2) + (n+1)*n*B*r0**(-n-2) # ligning (7)
ax1.plot(r,E,label='m={}, n={}'.format(m,n))
ax1.set_xlabel('$r$ [nm]')
ax1.set_ylabel('$E$ [eV]')
ax1.legend(loc='best')
ax1.set_ylim(E0*1.1,0)
ax1.grid(True)
ax2.plot(r,F ,label='k={:.0f}'.format(k))
ax2.set_xlabel('$r$ [nm]')
ax2.set_ylabel('$F$ [eV]')
ax2.legend(loc='best')
ax2.set_ylim(-5,16)
ax2.grid(True)
plt.tight_layout()
plt.draw()
Energipotensialet uttrykt som
\begin{equation} E = -\frac{A}{r^m} + \frac{B}{r^n} \tag{1} \end{equation}er en alternativ notasjon til Mie-potensialet som mer vanlig blir uttrykt på følgende form:
\begin{equation} E = C E_0 \bigg( \Big( \frac{r_c}{r} \Big)^m - \Big( \frac{r_c}{r} \Big)^n \bigg) \tag{9} \end{equation}hvor
$$C = \frac{n}{n-m}\Big( \frac{n}{m} \Big) ^{\frac{m}{n-m}} $$og parameteren $r_c$ angir hvor energien er null (til venstre for likevektsavstanden).
Sammenhengen mellom parametrene er illustrert i følgende kode:
E0 = -3
r0 = 0.25
m = 4
n = 8
A = -(E0*r0**m)/(1-m/n) # ligning (6)
B = (m/n)*A*r0**(n-m) # ligning (5)
r=np.linspace(0.15, 0.6, 50)
E = -A*r**(-m) + B*r**(-n) # ligning (1)
C = (n/(n-m))*(n/m)**(m/(n-m))
rc = (-A/(E0*C))**(1/m) # Oppgave: vis denne sammenhengen
E_Mie = C*E0*((rc/r)**m - (rc/r)**n) # ligning (9)
fig,ax=plt.subplots(nrows=1, ncols=1,figsize=(8,4))
ax.plot(r,E,color='red',label='(1)')
ax.plot(r,E_Mie,'.', color = 'blue', label='(9)')
ax.plot([rc,rc],[E0*1.1,0],'--', color = 'black') # illustrerer hvor rc er
ax.set_xlabel('$r$ [nm]')
ax.set_ylabel('$E$ [eV]')
ax.legend(loc='best')
ax.grid(True)
ax.set_xlim(0.1,max(r))
ax.set_ylim(E0*1.1,2)
ax.set_title('$C$={}, $r_c$={:.3f}, $A$={:.2E}, $B$={:.2E}'.format(C,rc,A,B), fontsize=10)
plt.tight_layout()
plt.draw()
Et annet energipotensial som er mye brukt er Morse-potensialet som har følgende form:
\begin{equation} E = -E_0 \big(\exp(-2a(r-r_0)-2\exp(-a(r-r_e)) \big) \tag{10} \end{equation}E0 = -3
r0 = 0.25
m = 4
n = 8
A = -(E0*r0**m)/(1-m/n) # ligning (6)
B = (m/n)*A*r0**(n-m) # ligning (5)
r=np.linspace(0.15, 0.6, 100)
E = -A*r**(-m) + B*r**(-n) # ligning (1)
a=16 # tilpasset på øyemål i forhold til (1)
E_Morse = -E0*(np.exp(-2*a*(r-r0)) -2*np.exp(-a*(r-r0)) ) # ligning (10)
fig,ax=plt.subplots(nrows=1, ncols=1,figsize=(8,4))
ax.plot(r,E,color='red',label='(1)')
ax.plot(r,E_Morse, color = 'blue', label='(10)')
ax.set_xlabel('$r$ [nm]')
ax.set_ylabel('$E$ [eV]')
ax.legend(loc='best')
ax.grid(True)
ax.set_xlim(0.1,max(r))
ax.set_ylim(E0*1.1,2)
plt.tight_layout()
plt.draw()
