Anharmonic Einstein model #271

SEckner · 2020-08-10T09:23:51Z

Cumulant formulas as given by Yokoyama (J. Synchrotron Rad. 6, 323-325, 1999) and Haug et al. (Phys. Rev. B 77, 184115, 2008, DOI: 10.1103/PhysRevB.77.184115).

I used it in my phd thesis (DOI: 10.22032/dbt.38705) and two publications (DOI: 10.1103/PhysRevB.97.195202, DOI: 10.1209/0295-5075/126/36002).

Perhaps someone finds it useful.

#!/usr/bin/env python
# anharmonic Einstein model for xafs 
import numpy as np
from larch import ValidateLarchPlugin, use_plugin_path

use_plugin_path('xray')

from xraydb_plugin import atomic_mass
import scipy.constants as consts
# EINS_FACTOR  = hbarc*hbarc/(2 * k_boltz * amu) = 24.254360157751783
#    k_boltz = 8.6173324e-5  # [eV / K]
#    amu     = 931.494061e6  # [eV / (c*c)]
#    hbarc   = 1973.26938    # [eV * A]
EINS_FACTOR = 1.e20*consts.hbar**2/(2*consts.k*consts.atomic_mass)
C3_FACTOR = 1.e40*consts.hbar**6/(consts.k**4*consts.atomic_mass**3) #for k3 to be in kg/(AA * s**2)

@ValidateLarchPlugin
def C3_eins(t, theta,k3, path=None, _larch=None):
    """calculate C3 for a Feff Path wih the einstein model

    third = C3_eins(t, theta, k3, path=None)

    Parameters:
    -----------
      t        sample temperature (in K)
      theta    Einstein temperature (in K)
      k3       cubic force constant
      path     FeffPath to caculate sigma2 for [None]

    if path is None, the 'current path'
    (_sys.paramGroup._feffdat) is used.

    """
    if path is None:
        try:
            path = _larch.symtable._sys.paramGroup
        except:
            pass
    try:
        fdat = path._feffdat
    except:
        return 0.00

    if theta < 1.e-5: theta = 1.e-5
    if t < 1.e-5: t = 1.e-5
    tx = theta/(2.0*t)
    return C3_FACTOR*k3/(theta**4 * fdat.rmass**3) * (3/(2*np.tanh(tx)**2)-1)
#enddef
    
@ValidateLarchPlugin
def C1_anharm(T,k0,k3,k4,terms=4,rmass=None,path=None,_larch=None):
    """calculate anharmonic first cumulant as in Yokoyama99 and Haug08
    
    Parameters:
    ------------
    t       sample temperature (in K)
    k0, k3,k4 force constants
    k0 in N/m, k3 in kg/(s^2 AA)
    path    FeffPath to calculate sigma2 for [None]
    
    if path is None, the 'current path'
    (_sys.paramGroup._feffdat) is used.

    """
    if rmass is None:
        if path is None:
            try:
                path = _larch.symtable._sys.paramGroup
            except:
                pass
        try:
            rmass = path._feffdat.rmass
        except:
            return 0.00
        
    if T < 1.e-5: T = 1.e-5
        
    omega = np.sqrt(k0/(rmass*consts.atomic_mass)) # 1/s
    sig02 = consts.hbar/(2*(rmass*consts.atomic_mass)*omega) #Js/(kg/s) = m2
    beta = 1/(consts.k*T)
    lnz = -beta*consts.hbar*omega
    z = np.exp(lnz)
    
    k3=k3*1.e10 #kg/(s^2 AA)=N/(m * AA) -> kg/(s^2 m)=N/(m^2)
    k4=k4*1.e20 #N/(m * AA^2) -> N/(m^3)
    
    C11 = 6*k3*sig02**2/(consts.hbar*omega)*(1+z)/(1-z)
    C12 = -k3*k4*(consts.hbar*omega)**2/(2*k0**4)*((31+z*(94+31*z))/((1-z)**2)-(18*z*(1+z)*lnz)/((1-z)**3))
    C13a= k3**3*3*(consts.hbar*omega)**2/(2*k0**5)*((11+z*(38+11*z))/((1-z)**2)-(15*z*(1+z)*lnz)/((1-z)**3))
    C13b= k3*k4**2*9*(consts.hbar*omega)**3/(8*k0**6)*(((1+z)*(187+z*(996+187*z)))/((1-z)**3)-(4*z*(64+z*(183+64*z))*lnz)/((1-z)**4)+(12*z*(1+z)*(1+z*(8+z))*lnz**2)/((1-z)**5))
    termlist = [C11,C12,C13a,C13b]
    return sum(termlist[0:terms])*1.e10 #m -> AA
