Well, thank you for your insightful comment!
I think you are right, for precise computation of log-softmax, a much lower breakpoint is sufficient.
What I did not mention (my bad) is that I want to apply gradient backprop through this computation.
This changes things a lot: If setting the breakpoint too low, the error in the gradient can be 7 orders of magnitude bigger than the error in the log-softplus output.
The following code is an extension of yours.
import torch
print (torch.__version__)

def log_softplusA (x):
    return torch.nn.functional.softplus (x).log()

def log_softplusB (x, brk = -15.0):
    return torch.where (x > brk, torch.nn.functional.softplus (x).log(), x)

x = torch.arange (-200, 100, 1).double() / 10.0
x.requires_grad = True
lsp = log_softplusA (x)   # double-precision "ground truth"
lsp.sum().backward()
lspg = x.grad

xA = torch.arange (-200, 100, 1) / 10.0
xA.requires_grad = True
lspA = log_softplusA (xA)
diffA = lspA - lsp
lspA.sum().backward()
lspgA = xA.grad
diffgA = lspgA - lspg

xB = torch.arange (-200, 100, 1) / 10.0
xB.requires_grad = True
lspB = log_softplusB (xB)
diffB = lspB - lsp
lspB.sum().backward()

lspgB = xB.grad
diffgB = lspgB - lspg


print ('diffA.abs().max() =', diffA.abs().max())
print ('diffB.abs().max() =', diffB.abs().max())
print ('diffB.abs().argmax() =', diffB.abs().argmax())
print ('x[diffB.abs().argmax()] =', x[diffB.abs().argmax()])

print ('diffgA.abs().max() =', diffgA.abs().max())
print ('diffgB.abs().max() =', diffgB.abs().max())
print ('diffgB.abs().argmax() =', diffgB.abs().argmax())
print ('x[diffgB.abs().argmax()] =', x[diffgB.abs().argmax()])
print ('lspgB[diffgB.abs().argmax()], lspg[diffgB.abs().argmax()] =', lspgB[diffgB.abs().argmax()], lspg[diffgB.abs().argmax()])

yields
1.8.1
diffA.abs().max() = tensor(8.4096e-07, dtype=torch.float64, grad_fn=<MaxBackward1>)
diffB.abs().max() = tensor(8.4096e-07, dtype=torch.float64, grad_fn=<MaxBackward1>)
diffB.abs().argmax() = tensor(61)
x[diffB.abs().argmax()] = tensor(-13.9000, dtype=torch.float64, grad_fn=<SelectBackward>)
diffgA.abs().max() = tensor(1.0000, dtype=torch.float64)
diffgB.abs().max() = tensor(0.1339, dtype=torch.float64)
diffgB.abs().argmax() = tensor(53)
x[diffgB.abs().argmax()] = tensor(-14.7000, dtype=torch.float64, grad_fn=<SelectBackward>)
lspgB[diffgB.abs().argmax()], lspg[diffgB.abs().argmax()] = tensor(0.8661) tensor(1.0000, dtype=torch.float64)

Upon variating the breakpoint and re-running the code, I empirically found the value -8.0 to perform best.