Convergence issue with Verilog-A model

I have developed Verilog-A model of a device that I am trying to simulate and every time it gives different value and at most of the occasions it just stops at a certain time of the simulation.

I have seen multiple posts around the same issue and it seems that this is due to a bad Verilog-A model or varying a bit the convergence parameters.

I would like to get your suggestions of either updating the model or changing the convergence parameters in a matter to get this module working properly.

The model is as follows;

//Start
`resetall
`include "constants.vams"
`include "disciplines.vams"


// check this for later updates//////////////////////////////
nature Magnetization
    access = Mgn;
    units = "SI";
    abstol = 1;
endnature

nature Ange
    access = Ang;
    units = "rad";
    abstol = 1;
endnature

discipline magnetization
  potential Magnetization;
enddiscipline

discipline ange
  potential Ange;
enddiscipline
/////////////////////////////////////////////////////////////

//Constants
//Define various physical constants in SI units
`define Pi     3.141592653589793        //rad     
//`define m_Bohr 9.27401E-24              // J/T       Bohr Magneton    
//`define h      6.626070040E-34          // J.s       Plank's constant
`define hbar   1.054571800e-34          // J.s       Reduced Planck's constant
//`define eV     1.60217653E-19           // C         Electron charge
`define KB     1.38064852E-23             // J/K       Boltzman constant
`define mu0    1.256637061435917E-06    // N/A^2     Permeability of free space
//`define ep0    8.854187817620389E-12    // F/m       Permativity of free space
`define e      1.60217653E-19            // C         Electron charge
//`define me     9.1093826E-31              // Kg        Electron rest mass
//`define a0     5.29177E-11                // m         Bohr radius

module asdw(T1,T2,R_out1,R_out2);
/*-------T1~2: Actual terminals coressponding to pinned layer and free layer--------*/
        inout T1,T2;        
        electrical T1,T2;
        /*-------Tmz: Virtual terminal for monitoring the Magnetization-------*/
        output R_out1,R_out2;
        
        electrical R_out1,R_out2;
//        magnetization Tmx,Tmy,Tmz;
        ange Th1, Ph1;
        ange Th2, Ph2;
        
///////////////////////////////////////////////////////////////////////////////////////////
/* Simulation time step bound                                                     */
///////////////////////////////////////////////////////////////////////////////////////////

    parameter real Ts = 2p; // Speed-accuracy tradeoff

    parameter integer TN    = 0;
///////////////////////////////////////////////////////////////////////////////////////////
/* TN is a model parameter to enable Thermal fluctuations                                                         */
/*      TN = 1    : Thermal fluctuations is activated                             */
/*      TN = 0   : Thermal fluctuations is NOT activated                         */      
///////////////////////////////////////////////////////////////////////////////////////////

//    parameter integer STT   = 0;

///////////////////////////////////////////////////////////////////////////////////////////
/* JH is a model parameter                                                               */
/*      JH = 1   : Joule heating mechanism is activated                                  */
/*      JH = 0   : Joule heating mechanism is NOT activated                              */      
///////////////////////////////////////////////////////////////////////////////////////////

//    parameter integer JH    = 0;

///////////////////////////////////////////////////////////////////////////////////////////
/* SHAPE is a model parameter                                                            */
/*      SHAPE = 1   : MTJ cross section is circular                                      */
/*      SHAPE = 0   : MTJ cross section is rectangular                                   */      
///////////////////////////////////////////////////////////////////////////////////////////

//    parameter integer SHAPE    = 1;

///////////////////////////////////////////////////////////////////////////////////////////
/* Device parameters                                                             */
///////////////////////////////////////////////////////////////////////////////////////////
//    parameter real Ms0        = 8E5;          
    parameter real Ms0        = 1.2E6;          
//    parameter real alpha1     = 0.01;           
    parameter real alpha1     = 0.02;       
    parameter real gamma      = 2.21276E5;      



