def force(box, r_cut, r):
    'Takes in box, cutoff range, and coordinate array, and calculates forces and potentials etc.'
    import numpy as np
    (n, d) = r.shape
    assert (d == 3), 'Dimension error in force'
    sr2_ovr = 1.77
    r_cut_box = (r_cut / box)
    r_cut_box_sq = (r_cut_box ** 2)
    box_sq = (box ** 2)
    sr2 = (1.0 / (r_cut ** 2))
    sr6 = (sr2 ** 3)
    sr12 = (sr6 ** 2)
    pot_cut = (sr12 - sr6)
    f = np.zeros_like(r)
    total = PotentialType(cut=0.0, pot=0.0, vir=0.0, lap=0.0, ovr=False)
    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)
            sr2 = np.where(in_range, (1.0 / rij_sq), 0.0)
            ovr = (sr2 > sr2_ovr)
            sr6 = (sr2 ** 3)
            sr12 = (sr6 ** 2)
            cut = (sr12 - sr6)
            vir = (cut + sr12)
            pot = np.where(in_range, (cut - pot_cut), 0.0)
            lap = (((22.0 * sr12) - (5.0 * sr6)) * sr2)
            fij = (vir * sr2)
            fij = (rij * fij[:, np.newaxis])
            total = (total + PotentialType(cut=np.sum(cut), pot=np.sum(pot), vir=np.sum(vir), lap=np.sum(lap), ovr=np.any(ovr)))
            f[i, :] = (f[i, :] + np.sum(fij, axis=0))
            f[(i + 1):, :] = (f[(i + 1):, :] - fij)
    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)
                    sr2 = (1.0 / rij_sq)
                    ovr = (sr2 > sr2_ovr)
                    sr6 = (sr2 ** 3)
                    sr12 = (sr6 ** 2)
                    cut = (sr12 - sr6)
                    vir = (cut + sr12)
                    pot = (cut - pot_cut)
                    lap = (((22.0 * sr12) - (5.0 * sr6)) * sr2)
                    fij = ((rij * vir) * sr2)
                    total = (total + PotentialType(cut=cut, pot=pot, vir=vir, lap=lap, ovr=ovr))
                    f[i, :] = (f[i, :] + fij)
                    f[j, :] = (f[j, :] - fij)
    f = (f * 24.0)
    total.cut = (total.cut * 4.0)
    total.pot = (total.pot * 4.0)
    total.vir = ((total.vir * 24.0) / 3.0)
    total.lap = ((total.lap * 24.0) * 2.0)
    return (total, f)
Takes in box, cutoff range, and coordinate array, and calculates forces and potentials etc.
def force(box, r_cut, r):
    
    import numpy as np
    (n, d) = r.shape
    assert (d == 3), 'Dimension error in force'
    sr2_ovr = 1.77
    r_cut_box = (r_cut / box)
    r_cut_box_sq = (r_cut_box ** 2)
    box_sq = (box ** 2)
    sr2 = (1.0 / (r_cut ** 2))
    sr6 = (sr2 ** 3)
    sr12 = (sr6 ** 2)
    pot_cut = (sr12 - sr6)
    f = np.zeros_like(r)
    total = PotentialType(cut=0.0, pot=0.0, vir=0.0, lap=0.0, ovr=False)
    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)
            sr2 = np.where(in_range, (1.0 / rij_sq), 0.0)
            ovr = (sr2 > sr2_ovr)
            sr6 = (sr2 ** 3)
            sr12 = (sr6 ** 2)
            cut = (sr12 - sr6)
            vir = (cut + sr12)
            pot = np.where(in_range, (cut - pot_cut), 0.0)
            lap = (((22.0 * sr12) - (5.0 * sr6)) * sr2)
            fij = (vir * sr2)
            fij = (rij * fij[:, np.newaxis])
            total = (total + PotentialType(cut=np.sum(cut), pot=np.sum(pot), vir=np.sum(vir), lap=np.sum(lap), ovr=np.any(ovr)))
            f[i, :] = (f[i, :] + np.sum(fij, axis=0))
            f[(i + 1):, :] = (f[(i + 1):, :] - fij)
    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)
                    sr2 = (1.0 / rij_sq)
                    ovr = (sr2 > sr2_ovr)
                    sr6 = (sr2 ** 3)
                    sr12 = (sr6 ** 2)
                    cut = (sr12 - sr6)
                    vir = (cut + sr12)
                    pot = (cut - pot_cut)
                    lap = (((22.0 * sr12) - (5.0 * sr6)) * sr2)
                    fij = ((rij * vir) * sr2)
                    total = (total + PotentialType(cut=cut, pot=pot, vir=vir, lap=lap, ovr=ovr))
                    f[i, :] = (f[i, :] + fij)
                    f[j, :] = (f[j, :] - fij)
    f = (f * 24.0)
    total.cut = (total.cut * 4.0)
    total.pot = (total.pot * 4.0)
    total.vir = ((total.vir * 24.0) / 3.0)
    total.lap = ((total.lap * 24.0) * 2.0)
    return (total, f)