It is shown how first-order nonadiabatic coupling matrix elements can be calculated using coupled cluster electronic structure methods. The formalism is consistent with the coupled cluster response theory approach for calculation of excitation energies and adiabatic transition properties. Expressions are derived that are in the limit of a complete coupled cluster expansion give results equivalent to the full configuration interaction results. Computational tractable expressions are given for the first-order nonadiabatic coupling matrix in coupled cluster theory. The final expressions are quite similar to those employed in the implementation of the analytical calculation of molecular gradients.