The current implementation uses the direct definition of fidelity, which works well for small dimensions. From my experience, I was getting values significantly greater than 1. For the fidelity between rho1 and rho2, I found that squaring the sum of the singular values of the product of square_root(rho1) and square_root(rho2) gives a more accurate calculation. This equation is equal to fidelity. The idea is not mine and came from qutip.

This is my version of the implementation that I use.

import scipy
def get_fidelity(rho1, rho2):
    rho1_sqrt=scipy.linalg.sqrtm(rho1)
    rho2_sqrt=scipy.linalg.sqrtm(rho2)
    return scipy.linalg.svdvals(rho1_sqrt @ rho2_sqrt).sum()**2