code: def folding_model_gradient(rvec, rcut):
      r
      rnorm = np.linalg.norm(rvec)
      if rnorm == 0.0:
          return np.zeros(rvec.shape)
      r = rnorm - rcut
      if r < 0.0:
          return -5.0 * r * rvec / rnorm
      return (1.5 * r - 2.0) * rvec / rnorm
docs: r"""computes the potential's gradient at point rvec