Lattice energy calculations using a model potential have been performed to model the crystal structures of cis-1,2,3,6-and 3,4,5,6-tetrahydrophthalic anhydrides. The optimized molecular models using the DFT method at the B3LYP/6-31G** level were found consistent with the available experimental evidence and allow to reproduce all differences observed in crystal packing between cis-1,2,3,6-and 3,4,5,6-tetrahydrophthalic anhydrides. Calculations provide evidence for the presence of dipole-dipole C=O…C=O intermolecular interactions and support the idea that the molecules distort from their ideal geometries, improving packing in both crystals. The search for minima in the lattice energy of both crystals amongst the more common space groups with Z'=1, using a simulated annealing crystal structure prediction procedure followed by lattice energy minimisation shows that the observed structure of 3,4,5,6-tetrahydrophthalic anhydride (Z'=2) is the thermodynamically most stable and allow us to justify why 3,4,5,6-tetrahydropthalic anhydride crystallises in such complex structure with 16 molecules in the unit cell. The computational model was successful to predict the second observed form at 173 K for cis-1,2,3,6-tetrahydropthalic anhydride as a polymorph and could predict several hypothetical structures with Z'=1 which appear competitive with the observed structures. The results of phonon estimates of zero point intermolecular vibrational energy and entropy suggest that crystal structures of cis-1,2,3,6-tetrahydropthalic anhydride cannot be predicted just on the basis of lattice energy, there are yet, other factors than thermodynamics favoring the observed structures.