Spring damping in Energy Function Supposing a behavior function $C(x_1, \ldots, x_n)$, then we have a scalar potential energy function  

$E = \frac{k_s}{2}C\cdot C$ where $k_S$ is stiffness constant.

Hence, the force is as follows $f_i = -\frac{\partial E}{\partial x_i} = -k_sC\frac{\partial C}{\partial x_i}$ 

I caught up so far, but couldn't understand when adding damping.

The paper read $f_i = (-k_sC - k_d\dot{C})\frac{\partial C}{\partial x_i}$ , where $k_d$ is damping constant and $\dot{C}$ is derivative of $C$.

How to derive the damping part?

I know damping equation is $f_d = -k_dv$, 

In detail, supposing there are two end of points in a spirng, called $p_a, p_b$.

then the damping force $f_d = -k_d\frac{\dot{l}\cdot l}{|l|}\frac{l}{|l|}$ in 3D sapce ($\dot{l}$ is $v_a - v_b$, $l = p_a - p_b$)

However, I could not correspond this equation to the above damping part.
The damping force you mentioned $f=-k \frac{\dot{l} \cdot l}{|l|} \frac{l}{|l|}$ is a special case of $f=-k \dot{C} \frac{\partial C}{\partial \mathbf{x}}$.

Let
$$
\begin{align}
C(\mathbf{x}) &amp;= \lVert\mathbf{l}\rVert-l_0 \\
\mathbf{l} &amp;= \mathbf{x}_1 - \mathbf{x}_2 \\
\mathbf{\dot{l}} &amp;= \mathbf{v}_1 - \mathbf{v}_2 \\
\mathbf{v} &amp;= \dot{\mathbf{x}} = \begin{pmatrix}\mathbf{v}_1 \\ \mathbf{v}_2\end{pmatrix}
\end{align}
$$
where $l_0$ is the rest length of the spring, and $\mathbf{x}_1$, $\mathbf{x}_2$ are the positions of two end points.
The damping force is
$$
\begin{align}
\mathbf{f} &amp;= -k \dot{C} \frac{\partial C}{\partial \mathbf{x}} \\
&amp;= -k (\frac{\partial C}{\partial \mathbf{x}})^T \dot{\mathbf{x}} \frac{\partial C}{\partial \mathbf{x}}\\
&amp;= -k (\mathbf{l}^T\mathbf{l})^{-1/2} \mathbf{l}^T \begin{pmatrix}\mathbf{I}&amp; -\mathbf{I} \end{pmatrix} \dot{\mathbf{x}} (\mathbf{l}^T\mathbf{l})^{-1/2} \begin{pmatrix}\mathbf{I}&amp; -\mathbf{I} \end{pmatrix}^T \mathbf{l} \\
&amp;= -k (\mathbf{l}^T\mathbf{l})^{-1} \begin{pmatrix}\mathbf{l}^T&amp; -\mathbf{l}^T \end{pmatrix} \mathbf{v} \begin{pmatrix}\mathbf{l} \\ -\mathbf{l} \end{pmatrix}\\
&amp;= -k \frac{(\mathbf{l}\cdot \mathbf{\dot{l}}) }{\lVert \mathbf{l} \rVert^2}
\begin{pmatrix}\mathbf{l} \\ -\mathbf{l} \end{pmatrix}
\end{align}
$$

This matches your damping force equation. The $-\mathbf{l}$ part is just the force for the other particle, which has the opposite direction as expected. Bold letters are used to denote vector quantities.

As mentioned in [baraff98], this damping force is along the direction $\frac{\partial C}{\partial \mathbf{x}}$,
and its magnitude is the projected component of velocity $\frac{\partial C}{\partial \mathbf{x}}^T \mathbf{v}=\frac{\partial C}{\partial \mathbf{x}}^T \mathbf{\dot{x}}=\dot{C}$.