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 si\mulator’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 .

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 , 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 (temperature) is also a parameter hidden in .

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

If the probabilities and are defined as:

and we define a function as a Gaussian sample, such as the Box-\muller Transform:

, then

. If , , ,

, then

. Therefore plugging everything in (really ugly!), also setting :

.

The problem here, as you can see in the above mess, is that there is no way to refor\mulate that in terms of . 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 si\mulator 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 si\mulation engines, we can refor\mulate 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 equation now looks like this:

The change in the number of switches, scaled between 0 and 1, is:

If the probabilities and are defined as:

. The number of switches switching state is thus:

. Therefore plugging everything in, also setting :

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 and . 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 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

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 ,

Finally, by Ohm’s Law the current is

Preliminary Results in JSPICE

We call this new version of the MSS model the MMSS model, which stands for Mean Metastable Switch. Si\mulations comparing the original MSS and the new MMSS model in our own custom circuit si\mulator 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

Master::updateState 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; } 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.

Sacado 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; } 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 Si\mulation Results

The Xyce si\mulation file used to test the new MSS device follows.

* 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 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 , integration into Xyce was possible and pretty straight forward. The following screenshot shows the matching results from the above custom Java si\mulation.

Generalized Metastable Switch Model

In the Generalized model, the total current through the device comes from both a memory-dependent current component (MSS), , and a Schottky diode current, in parallel:

, where . A value of 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, , 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:

, where , , , and 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:

Sacado Templates template <typename ScalarT> ScalarT SchottkyCurrent( const ScalarT & V1, const ScalarT & V2, double SchottkyForward\alpha, double SchottkyForward\beta, double SchottkyReverse\alpha, double SchottkyReverse\beta ) { return SchottkyReverse\alpha * (-1 * exp(-1 * SchottkyReverse\beta * (V1-V2))) + SchottkyForward\alpha * (exp(SchottkyForward\beta * (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 SchottkyForward\alpha, double SchottkyForward\beta, double SchottkyReverse\alpha, double SchottkyReverse\beta ){ ScalarT Gval= Geff( X, RON, ROFF ); ScalarT MSSCurrentval = (V1-V2)*Gval; ScalarT SchottkyCurrentval = SchottkyCurrent(V1, V2, SchottkyForward\alpha, SchottkyForward\beta, SchottkyReverse\alpha, SchottkyReverse\beta); ScalarT fval = PHI * MSSCurrentval + (1 - PHI) * SchottkyCurrentval; return fval; } 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 SchottkyForward \ alpha , double SchottkyForward \ beta , double SchottkyReverse \ alpha , double SchottkyReverse \ beta ) { return SchottkyReverse \ alpha * ( - 1 * exp ( - 1 * SchottkyReverse \ beta * ( V1 - V2 ) ) ) + SchottkyForward \ alpha * ( exp ( SchottkyForward \ beta * ( 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 SchottkyForward \ alpha , double SchottkyForward \ beta , double SchottkyReverse \ alpha , double SchottkyReverse \ beta ) { ScalarT Gval = Geff ( X , RON , ROFF ) ; ScalarT MSSCurrentval = ( V1 - V2 ) * Gval ; ScalarT SchottkyCurrentval = SchottkyCurrent ( V1 , V2 , SchottkyForward \ alpha , SchottkyForward \ beta , SchottkyReverse \ alpha , SchottkyReverse \ beta ) ; ScalarT fval = PHI * MSSCurrentval + ( 1 - PHI ) * SchottkyCurrentval ; return fval ; }

After recompiling and running the previous si\mulation with the modified model:

.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) 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