  We consider harmonic maps$u(z): \mathcal{X}_z\to N$ in a fixed homotopy class
from Riemann surfaces $\mathcal{X}_z$ of genus $g\geq 2$ varying in the
Teichm\"u{}ller space $\mathcal T$ to a Riemannian manifold $N$ with
non-positive Hermitian sectional curvature.
  The energy function $E(z)=E(u(z))$ can be viewed as a function on $\mathcal
T$ and we study its first and the second variations. We prove that the
reciprocal energy function $E(z)^{-1}$ is plurisuperharmonic on Teichm\"uller
space. We also obtain the (strict) plurisubharmonicity of $\log E(z)$ and
$E(z)$. As an application, we get the following relationship between the second
variation of logarithmic energy function and the Weil-Petersson metric if the
harmonic map $u(z)$ is holomorphic or anti-holomorphic and totally geodesic,
i.e., $$ \sqrt{-1}\p\b{\p}\log E(z)=\frac{\omega_{WP}}{2\pi(g-1)}. $$ We
consider also the energy function $E(z)$ associated to the harmonic maps from a
fixed compact K\"ahler manifold $M$ to Riemann surfaces
${\mathcal{X}_z\}_{z\in\mathcal{T}}$ in a fixed homotopy class. If $u(z)$ is
holomorphic or anti-holomorphic, then the above equation is also proved.
