def S(r, rc1, rc2, derivative=False):
    """
    Calculate the switching function S(r) which decays continuously
    between 1 and 0 in the range from rc1 to rc2 (rc2>rc1):
       S(r) = (rc2^2 - r^2)^2 * (rc2^2 + 2*r^2 - 3*rc1^2) / (rc2^2-rc1^2)^3
    I'm using the same smoothing/switching cutoff function used by the CHARMM
    force-fields.  (I'm even using the same code to implement it, taken
    from lammps charmm/coul/charmm pair style, rewritten in python.)

    """
    assert(rc2>rc1)
    rsq           = r*r
    rc1sq         = rc1*rc1
    rc2sq         = rc2*rc2
    denom_lj_inv  = (1.0 / ((rc2sq-rc1sq)*
                            (rc2sq-rc1sq)*
                            (rc2sq-rc1sq)))
    if rsq > rc2sq:
        return 0.0
    elif rsq < rc1sq:
        if derivative:
            return 0.0
        else:
            return 1.0
    else:
        rc2sq_minus_rsq    = (rc2sq - rsq)
        rc2sq_minus_rsq_sq = rc2sq_minus_rsq * rc2sq_minus_rsq
        if derivative:
            return (12.0 * rsq * rc2sq_minus_rsq * (rsq-rc1sq) * denom_lj_inv)
        else:
            return (rc2sq_minus_rsq_sq *
                    (rc2sq + 2.0*rsq - 3.0*rc1sq) * denom_lj_inv)