Incorporation of the Hubbard U potential into the Kohn-Sham scheme has extended the realm of the density functional theory (DFT) to correlated systems, such as Mott-insulator state of transition metal oxides. In this scheme, the Coulomb and exchange energies are parametrized in the projected subspace. In a sense, the presence of the Hubbard U potential restores the derivative discontinuity otherwise missed in the continuous exchange-correlation potential. The DFT plus Hubbard U potential scheme is also effective for the surface adsorption problem in which a long-range charge transfer occurs between a solid surface and isolated atomic or molecular species. For the projection, localized atomic orbitals centered on each atom have been used. However, for the adsorption of a highly symmetric molecule, the implementation of the Hubbard U method on such atomic orbitals does not work because of the molecular symmetry: for example, in diatomic molecules like O2 and F2, the molecular orbitals comprise the bonding and antibonding state of atomic orbitals. Thus the density matrices with various occupations should not break such molecular symmetry, and the molecular orbital, not atomic orbitals, should be used as a unit for the correction of correlation. Here, we propose to add the Hubbard U potentials on molecular orbitals. The method for the linear combination of atomic orbitals is explained in which the molecular symmetry is preserved in the density matrices for various occupations. We derived the formula for the atomic force as gradient with respect to atom positions over the total energy functional with the Hubbard U potential on the molecular orbitals. Bench marking test with other method, including wavefunction theories, are also presented.