///////////////////////////////////////////////////////////////////////////////////////
/* Device size                                                                       */
///////////////////////////////////////////////////////////////////////////////////////    

//    parameter real t_FL       = 2e-9;               /
    parameter real t_FL       = 1.1E-9;               
//    parameter real w_FL       = 70e-9;                //
    parameter real w_FL1       = 170E-9;                //
//    parameter real l_FL       = 140e-9;               // (
    parameter real l_FL1      = 60E-9;               // (

    parameter real w_FL2       = 170E-9;                //
//    parameter real l_FL       = 140e-9;               //
    parameter real l_FL2       = 60E-9;               // (m)

//    parameter real t_ox       = 0.8e-9;               // (m)
    parameter real t_ox       = 1.1E-9;               // (m)

///////////////////////////////////////////////////////////////////////////////////////
/* Magnetoresistance and electrical parameters                                       */
///////////////////////////////////////////////////////////////////////////////////////

//    parameter real     b_h            =     3.28;          //ev       
    parameter real     b_h            =     0.39;          //ev       
//    parameter real     v_h            =     0.75;          //V       
    parameter real     v_h            =     0.75;          //V       

    parameter real     alpha_sd    =     1.4e-5;         // K^(-3/2)
    parameter real    P01        =     0.591;         //59.1%
    parameter real    P02        =    0.591;         //59.1%
    parameter real    S        =    1.5e-12;       //(ohm.(um)^2)^-1 . K^(-4/3)


    parameter real TMR0      = 1.5;            // TMR @ zero voltage

///////////////////////////////////////////////////////////////////////////////////////
/* Effective field parameters                                       */
///////////////////////////////////////////////////////////////////////////////////////    

//    parameter real Ki0      = 2.6E5 * t_FL * 0.1; //0.32e-3         //
    parameter real Ki0      = 0; //0.32e-3         //(J/m^2)  v=0

    parameter real eta_VCMA      = 0;                         // J/(V.m)
//External Field
    parameter real hx_external=    0; //-31830 A/m               external magnetic field along x axis
    parameter real hy_external=    0;                          // A/m               mexternal magnetic field along y axis
    parameter real hz_external=    0;                          // A/m               external magnetic field along z axis

//Demagnetization Fields, Shape anistropy field
//    parameter real Nx=         0.018;//0.0168            X-axis component of the demagnetization facator
//    parameter real Ny=         0.038;//0.0168;         Y-axis component of the demagnetization facator
//    parameter real Nz=         0.94;//0.9664;         Z-axis component of the demagnetization facator

    parameter real Nx=         0.0045;//0.0168            X-axis component of the demagnetization facator
    parameter real Ny=         0.0152;//0.0168;         Y-axis component of the demagnetization facator
    parameter real Nz=         0.9803;//0.9664;         Z-axis component of the demagnetization facator

//////////////////////////////////////////////////////////////////////////////////////////
/* Spin Transfer Torque Parameters                                                                 */
///////////////////////////////////////////////////////////////////////////////////////
    parameter real th_initial1= `Pi/2;                                // Azimuth initial angle (with Z direction)
//    parameter real ph_initial1= 0.001;                                // Elevation angle  (with X direction in X-Y plane)
    parameter real ph_initial1= `Pi-0.1;                                // Elevation angle  (with X direction in X-Y plane)
    parameter real th_initial2= `Pi/2;                                // Azimuth initial angle (with Z direction)
//    parameter real ph_initia2= 0.001;                                // Elevation angle  (with X direction in X-Y plane)
    parameter real ph_initial2= `Pi-0.1;                                // Elevation angle  (with X direction in X-Y plane)

//    parameter real pin_x= 0;                                // pinned layer orientation
//    parameter real pin_y=0;                                // pinned layer orientation
//    parameter real pin_z=1;                                // pinned layer orientation

    parameter real pin_x= 1;                                // pinned layer orientation
    parameter real pin_y=0;                                // pinned layer orientation
    parameter real pin_z=0;                                // pinned layer orientation

