In this paper, we introduce and discuss an exactly energy-conserving particle-in-cell method for arbitrary curvilinear coordinates. The flexibility provided by curvilinear coordinates enables the study of plasmas in complex-shaped domains by aligning the grid to the given geometry or by focusing grid resolution on regions of interest without overresolving the surrounding, potentially uninteresting domain. We have achieved this through the introduction of the metric tensor, the Jacobian matrix, and contravariant operators combined with an energy-conserving fully implicit solver. We demonstrate the method’s capabilities using a Python implementation to study several one- and two-dimensional test cases: the electrostatic two-stream instability, the electromagnetic Weibel instability, and the geomagnetic environment modeling reconnection challenge. The test results confirm the capability of our new method to reproduce theoretical expectations (e.g., instability growth rates) and the corresponding results obtained with a Cartesian uniform grid when using curvilinear grids. Simultaneously, we show that the method conserves energy to machine precision in all cases.