Metamaterial antennas are appealing for next-generation wireless networks due to their simplified hardware and much-reduced size, power, and cost. This paper investigates the holographic multiple-input multiple-output (HMIMO)-aided multi-cell systems with practical per-radio frequency (RF) chain power constraints. With multiple antennas at both base stations (BSs) and users, we design the baseband digital precoder and the tuning response of HMIMO metamaterial elements to maximize the weighted sum user rate. Specifically, under the framework of block coordinate descent (BCD) and weighted minimum mean square error (WMMSE) techniques, we derive the low-complexity closed-form solution for baseband precoder without requiring bisection search and matrix inversion. Then, for the design of HMIMO metamaterial elements under binary tuning constraints, we first propose a low-complexity suboptimal algorithm with closed-form solutions by exploiting the hidden convexity (HC) in the quadratic problem and then further propose an accelerated sphere decoding (SD)-based algorithm which yields global optimal solution in the iteration. For HMIMO metamaterial element design under the Lorentzian-constrained phase model, we propose a maximization-minorization (MM) algorithm with closed-form solutions at each iteration step. Furthermore, in a simplified multiple-input single-output (MISO) scenario, we derive the scaling law of downlink single-to-noise (SNR) for HMIMO with binary and Lorentzian tuning constraints and theoretically compare it with conventional fully digital/hybrid arrays. Simulation results demonstrate the effectiveness of our algorithms compared to benchmarks and the benefits of HMIMO compared to conventional arrays.