The random phase approximation (RPA) is an increasingly popular method for computing molecular ground-state correlation energies within the adiabatic connection fluctuation-dissipation theorem framework of density functional theory. We present an efficient analytical implementation of first-order RPA molecular properties and nuclear forces using the resolution-of-the-identity (RI) approximation and imaginary frequency integration. The centerpiece of our approach is a variational RPA energy Lagrangian invariant under unitary transformations of occupied and virtual reference orbitals, respectively. Its construction requires the solution of a single coupled-perturbed Kohn-Sham equation independent of the number of perturbations. Energy gradients with respect to nuclear displacements and other first-order properties such as one-particle densities or dipole moments are obtained from partial derivatives of the Lagrangian. Our RPA energy gradient implementation exhibits the same O(N4logâ¡N) scaling with system size N as a single-point RPA energy calculation. In typical applications, the cost for computing the entire gradient vector with respect to nuclear displacements is ∼5 times that of a single-point RPA energy calculation. Derivatives of the quadrature nodes and weights used for frequency integration are essential for RPA gradients with an accuracy consistent with RPA energies and can be included in our approach. The quality of RPA equilibrium structures is assessed by comparison to accurate theoretical and experimental data for covalent main group compounds, weakly bonded dimers, and transition metal complexes. RPA outperforms semilocal functionals as well as second-order Møller-Plesset (MP2) theory, which fails badly for the transition metal compounds. Dipole moments of polarizable molecules and weakly bound dimers show a similar trend. RPA harmonic vibrational frequencies are nearly of coupled cluster singles, doubles, and perturbative triples quality for a set of main group compounds. Compared to the ring-coupled cluster based implementation of Rekkedal et al. [J. Chem. Phys. 2013, 139, 081101.], our method scales better by two powers of N and supports a semilocal Kohn-Sham reference. The latter is essential for the good performance of RPA in small-gap systems.