The paper deals with the reliability–based design optimization (RBDO) of concrete gravity dams subjected to earthquake load using subset simulation. The optimization problem is formulated such that the optimal shape of concrete gravity dam described by a number of variables is found by minimizing the total cost of concrete gravity dam for the given target reliability. In order to achieve this purpose, a framework is presented whereby subset simulation is integrated with a hybrid optimization method to solve the RBDO approach of concrete gravity dam. Subset simulation with Markov Chain Monte Carlo (MCMC) sampling is utilized to estimate accurately the failure probability of dams with a minimum number of samples. In this study, the concrete gravity dam is treated as a two–dimensional structure involving the material nonlinearity effects and dam–reservoir–foundation interaction. An efficient metamodel in conjunction with subset simulation–MCMC is provided to reduce the computational cost of dynamic analysis of dam–reservoir–foundation system. The results demonstrate that the RBDO approach is more appropriate than the deterministic optimum approach for the optimal shape design of concrete gravity dams.