Perform Monte Carlo simulation in ngspice
* 25 stage Ring-Osc. BSIM3 or 4 with statistical variation of model parameters
* Model parameters are varied according to the PDK selection.
* Tested with 3 different commercial HSPICE libraries from 2 vendors.
* Add your library to mc_ring_circ.net and choose transistors accordingly.
* Add the library path to the .LIB statement.
* A simple BSIM3 inverter R.O. serves as an MC example.

.options noacct

vin in out dc 0.5 pulse 0.5 0 0.1n 5n 1 1 1
vdd dd 0 dc 3.3
vss ss 0 dc 0
ve  sub  0 dc 0
vpe well 0 dc 3.3

* transistors to be selected according to the library (here: p33ll and n33ll or pch_5_mac and nch_5_mac
* or pe3 and ne3 or p1 and n1 (these models see below))
.subckt inv1 dd ss sub well in out
*XMP1 out in  dd  well p33ll w=5u l=800n m=3 nf=1 ad=1.35p as=1.35p pd=9.6u ps=9.6u mosmis_mod=1
*XMN1 out in  ss  sub n33ll w=5u l=800n m=1 nf=3 ad=0.9p as=0.9p pd=6.6u ps=6.6u mosmis_mod=1
*XMP1 out in  dd  well pch_5_mac w=5u l=800n m=3 nf=1 ad=1.35p as=1.35p pd=9.6u ps=9.6u mosmis_mod=1
*XMN1 out in  ss  sub nch_5_mac w=5u l=800n m=1 nf=3 ad=0.9p as=0.9p pd=6.6u ps=6.6u mosmis_mod=1
*XMP1 out in  dd  well pe3 w=5u l=800n m=3 nf=1 ad=1.35p as=1.35p pd=9.6u ps=9.6u mosmis_mod=1
*XMN1 out in  ss  sub ne3 w=5u l=800n m=1 nf=3 ad=0.9p as=0.9p pd=6.6u ps=6.6u mosmis_mod=1
MP1 out in  dd  well p1 w=5u l=800n m=3 ad=1.35p as=1.35p pd=9.6u ps=9.6u
MN1 out in  ss  sub n1 w=5u l=800n m=1 ad=0.9p as=0.9p pd=6.6u ps=6.6u
.ends inv1

.subckt inv5 dd ss sub well in out
xinv1 dd ss sub well in 1 inv1
xinv2 dd ss sub well 1  2 inv1
xinv3 dd ss sub well 2  3 inv1
xinv4 dd ss sub well 3  4 inv1
xinv5 dd ss sub well 4 out inv1
.ends inv5

xinv1 dd ss sub well in out5 inv5
xinv2 dd ss sub well out5 out10 inv5
xinv3 dd ss sub well out10 out15 inv5
xinv4 dd ss sub well out15 out20 inv5
xinv5 dd ss sub well out20 out inv5
xinv11 dd 0 sub well out buf inv1
cout  buf ss 0.2pF

 *** Model library files.
* Add your library here (full path required, or path relative to path
* of ngspice executable (interactive mode), or relative to path of
* input file (batch mode))
* Chose the transistors for XMP1 and XMN1 according to the library
*.lib "jc_usage.l" MC_LIB
*.lib "../../../various/lib-test/my_usage.l" MC_LIB
*.lib "D:\Spice_general\tests\lib-test\ts14\my_ts_usage.l" MC_LIB
*.lib "x_usage.l" MC_LIB

* or use the BSIM3 model with internal parameters except Vth0
* that varies the threshold voltage +-3 sigma around a mean of +-0.6V
.model p1 PMOS version=3.3.0 Level=8 Vth0=agauss(-0.6, 0.1, 3)
.model n1 NMOS version=3.3.0 Level=8 Vth0=agauss(0.6, 0.1, 3)

.control
  let mc_runs = 10               $ number of runs for monte carlo
  let run = 0                     $ number of actual run
  set curplot = new               $ create a new plot
  set curplottitle = "Transient outputs"
  set plot_out = $curplot         $ store its name to 'plot_out'
  set curplot = new               $ create a new plot
  set curplottitle = "FFT outputs"
  set plot_fft = $curplot         $ store its name to 'plot_fft'
  set curplot = new               $ create a new plot
  set curplottitle = "Oscillation frequency"
  set max_fft = $curplot          $ store its name to 'max_fft'
  let mc_runsp = mc_runs + 1
  let maxffts = unitvec(mc_runsp) $ vector for storing max measure results
  let halfffts = unitvec(mc_runsp)$ vector for storing measure results at -40dB rising
  unlet mc_runsp

  set mc_runs = $&mc_runs        $ create a variable from the vector
  let seeds = mc_runs + 2
  setseed $&seeds
  unlet seeds

  save buf                        $ we just need buf, save memory by more than 10x

