Computational chemistry simulates chemical structures and reactions numerically, based on the fundamental laws of physics using two different methods including electronic structure and molecular mechanics. In this work, a four dimension intermolecular potential energy surface, for CS2-CS2 system, is obtained by ab initio calculation of the potential energies for a range of configurations and intermolecular distances. The calculations were performed using the supermolecular approach on the MP2 level of theory with the augmented correlation consistent (aug-cc-pVXZ(X=D,T) basis sets using rigid rotor approximation. The number of configuration points were selected for calculations is 1387. The calculated energy points were corrected for the basis set superposition error using the full counterpoise correction method suggested by Boys and Bernardi. A two point extrapolation method was used to extrapolate the calculated correlation energy points to the complete basis set limit and the Hartree-Fock energy part has not been extrapolated but taken from the aug-cc-pVTZ calculations. The extrapolated intermolecular energy points were fitted to the exp-6 potential function to represent the potential energy surface analytically. The results show that the most stable configuration of CS2-CS2 system is cross configuration ( q 1 =90°, q 2= 90° and f =90° at R=6.23 Bohr), with a well depth of 3.98±0.018 kcal/mol and the most unstable configuration is linear ( q 1 =0°, q 2= 0° and f =0° at R=12.548 Bohr) , with a well depth of 1.48±0.018 kcal/mol. The second virial coefficients have been calculated using the scaled potential energy surface over a range of temperature (280-470 K) and compared with the experimental values to test the quality of the presented potential energy surface. Also, we used global fitting to have a simple analytical spherical one-center intermolecular potential using exp-6 and lennard-jones potential, separately. The parameters of the potential were obtained by radial fitting of the intermolecular interaction potential of all configurations with exp-6 and lennardjones simultaneously. The obtained potential parameters for the exp-6 and lennard-jones potentials are ( e = .00265±.0000943 hartree, Rm = 6.45 ±.0258Bohr, a =10.0 ±.231 )and( e = .00304028±.0000972 hartree, s = 5.478±.0122 ), respectively. .