In a previous post titled Native Memristor Device Development in Xyce, I showed how to create the Joglekar memristor model natively in Xyce, by first cloning an existing Yakopcic device model followed by implementing the Joglekar’s `dX/dt`

behavior. In this post, I show how we first modified our current memristor model – the Metastable Switch Model – in order to fit into the Xyce simulator’s constraints followed by the code showing the implementation. Full source code for the model can be found at the memristor-models-4-all Project on Github.

## Memristor Modeling in Xyce

In order to be more in-line with conventional memristor models, and to tie into the “Xyce Way”, we needed to adapt our current MSS model in 4 ways:

- Define on and off conductances as resistances
- Define the two main resistance levels as On and Off rather than A and B
- Define the main state of the device as
`X [0,1]`

, rather than`N_off`

. It’s the same but it’s scaled to one and is continuous. - Define the state, X, in terms of [latex] dX/dt [/latex].

Defining On and Off aren’t strictly necessary for Xyce, but much easier for everyone else to understand the model at first sight, as it’s become the standard convention in memristor modeling. For a different reference to see how other groups define [latex] dX/dt [/latex], see the Yakopcic Model, equation 7.

### Model Conversion to dX/dt

#### Parameters

First off, here is a summary of the model parameters. Note that for now, I removed the two oppositely connected Schottky diodes in parallel with the collection of metastable switches.

Name | Units | Description |
---|---|---|

ROFF | Ohm | Off resistance |

RON | Ohm | On resistance |

VOFF | Volts | Threshold voltage to turn device off |

VON | Volts | Threshold voltage to turn device on |

TAU | Seconds | Time constant |

Note that [latex] T [/latex] (temperature) is also a parameter hidden in [latex] \beta = \frac{q}{{kT}} = {({V_{T}})^{-1}} [/latex].

Starting with what we already had, the change in the number of switches, scaled between 0 and 1, [latex] dX [/latex] is:

[latex isblock=true] dX = \frac{N_{OFF \rightarrow ON} – N_{ON \rightarrow OFF} }{N} [/latex]

If the probabilities [latex] P_{OFF \rightarrow ON} [/latex] and [latex] P_{ON \rightarrow OFF} [/latex] are defined as:

[latex isblock=true] P_{OFF \rightarrow ON} = \alpha \frac{1}{{1 + {e^{ \beta \left( {V – {V_{ON}}} \right)}}}} [/latex]

[latex isblock=true] P_{ON \rightarrow OFF} = \alpha \left(1 – \left(\frac{1}{{1 + {e^{ \beta \left( {V + {V_{OFF}}} \right)}}}}\right)\right) [/latex]

and we define a function [latex] N’ [/latex] as a Gaussian sample, such as the Box-Muller Transform:

[latex isblock=true] N’\left( {\mu ,{\sigma ^2}} \right) = Some algo returning a random sample [/latex]

, then

[latex isblock=true] dX = \frac{N’\left( {\mu_{OFF} ,{\sigma_{OFF \rightarrow ON} ^2}} \right) – N’\left( {\mu_{ON} ,{\sigma_{ON \rightarrow OFF} ^2}} \right) }{N} [/latex]

. If [latex] \mu_{OFF} = N-X N [/latex], [latex] \mu_{ON} = X N [/latex], [latex] \sigma_{OFF \rightarrow ON} = N P_{ON \rightarrow OFF}(1-P_{OFF \rightarrow ON}) [/latex], [latex] \sigma_{ON \rightarrow OFF} = N P_{ON \rightarrow OFF}(1-P_{ON \rightarrow OFF}) [/latex]

, then

