How to solve $\binom{n}{1}^2+2\binom{n}{2}^2 + 3\binom{n}{3}^2 + 4\binom{n}{4}^2+\cdots + n\binom{n}{n}^2$?
First recall that the coefficient of $x^n$ in $(1+x)^n(1+x)^n=(1+x)^{2n}$ implies $$ \begin{align} \sum_{k=0}^n\binom{n}{k}^2 &=\sum_{k=0}^n\binom{n}{k}\binom{n}{n-k}\\ &=\binom{2n}{n}\tag{1} \end{align} $$ and then note that $$ \begin{align} \sum_{k=0}^nk\binom{n}{k}^2 &=\sum_{k=0}^nk\binom{n}{n-k}^2\\ &=\sum_{k=0}^n(n-k)\binom{n}{k}^2\tag{2} \end{align} $$ Adding the first and last parts of $(2)$ yields $$ \begin{align} 2\sum_{k=0}^nk\binom{n}{k}^2 &=n\sum_{k=0}^n\binom{n}{k}^2\\ &=n\binom{2n}{n}\tag{3} \end{align} $$ Therefore, $$ \sum_{k=0}^nk\binom{n}{k}^2=\frac{n}{2}\binom{2n}{n}\tag{4} $$
This kind of looks like you want to appeal to Vandermonde's convolution, or at least the method you'd use to prove it. It can be applied directly as follows:
Let $S = \sum\limits_{k=0}^n k\binom{n}{k}^2$ be the sum we want to compute. Note that $S = \sum\limits_{k=0}^n (n-k) \binom{n}{n-k}^2 = \sum\limits_{k=0}^n (n-k) \binom{n}{k}^2$. Therefore $2S = n\sum\limits_{k=0}^n \binom{n}{k}^2 = n \binom{2n}{n}$. Then
$$S = \frac{n}{2}\binom{2n}{n}.$$
Building on the excellent work of @thomas-belulovich & @robjohn, but addressing your subsequently revealed goal....
You want $$ T(n) =S(n)-m\left\lfloor\frac{S(n)}m\right\rfloor =S(n)\text{ mod }m \qquad\text{for}\qquad S(n) =\frac{n}{2}{2n\choose n}. $$ Note that $$ \frac{S(n)}{S(n-1)} =\frac{n}{n-1}\cdot\frac{2n(2n-1)}{n^2} =2\frac{2n-1}{n-1}. $$ A brute force approach to this might be to calculate $$ R(k)=2(2k-1)(k-1)^{-1}\pmod m $$ for each $1<k\le n$ and then multiply modulo $m$ termwise to obtain $$ T(n)=S(1)\prod_{k=2}^nR(k). $$ Calculating each $R(k)=(4k-2)x$ requires $O(\log k)$ divisions using the extended Euclidean algorithm to find an $x$ so that $$ (k-1)x+my=1. $$
Extended Euclidean Algorithm: Given nonzero $a,b\in\mathbb{Z}$, find $x,y\in\mathbb{Z}$ so that $ax+by=\text{gcd}(a,b)$. (Adapted from David Joyner's book.)
Set $\overline{u}=\langle{u_0,u_1,u_2}\rangle$ $\leftarrow\langle{a,1,0}\rangle$ and $\overline{v}\leftarrow\langle{b,0,1}\rangle$ (vectors in $\mathbb{Z}^3$). Set $\overline{v}\leftarrow-\overline{v}$ if $b<0$.
While $v_0 \ne 0$, repeat:
$\qquad$ Calculate the quotient $q=\text{quo}(u_0,v_0)=u_0-v_0\left\lfloor\frac{u_0}{v_0}\right\rfloor$.
$\qquad$ Set $\quad\overline{w}=\overline{u}-q\overline{v},\quad$ then $\quad\overline{u}=\overline{v},\quad$ and then $\quad\overline{v}=\overline{w}.\quad$
Return $a=u_0$, $x=u_1$, and $y=u_2$.
This is doable in not too many lines of C code, and it works as long as $m=10^9+7$ has no (prime) factors $p\le n$ (in fact, this $m$ is prime). If you needed a more efficient algorithm and $m$ were a composite product of distinct primes (which it isn't), you could use the prime factorization of $m$ and a nice fact about binomial coefficients modulo primes [Lukas 1878] to separately calculate residues $$ {a\choose b}\equiv{[a/p]\choose[b/p]}{a\mod p\choose b\mod p}\pmod p $$ modulo each factor $p$, and then recover $T(n)$ using the Chinese Remainder Theorem.
Here's a few pre-computed values: $$\matrix{ k& n=10^k &{2n\choose n}\text{ mod }m &T(n) \\0& 1 &2 &1 \\1& 10 &184756 &923780 \\2& 100 &407336795 &366839610 \\3& 1000 &72475738 &237868748 \\4& 10000 &703593270 &966325381 \\5& 100000 &879467333 &366342189 \\6& 1000000 &192151600 &799327475 }$$
If you want to find an efficient method to solve this problem, you'll probably want to look at the works of Donald Knuth, Andrew Granville, and Richard Stanley, who also has compiled outstanding lists of bijective proof problems and characterizations of the Catalan numbers $C_n=\frac1{n+1}{2n\choose n}$, to which our $S(n)$ are closely related since $S(n)={n+1\choose2}C_n$.
One might be tempted to try to shorten the computation using Wilson's theorem, $$ p\text{ prime} \quad\implies\quad (p-1)!\equiv-1\pmod p $$ Morley's (1895) congruence, $$ p=2n+1\text{ prime} \quad\implies\quad {2n\choose n}\equiv(-1)^n4^{2n}\pmod{p^3} $$ and/or Kummer's theorem using enough "landmark" primes near $n$ and $2n$, with extended Euclidean algorithm for inverses and the CRT (there's the catch!) to get the final result. For example, here are some prime pairs $q_i=2p_i+1$ near $n=10^6$ and $2n$: $$ \matrix{ i & p_i-n & q_i-2n \\ 1 & -251 & -501 \\ 2 & -191 & -381 \\ 3 & 151 & 303 \\ 4 & 193 & 387 \\ } $$ From Wilson's theorem for odd primes $q$, grouping $(q-1)!=(1\cdots\frac{q-1}2)(\frac{q+1}2\cdots q-1)$ and noting that the second term is $(-1)^\frac{q-1}2$ times the first term modulo $q$, we find that $$ \left(\left(\tfrac{q-1}{2}\right)!\right)^2 \equiv(-1)^{\frac{q+1}2} \pmod q $$ so that for prime pairs $q_i=2p_i+1$, $$ {2p_i\choose p_i}\equiv(-1)^{p_i}=-1\pmod{q_i}. $$ Thus we can compute (using the extended Euclidean algorithm for the inverses of $k$ modulo $q_i$) $$ {2n\choose n}\equiv-\prod_{k=p_i+1}^n\frac{4k-2}{k} \pmod{q_i} $$ for $i=1,2$ above. However, we cannot use the CRT to get $T(n)$. We would have to have enough prime pairs to reconstruct $S(n)$, and then compute its residue at the end. Since the central binomial coefficient is approximately $${2n\choose n}\approx\frac{4^n}{\sqrt{\pi n}}\left(1-\frac1{nc_n}\right),\qquad c_n\in(8,9)$$ we would need about 96 thousand pairs (p,q), making this approach infeasible.