We propose an approach that is under the framework of Gutzwiller wave function but goes beyond the commonly adopted Gutzwiller approximation to improve the accuracy and flexibility in treating the correlation effects. Detailed formalism is described for a dimer which is straightforwardly generalized later to more complicated periodic bulk systems. The accuracy of the approach is demonstrated by evaluating the potential energy curves of spin-singlet N2 dimer, spin-triplet O2 dimer, and 1D hydrogen chain. The computational workload of the approach can be easily handled by efficient parallel computing.