Handling underflows

Question

I am trying to implement an algorithm described for example in this paper, Section III, or referenced here, from where numerical parameters have been taken. I added the snippet I use to compute the matrix from Eq. 3.1 in the first reference at the end, if it could be of any benefit. The need arises to compute very small numbers. Concretely, set

    Km = 0.025;
    Dm = 0.063;
    am = 4.2;
    kb = 8.617*^-5;
    V[delta_] := Dm*(1 - Exp[-delta])^2;

numbers such as

   Exp[-V[-40]]

must be calculated. I read these answers How to check for underflow, and find a constant to correct it? and here New General::munfl error and loss of precision but I am not getting anywhere, probably because I do not have yet any understanding of the underlying theory.

I tried

           Activate@SetPrecision[Inactivate[Exp[-V[-40]]], 15]

as well as

           Exp[SetPrecision[-V[-40], 15]]

or

            Quiet[Exp[-V[-40]], General::unfl]

but I still get underflows errors. The question whether these figures are to be computed at all is also interesting, but I am trying to replicate the results of the paper.

How should these overflows be handled?

As a layman, I would not in principle mind something simple as Julia does

 julia> dm = 0.063
 0.063

 julia> morse(x)=dm*(1-exp(-x))^2
 morse (generic function with 1 method)

 julia> morse(-40)
 3.4905921021679106e33

 julia> exp(-morse(-40))
 0.0

although I understand this could be insufficient for more complex applications. I assume underflow messages a warning about the achieved precision, and automatic truncation to 0 might be undesirable, but I do not know.

** Matrix to be diagonalised**

   lmin = -40;
   lmax = 5000;

   ndiv = 20;
   divis = (lmax - lmin)/ndiv;
   W[delta1_, delta2_] := (Km/(2*am^2))*(delta1 - delta2)^2;

   mater[tempe_] := 
   Table[(divis)*
   Exp[ -(1/(kb*tempe))*(0.5*V[lmin + divis*i] + 
    0.5*V[lmin + divis*j] + 
    W[lmin + divis*i, lmin + divis*j])] , {i, 0, ndiv}, {j, 0, 
    ndiv}];