* run the simulation loop

* We have to figure out what to do if a single simulation will not converge.
* There is now the variable sim_status, that is 0 if simulation ended regularly,
* and 1 if the simulation has been aborted with error message '...simulation(s) aborted'.
* Then we skip the rest of the run and continue with a new run.

  dowhile run <= mc_runs

    set run = $&run              $ create a variable from the vector

    * run=0 simulates with nominal parameters
    if run > 0
      echo
      echo * * * * * *
      echo Source the circuit again internally for run no. $run
      echo * * * * * *
      setseed $run
      mc_source  $ re-source the input file
    else
      echo run no. $run
    end
    echo simulation run no. $run of $mc_runs
    tran 100p 1000n 0
    echo Simulation status $sim_status
    let simstat = $sim_status
    if simstat = 1
      if run = mc_runs
        echo go to end
      else
        echo go to next run
      end
      destroy $curplot
      goto next
    end

* select stop and step so that number of data points after linearization is not too
* close to 8192, which would yield varying number of line length and thus scale for fft.
*
    set dt0 = $curplot
    * save the linearized data for having equal time scales for all runs
    linearize buf                 $ linearize only buf, no other vectors needed
    set dt1 = $curplot             $ store the current plot to dt (tran i+1)
    setplot $plot_out             $ make 'plt_out' the active plot
    * firstly save the time scale once to become the default scale
    if run=0
       let time={$dt1}.time
    end
    let vout{$run}={$dt1}.buf     $ store the output vector to plot 'plot_out'
    setplot $dt1                  $ go back to the previous plot (tran i+1)
    fft buf $ run fft on vector buf
    let buf2=db(mag(buf))
    * find the frequency where buf has its maximum of the fft signal
    meas sp fft_max MAX_AT buf2 from=0.05G to=0.7G
    * find the frequency where buf is -40dB at rising fft signal
    meas sp fft_40 WHEN buf2=-40 RISE=1 from=0.05G to=0.7G
    * store the fft vector
    set dt2 = $curplot            $ store the current plot to dt (spec i)
    setplot $plot_fft             $ make 'plot_fft' the active plot
    if run=0
       let frequency={$dt2}.frequency
    end
    let fft{$run}={$dt2}.buf      $ store the output vector to plot 'plot_fft'
    * store the measured value
    setplot $max_fft              $ make 'max_fft' the active plot
    let maxffts[{$run}]={$dt2}.fft_max
    let halfffts[{$run}]={$dt2}.fft_40
    destroy $dt0  $dt1  $dt2      $ save memory, we don't need this plot (spec) any more

    label next
    remcirc
    let run = run + 1
  end
***** plotting **********************************************************
if $?batchmode
  echo
  echo Plotting not available in batch mode
  echo Write linearized vout0 to vout{$mc_runs} to rawfile $rawfile
  echo
  write $rawfile {$plot_out}.allv
  rusage
  quit
else
  plot {$plot_out}.vout0          $ just plot the tran output with run 0 parameters
  setplot $plot_fft
  plot db(mag(ally)) xlimit 0 1G ylimit -80 10
*
* create a histogram from vector maxffts
  setplot $max_fft                $ make 'max_fft' the active plot
  set startfreq=50MEG
  set bin_size=1MEG
  set bin_count=100
  compose osc_frequ start=$startfreq step=$bin_size lin=$bin_count $ requires variables as parameters
  settype frequency osc_frequ
  let bin_count=$bin_count        $ create a vector from the variable
  let yvec=unitvec(bin_count)     $ requires vector as parameter
  let startfreq=$startfreq
  let bin_size=$bin_size
  * put data into the correct bins
  let run = 0
  dowhile run < mc_runs
    set run = $&run             $ create a variable from the vector
    let val = maxffts[{$run}]
    let part = 0
    * Check if val fits into a bin. If yes, raise bin by 1
    dowhile part < bin_count
      if ((val < (startfreq + (part+1)*bin_size)) & (val > (startfreq + part*bin_size)))
        let yvec[part] = yvec[part] + 1
                break
      end
      let part = part + 1
    end
    let run = run + 1
  end
  * plot the histogram
  set plotstyle=combplot
  let counts = yvec - 1             $ subtract 1 because we started with unitvec containing ones
  plot counts vs osc_frequ
* calculate jitter
  let diff40 = (vecmax(halfffts) - vecmin(halfffts))*1e-6
  echo
  echo Max. jitter is "$&diff40" MHz
end
  rusage
*  quit
.endc

.end