#enddef

@ValidateLarchPlugin
def C2_anharm(T,k0,k3,k4,terms=6,rmass=None,path=None,_larch=None):
    """calculate anharmonic second cumulant as in Yokoyama99 and Haug08
    
    Parameters:
    ------------
    t       sample temperature (in K)
    k0, k3,k4 force constants
    path    FeffPath to calculate sigma2 for [None]
    
    if path is None, the 'current path'
    (_sys.paramGroup._feffdat) is used.

    """
    if rmass is None:
        if path is None:
            try:
                path = _larch.symtable._sys.paramGroup
            except:
                pass
        try:
            rmass = path._feffdat.rmass
        except:
            return 0.00
        
    if T < 1.e-5: T = 1.e-5
        
    omega = np.sqrt(k0/(rmass*consts.atomic_mass)) # 1/s
    sig02 = consts.hbar/(2*(rmass*consts.atomic_mass)*omega) #Js/(kg/s) = m2
    beta = 1/(consts.k*T)
    lnz = -beta*consts.hbar*omega
    z = np.exp(lnz)
    
    k3=k3*1.e10 #kg/(s^2 AA) -> kg/(s^2 m)
    k4=k4*1.e20 #N/(m * AA^2) -> N/(m^3)
    
    C20 = sig02*(1+z)/(1-z)
    C21 = -k4*sig02**3/(consts.hbar*omega)*12*(1+z)**2/((1-z)**2) - k4*sig02**3/(consts.k*T)*24*z*(1+z)/((1-z)**3)
    C22a = k3**2*sig02**4/((consts.hbar*omega)**2)*4*(13*z**2+58*z+13)/((1-z)**2) + k3**2*sig02**4/(consts.hbar*omega*consts.k*T)*120*z*(1+z)/((1-z)**3)
    C22b= k4**2*3*(consts.hbar**3*omega**3)/(8*k0**5)*((5*(1+z)*(7+z*(20+7*z)))/((1-z)**3)-(12*z*(8+z*(21+8*z))*lnz)/((1-z)**4)+(12*z*(1+z)*(1+z*(8+z))*lnz**2)/((1-z)**5))
    C23a= -k3**2*k4*3*(consts.hbar*omega)**3/(8*k0**6)*(((1+z)*(275+z*(1976+275*z)))/((1-z)**3)-(12*z*(71+z*(209+71*z))*lnz)/((1-z)**4)+(60*z*(1+z)*(1+z*(8+z))*lnz**2)/((1-z)**5))
    C23b= -k4**3*9*(consts.hbar*omega)**4/(8*k0**7)*((4*(1+z)**2*(37+z*(176+37*z)))/((1-z)**4)-(z*(1-z)*(535+z*(2072+535*z))*lnz)/((1-z)**5)+(12*z*(6+z*(93+z*(205+93*z+6*z**2)))*lnz**2)/((1-z)**6)-(4*z*(1+z)*(1+z*(38+z*(144+z*(38+z))))*lnz**3)/((1-z)**7))
    termlist = [C20,C21,C22a,C22b,C23a,C23b]
    return sum(termlist[0:terms])*1.e20 #m2 -> AA2
#enddef

