Please look at this discussion on how to compute the intermolecular potential energy. For a molecular system described by a type-I force field, the cohesive energy corresponds to the pairwise intermolecular energy [citation needed]. Including the contribution of long-range electrostatics is tricky, and, to me, contradictory advice has been given on the LAMMPS mailing list.
From the discussion I referred to, it seems that you can compute the intermolecular energy as a variable:
compute intra all group/group all pair yes kspace yes molecule intra
variable enb equal evdwl+ecoul+elong
variable inter equal v_enb-c_intra
Then, you normalise it by volume, molecules, or atoms, depending on your needs.
Could you please report your findings and any comparison with reference data? That would be very useful for the community (and for me to finish this work!).