A four-dimensional intermolecular potential energy surface for the carbon dioxide dimer has been computed using the many-body symmetry-adapted perturbation theory (SAPT) and a large 5s3p2d1f basis set including bond functions. The SAPT level applied is approximately equivalent to the supermolecular many-body perturbation theory at the second-order level. An accurate fit to the computed data has been obtained in a form of an angular expansion incorporating the asymptotic coefficients computed ab initio at the level consistent with the applied SAPT theory. A simpler site-site fit has also been developed to facilitate the use of the potential in molecular dynamics and Monte Carlo simulations. The quality of the new potential has been tested by computing the values of the second virial coefficient which agree very well with the experimental data over a wide range of temperatures. Our potential energy surface turns out to be substantially deeper than previous ab initio potentials. The minimum of -484 cm(-1) has been found for the slipped parallel geometry at the intermolecular separation R = 3.54 Angstrom and a saddle point at -412 cm(-1) for the T-shaped configuration and R = 4.14 Angstrom. Three minima and two first-order saddle points have been located on the pairwise-additive potential energy surface of the CO2 trimer. The nonplanar structure of C-2 symmetry has been found to be 48.8 cm(-1) more stable than the cyclic planar form of C-3h symmetry, in disagreement with experimental observation. It is suggested that the relative stability of the two isomers cannot be reliably determined by pairwise-additive potential and inclusion of three-body forces is necessary for this purpose.