This information is written into the .rwf file generated by Gaussian calculations. To save the .rwf file beyond the normal termination of the run, you have to explicitly name it in the input file. I recommend you send a message to Pedro A Derosa and ask him to provide more information. It is possible that MOPAC, GAMESS, etc., also provide this info.
I haven't done this in a while, but I know you can extract this information from a Gaussian calculation. You need to save the rwf file as Ramu Ramachandran said . Then use the rwfdump utility, the syntax is:
rwfdump rwf_file_name.rwf yourfile code.
Where code is the code for the matrix you are looking for. For instance 514r is the code for the overlap matrix 536r and 538r will dump the Fock matrices (alpha and beta respectively). "r" is a format instruction. It'll save it as a lower triangular matrix. It is formatted in 5 columns so the way you read it is, (from left to right and from top to bottom
You can see the full list in the programmer's reference manual, chapter 3 I believe.
now the self-energy terms, you will have to calculate those with some input from the electrodes, it won't come from a molecular calculation.
You mention ATK, I don't know if you can extract the matrices you mention, but ATK will do the NEGF calculation, it is actually the original purpose of that software. The same for QuantumExpresso.