Abstract We have performed N-body simulation on final accretion stage of terrestrial planets, including the effect of damping of eccentricity and inclination caused by tidal interaction with a remnant gas disk. As a result of runway and oligarchic accretion, about 20 Mars-sized protoplanets would be formed in nearly circular orbits with orbital separation of several to ten Hill radius. The orbits of the protoplanets would be eventually destabilized by long-term mutual gravity and/or secular resonance of giant gaseous planets. The protoplanets would coalesce with each other to form terrestrial planets through the orbital crossing. Previous N-body simulations, however, showed that the final eccentricities of planets are around 0.1, which are about 10 times higher than the present eccentricities of Earth and Venus. The obtained high eccentricities are the remnant of orbital crossing. We included the effect of eccentricity damping caused by gravitational interaction with disk gas as a drag force (“gravitational drag”) and carried out N-body simulation of accretion of protoplanets. We start with 15 protoplanets with 0.2 M⊕ and integrate the orbits for 10 7 years, which is consistent with the observationally inferred disk lifetime (in some runs, we start with 30 protoplanets with 0.1 M⊕). In most runs, the damping time scale, which is equivalent to the strength of the drag force, is kept constant throughout each run in order to clarify the effects of the damping. We found that the planets' final mass, spatial distribution, and eccentricities depend on the damping time scale. If the damping time scale for a 0.2 M⊕ mass planet at 1 AU is longer than 10 8 years, planets grow to Earth's size, but the final eccentricities are too high as in gas-free cases. If it is shorter than 10 6 years, the eccentricities of the protoplanets cannot be pumped up, resulting in not enough orbital crossing to make Earth-sized planets. Small planets with low eccentricities are formed with small orbital separation. On the other hand, if it is between 10 6 and 10 8 years, which may correspond to a mostly depleted disk (0.01–0.1% of surface density of the minimum mass model), some protoplanets can grow to about the size of Earth and Venus, and the eccentricities of such surviving planets can be diminished within the disk lifetime. Furthermore, in innermost and outermost regions in the same system, we often find planets with smaller size and larger eccentricities too, which could be analogous to Mars and Mercury. This is partly because the gravitational drag is less effective for smaller mass planets, and partly due to the “edge effect,” which means the innermost and outermost planets tend to remain without collision. We also carried out several runs with time-dependent drag force according to depletion of a gas disk. In these runs, we used exponential decay model with e-folding time of 3×10 6 years. The orbits of protoplanets are stablized by the eccentricity damping in the early time. When disk surface density decays to ⋦1% of the minimum mass disk model, the damping force is no longer strong enough to inhibit the increase of the eccentricity by distant perturbations among protoplanets so that the orbital crossing starts. In this disk decay model, a gas disk with 10 −4–10 −3 times the minimum mass model still remains after the orbital crossing and accretional events, which is enough to damp the eccentricities of the Earth-sized planets to the order of 0.01. Using these results, we discuss a possible scenario for the last stage of terrestrial planet formation.