This paper proposes and analyzes a fully discrete scheme for nonlinear biharmonic Schrödinger equations. We first write the single equation into a system of problems with second-order spatial derivatives and then discretize the space variable with an ultraweak discontinuous Galerkin scheme and the time variable with the Crank–Nicolson method. The proposed scheme proves to be computationally more efficient compared to the local discontinuous Galerkin method in terms of the number of equations needed to be solved at each single time step, and it is unconditionally stable without imposing any penalty terms. It also achieves optimal error convergence in L2 norm both in the solution and in the auxiliary variable with general nonlinear terms. We also prove several physically relevant properties of the discrete schemes, such as the conservation of mass and the Hamiltonian for the nonlinear biharmonic Schrödinger equations. Several numerical studies demonstrate and support our theoretical results.