This post mainly covers different famous memristor models, window functions along with their SPICE simulations(done on LTSPICE). The window function is used to add non-linearity at the boundaries. In contrast, Models are used to establish a linear/non-linear relationship between the rate of state change and the stimulus (current/voltage).

Content:

• Models
  • Linear
  • Non-linear
  • Threshold
  • Macro-model
• Window functions
  • Properties of window functions
  • Types of window functions:
    • Struckov
    • Benderli
    • Joglekar
    • Biolek
    • Prodomakis
    • Jinxiang
• SPICE simulations

Memristor Models

Linear Ion drift Model (Ideal Model)

In this model, there exists a linear relationship between the state derivative and the stimuli(Voltage/Current). It assumes that vacancies have the freedom to move around the entire length of the device. It is not true since vacancies slow down a lot at boundaries. It’s one of the reasons, along with the tunneling effect, why the linear model can’t be used at higher voltages.

Further, it predicts an inverse relationship between switching time(T_o) and applied voltage(V_o),  stimuli are voltage here. However, actual experiments show a logarithmic inverse relationship.

More details about the linear model can be found here: http://www.iue.tuwien.ac.at/phd/mahmoudi/dissse9.html

Non-Linear Ionic drift model

It is the most accurate model which predicts both static and switching dynamic behavior of a memristor. It proposes a non-linear relationship between stimuli and state derivative with asymmetric switching behavior.f(w) here, is the window function.  i(t) is made up of two components sinh() part that is due to tunneling effect and exp() part that is due to PN junction behavior.

Fig 1 depicts the I-V relationship for a non-linear memristor model when input is a sinusoidal wave and rectangular pulse. This model properly describes the non-linear switching dynamics arising from the ionic motion, which modulates the effective width of the insulating region as well as the electron tunneling effect through the insulating region, which is a function of the width of the insulating region and applied voltage.

Figure 2 shows the variation of memristance with the transition in the state for linear and non-linear models.

Visit here for detailed information on non-linear memristor model: http://www.iue.tuwien.ac.at/phd/mahmoudi/dissse9.html#x21-240003.3.2

Threshold Adaptive memristor (TEAM) Model

It’s a more straightforward and generalized model that approximates the non-linear drift model of the memristor to improve computational efficiency. The crux point here is that below a certain threshold, the state variable doesn’t change.
Instead of exponential dependency, there is a polynomial dependence involved between the memristor current and the internal state drift derivative.

where k_off > 0 and k_on < 0. alpha_off and alpha_on being constants. i_on and i_off are current thresholds.

TEAM is accurate enough with a mean error of 0.2% and can boost simulation run-time by 47.5%. It also satisfies convergence conditions, computational efficiency, and it can fit any of the existing models like the linear ion drift model.

VTEAM is another variation in the TEAM model with some changes in threshold and thus a bit change in performance.

MACRO Model

Macro model of memristor, unlike another previous model, is made up of a circuit consisting of various discrete electrical components that altogether give characteristics like memristor.

Window functions

It’s used to prevent the state variable from getting out of the bounds [0, D] and also to add more non-linearity near the bounds. The derivative of the state variable is multiplied by a window function to force it to reach 0 when the state variable is at the bounds.

Properties of window functions:

1. Boundary conditions – whether the boundary conditions at the top or the bottom electrode of the device are taken into considerations or not.
2. Boundary lock – whether a window function is an able account for boundary lock problem*.
3. The linearity of drift – whether the drift at boundaries non-linear or not.
4. The linkage between linear and non-linear model – whether at low voltages, the non-linear model performs the same as a linear model both with the same window functions.
5. Scalability – whether it provides full scalability or not
6. Flexibility – Whether it provides control parameters for several fine-tuning properties.

* Boundary lock problem – When a model is dependent only on state function and not on stimuli, usually it states can’t be changed from one extreme to another. The blue line in the below images is the memristance value, and the green line is the input voltage waveform. Note that in case of the boundary-lock problem after reaching one extreme, memristance value doesn’t change further irrespective of input.

