Good day!

May I please kindly ask how can I derive the inverse Mills ratio from a heteroscedastic probit model (estimated using the command hetprobit)? I will be using it for the second stage estimation of a wage model.
Code:
hetprobit h800job_2 g702age g702age_2 i.b202sex_2 i.partner i.h812recv_2 child_bmov i.relgn_christ, het(g702age i.b202sex_2 i.partner) vce(robust)
Shown below is the results of the command:
HTML Code:
Heteroskedastic probit model                    Number of obs     =      3,167
                                                Zero outcomes     =      1,489
                                                Nonzero outcomes  =      1,678

                                                Wald chi2(7)      =      36.35
Log pseudolikelihood = -1756.303                Prob > chi2       =     0.0000

--------------------------------------------------------------------------------
               |               Robust
     h800job_2 |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
---------------+----------------------------------------------------------------
h800job_2      |
       g702age |    .334209   .0642149     5.20   0.000     .2083502    .4600679
     g702age_2 |  -.0041946    .000836    -5.02   0.000    -.0058331   -.0025562
   1.b202sex_2 |  -1.954624   .3479778    -5.62   0.000    -2.636648     -1.2726
     1.partner |  -.0644116   .1065569    -0.60   0.546    -.2732594    .1444361
  1.h812recv_2 |  -.4597904   .1403239    -3.28   0.001    -.7348202   -.1847605
    child_bmov |   .0350233   .0520128     0.67   0.501    -.0669199    .1369665
1.relgn_christ |   .2872648   .1473732     1.95   0.051    -.0015813    .5761108
         _cons |  -4.805648   .9846948    -4.88   0.000    -6.735615   -2.875682
---------------+----------------------------------------------------------------
lnsigma2       |
       g702age |   .0170492   .0040073     4.25   0.000      .009195    .0249034
   1.b202sex_2 |   .6008698   .1568638     3.83   0.000     .2934223    .9083172
     1.partner |  -.3166569   .0864646    -3.66   0.000    -.4861245   -.1471893
--------------------------------------------------------------------------------
Wald test of lnsigma2=0: chi2(3) = 42.92                  Prob > chi2 = 0.0000
Initially, I used the following commands but later realized that I might be treating it as if it's a homoscedastic probit model:
Code:
predict phat, xb
gen imr = normalden(phat)/normal(phat)
Thank you very much!