We investigate a method of combining the techniques of quantum chemistry with QED. In our theory, we treat the N-electron system and the Dirac sea on an equal footing; we regard both of them as the dynamical degrees of freedom of a many-body system. After the introduction of our QED-SCF method, the QED-SCF solutions are classified on the basis of group-theoretical operations such as time-reversal, parity and O(3) rotational symmetry. The natural orbitals of general QED-SCF solutions are determined by diagonalizing the first-order density matrix. Thus, we obtain the possibility of treating the system under strong QED effect by the methods of quantum chemistry, such as QED-MCSCF and QED coupled-cluster approaches. (C) 2001 John Wiley & Sons, Inc.