Well-posed Memristor Modeling with Xyce and Verilog-A

Are you interested in nanotechnology based AI processors?

Sign up for our newsletter for progress updates, product announcements and news.

Knowm Memristor Array

In this post, I talk about the current state of memristor modeling after having rolled up my sleeves for the last several weeks, getting my hand dirty with as many open source circuit simulators as possible, and simulating memristors using various techniques and approaches. I’ve gained a lot of insight into the current state of the art of memristor modeling that I thought would be worth documenting.

“Hacked” SPICE Subcircuit Memristor Models

Around 2008, the first memristor models designed for SPICE started showing up. A paper by Biolek et al. showed how to implement a memristor model in SPICE as a .subckt. I call this a “hacked” model because in order to implement it you are constrained in using only existing SPICE elements such as resistors, capacitors, and various specialized voltage and current sources, hacking together a model from these pieces. In a previous post of mine The Joglekar Resistance Switch Memristor Model in LTSpice, I explain in great detail how it works, and I demonstrate a working model and simulation in LTSpice.

The following is the SPICE model from the LTSpice link above:

While the simulation was successful, there are some disadvantages of this type of model:

1. The model is unintuitive and hard to understand because it is hacked together.
2. It only works for transient analysis, not DC sweeps.
3. It suffers from convergence issues at initial conditions and edge cases (See Note)
4. Netlist size grows quickly for large circuits with memristors.

Note: Simulating this same mode in NGSpice and JSpice (our internal Java-based SPICE simulator) yielded convergence issues at the first time step. It did however work just fine in Xyce.

Ill-posed Models

In another class of models, the memristor is usually defined in a more flexible and less hacked manner in a language such as Verilog-A, MATLAB or ModSPEC. In such simulation environments, you are given a blank slate to formulate your model using the language’s full API and methods. Inside your model implementation, you are given access to the global voltage drop across the device and from there you define the current. You are also given the simulation time, and you can use this to calculate how long a certain voltage was applied since the last time-step–important in modeling memristors given their memory/flux behavior.

In order to demonstrate this, here is the Joglekar model again, but implemented in Verilog-A and simulated in Xyce. For anyone interested in running this yourself you may want to check out Build Xyce from Source for ADMS Verilog-A Model Integration.

Joglekar Resistance Switch Model in Verilog-A (Xyce Flavor) memristor.va

The following code can be found on the memristor-models-4-all project.

Simulation File (Xyce) memristor_sim.cir

Simulation Commands

Result

Xyce Joglekar Memristor Model

While the simulation was successful, there are some disadvantages of this type of model:

1. It only works for transient analysis, not DC sweeps
2. Need to calculate dt as time_delta in order to perform an integration
3. Must declare a real as the general variable to hold the state of the device
4. Explicitly guarding against non-physical conditions for w

In fact, these issues are only the tip of the iceberg. If you understand the internal workings of the simulator, there are more issues at play here.

Well-posed Models

In May, 2016, Wang and Roychowdhury published an excellent paper titled Well-Posed Models of Memristive Devices, in which they plainly list out all of the problems with the current approaches and implementations of memristor models. In their words:

In spite of seemingly significant efforts to devise memristor models, every model we are aware of in the literature [25–30] suffers from one or more of the above-mentioned types of ill-posedness.

The types are:

Common ill-posedness mechanisms in models include division-by-zero errors, often due to expressions like 1/(x−a), which become unbounded (and feature a “doubly infinite” jump) at x=a; the use of log(), 1/(x−a) or sqrt() without ensuring that their arguments are always positive, regardless of bias; the fundamental misconception that non-real (i.e., complex) numbers or infinity are legal values for device models (they are not!); and “sharp”/“pointy” functions like |x|, whose derivatives are not continuous. Another key aspect of well-posedness is that the model’s equations must produce mathematically valid outputs for any mathematically valid input to the device. Possibly the most basic kind of input is one that is constant (“DC”) for all time. DC solutions (“operating points”) are fundamental in circuit design; they are typically the kind of solution a designer seeks first of all, and are used as starting points for other analyses like transient, small signal AC, etc. If a model’s equations do not produce valid DC outputs given DC inputs, it fails a very fundamental well-posedness requirement. For example, the equation d/dt o(t) = i(t) is ill posed, since no DC (constant) solution for the output o(t) is possible if the input i(t) is any non-zero constant. Such ill posedness is typically indicative of some fundamental physical principle being violated; for example, in the case of d/dt o(t) = i(t), the system is not strictly stable [24]. Indeed, a well-posed model that is properly written and implemented should work consistently in every analysis (including DC, transient, AC, etc.).

