A simple formula is obtained for coupling electrons in a complex system to the electromagnetic field. It includes the effect of intra-atomic excitations and nuclear motion, and can be applied in, e.g., first-principlesbased simulations of the coupled dynamics of electrons and nuclei in materials and molecules responding to ultrashort laser pulses. Some additional aspects of nonadiabatic dynamical simulations are also discussed, including the potential of "reduced Ehrenfest" simulations for treating problems where standard Ehrenfest simulations will fail. DOI: 10.1103/PhysRevB.78.064305 PACS number͑s͒: 78.20.Bh, 71.15.Pd It is now possible to perform first-principles simulations of the coupled dynamics of electrons and nuclei with all nuclear coordinates included 1-4 rather than a subset of nominal reaction coordinates. For very large systems or when many trajectories are necessary, it is convenient to use a first-principles-based scheme 5-8 with a valence-electron Hamiltonian and ion-ion repulsive potential derived from calculations using density functional or other first-principles techniques. Here we are mainly concerned with the issue of how one can efficiently and accurately couple electrons to the electromagnetic field in such an approach, where matrix elements of various operators between localized basis functions ͑or "atomic orbitals"͒ can be calculated from first principles, and then used in large-scale calculations for complex systems, such as materials and molecules, responding to applied fields, such as ultrashort laser pulses. [9][10][11][12][13][14][15][16][17][18][19] Our starting point is, of course, the time-dependent Schrödinger equation, iប ץ ץt ͑x,t͒ = Ĥ ͑x,t͒, ͑1͒Some time ago, Graf and Vogl 20 obtained a result, used in Refs. 13-19, which is the time-dependent version of the Peierls substitution: If H 0 is the Hamiltonian matrix in a localized basis with no applied field,and H is the approximate Hamiltonian when there is an applied field with vector potential A͑x , t͒, then they are related bywith A͑t͒ = ͓A͑XЈ,t͒ + A͑X,t͔͒/2. ͑5͒Here ᐉ labels a localized basis function centered on a nucleus whose instantaneous position is X͑ᐉ , t͒, and we adopt the convention of normally suppressing the indices ᐉ and ᐉЈ as well as the time t by just writing X and XЈ. We will ignore any applied scalar potential A 0 , any B · B spin interactions, and the coupling of ion cores or nuclei to the applied fields since these effects can be easily included when necessary.With the prescription of Eq. ͑4͒, one does not need any new parameters in a calculation that employs either a semiempirical 13,14 or a first-principles-based [15][16][17][18][19] Hamiltonian H 0 whose elements are known as a function of ͑X − XЈ͒. On the other hand, this prescription is in one respect a rather crude approximation: It omits intra-atomic excitations and would therefore give no excitation at all for isolated atoms.Here a more general version of the result of Ref. 20 will be obtained in a form that is almost equally convenient for...