With Schaefer and Harris method, i.e., constructing LS eigenfunctions by diagonlization of L²+λS², a FORTRAN code with two modes is developed to calculate the most concise LS eigenfunctions under certain symmetry and all LS coupling atomic states (both non-equivalent and equivalent electronic configurations) with given orbital occupancy. As an example, LS coupling atomic states and the simplest LS eigenfunctions for ⁶I symmetry of Np (Z=93) 5f⁴6d electronic occupancy are calculated.