def hessian(box, r_cut, r, f):
    'Calculates Hessian function (for 1/N correction to config temp).'
    import numpy as np
    (n, d) = r.shape
    assert (d == 3), 'Dimension error in hessian'
    assert np.all((r.shape == f.shape)), 'Dimension mismatch in hessian'
    r_cut_box = (r_cut / box)
    r_cut_box_sq = (r_cut_box ** 2)
    box_sq = (box ** 2)
    hes = 0.0
    if fast:
        for i in range((n - 1)):
            rij = (r[i, :] - r[(i + 1):, :])
            rij = (rij - np.rint(rij))
            rij_sq = np.sum((rij ** 2), axis=1)
            in_range = (rij_sq < r_cut_box_sq)
            rij_sq = (rij_sq * box_sq)
            rij = (rij * box)
            fij = (f[i, :] - f[(i + 1):, :])
            ff = np.sum((fij * fij), axis=1)
            rf = np.sum((rij * fij), axis=1)
            sr2 = np.where(in_range, (1.0 / rij_sq), 0.0)
            sr6 = (sr2 ** 3)
            sr8 = (sr6 * sr2)
            sr10 = (sr8 * sr2)
            v1 = ((24.0 * (1.0 - (2.0 * sr6))) * sr8)
            v2 = ((96.0 * ((7.0 * sr6) - 2.0)) * sr10)
            hes = ((hes + np.sum((v1 * ff))) + np.sum((v2 * (rf ** 2))))
    else:
        for i in range((n - 1)):
            for j in range((i + 1), n):
                rij = (r[i, :] - r[j, :])
                rij = (rij - np.rint(rij))
                rij_sq = np.sum((rij ** 2))
                if (rij_sq < r_cut_box_sq):
                    rij_sq = (rij_sq * box_sq)
                    rij = (rij * box)
                    fij = (f[i, :] - f[j, :])
                    ff = np.dot(fij, fij)
                    rf = np.dot(rij, fij)
                    sr2 = (1.0 / rij_sq)
                    sr6 = (sr2 ** 3)
                    sr8 = (sr6 * sr2)
                    sr10 = (sr8 * sr2)
                    v1 = ((24.0 * (1.0 - (2.0 * sr6))) * sr8)
                    v2 = ((96.0 * ((7.0 * sr6) - 2.0)) * sr10)
                    hes = ((hes + (v1 * ff)) + (v2 * (rf ** 2)))
    return hes
Calculates Hessian function (for 1/N correction to config temp).
def hessian(box, r_cut, r, f):
    
    import numpy as np
    (n, d) = r.shape
    assert (d == 3), 'Dimension error in hessian'
    assert np.all((r.shape == f.shape)), 'Dimension mismatch in hessian'
    r_cut_box = (r_cut / box)
    r_cut_box_sq = (r_cut_box ** 2)
    box_sq = (box ** 2)
    hes = 0.0
    if fast:
        for i in range((n - 1)):
            rij = (r[i, :] - r[(i + 1):, :])
            rij = (rij - np.rint(rij))
            rij_sq = np.sum((rij ** 2), axis=1)
            in_range = (rij_sq < r_cut_box_sq)
            rij_sq = (rij_sq * box_sq)
            rij = (rij * box)
            fij = (f[i, :] - f[(i + 1):, :])
            ff = np.sum((fij * fij), axis=1)
            rf = np.sum((rij * fij), axis=1)
            sr2 = np.where(in_range, (1.0 / rij_sq), 0.0)
            sr6 = (sr2 ** 3)
            sr8 = (sr6 * sr2)
            sr10 = (sr8 * sr2)
            v1 = ((24.0 * (1.0 - (2.0 * sr6))) * sr8)
            v2 = ((96.0 * ((7.0 * sr6) - 2.0)) * sr10)
            hes = ((hes + np.sum((v1 * ff))) + np.sum((v2 * (rf ** 2))))
    else:
        for i in range((n - 1)):
            for j in range((i + 1), n):
                rij = (r[i, :] - r[j, :])
                rij = (rij - np.rint(rij))
                rij_sq = np.sum((rij ** 2))
                if (rij_sq < r_cut_box_sq):
                    rij_sq = (rij_sq * box_sq)
                    rij = (rij * box)
                    fij = (f[i, :] - f[j, :])
                    ff = np.dot(fij, fij)
                    rf = np.dot(rij, fij)
                    sr2 = (1.0 / rij_sq)
                    sr6 = (sr2 ** 3)
                    sr8 = (sr6 * sr2)
                    sr10 = (sr8 * sr2)
                    v1 = ((24.0 * (1.0 - (2.0 * sr6))) * sr8)
                    v2 = ((96.0 * ((7.0 * sr6) - 2.0)) * sr10)
                    hes = ((hes + (v1 * ff)) + (v2 * (rf ** 2)))
    return hes