Molcas Forum

lijingbai2009

Member

Registered: 2019-08-28

Posts: 29

Computing surface hopping probability

Dear OpenMolcas developers,

I have been looking at the source code (tully.f) of surface hopping in OpenMolcas. In the source code, the B matrix is always re-computed with the state population (A matrix) and derivative coupling (D matrix) at each substep. So the surface hopping probability from state 2 to state 1 at time t is g(t) = B12(t)*dT/(nsubstep*A22(t)), where dT is the MD time step (e.g., 20 au), and nsubstep is the number of substeps (e.g., 200).

As you can see, in each substep, g is usually small as it was divided by nsubstep. I noticed that even at some surface hopping events, g was smaller than 0.01, suggesting a sudden state transition.

My question is why the B matrix is re-computed instead of being accumulated throughout the substeps. My understanding is that B12*dT/nsubstep approximates the population of the state transition in a tiny time step. The numerical integration over nsubstep in an MD time step dT would approximately be the sum of B12(t)*dT/nsubstep. The probability g(t) = [(B12(t0)+...+B12(t))*dt/nsubstep]/A22(t). This can give a gradually increasing probability near the surface hopping event.

Do I have any misunderstandings? Thanks!