The highly polymorphic human cytochrome P450 2D6 (CYP2D6) metabolizes about 25% of currently used drugs. In this study, we have explored the interaction of a large number of substrates (n = 120) with wild-type and mutated CYP2D6 by molecular docking using the CDOCKER module. Before we conducted the molecular docking and virtual mutations, the pharmacophore and QSAR models of CYP2D6 substrates were developed and validated. Finally, we explored the interaction of a traditional Chinese herbal formula, Fangjifuling decoction, with CYP2D6 by virtual screening. The optimized pharmacophore model derived from 20 substrates of CYP2D6 contained two hydrophobic features and one hydrogen bond acceptor feature, giving a relevance ratio of 76% when a validation set of substrates were tested. However, our QSAR models gave poor prediction of the binding affinity of substrates. Our docking study demonstrated that 117 out of 120 substrates could be docked into the active site of CYP2D6. Forty one out of 117 substrates (35.04%) formed hydrogen bonds with various active site residues of CYP2D6 and 53 (45.30%) substrates formed a strong π-π interaction with Phe120 (53/54), with only carvedilol showing π-π interaction with Phe483. The active site residues involving hydrogen bond formation with substrates included Leu213, Lys214, Glu216, Ser217, Gln244, Asp301, Ser304, Ala305, Phe483, and Phe484. Furthermore, the CDOCKER algorithm was further applied to study the impact of mutations of 28 active site residues (mostly non-conserved) of CYP2D6 on substrate binding modes using five probe substrates including bufuralol, debrisoquine, dextromethorphan, sparteine, and tramadol. All mutations of the residues examined altered the hydrogen bond formation and/or aromatic interactions, depending on the probe used in molecular docking. Apparent changes of the binding modes have been observed with the Glu216Asp and Asp301Glu mutants. Overall, 60 compounds out of 130 from Fangjifuling decoction matched our pharmacophore model for CYP2D6 substrates. Fifty four out of these 60 compounds could be docked into the active site of CYP2D6 and 24 of 54 compounds formed hydrogen bonds with Glu216, Asp301, Ser304, and Ala305 in CYP2D6. These results have provided further insights into the factors that determining the binding modes of substrates to CYP2D6. Screening of high-affinity ligands for CYP2D6 from herbal formula using computational models is a useful approach to identify potential herb-drug interactions.