[latex isblock=true] dX = \frac{N’\left( {N-X N ,{N P_{OFF \rightarrow ON}(1-P_{OFF \rightarrow ON}) ^2}} \right) – N’\left( {X N ,{N P_{ON \rightarrow OFF}(1-P_{ON \rightarrow OFF} ^2}} \right) }{N} [/latex]

. Therefore plugging everything in (really ugly!), also setting [latex] \alpha = \frac{dT}{\tau} [/latex]:

[latex isblock=true] dX = \frac{N’\left( {N-X N ,{N \frac{dT}{\tau} \frac{1}{{1 + {e^{ \beta \left( {V – {V_{ON}}} \right)}}}}(1-\frac{dT}{\tau} \frac{1}{{1 + {e^{ \beta \left( {V – {V_{ON}}} \right)}}}}) ^2}} \right) – N’\left( {X N ,{N \frac{dT}{\tau} \left(1 – \left(\frac{1}{{1 + {e^{ \beta \left( {V + {V_{OFF}}} \right)}}}}\right)\right)(1-\frac{dT}{\tau} \left(1 – \left(\frac{1}{{1 + {e^{ \beta \left( {V + {V_{OFF}}} \right)}}}}\right)\right) ^2}} \right) }{N} [/latex]

.

The problem here, as you can see in the above mess, is that there is no way to reformulate that in terms of [latex] dX/dt [/latex]. The MSS model was developed specifically to investigate the behavior of our learning circuits as memristors transitioned from well-behaved incremental devices to low-precision stochastic or binary devices. It can do this because of the calls the the random number generators. Unfortunately, this means that the rate of change of the device conductance is non-deterministic and the ‘dt’ is wrapped up in the function call to the random number generator. While this is a more realistic and flexible device model, it is not compatible with some circuit simulator engines. Fortunately, the derivation of the MSS model gives us an easy way to fix this problem.

## The Mean Metastable Switch Memristor Model

Given the complications introduced from the Gaussian distribution sampling in regards to the constraints of the Xyce simulation engines, we can reformulate the model by using the *mean* of the Gaussian sample rather than the sample directly. That is, we take the limit of the MSS model as the number of switches goes to infinity. The [latex] dX/dt [/latex] equation now looks like this:

The change in the number of switches, scaled between 0 and 1, [latex] dX [/latex] is:

[latex isblock=true] dX = N_{OFF \rightarrow ON} – N_{ON \rightarrow OFF} [/latex]

If the probabilities [latex] P_{OFF \rightarrow ON} [/latex] and [latex] P_{ON \rightarrow OFF} [/latex] are defined as:

[latex isblock=true] P_{OFF \rightarrow ON} = \alpha \frac{1}{{1 + {e^{ -\beta \left( {V – {V_{ON}}} \right)}}}} [/latex]

[latex isblock=true] P_{ON \rightarrow OFF} = \alpha \left(1 – \left(\frac{1}{{1 + {e^{ -\beta \left( {V + {V_{OFF}}} \right)}}}}\right)\right) [/latex]

. The number of switches switching state is thus:

[latex isblock=true] N_{OFF \rightarrow ON} = P_{OFF \rightarrow ON} \cdot (1-X) [/latex]

[latex isblock=true] N_{ON \rightarrow OFF} = P_{ON \rightarrow OFF} \cdot (X) [/latex]

. Therefore plugging everything in, also setting [latex] \alpha = \frac{dT}{\tau} [/latex]:

[latex isblock=true] dX = P_{OFF \rightarrow ON} \cdot (1-X) – P_{ON \rightarrow OFF} \cdot X [/latex]

[latex isblock=true] dX = \frac{dt}{\tau} \frac{1}{{1 + {e^{ -\beta \left( {V – {V_{ON}}} \right)}}}} \cdot (1-X) – \frac{dt}{\tau} \left(1 – \left(\frac{1}{{1 + {e^{ -\beta \left( {V + {V_{OFF}}} \right)}}}}\right)\right) \cdot X [/latex]

[latex isblock=true] dX/dt = 1/\tau \left[ \frac{1}{{1 + {e^{ -\beta \left( {V – {V_{ON}}} \right)}}}} \cdot (1-X) – \left(1 – \left(\frac{1}{{1 + {e^{ -\beta \left( {V + {V_{OFF}}} \right)}}}}\right)\right) \cdot X \right] [/latex]

This is now in a form that Xyce can process.

## Visualization and Sanity Check

Before continuing with the implementation, it’s worth taking a look at some of the equations above to get a better intuitive sense of what’s going on.

### Probability as a Function of Voltage

The following plot shows the Logistics Function (standard) and the two scaled versions of it adapted to [latex] P_{ON \rightarrow OFF} [/latex] and [latex] P_{OFF \rightarrow ON} [/latex]. We can see that the probabilities are between 0 and 1 and that as the voltages approach the threshold voltages (+0.27 and -0.27 V) the probabilities are 0.5. Going beyond the threshold voltages in either direction drive the probabilities towards 1.

Note that this plot looks very similar to the window function you will see in descriptions of other memristor models.

### DX/dt as a Function of X at Different Applied Voltages

Here we plot the [latex] dX/dt [/latex] equation as a function of x at different applied voltages. The following plot shows a the correct behavior desired from the model:

- If V=0, no change in X occurs
- If V=+0.5V and X = 1, there will no more change (it’s already max.)
- If V=+0.5V and X = 0, X will go to 1
- If V=-0.5V and X = 1, X will go to 0 (dX=-1)
- If V=-0.5V and X = 0, there will no more change (it’s already min.)

The same behavior occurs at the threshold voltages as well except that instead of at +/-1 change it’s at +/-0.5

One nice side effect of this model is that there are no extra checks on the bounds of X necessary like some other models. This is due to the Logistics function being bounded between 0 and 1.

## Current as a Function of X

The final step to complete the MMSS model description is to define the current as a function of X. We can break that definition down into two steps, first the conductance, followed by the current. The conductance as a function of X is

[latex isblock=true] G = \frac{X}{ R_{ON} }+ \frac{1 – X}{R_{OFF}} [/latex]

Note that this describes a summation of conductances. Relating this to the real-world it tells us that our model is in fact 2 parallel resistors whose resistance values are coupled to each other via X.

For completeness, given a memristance value, for example [latex]R_{init}[/latex],

[latex isblock=true] X = \frac{R_{ON} (R_{init} – R_{OFF} )}{ R_{init} (R_{ON} – R_{OFF})} [/latex]

Finally, by Ohm’s Law the current is

[latex isblock=true] I = V \cdot G [/latex]

## Preliminary Results in JSPICE

We call this new version of the MSS model the MMSS model, which stands for Mean Metastable Switch. Simulations comparing the original MSS and the new MMSS model in our own custom circuit simulator written is Java look the same, which is to be expected for large values of switches.

Given the results in JSPICE, I was confident that we could get this to work in Xyce.

## Xyce Implementation

Most of the changes that were necessary for implementing this device in Xyce are exactly the same as in the previous post: Native Memristor Device Development in Xyce. The interesting portions of this device implementation are in the function `Master::updateState`

and in the Sacdo templates.

### Master::updateState

In this method I setup the two method calls to both `I_V`

and `dXdt`

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 |
bool Master::updateState(double * solVec, double * staVec, double * stoVec) { for (InstanceVector::const_iterator it = getInstanceBegin(); it != getInstanceEnd(); ++it) { Instance & ri = *(*it); double v_pos = solVec[ri.li_Pos]; double v_neg = solVec[ri.li_Neg]; double x = solVec[ri.li_x]; if (DEBUG_DEVICE){ Xyce::dout() << " x_before = " << x << std::endl; } { Sacado::Fad::SFad<double,3> varV1( 3, 0, v_pos ); Sacado::Fad::SFad<double,3> varV2( 3, 1, v_neg ); Sacado::Fad::SFad<double,3> varX( 3, 2, x ); Sacado::Fad::SFad<double,3> resultFad; resultFad = I_V( varV1, varV2, varX, ri.model_.Ron_, ri.model_.Roff_); ri.i0 = resultFad.val(); // current ri.G = resultFad.dx(0); // di/dv = conductance ri.dIdx = resultFad.dx(2); // di/dx } { // evaluate the state variable equation Sacado::Fad::SFad<double,3> varV1( 3, 0, v_pos ); Sacado::Fad::SFad<double,3> varV2( 3, 1, v_neg ); Sacado::Fad::SFad<double,3> varX( 3, 2, x ); Sacado::Fad::SFad<double,3> resultFad; // TODO make thermal voltage dynamic double vT = 0.026; resultFad = dXdt( varV1, varV2, varX, ri.model_.Ron_, ri.model_.Roff_, ri.model_.Von_, ri.model_.Voff_, ri.model_.Tau_, vT ); ri.xVarFContribution = resultFad.val(); if( getSolverState().dcopFlag ) { ri.xVarFContribution = 0; } ri.dxFEqdVpos = resultFad.dx(0); ri.dxFEqdVneg = resultFad.dx(1); ri.dxFEqdx = resultFad.dx(2); } } return true; } |

### Sacado Templates

Here are the Sacado templates used, which implement the bulk of the model.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 |
template <typename ScalarT> ScalarT p0ff2on( const ScalarT & V1, const ScalarT & V2, double VON, double TC, double VT ) { ScalarT exponent = -1 * ((V1-V2) - VON) / VT; double alpha = 1 / TC; ScalarT fval = alpha / (1.0 + exp(exponent)); return fval; } template <typename ScalarT> ScalarT pOn2Off( const ScalarT & V1, const ScalarT & V2, double VOFF, double TC, double VT ) { ScalarT exponent = -1 * ((V1-V2) + VOFF) / VT; double alpha = 1 / TC; ScalarT fval = alpha * (1.0 - 1.0 / (1.0 + exp(exponent))); return fval; } template <typename ScalarT> ScalarT dXdt( const ScalarT & V1, const ScalarT & V2, const ScalarT & X, double RON, double ROFF, double VON, double VOFF, double TC, double VT) { // Probabilities ScalarT p0ff2onVal = p0ff2on(V1, V2, VON, TC, VT); ScalarT pOn2OffVal = pOn2Off(V1, V2, VOFF, TC, VT); // Number of switches making a transition ScalarT n0ff2on = (1 - X) * p0ff2onVal; ScalarT nOn2Off = X * pOn2OffVal; ScalarT fval = n0ff2on - nOn2Off; return fval; } |

## Xyce Simulation Results

The Xyce simulation file used to test the new MSS device follows.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 |
* Voltage Sources V1 1 0 SIN(0V 0.5V 100Hz) * Memristors YKNOWM mr1 1 0 mrm1 Rinit=500 .MODEL mrm1 knowm (level=1 Roff=1500 Ron=500 Voff=0.27 Von=0.27 TAU=0.0001) * Analysis Command .TRAN .1ms .04s * Output .PRINT TRAN V(1) I(V1) N(YKNOWM!mr1:R) .END |

Once rephrasing the MSS model in terms of [latex] dX/dt [/latex], integration into Xyce was possible and pretty straight forward. The following screenshot shows the matching results from the above custom Java simulation.

## Generalized Metastable Switch Model

In the **Generalized** model, the total current through the device comes from both a memory-dependent current component (MSS), [latex] I_{m} [/latex], and a Schottky diode current, [latex] I_{s} [/latex] in parallel:

[latex isblock=true] I=\phi I_{m}(V,t)+(1-\phi)I_{s}(V) [/latex]

, where [latex] \phi\in{[0,1]} [/latex]. A value of [latex] \phi=1 [/latex] represents a device that contains no Schottky diode effects.

We added the diode current component after realizing that the MSS component alone didn’t allow enough flexibility to fit a wide range of existing devices. After adding the diode component, we were able to fit data of a diverse range of memristors to the general model.

The Schottky component, [latex] I_{s}(V) [/latex], follows from the fact that many memristive devices contain a Schottky barrier formed at a metal–semiconductor junction. The Schottky component is modeled by forward bias and reverse bias components as follows:

[latex isblock=true] I_{s}=\alpha_{f}e^{\beta_{f}V}-\alpha_{r}e^{-\beta_{r}V} [/latex]

, where [latex] \alpha_{f} [/latex], [latex] \beta_{f} [/latex], [latex] \alpha_{r} [/latex], and [latex] \beta_{r} [/latex] are positive valued parameters setting the exponential behavior of the forward and reverse biases exponential current flow across the Schottky barrier.

Adding the Schottky component to our device model required minor effort including adding 5 new parameters, and modifying the `I_V`

Sacado template as shown here:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 |
template <typename ScalarT> ScalarT SchottkyCurrent( const ScalarT & V1, const ScalarT & V2, double SchottkyForwardAlpha, double SchottkyForwardBeta, double SchottkyReverseAlpha, double SchottkyReverseBeta ) { return SchottkyReverseAlpha * (-1 * exp(-1 * SchottkyReverseBeta * (V1-V2))) + SchottkyForwardAlpha * (exp(SchottkyForwardBeta * (V1-V2))); } template <typename ScalarT> ScalarT Geff( const ScalarT & X, double RON, double ROFF ) { return X / RON + (1 - X) / ROFF; } template <typename ScalarT> ScalarT I_V( const ScalarT & V1, const ScalarT & V2, const ScalarT & X, double RON, double ROFF, double PHI, double SchottkyForwardAlpha, double SchottkyForwardBeta, double SchottkyReverseAlpha, double SchottkyReverseBeta ){ ScalarT Gval= Geff( X, RON, ROFF ); ScalarT MSSCurrentval = (V1-V2)*Gval; ScalarT SchottkyCurrentval = SchottkyCurrent(V1, V2, SchottkyForwardAlpha, SchottkyForwardBeta, SchottkyReverseAlpha, SchottkyReverseBeta); ScalarT fval = PHI * MSSCurrentval + (1 - PHI) * SchottkyCurrentval; return fval; } |

After recompiling and running the previous simulation with the modified model:

1 2 |
.MODEL mrm1 knowm (level=1 Roff=1500 Ron=500 Voff=0.27 Von=0.27 TAU=0.0001 PHI=0.8 SFA=0.0008 SFB=4 SRA=0.0008 SRB=4) |

, the following hysteresis plot emerges.

We can see that the addition of the forward and reverse Schottky diodes causes the shape of the forward and reverse current to take on a more exponential behavior. With this addition, fitting a very wide range of devices is possible, making it a good general model.

## Conclusion

In this post I introduced the mathematics and implementation as a native device in Xyce of our generalized mean metastable switch (MMSS) memristor model. As a reminder, the model source files can be found at the memristor-models-4-all project on Github. Be sure to stay tuned by signing up for our newsletter and/or subscribing to our RSS feed to get the latest information about our memristor modeling progress. In the future, we will be blogging about more SPICE tools, other model implementations, and other memristor modeling advancements!

## Further Resources

- Knowm Memristors
- The Generalized Metastable Switch Memristor Model
- The Problem is Not HP’s Memristor–It’s How They Want To Use It
- The Joglekar Resistance Switch Memristor Model in LTSpice
- Build Xyce from Source for ADMS Verilog-A Model Integration
- The Pershin Voltage Threshold Memristor Model in NGSpice
- memristor-models-4-all Project on Github
- Well-posed Memristor Modeling with Xyce and Verilog-A
- Native Memristor Device Development in Xyce

## Subscribe To Our Newsletter

Join our low volume mailing list to receive the latest news and updates from our team.