A computer docking study has been carried out on the crystal surfaces of cellulose Ialpha crystal models for the carbohydrate binding module (CBM) protein of the cellobiohydrolase Cel7A produced by Trichoderma reesei. Binding free energy maps between the CBM and the crystal surface were obtained by calculating the noncovalent interactions and the solvation free energy at grid points covering the area of the unit cell dimensions at the crystal surface. The potential maps obtained from grid searches of the hydrophobic (110) crystal surface exhibited two distinct potential wells. These reflected the 2-fold helical symmetry of the cellulose chain and had lower binding energies at the minimum positions than those for the hydrophilic (100) and (010) crystal surfaces. The CBM-cellulose crystal complex models derived from the minimum positions were then subjected to molecular dynamics (MD) simulation under an explicit solvent system. The (110) complex models exhibited larger affinities at the interface than the (100) and (010) ones. The CBM was more stably bound to the (110) surface when it was placed in an antiparallel orientation with respect to the cellulose fiber axis. In the solvated dynamics state, the curved (110) surface resulting from the fiber twist somewhat assisted a complementary fit with the CBM at the interface. In addition to the conventional Generalized Born (GB) method, the three-dimensional reference interaction site model (3D-RISM) theory was adopted to assess a solvent effect for the solvated MD trajectories. Large exothermic values for the noncovalent interactions appeared correlated to and were mostly compensated by endothermic values for the solvation free energy. These gave total binding free energies of -13 to -28 kcal/mol. Results also suggested that the hydrogen bonding scheme was not essential for substrate specificity.