Several window functions(w.f.)

1.  Strukov’s w.f.
   • 
   • where w is the state variable, D be the low conductivity portion length. w = W/D
2.  Benderli’s w.f.
   • 
   • Where w is the state variable, D be the low conductivity portion length. w = W/D, W be the length of the high-conductivity portion
3.  Joglekar w.f.
   • 
   • Where w is the state variable, p is the parameter that controls the amount of non-linearity.
   • Fig-4(a) shows Joglekar w.f. Characteristics. Here, Boundary conditions are satisfied with a maximum value at the center, with the only problem being boundary lock.
4.  Biolek w.f.
   • 
   • Where w is the state variable, p is the parameter that controls the amount of non-linearity, i is the current flowing.
   • Fig-4(b) shows biolek window function characteristics. Problem with biolek’s w.f. is a discontinuity at boundaries, and it’s not scalable.

     

5.  Prodomakis w.f.
   • 
   • Where w is the state variable, p is the parameter that controls the amount of non-linearity, j is the scalability factor.
   • The main issue with prodomain w.f. is the boundary lock problem.
6.  Jinxiang w.f.
   • 
   • Where w is the state variable, p is the parameter that controls the amount of non-linearity, j is the scalability factor, 0<a<1 a controls the resolution.
   • It solves most of the issues, taking into care boundary lock problems. However, it suffers a disadvantage at high voltages.

One of the recent w.f. As mentioned here, https://doi.org/10.1016/j.aeue.2018.04.020 is by far( as of Jan 2019) the best.

SPICE SIMULATION CODES

Biolek window function

* HP Memristor SPICE Model Using Biolek Window

* Connections:
* TE:  Top electrode
* BE:  Bottom electrode
* XSV: External connection to plot state variable
*      that is not used otherwise

.SUBCKT MEM_BIOLEK TE BE XSV 

* Ron:  Minimum device resistance
* Roff: Maximum device resistance
* D:    Width of the thin film
* uv:   Dopant mobility
* p:    Parameter for window function
* x0:   State varaible initial value

.params Ron=100 Roff=16K x0=.076 D=16N uv=40F p=7

* Biolek Window Function
.func f(V1,I1)={1-pow((V1-stp(-I1)),(2*p))}

* Memristor I-V Relationship
.func IVRel(V1,V2) = V1/(Ron*V2 + Roff*(1-V2))

* Circuit to determine state variable
Gx 0 XSV value={I(Gmem)*Ron*uv*f(V(XSV,0),I(Gmem))/pow(D,2)}
Cx XSV 0 {1}
.ic V(XSV) = x0

* Current source representing memristor
Gmem TE BE value={IVRel(V(TE,BE),V(XSV,0))}

.ENDS MEM_BIOLEK

Joglekar window function

***********************************************
* HP Memristor SPICE Model
* For Transient Analysis only
* created by Zdenek and Dalibor Biolek
***********************************************
* Ron, Roff - Resistance in ON / OFF States
*
* Rinit - Resistance at T=0
*
* D - Width of the thin film
*
* uv - Migration coefficient
*
* p - Parameter of the WINDOW-function for
*     modeling nonlinear boundary conditions
*
* x - W/D Ratio, W is the actual width
*     of the doped area (from 0 to D)
*
***********************************************
.SUBCKT memristor plus minus PARAMS:
+ Ron=100 Roff=16K Rinit=80K D=10N uv=10F p=1
 
***********************************************
* DIFFERENTIAL EQUATION MODELING *
***********************************************
Gx 0 x value={I(Emem)*uv*Ron/D**2*f(V(x),p)}
Cx x 0 1 IC={(Roff-Rinit)/(Roff-Ron)}
Raux x 0 1000000
 
***********************************************
* RESISTIVE PORT OF THE MEMRISTOR *
***********************************************
Emem plus aux value={-I(Emem)*V(x)*(Roff-Ron)}
Roff aux minus {Roff}
 
***********************************************
* FLUX COMPUTATION *
***********************************************
Eflux flux 0 value={SDT(V(plus,minus))}
 
