Try the following route line (where you change B3LYP/6-31G* to your method of choice):

#P scf(conventional) B3LYP/6-31G* Iop(3/33=6) Extralinks=L316 Noraff Symm=Noint Iop(3/33=1) pop(full)

In the log file, look for the section starting:

*** Dumping Two-Electron integrals ***

The integrals are given according to I, J, K and L indices (rather than r,s,t,u)

Thank you deeply for sharing your valuable advices.

I mistakenly thought that the two-electron integral was also provided in a triangular matrix, like *** Overlap *** and *** Kinetic Energy ***.

Perhaps you can try: https://joshuagoings.com/2013/09/02/gaussian-matrix-parser/

Dear Prof. Paton,

Thank you for your valuable information.

I will try it immediately.

For DFT calculation, 1/2 , Ex and Ec are printed in Hartree units to the log file. Isn't there a way to have each of these output as an nxn matrix?