Analytic gradient expressions for the spin-conserving and spin-flipping equation-of-motion coupled-cluster models with single and double substitutions are derived using a Lagrangian approach for the restricted and unrestricted Hartree-Fock references, both for the case of all orbitals being active in correlated calculations and for the frozen core and/or virtual orbitals. Details of the implementation within the Q-CHEM electronic structure package are discussed. The capabilities of the new code are demonstrated by application to cyclobutadiene.