//    parameter real zeta= 1;//%0.58;                       // Spin polarized percent STT Coefficent
//    parameter real zeta_FL=    0;                        // Field-like STT torque coefficient

    parameter real zeta_FL= 1;//%0.58;                       // Spin polarized percent STT Coefficent
    parameter real zeta_DL=    0;                        // Field-like STT torque coefficient


    integer seed;                        
    parameter real Ta      = 300;       // Ambient temperature (K)
    
///////////////////////////////////////////////////////////////////////////////////////
// Model computing parameters:                                 //
///////////////////////////////////////////////////////////////////////////////////////
    real gamma_llg, surface1, surface2,volume1,volume2;
    real H_Demx1, H_Demy1, H_Demz1,H_Demx2, H_Demy2, H_Demz2;
    real Hkp1,Hkp2;
    real H_effx1, H_effy1, H_effz1,H_effx2, H_effy2, H_effz2;
    real Hlx1, Hly1, Hlz1,Hlx2, Hly2, Hlz2, Temp, Q1,Q2;
    real AJ1, BJ1,AJ2, BJ2;
    // real AJ_polarizer, BJ_polarizer;
    real st1, ct1, sp1, cp1,st2, ct2, sp2, cp2;
    real Heff_spherical_theta1, Heff_spherical_phi1, dMdt_eff_theta1, dMdt_eff_phi1;
    real Heff_spherical_theta2, Heff_spherical_phi2, dMdt_eff_theta2, dMdt_eff_phi2;
    real T_STT_spherical_pinned_theta1, T_STT_spherical_pinned_phi1, dMdt_STT_pinned_theta1, dMdt_STT_pinned_phi1;
    real T_STT_spherical_pinned_theta2, T_STT_spherical_pinned_phi2, dMdt_STT_pinned_theta2, dMdt_STT_pinned_phi2;
    
// real T_STT_spherical_polarizer_theta, T_STT_spherical_polarizer_phi, dMdt_STT_polarizer_theta, dMdt_STT_polarizer_phi;
    real dMdt_theta1, dMdt_phi1,dMdt_theta2, dMdt_phi2;
    real M_theta1, M_phi1, M_theta2, M_phi2;
    real I_STT,R1,R2;
    real TD,Time1,Time2;
    real C,sin_temp,G_T1,G_T2, G_V, G_SI, P1_T, P2_T, G1,G2;
    real t_ox_angestrom,C0,C1,C2,A_cm2_1,A_cm2_2,G01,G02,TMR_mod,F;
    
///////////////////////////////////////////////////////////////////////////////////////
/* Branches                                                                     */
///////////////////////////////////////////////////////////////////////////////////////
//    branch (T1,T2) res_mtj;
    
///////////////////////////////////////////////////////////////////////////////////////
/* Main                                                                              */
///////////////////////////////////////////////////////////////////////////////////////

analog begin
    
    $bound_step(Ts); // define simuation time step

    @(initial_step) begin
        M_theta1= th_initial1;
        M_phi1=ph_initial1;
        M_theta2=th_initial2;
        M_phi2=ph_initial2;

    
        seed = 25;
        //T1   = 0;
        F=1;
        t_ox_angestrom  =t_ox*1e10;
        surface1          = `Pi*l_FL1*w_FL1/4;                    // (m^2) Surface area
        surface2          = `Pi*l_FL2*w_FL2/4;                    // (m^2) Surface area
                
        volume1          = surface1*t_FL;                        // (m^3) Surface area
        volume2          = surface2*t_FL;                        // (m^3) Surface area
        gamma_llg      = gamma/(1+alpha1*alpha1);          // Hz/T (1/T.m)      Reduced gyromagnetic ratio
        A_cm2_1           = surface1*1e4;
        A_cm2_2           = surface2*1e4;

        G01 = ((3.16*1e10)*sqrt(b_h)/t_ox_angestrom) * exp(-1.025 * t_ox_angestrom * sqrt(b_h)) * A_cm2_1 * F; // 1/ohm
        G02 = ((3.16*1e10)*sqrt(b_h)/t_ox_angestrom) * exp(-1.025 * t_ox_angestrom * sqrt(b_h)) * A_cm2_2 * F; // 1/ohm