@ValidateLarchPlugin
def C3_anharm(T,k0,k3,k4,terms=6,rmass=None,path=None,_larch=None):
    """calculate anharmonic third cumulant as in Yokoyama99 and Haug08
    
    Parameters:
    ------------
    t       sample temperature (in K)
    k0, k3 force constants
    path    FeffPath to calculate C3 for [None]
    
    if path is None, the 'current path'
    (_sys.paramGroup._feffdat) is used.

    """
    if rmass is None:
        if path is None:
            try:
                path = _larch.symtable._sys.paramGroup
            except:
                pass
        try:
            rmass = path._feffdat.rmass
        except:
            return 0.00
        
    if T < 1.e-5: T = 1.e-5
        
    omega = np.sqrt(k0/(rmass*consts.atomic_mass))
    sig02 = consts.hbar/(2*(rmass*consts.atomic_mass)*omega)
    beta = 1/(consts.k*T)
    lnz = -beta*consts.hbar*omega
    z = np.exp(lnz)
    
    k3=k3*1.e10 #kg/(s^2 AA) -> kg/(s^2 m)
    k4=k4*1.e20 #N/(m * AA^2) -> N/(m^3)

    C31 = k3*sig02**3/(consts.hbar*omega)*4*(z**2+10*z+1)/((1-z)**2)
    C32 = -k3*k4*3*(consts.hbar*omega)**3/(4*k0**5)*(((1+z)*(17+z*(224+17*z)))/((1-z)**3)-(12*z*(2+z*(9+2*z))*lnz)/((1-z)**4))
    C33a= k3**3*3*(consts.hbar*omega)**3/(8*k0**6)*(((1+z)*(41+z*(680+41*z)))/((1-z)**3)-(60*z*(2+z*(9+2*z))*lnz)/((1-z)**4))
    C33b= k3*k4**2*(consts.hbar*omega)**4/(16*k0**7)*((4717+z*(92234+z*(215994+z*(92234+4717*z))))/((1-z)**4)-(180*z*(1+z)*(68+z*(497+68*z))*lnz)/((1-z)**5)+(216*z*(2+z*(49+z*(124+z*(49+2*z))))*lnz**2)/((1-z)**6))
    C34a= -k3**3*k4*3*(consts.hbar*omega)**4/(16*k0**8)*((5483+z*(134182+z*(333174+z*(134182+5483*z))))/((1-z)**4)-(36*z*(1+z)*(572+z*(4483+572*z))*lnz)/((1-z)**5)+(360*z*(2+z*(49+z*(124+z*(49+2*z))))*lnz**2)/((1-z)**6))
    C34b= -k3*k4**3*3*(consts.hbar*omega)**5/(32*k0**9)*((5*(1+z)*(15035+z*(378526+z*(1265526+z*(378526+15035*z)))))/((1-z)**5)-(12*z*(23215+z*(310531+z*(668432+z*(310531+23215*z))))*lnz)/((1-z)**6)+(144*z*(1+z)*(97+z*(3564+z*(13435+z*(3564+97*z))))*lnz**2)/((1-z)**7)-(144*z*(2+z*(187+z*(1712+z*(3404+z*(1712+z*(187+2*z))))))*lnz**3)/((1-z)**8))
    termlist = [C31,C32,C33a,C33b,C34a,C34b]
    return sum(termlist[0:terms])*1.e30 #m3 -> AA3
#enddef

