There are many ways to realize this aim. For example, in your subroutine, you can copy the current orbital coefficients (e.g. CO or CObasa matrix) to a temporary array, and then call for example "readinfile" subroutine to read wavefunction information from another file.

subroutine "orbcorres" (subfunction 6 of main function 200) utilizes orbitals from two wavefunction files, you can take it as instance.

Dear Multiwfn programmer.

I am interested in developing a real-space function that depends on two different wavefunctions (with the same nuclear configuration). For this purpose, I want to use Multiwfn. Could you explain me or give me some insight of how such a function may be declared into Multiwfn code?