//        G0=2.77414e-4;                             //(APPLY kAZ PARAMETERS)

    end
        

    Time2         =     $abstime;
    TD             =     Time2 - Time1;            //Time duration
    Time1         =     Time2;
    
    Temp            =     Ta;


        if (TN == 1)
            begin
                Q1 = sqrt((2 * alpha1 * `KB * Temp / (`mu0*gamma_llg * Ms0 * volume1 * TD)));
                Q2 = sqrt((2 * alpha1 * `KB * Temp / (`mu0*gamma_llg * Ms0 * volume2 * TD)));
            end
        else
            begin
                Q1 = 0;
                Q2 = 0;
            end
    
        Hlx1  = $rdist_normal(seed,0,Q1);        //$dist_normal (seed, mean, standard_deviation) ; seed const means get const random value
        Hly1  = $rdist_normal(seed,0,Q1);
        Hlz1  = $rdist_normal(seed,0,Q1);
        Hlx2  = $rdist_normal(seed,0,Q2);        //$dist_normal (seed, mean, standard_deviation) ; seed const means get const random value
        Hly2  = $rdist_normal(seed,0,Q2);
        Hlz2  = $rdist_normal(seed,0,Q2);

        st1   = sin(Ang(Th1));
        ct1   = cos(Ang(Th1));
        sp1   = sin(Ang(Ph1));
        cp1   = cos(Ang(Ph1));    

        st2   = sin(Ang(Th2));
        ct2   = cos(Ang(Th2));
        sp2   = sin(Ang(Ph2));
        cp2   = cos(Ang(Ph2));    
    
// as H_dem = - Ms (N(vector).M(vector) (Nx,Ny,Nz). (M_theta,Mphi) cross product
        H_Demx1=         Nx*Ms0 * st1 * cp1 ;
        H_Demy1=         Ny*Ms0 * st1 * sp1;
        H_Demz1=         Nz*Ms0 * ct1;
        H_Demx2=         Nx*Ms0 * st1 * cp1 ;
        H_Demy2=         Ny*Ms0 * st1 * sp1;
        H_Demz2=         Nz*Ms0 * ct1;
//Anisotropy Field, Interfacial
        Hkp1=            (2*Ki0/(t_FL*`mu0*Ms0)-2*eta_VCMA*V(T1,T2)/(t_FL*`mu0*Ms0*t_ox))*ct1;     // The effective anistropy with VCMA effect
        Hkp2=            (2*Ki0/(t_FL*`mu0*Ms0)-2*eta_VCMA*V(T1,T2)/(t_FL*`mu0*Ms0*t_ox))*ct2;     // The effective anistropy with VCMA effect
//Effective anistropy field with VCMA
        H_effx1=       Hlx1 + hx_external- H_Demx1;
        H_effy1=       Hly1 + hy_external- H_Demy1;
        H_effz1=       Hlz1 + hz_external- H_Demz1 + Hkp1;

        H_effx2=       Hlx2 + hx_external- H_Demx2;
        H_effy2=       Hly2 + hy_external- H_Demy2;
        H_effz2=       Hlz2 + hz_external- H_Demz2 + Hkp2;

