The metabolic mechanisms for imidazolidine oxidation of imidacloprid (IMI) by cytochrome P450 3A4 (CYP3A4) have been investigated using quantum mechanical/molecular mechanical (QM/MM) calculations. The binding mode of CYP3A4 with IMI is examined by molecular docking in collaboration with molecular dynamics (MD) simulations. The results show that there are six amino acid residues, involving Arg192, Phe195, Ile349, Ala285, Phe284 and Phe88, closely distributed around the IMI. The binding free energy analysis exhibits that the CYP3A4-IMI binding structure is stabilized by electrostatic interaction and van der Waals interaction. Arg192 plays a major role in the binding of CYP3A4 with IMI based on its polarity and the hydrogen bond between the H atom in Arg192 side chain and the nitryl O atom of IMI. Two possible pathways, pathway 1 and pathway 2, are evaluated. Two spin states of the Fe (III) center, quartet and doublet, are considered. The free energy calculations are done using QM/MM steered molecular dynamics (SMD) simulation at the B3LYP/6-31 + G(d):ff14SB level for two pathways. The ONIOM QM/MM single-point calculations at the B3LYP/6-311 + G(2d,2p):ff99SB//B3LYP/6-31 + G(d): ff14SB and M06-2X/6-311 + G(2d,2p):ff99SB//B3LYP/6-31 + G(d):ff14SB levels are carried out to obtain more credible energy information. The results indicate that for both pathways, the free energy barriers on the low-spin doublet state are lower than those on the high-spin quartet state. Both pathways are the stepwise processes. Pathway 1 has higher possibility to occur with the free energy barriers being lower by 10-15 kcal·mol-1 compared with pathway 2, which gives rise to trans-5'-hydroxyl-IMI as the final product. The first proton-transfer is the rate-limiting step and the calculated activation free energy is consistent with the experimental conclusion.