The quantum dynamics of the primary photoisomerization event in bacteriorhodopsin is studied by a semiclassical trajectory approach. The relevant surface crossing probability is evaluated from the wave functions and potential surfaces of a hybrid quantum mechanical/molecular mechanics (QM/MM) Hamiltonian of the complete chromophore-protein-solvent system. The QM/MM model combines consistently the quantum mechanical Hamiltonian of the chromophore with the microscopic electric field of the ionized groups and induced dipoles of the protein-solvent system. The QCFF/PI Hamiltonian of the chromophore is adjusted to reproduce relevant ab initio results. The nonadiabatic coupling term < psi (1)/partial derivative psi (0)/partial derivativet > calculated numerically from the corresponding wave functions. The simulations are performed by combining the ENZYMIX and QCFF/PI molecular modeling programs. The effect of the protein on the absorption spectrum of the chromophore is examined. It is found that this spectrum reflects the effect of the protein permanent dipoles, ionized residues, water molecules (in and around the protein), and the induced dipoles of the protein plus water system. Next, we probe the motion along the excited state surface. It is demonstrated, in agreement with our early study and more recent works, that the motion starts with bond vibrations and evolves to a torsional motion. It is also found that we are dealing with an overdumped motion. Major emphasis is placed on the nature of the surface crossing process. In particular, we try to examine the origin of the very large probability of crossing in the pi /2 region. A large crossing probability was obtained first in our early simulation (Warshel, A. Nature 1976, 260, 679), but its origin was not explored in details. Such large crossing probabilities can be obtained by passing through strict conical intersections (where the, two surfaces "touch" each other) or by passing through regions with large nonadiabatic coupling and small energy gap (such regions are usually close to conical intersections). It is found that some trajectories pass through strict conical intersections whereas others cross through regions with nonzero energy gap and a large nonadiabatic coupling. This feature helps probably to ensure the stability of the photobiological process with regards to various mutations. The average surface crossing probability and our previously derived expression (Weiss, R. M.; Warshel, A. J. Am. Chem. Soc. 1979, 101, 6131) appear to provide an excellent approximation for the calculated quantum yield. Furthermore, the calculated quantum yield reproduces the corresponding observed value. finally, we examine the behavior of trajectories that cross to the ground state before the pi /2 region. Our finding that these trajectories are deflected backward allow us to exclude models where the surface crossing occurs before the pi /2 region.