***********************************************
* CHARGE COMPUTATION *
***********************************************
Echarge charge 0 value={SDT(I(Emem))}
 
***********************************************
* WINDOW FUNCTIONS
* FOR NONLINEAR DRIFT MODELING *
***********************************************
* window function, according to Joglekar
.func f(x,p)={1-(2*x-1)**(2*p)}
.ENDS memristor

Prodomakis window function

* HP Memristor SPICE Model Using Prodomakis Window

* Connections:
* TE:  Top electrode
* BE:  Bottom electrode
* XSV: External connection to plot state variable
*      that is not used otherwise

.SUBCKT MEM_PRODOMAKIS TE BE XSV 

* Ron:  Minimum device resistance
* Roff: Maximum device resistance
* D:    Width of the thin film
* uv:   Dopant mobility
* p:    Parameter for window function
* x0:   State varaible initial value

.params Ron=100 Roff=16K x0=.076 D=16N uv=40F p=7 j=1

* prodomakis Window Function
.func f(V1)=  j*{1-pow((pow((V1 - 0.5),2) + 3/4 ),p)}

* Memristor I-V Relationship
.func IVRel(V1,V2) = V1/(Ron*V2 + Roff*(1-V2))

* Circuit to determine state variable
Gx 0 XSV value={I(Gmem)*Ron*uv*f(V(XSV,0))/pow(D,2)}
Cx XSV 0 {1}
.ic V(XSV) = x0

* Current source representing memristor
Gmem TE BE value={IVRel(V(TE,BE),V(XSV,0))}

.ENDS MEM_PRODOMAKIS

TEAM memristor model

*threshold model*
.subckt MEMTHRESH pl mn PARAMS: Ron=1K Roff=16K Rinit=5K beta=1E13 Vtp=4.6 Vtm=4.6 nu1=0.0001 nu2=0.1
Gx 0 x value={f1(V(pl)-V(mn))*(f2(f1(V(pl)-V(mn)))*f3(Roff-V(x))+f2(-f1(V(pl)-V(mn)))*f3(V(x)-Ron))}
Raux x 0 1E12
Cx x 0 1 IC={Rinit}
Gpm pl mn value={(V(pl)-V(mn))/V(x)}
.func f1(y)={beta*(y-Vtp)/(exp(-(y-Vtp)/nu1)+1)+beta*(y+Vtm)/(exp(-(-y-Vtm)/nu1)+1)}
.func f2(y1)={1/(exp(-y1/nu1)+1)}
.func f3(y)={1/(exp(-y/nu2)+1)}
.ends MEMTHRESH

Jinxiang memristor model

.SUBCKT mem_jinxiang plus minus
.param ron = 100
.param roff = 16000
.param rinit = 18000
.param d = 40E-9
.param uv = 10E-15
.param p = 7
.param x0 = 0.076
.param j = 3
.param a = 0.5
 
***********************************************
* DIFFERENTIAL EQUATION MODELING *
***********************************************
Gx 0 x value={I(Emem)*uv*ron/d**2*f(V(x),I(Emem),p)}
Cx x 0 1 IC={(roff-rinit)/(roff-ron)}
Raux x 0 1000000
 
***********************************************
* RESISTIVE PORT OF THE MEMRISTOR *
***********************************************
Emem plus aux value={-I(Emem)*V(x)*(roff-ron)}
roff aux minus {roff}
 
***********************************************
* FLUX COMPUTATION *
***********************************************
Eflux flux 0 value={SDT(V(plus,minus))}
 
***********************************************
* CHARGE COMPUTATION *
***********************************************
Echarge charge 0 value={SDT(I(Emem))}
 
***********************************************
* WINDOW FUNCTIONS
* FOR NONLINEAR DRIFT MODELING *
***********************************************
* window function, according to Jinxiang
*.func f(x,p)={1-pow((V1-stp(-I1)),(2*p))}
.func f(x,i,p)=j*{1-pow((a*pow((x - stp(-i)),(2*p)) + (1 - a)),(p))}
.ENDS mem_jinxiang