The analytical total energy gradient scheme for the Hartree-Fock and density functional crystalline orbital theory is formulated for infinitely extended periodic systems of general dimensions and implemented for those of two dimensions. Two major differences between the analytical gradient scheme for extended systems and that for molecular systems are described in detail. The first is the treatment of the long-range Coulomb interactions, which arise due to the infinite nature of the system size. The long-range effect is efficiently included by the multipole expansion technique and its extension, the fast multipole method. The use of the fast multipole method enables us to include the long-range effect up to the order of micrometer to millimeter region around the reference unit cell by virtue of the logarithmic cost scaling of the algorithm achieved by regrouping distant multipoles together and reducing the number of pairwise interactions. The second is the formulation of analytical gradient expressions with respect to unit cell parameters. In HF theory they can be calculated by accumulating forces acting on atoms multiplied by some appropriate factors, while there is an extra term which requires a special numerical treatment in grid-based density functional theory. Specifically, it is shown that the quadrature weight derivatives do not vanish even in the limit of infinitely fine grid when calculating the gradients with respect to unit cell parameters, and are essential in evaluating those gradients. Combining the analytical gradient scheme and an efficient inclusion of the long-range interaction makes it feasible to perform a full geometry optimization of extended systems at ab initio levels. As an illustration, the long-range interaction energies are computed for a two-dimensional sheet of hydrogen-fluoride. The CPU time reduction on going from the explicit evaluation of the two-electron integrals to multipole expansion, and from the multipole expansion to fast multipole method is significant. Geometry optimizations are performed on an infinite two-dimensional hexagonal boron-nitride sheet and the dependence of the gradients on various parameters are investigated. (C) 2003 American Institute of Physics.