We present a general method suitable for automatic generation of trial wave functions for generalized valence bond (GVB) descriptions of large molecules. This method uses pseudo-Haruee-Fock (P-HF) molecular orbitals formed from HF atomic orbitals but without Fock matrix diagonalization. The occupied P-HF orbitals are projected onto atomic basis functions to obtain GVB first natural orbitals, and the unoccupied HF orbitals are projected to obtain GVB second natural orbitals. This method (denoted GVB-INIT) is fast because no HF wave functions need be calculated and because the localization is piecewise atomic. In conjunction with the recently developed GVB-DIIS method for converging GVB wave functions and the new pseudospectral programs (PS-GVB) for the Fock matrix elements, GVB-INIT makes calculation of highly correlated GVB wave functions quite practical. The efficacy of GVB-INIT is illustrated by application to several cases including GVB wave functions with up to 26 correlated pairs.