This is possible in a couple different ways. The easiest way (i.e. the way that doesn't require modifying the source code) is to concatenate your two geometries into one geometry and then run a calculation on the combined geometry. For example, if you were interested in the AO overlap matrix between H2 with a bond length of 1.0 and a bond length of 1.1, you could set up the combined geometry as
set geom:dont_verify t
geometry noautosym nocenter
H 0.0 0.0 0.5
H 0.0 0.0 0.5
H 0.0 0.0 0.5
H 0.0 0.0 0.6
end
Important things to note, you will probably need to tell the code not to check if you are using a reasonable geometry since I imagine you will have overlapping (or nearly overlapping) atoms. This is accomplished with the "set geom:dont_verify t" directive. Also it is probably a good idea, if not required, to turn off symmetry and centering of the system in order to make sure you are getting the overlap that you want.
Then you just need to add
print "ao overlap"
to your DFT/SCF input block to get the code to print the overlap matrix to the output file. The ordering of the elements in the overlap matches the ordering of the atoms in the geometry input. So the printed overlap matrix will have a structure like
AA AB
BA BB
where AA corresponds to the AO overlap matrix for the first geometry, BB is for the second geometry, and AB/BA is for the overlap of the first geometry with the second geometry. The dimensions of each block are the number of basis functions for each geometry. For example, if we used the STO3G basis with the H2 example above (2 basis functions per molecule), S(1:2,1:2) would correspond to the AO overlap matrix for the H2 molecule with bond length 1.0; S(3:4,3:4) would correspond to the AO overlap matrix for the H2 molecule with bond length 1.1; and S(1:2,3:4)=S(3:4,1:2) would be the AO overlap matrix you are interested in.
