I do no think this is accurate, GAMESS-US does print the molecular orbitals block after the optimization being finished. Maybe I am missing something ? In all my optimization files, after the "***** EQUILIBRIUM GEOMETRY LOCATED *****" it prints the coordinates of all atoms followed by a "MOLECULAR ORBITALS" block. What is missing in that block ? For the 1st step in the geometry optimization this block appears as "EIGENVECTORS" which it is what you are using but the name changes after geometry convergence. I have surpassed this difficult by reading the eigenvectors and doing a PROP type calculation and writing the initial orbitals (PRTMO=.TRUE. in $GUESS group) and then changing the part where it read "INITIAL GUESS ORBITALS" to read "EIGENVECTORS". Hope this information is useful.

What I meant is all orbitals.
This is an example of optimization task: opt_RHF.zip
It can be seen that at initial stage, coefficents of all 24 orbitals are given; however, at the end of optimization, only 15 orbitals are printed. (It seems that GAMESS only prints HOMO+10 at final stage... I didn't check its source code)

Do you want to let Multiwfn automatically load the wavefunction at optimized geometry?
If yes, perhaps this cannot be realized, because I found GAMESS doesn't print expansion coefficients of all orbitals for final geometry. If you have a way to make GAMESS output all orbital coefficients of final geometry in output file, I can modify Multiwfn to load them instead of loading the EIGENVECTORS of the initial geometry.

I do no think this is accurate, GAMESS-US does print the molecular orbitals block after the optimization being finished. Maybe I am missing something ? In all my optimization files, after the "***** EQUILIBRIUM GEOMETRY LOCATED *****" it prints the coordinates of all atoms followed by a "MOLECULAR ORBITALS" block. What is missing in that block ? For the 1st step in the geometry optimization this block appears as "EIGENVECTORS" which it is what you are using but the name changes after geometry convergence. I have surpassed this difficult by reading the eigenvectors and doing a PROP type calculation and writing the initial orbitals (PRTMO=.TRUE. in $GUESS group) and then changing the part where it read "INITIAL GUESS ORBITALS" to read "EIGENVECTORS". Hope this information is useful.

Do you want to let Multiwfn automatically load the wavefunction at optimized geometry?

Yes.

Well. Thanks!
I will check GAMESS documentation and ask for help in GAMESS  google groups. In case a solution is found, I will make an update in this thread to inform you, Professor.
Meanwhile, I will keep using the usual method, i.e., carrying out another GAMESS job using the optimized geometry then loading the output into Multiwfn as input.

I use GAMESS for geometry optimization and get .out files. In order to use Multiwfn, I have to rerun single point energy calculation for the new optimized geometry to use the new .out files as input into Multiwfn  (after changing the extension to .gms, of course)

Whenever I use the optimization .out files in Multiwfn I get this message:

"Warning: This is an optimization task, only wavefunction corresponding to initial geometry will be loaded
Press ENTER button to continue"

Is there a way to get around this problem?

Thanks,