// Copyright 2020 The SkyWater PDK Authors
//
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
//     https://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.
//
// SPDX-License-Identifier: Apache-2.0

`include "constants.vams"
`include "disciplines.vams"

module rram2(TE, BE);
    inout TE; // top electrode
    inout BE; // bottom electrode
    electrical TE, BE;

    // ==================================================
    // input parameters
    // --------------------------------------------------
    // physical area and thickness
    parameter real area_ox             = 0.1024e-12 from (0:inf);                         // area of TE/BE overlap [meters^2]
    parameter real Tox                 = 5.0e-9     from (0:inf);                         // thickness of oxide between TE and BE [meters]
    parameter real Tfilament_max       = 4.9e-9     from (0:inf);                         // maximum thickness of conductive filament (for minimum resistance) [meters]
    parameter real Tfilament_min       = 3.3e-9     from (0:inf);                         // minimum thickness of conductive filament (for maximum resistance) [meters]
    parameter real Tfilament_0         = 3.3e-9     from [Tfilament_min:Tfilament_max];   // initial thickness of conductive filament (at t=0 for transient simulation) [meters]
    // activation energy
    parameter real Eact_generation     = 1.501      from (0:inf);                         // activation energy for vacancy generation [eV]
    parameter real Eact_recombination  = 1.500      from (0:inf);                         // activation energy for vacancy recombination [eV]
    // calibration parameters: I-V
    parameter real I_k1                = 6.140e-5   from (0:inf);                         // current calibration parameter [Amps]
    parameter real Tfilament_ref       = 4.7249e-9  from (0:inf);                         // filament thickness calibration parameter [meters]
    parameter real V_ref               = 0.430      from (0:inf);                         // voltage calibration parameter [Volts]
    // calibration parameters: filament growth
    parameter real velocity_k1         = 150        from (0:inf);                         // velocity calibration parameter [meters/second]
    parameter real gamma_k0            = 16.5       from (0:inf);                         // enhancement factor calibration parameter [unitless]
    parameter real gamma_k1            = -1.25      from (-inf:inf);                      // enhancement factor calibration parameter [unitless]
    // calibration parameters: temperature
    parameter real Temperature_0       = 300        from (0:inf);                         // intial temperature [degrees Kelvin]
    parameter real C_thermal           = 3.1825e-16 from (0:inf);                         // effective thermal capacitance [Joules/Kelvin]
    parameter real tau_thermal         = 0.23e-9    from (0:inf);                         // effective thermal time constant [seconds]
    // simulation control
    parameter real t_step              = 1.0e-9     from (0:inf);                         // maximum time step [seconds]
    // ==================================================

    // internal parameters
    real a0 = 0.25e-9;           // atomic distance [m]
    real Tfilament_current;      // current filament thickness [m]
    real Tfilament_dTdt;         // current filament thickness, derivative w.r.t. time [m/s]
    real gamma;                  // local enhancement factor
    real Temperature_current;    // current temperature
    real kT_over_q;              // e.g., 0.0259 at 300 degrees K [eV]

    // transient simulation parameters
    real t_current;              // current time step [s]
    real t_previous;             // previous time step [s]
    real t_delta;                // difference between current vs. previous time step

    // local functions
    analog function real soft_minmax;
        input x, xmin, xmax;
        real x, xmin, xmax;
        real p = 10;
        begin
            if (x > xmax) begin
                soft_minmax = xmax;
            end else if (x < xmin) begin
                soft_minmax = xmin;
            end else begin
                soft_minmax = x;
            end
        end
    endfunction // soft_minmax

    // core equations
    analog begin
        @(initial_step) begin
            // initial condition
            Temperature_current = Temperature_0;
            Tfilament_current = Tfilament_0;
	end
        $bound_step(t_step); // bound maximum time step
        t_current = $abstime; // current time
        t_delta = t_current - t_previous;
        gamma = gamma_k0 + gamma_k1 * pow((Tox - Tfilament_current)/1.0e-9, 3);
	kT_over_q = (`P_K * Temperature_current) / `P_Q;
        Tfilament_dTdt = velocity_k1 * (exp(-Eact_generation    / kT_over_q) * exp( gamma * a0/Tox * V(TE,BE) / kT_over_q) -
                                        exp(-Eact_recombination / kT_over_q) * exp(-gamma * a0/Tox * V(TE,BE) / kT_over_q));
        Tfilament_current = Tfilament_current + Tfilament_dTdt * t_delta; // 1st-order update to filament thickness
        Tfilament_current = soft_minmax(Tfilament_current, Tfilament_min, Tfilament_max); // bound filament thickness
        I(TE,BE) <+ I_k1 * exp(-(Tox - Tfilament_current)/(Tox - Tfilament_ref)) * sinh( V(TE,BE)/V_ref );
        Temperature_current = (Temperature_current + t_delta * (abs(V(TE,BE)*I(TE,BE)) / C_thermal + Temperature_0/tau_thermal))
                              / (1 + t_delta/tau_thermal); // 1st-order update to temperature
	t_previous = $abstime; // current time step is previous time step for next iteration
    end
endmodule