How do I combine z-matrices of multiple molecules into one input stack?
The Setup

The first thing you need to do is to get the difference between a Z-matrix (sometimes also called internal coordinates) and cartesian coordinates straight. They are truly different things.

Secondly, take the Z-matrices of a sufficiently large molecule. The basic structure is:

C1 .  .  .
C2 r2 .  .
C3 r3 a3 .
C4 r4 a4 d4
C5 r5 a5 d5
(and so on from this point)


where in this simplified notation, I have omitted the anchor centers for the coordinates. One notices that there are 6 coordinates left out, marked by dots. They can be omitted because three degrees of freedom for translation and three degrees of freedom for rotation of the entire molecule are irrelevant for the calculation of properties.

Now imagine a second molecule tacked onto the other one:

C1 .  .  .
C2 r2 .  .
C3 r3 a3 .
C4 r4 a4 d4
C5 r5 a5 d5
C6 *  *  *
C7 r7 *  *
C8 r8 a8 *
C9 r9 a9 d9


This Z-matrix is invalid. One needs to fill the missing coordinates marked by the asterisks. In a sense, the second molecule does not have the 6 degrees of freedom mentioned above, because any change in them would change the orientation relative to the first one, impacting all properties.

So, what can you do?

Filling in the missing coordinates is hard, at the very least by hand or some other automated fashion. Doing it in an editor is basically the same amount of work as the solution I actually suggest: Create the system of both molecules in one Z-matrix and split that Z-matrix. Remove the superfluous coordinates in the resulting matrices and run calculations as necessary.