@ValidateLarchPlugin
def C4_anharm(T,k0,k3,k4,terms=4,rmass=None,path=None,_larch=None):
    """calculate anharmonic fourth cumulant as in Yokoyama99 and Haug08
    
    Parameters:
    ------------
    t       sample temperature (in K)
    k0, k3, k4 force constants
    path    FeffPath to calculate C3 for [None]
    
    if path is None, the 'current path'
    (_sys.paramGroup._feffdat) is used.

    """
    if rmass is None:
        if path is None:
            try:
                path = _larch.symtable._sys.paramGroup
            except:
                pass
        try:
            rmass = path._feffdat.rmass
        except:
            return 0.00
        
    if T < 1.e-5: T = 1.e-5
        
    omega = np.sqrt(k0/(rmass*consts.atomic_mass))
    sig02 = consts.hbar/(2*(rmass*consts.atomic_mass)*omega)
    beta = 1/(consts.k*T)
    lnz = -beta*consts.hbar*omega
    z = np.exp(lnz)

    k3=k3*1.e10 #kg/(s^2 AA) -> kg/(s^2 m)
    k4=k4*1.e20 #N/(m * AA^2) -> N/(m^3)
    
    C41 = -k4*sig02**4/(consts.hbar*omega)*12*(z**3+9*z**2+9*z+1)/((1-z)**3) - k4*sig02**4/(consts.k*T)*144*z**2/((1-z)**4)
    C42 = k3**2*sig02**5/(consts.hbar**2*omega**2)*12*(5*z**3+109*z**2+109*z+5)/((1-z)**3) + k3**2*sig02**5/(consts.hbar*omega*consts.k*T)*720*z**2/((1-z)**4)
    C43a= -k3**2*k4*9*(consts.hbar*omega)**4/(16*k0**7)*((237+z*(5930+z*(15482+z*(5930+237*z))))/((1-z)**4)-(4*z*(1+z)*(92+z*(1613+92*z))*lnz**2)/((1-z)**5)+(120*z**2*(5+z*(16+5*z))*lnz**2)/((1-z)**6))
    C43b= -k4**3*3*(consts.hbar*omega)**5/(32*k0**8)*(((1+z)*(3707+z*(59310+z*(133046+z*(59310+3707*z)))))/((1-z)**5)-(12*z*(861+z*(13803+z*(29524+3*z*(4601+287*z))))*lnz)/((1-z)**6)+(432*z*(1+z)*(1+z*(82+z*(371+z*(82+z))))*lnz**2)/((1-z)**7)-(144*z**2*(19+z*(230+z*(492+z*(230+19*z))))*lnz**3)/((1-z)**8))
    termlist=[C41,C42,C43a,C43b]
    return sum(termlist[0:terms])*1.e40 #m4 ->AA4
#enddef
    
def registerLarchPlugin():
    return ('_xafs', {  'C3_eins':C3_eins,
                        'C1_anharm':C1_anharm,
                        'C2_anharm':C2_anharm,
                        'C3_anharm':C3_anharm,
                        'C4_anharm':C4_anharm})

The text was updated successfully, but these errors were encountered:

newville · 2020-08-11T14:49:30Z

@SEckner Thanks -- that looks very useful!

Well, OK the math looks a bit challenging, and I did not check it ;). It would be nice to have a realistic check on these calculations so that future changes don't mess up the results. Do you have a good example of these calculations from your dissertation or papers that you could provide?

Is it OK if I make some cleanups and add these functions? I think the code should have your name and links to your papers in it too. Is that acceptable to you?

SEckner · 2020-08-12T07:05:39Z

Please feel free to clean up and to add the functions. Adding my name and a reference to the papers would be nice. Certainly a reference to the papers stating the formulas (Yokoyama + Haug et al.) should go there, too.

I can provide a testcase. But I will need some days since I used a bunch of scripts to comfortably handle temperature and composition series and now have to extract a suitable set of spectra by hand.

newville · 2020-08-12T12:34:33Z

@SEckner Great! I'm in no rush for a test case, it's just that it would be nice as part of the code documentation so that in years to come (when you are showing it to your students, etc), we can all be sure that it is still working correctly.

maurov · 2021-06-14T05:11:18Z

@SEckner any news on this? Are you still willing to provide a testcase?

newville · 2021-06-14T15:32:07Z

@maurov @SEckner perhaps we should just merge this -- tests and example would be nice, but I think it is also important to not lose this.

maurov added the enhancement label Jun 14, 2021

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Anharmonic Einstein model #271

Anharmonic Einstein model #271

SEckner commented Aug 10, 2020

newville commented Aug 11, 2020

SEckner commented Aug 12, 2020

newville commented Aug 12, 2020

maurov commented Jun 14, 2021

newville commented Jun 14, 2021

Anharmonic Einstein model #271

Anharmonic Einstein model #271

Comments

SEckner commented Aug 10, 2020

newville commented Aug 11, 2020

SEckner commented Aug 12, 2020

newville commented Aug 12, 2020

maurov commented Jun 14, 2021

newville commented Jun 14, 2021