In the present study, we carried out thermodynamic integration molecular dynamics simulation for a pair of analogous inhibitors binding with Erk kinase to investigate how computation performs in reproducing the relative binding free energy. The computation with BCC-AM1 charges for ligands gave −1.1 kcal/mol, deviated from experimental value of −2.3 kcal/mol by 1.2 kcal/mol, in good agreement with experimental result. The error of computed value was estimated to be 0.5 kcal/mol. To obtain convergence, switching vdw interaction on and off required approximately 10 times more CPU time than switching charges. Residue-based contributions and hydrogen bonding were analyzed and discussed. Furthermore, subsequent simulation using RESP charge for ligand gave ΔΔG of −1.6 kcal/mol. The computed results are better than the result of −5.6 kcal/mol estimated using PBSA method in a previous study. Based on these results, we further carried out computations to predict ΔΔG for five new analogs, focusing on placing polar and nonpolar functional groups at the meta site of benzene ring shown in the Fig. 1, to see if these ligands have better binding affinity than the above ligands. The computations resulted that a ligand with polar –OH group has better binding affinity than the previous examined ligand by ~2.0 kcal/mol and two other ligands have better affinity by ~1.0 kcal/mol. The predicted better inhibitors of this kind should be of interest to experimentalist for future experimental enzyme and/or cell assays.