Introduction: Since the discovery of monolayer graphene in the year 2004 [2], this two-dimensional material remains in the focus of active research motivated by its novel physical properties and a promising potential for applications[3–5]. Experimental groups have enabled preparation and study of systems with one or a small number of graphene layers [6]. Bilayer graphene, which is made of two stacked graphene layers, is considered to be particularly importance for electronics applications because of its special band structure [7]. Coupling of the two monolayer graphene sheets in the usual A-B stacking of bilayer graphene yields pairs of hyperbolically dispersing valence and conduction bands that are split from one another by the interlayer interaction[8]. The band gap of the bilayer can be easily opened when a difference between the electrostatic potential of the two layer is introduced, either by chemical doping or by applying gate voltage [9–17], which makes graphene channels have a high resistance for the OFF state. Angle-resolved photoemission (ARPES) measurements indicate such a gap in potassium doped bilayer graphene epitaxially grown on SiC [14], and infrared spectroscopy measurements also detect the similar gap in the electrostatically gated bilayer graphene [10–12]. In this work, we carried out first-principle calculations and theoretical analysis to explore the electronic and magnetic properties of bilayer graphene with intercalated atoms C, N and O. The calculations were performed using the projector augmented wave (PAW) formalism of density functional theory (DFT) as implemented in the Vienna ab initio simulation pack- age (VASP) [18]. Because generalized gradient approximation (GGA) [19] gives essentially no bonding between graphene planes and leads to excessively large values of bilayer distance [13], we performed the calculations within the localized density approximation (LDA). We find that LDA gives rise to a bilayer distance of 3.34 ˚ A , in good agreement with the ex- periment value. Some previous investigations used LDA to optimize the structure to get a reasonable bilayer distance, and sequently used GGA to calculate the electronic structure [9]. To keep consistent, we use LDA in all our calculations for the investigated system. An energy cutoff of 400 eV for plane-wave expansion of the PAWs is used. The model system here consists of a 4 × 4 supercell with a foreign atom intercalated between the two coupling graphene layers. The supercell parameters are set to be the same as a = b = 9 . 84 ˚ A in the xy plane ( a and b indicate the crystal lattice constants). The Brillouin zone is sampled using a 11 × 11 × 1 Γ centered k-point grid. For geometry optimization, all the internal coordinates are relaxed until the Hellmann-Feynman forces are less than 0.01 eV/ ˚ A. The 2