// write disturbance
//Spin Transfer Torque components
             
        AJ1=                   (`hbar*zeta_FL*I(T1,T2)/(2*`e*`mu0*Ms0*surface1*t_FL));  //0; %0.3e4;                       //A/m STT coefficient Normalized
        BJ1=                   (`hbar*zeta_DL*I(T1,T2)/(2*`e*`mu0*Ms0*surface1*t_FL));                            //  Normalized     
        AJ2=                   (`hbar*zeta_FL*I(T1,T2)/(2*`e*`mu0*Ms0*surface2*t_FL));  //0; %0.3e4;                       //A/m STT coefficient Normalized
        BJ2=                   (`hbar*zeta_DL*I(T1,T2)/(2*`e*`mu0*Ms0*surface2*t_FL));                            //  Normalized    

    //

    T_STT_spherical_pinned_theta1=        pin_x * ct1 * cp1 + pin_y * ct1 * sp1 - pin_z * st1;
    T_STT_spherical_pinned_phi1=      -pin_x * sp1 + pin_y * cp1;
    T_STT_spherical_pinned_theta2=        pin_x * ct2 * cp2 + pin_y * ct2 * sp2 - pin_z * st2;
    T_STT_spherical_pinned_phi2=      -pin_x * sp2 + pin_y * cp2;


// Heff_spherical =   @(t,M) [cos(M(1)).* cos(M(2)), cos(M(1)).*sin(M(2)),    -sin(M(1)); -sin(M(2)), cos(M(2)), 0] * H(t,M);    %A/m Magnetic Field Strength normalized
    Heff_spherical_theta1=               H_effx1 * ct1 * cp1 + H_effy1 * ct1 * sp1 - H_effz1 * st1;
    Heff_spherical_phi1=                -H_effx1 * sp1 + H_effy1 * cp1;
    Heff_spherical_theta2=               H_effx2 * ct2 * cp2 + H_effy2 * ct2 * sp2 - H_effz2 * st2;
    Heff_spherical_phi2=                -H_effx2 * sp2 + H_effy2 * cp2;

    dMdt_eff_theta1=             alpha1 * Heff_spherical_theta1 + Heff_spherical_phi1;
    dMdt_eff_phi1=                -Heff_spherical_theta1/st1 + alpha1 * Heff_spherical_phi1 / st1;
    dMdt_eff_theta2=             alpha1 * Heff_spherical_theta2 + Heff_spherical_phi2;
    dMdt_eff_phi2=                -Heff_spherical_theta2/st2 + alpha1 * Heff_spherical_phi2 / st2;    

    dMdt_STT_pinned_theta1=          (-AJ1 - alpha1 * BJ1) * T_STT_spherical_pinned_theta1 + (alpha1 * AJ1 - BJ1) * T_STT_spherical_pinned_phi1;
    dMdt_STT_pinned_phi1=         (BJ1 - alpha1 * AJ1) * T_STT_spherical_pinned_theta1/st1 + (- AJ1 - alpha1 * BJ1) * T_STT_spherical_pinned_phi1 / st1;
    dMdt_STT_pinned_theta2=          (-AJ2 - alpha1 * BJ2) * T_STT_spherical_pinned_theta2 + (alpha1 * AJ2 - BJ2) * T_STT_spherical_pinned_phi2;
    dMdt_STT_pinned_phi2=         (BJ2 - alpha1 * AJ2) * T_STT_spherical_pinned_theta2/st2 + (- AJ2 - alpha1 * BJ2) * T_STT_spherical_pinned_phi2 / st2;

    dMdt_theta1=         gamma_llg * (dMdt_eff_theta1 + dMdt_STT_pinned_theta1 );
    dMdt_phi1=         gamma_llg * (dMdt_eff_phi1 + dMdt_STT_pinned_phi1 );
    dMdt_theta2=         gamma_llg * (dMdt_eff_theta2 + dMdt_STT_pinned_theta2 );
    dMdt_phi2=         gamma_llg * (dMdt_eff_phi2 + dMdt_STT_pinned_phi2 );    
    
    