They do a good job in explaining the problems in detail, while also proposing a better way to approach memristor model development at a fundamental level, which I won’t repeat here. Rather, I take their Verilog-A/Xyce example for a generic memristor directly from their paper, add some necessary tweaks, and show the results.

Generic Model in Verilog-A (Xyce Flavor) hys.va

The following code can be found on the memristor-models-4-all project.

Forward DC Sweep test_hys_forward.cir

Reverse DC Sweep test_hys_reverse.cir

Transient test_hys_transient.cir

Simulation Commands

Xyce Well-posed Memristor Model

The currents are negative, and the reverse DC sweep X-axis and Y-Axis should be flipped, but the results are the same as shown in the above referenced paper by Wang et al. You can see here that the results are valid for both transient and DC analysis!

The strategy implemented in the above model is explained as:

To circumvent these issues and write a robust Verilog-A model for hys_example that should work consistently in all simulators and all simulation algorithms, we model state variable s as a voltage. We declare an internal branch, whose voltage represents s. One end of the branch is an internal node that doesn’t connect to any other branches. In this way, by contributing V − s^3 + s and ddt(-tau * s) both to this same branch, the KCL at the internal node will enforce the implicit differential. Declaring s as a voltage is not the only way to model hys_example in Verilog-A. Depending on the physical nature of s, one can also use Verilog-A’s multiphysics support and model it as a mechanical property, such as a position from the kinematic discipline. This may be closer to the actual meaning of s for MIM-structured RAM devices. Alternatively, we can also use the property for potential from the thermal or magnetic discipline. One can also switch potential and flow by defining s as a flow instead. These alternatives may make the model look more physical, but they do not make a difference mathematically, except from the scale of tolerances in each discipline, which we will discuss in more detail in Sec. IV-C. The essence of our approach is to recognize that state variable s is a circuit unknown, and thus should be modeled as a potential or flow in Verilog-A, for the consistent support from different simulators in various circuit analyses.

Again, you should read the full paper. Hats off to Wang and Roychowdhury, as this is a push in the right direction that the memristor modeling community needs.

Homotopy Analysis

Continuation algorithms, also known as homotopy algorithms are normally used to sweep the magnitudes of independent voltage and current sources when a DC operating point calculation fails to converge. Homotopy analysis can track the DC solution curve in the state space. Xyce will apply source stepping automatically during the DC operating point calculation if both the standard Newton’s method and GMIN stepping fail. So, normally, this technique of manually forcing stepping of a single source is unnecessary. But we can actually use it to take a look at some of the parameters in our memristor model as pointed out in the Wang and Roychowdhury paper. For completeness, I demonstrate this as well. As you will see, the results show that as all DC points, the model is continuous.

From the Xyce 6.6 User Manual:

Another basic type of continuation in Xyce is accessed by setting .options nonlin continuation=1. This allows the user to sweep existing device parameters (models and instances), as well as a few reserved artificial parameter cases. The most obvious natural parameter to use is the magnitude(s) of independent voltage or current sources, the choice of which is equivalent to “source stepping” in SPICE. Section 8.2 provides a Xyce source-stepping example. For some circuits (as in the aforementioned Schmidt trigger), source stepping leads to turning points in the continuation.

Homotopy test_hys_homotopy.cir

Simulation Commands

Results

Current	State

Indeed the DC analysis forms a smooth curve in the state space. The corresponding curve in current consists of all the unstable solutions in the i–v relationship (wherever the slope of S is negative).

This curve was previously missing in DC and transient sweep results, and now displayed by the homotopy analysis. These results from homotopy analysis provide us with important insights into the model. They reveal that there is a single smooth and continuous DC solution curve in the state space, which is an indicator of the well-posedness of the model. They also illustrate that it is the folding in the smooth DC solution curve that has created the discontinuities in DC sweep results. These insights are important for the proper modelling of hysteresis. Moreover, the top view explains the use of internal state s for modeling hysteresis from another angle. Without the internal state, it would be difficult if not impossible to write a single equation describing the i–v relationship. With the help of s, we can easily choose two simple model equations, and the complex i–v relationship forms naturally.

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