//idtmod(integrand x, initial condition y0, modulus m, offset b, tolerance) preferable than just idt
    M_theta1 =  idtmod(dMdt_theta1,th_initial1,2*`Pi);    
    M_phi1   =  idtmod(dMdt_phi1,ph_initial1,2*`Pi);
    M_theta2 =  idtmod(dMdt_theta2,th_initial2,2*`Pi);    
    M_phi2   =  idtmod(dMdt_phi2,ph_initial2,2*`Pi);


// MTJ resistance

    C          =     1.387e-4 * t_ox_angestrom/sqrt(b_h); //Check the parameters later
    sin_temp =      sin(C*Temp);
    G_T1      =     G01* (C*Temp)/(sin_temp);
    G_T2      =     G02* (C*Temp)/(sin_temp);
    G_SI        =     S*pow(Temp,(4/3));            //
//Method1    
    P1_T    =     P01*(1-alpha_sd*pow(Temp,1.5));
    P2_T    =     P02*(1-alpha_sd*pow(Temp,1.5));

    G1    =     G_T1*(1+P1_T*P2_T*cp1 + G_SI);    //
    G2    =     G_T2*(1+P1_T*P2_T*cp2 + G_SI);    //
    R1     = 1 / (G1);
    R2     = 1 / (G2);  



//    I(T1,T2) <+ V(T1,T2)/R;
     V(T1,T2) <+ I(T1,T2) * (R1+R2);
    V(R_out1) <+  R1;
    V(R_out2) <+  R2;
    Ang(Th1) <+ M_theta1;
    Ang(Ph1) <+ M_phi1;
    Ang(Th2) <+ M_theta2;
    Ang(Ph2) <+ M_phi2;

    //

    
    end // end of analog

endmodule   // end of module

the input.scs

simulator lang=spectre
global 0
parameters I=600u

// View name: schematic
I1 (0 Vout) isource dc=I type=dc
I5 (Vout 0 Rout1 Rout2) asdw Ts=2e-12 TN=0 Ms0=1.2e+06 alpha1=0.02 \
        gamma=221276 t_FL=1.1e-09 w_FL1=1.7e-07 l_FL1=6e-08 w_FL2=1.5e-07 \
        l_FL2=5e-08 t_ox=1.1e-09 b_h=0.39 v_h=0.75 alpha_sd=1.4e-05 \
        P01=0.591 P02=0.591 S=1.5e-12 TMR0=1.5 Ki0=0 eta_VCMA=0 \
        hx_external=0 hy_external=0 hz_external=0 Nx=0.0045 Ny=0.0152 \
        Nz=0.9803 th_initial1=1.5708 ph_initial1=0.1 th_initial2=1.5708 \
        ph_initial2=0.1 pin_x=1 pin_y=0 pin_z=0 zeta_FL=1 zeta_DL=0 Ta=300
simulatorOptions options reltol=1e-3 vabstol=1e-6 iabstol=1e-12 temp=27 \
    tnom=27 scalem=1.0 scale=1.0 gmin=1e-12 rforce=1 maxnotes=5 maxwarns=5 \
    digits=5 cols=80 pivrel=1e-3 sensfile="../psf/sens.output" \
    checklimitdest=psf
tran tran stop=800n param=I param_vec=[0 0 30n 450u 50n 0 70n 600u 100n 0 \
    150n -620u 170n 0 200n -630u 250n 0 280n 450u 300n 0 330n -630u 360n \
    0] step=100000 write="spectre.ic" writefinal="spectre.fc" \
    annotate=status maxiters=5
finalTimeOP info what=oppoint where=rawfile
modelParameter info what=models where=rawfile
element info what=inst where=rawfile
outputParameter info what=output where=rawfile
designParamVals info what=parameters where=rawfile
primitives info what=primitives where=rawfile
subckts info what=subckts where=rawfile
save I5:T1
saveOptions options save=allpub
ahdl_include "/Directory/asdw/veriloga/veriloga.va"

The circuit