@@ -106,23 +106,23 @@ subroutine lz_init(pot)
106
106
end subroutine lz_init
107
107
108
108
subroutine lz_init_terash ()
109
- use mod_general, only: natom
110
- use mod_sh, only: en_array, tocalc, nacx, nacy, nacz
111
-
112
- nstate = nstate_lz ! Needed in init_terash
113
- inac = 2 ! Turns off couplings calculation
114
-
115
- ! Based on sh_init() routine, sharing most of the functions
116
- ! TODO: Break dependence - separation of MPI interface needed
117
- allocate (en_array(nstate_lz))
118
- allocate (nacx(natom, nstate_lz, nstate_lz))
119
- allocate (nacy(natom, nstate_lz, nstate_lz))
120
- allocate (nacz(natom, nstate_lz, nstate_lz))
121
- allocate (tocalc(nstate, nstate))
122
- en_array = 0.0D0
123
- tocalc = 0
124
- tocalc(istate_lz, istate_lz) = 1
125
- call set_current_state(istate_lz)
109
+ use mod_general, only: natom
110
+ use mod_sh, only: en_array, tocalc, nacx, nacy, nacz
111
+
112
+ nstate = nstate_lz ! Needed in init_terash
113
+ inac = 2 ! Turns off couplings calculation
114
+
115
+ ! Based on sh_init() routine, sharing most of the functions
116
+ ! TODO: Break dependence - separation of MPI interface needed
117
+ allocate (en_array(nstate_lz))
118
+ allocate (nacx(natom, nstate_lz, nstate_lz))
119
+ allocate (nacy(natom, nstate_lz, nstate_lz))
120
+ allocate (nacz(natom, nstate_lz, nstate_lz))
121
+ allocate (tocalc(nstate, nstate))
122
+ en_array = 0.0D0
123
+ tocalc = 0
124
+ tocalc(istate_lz, istate_lz) = 1
125
+ call set_current_state(istate_lz)
126
126
end subroutine lz_init_terash
127
127
128
128
! LZ singlets hop
@@ -176,9 +176,9 @@ subroutine lz_hop(x, y, z, vx, vy, vz, fxc, fyc, fzc, amt, dt, eclas, chpot)
176
176
end if
177
177
178
178
do ist1 = ibeg, iend
179
- if (ist1 == ist) cycle
179
+ if (ist1 == ist) cycle
180
180
! only closest states are considered for hopping
181
- if (ist1 > (ist + 1 ) .or. ist1 < (ist - 1 )) cycle
181
+ if (ist1 > (ist + 1 ) .or. ist1 < (ist - 1 )) cycle
182
182
183
183
do ihist = 1 , 4
184
184
en_diff(ihist) = abs (en_array_lz(ist, ihist) - en_array_lz(ist1, ihist))
@@ -204,31 +204,31 @@ subroutine lz_hop(x, y, z, vx, vy, vz, fxc, fyc, fzc, amt, dt, eclas, chpot)
204
204
! only for significant probabilities, for tiny probabilities this does not make sense
205
205
! e.g. for parallel states (the whole LZ does not make sense for this case)
206
206
if (prob(ist1) > 0.01 ) then
207
- ! Calculating backward second derivative formula. We are not using the newest energy
208
- ! but the previous three. Discontinuity would not affect this formula.
209
- second_der_back = ((en_diff(2 ) - 2 * en_diff(3 ) + en_diff(4 )) / dt** 2 )
210
- ! We have central and backward second derivative formulas and compare them.
211
- ! If the change is too large, we either almost hit the CI or we have discontinuity.
212
- der_check = abs ((second_der - second_der_back) / second_der)
213
- ! If they differ by more then 130%, we have aa unphysical change of curvature --> certain discontinuity.
214
- if (der_check > 1.3 ) then
215
- write (stdout,* ) " ERROR: Change of curvature --> discontinuity in PES!"
216
- write (stdout,* ) " Probability set to 0!"
217
- prob(ist1) = 0.0D0
218
- ! 30% threshold was set empirically and should capture most discontinuities
219
- ! yet it can also be a conical intersection. Thus, we just issue an warning
220
- ! and let the user to evaluate on his own.
221
- else if (der_check > 0.3 ) then
222
- write (stdout,* ) " WARNING: Possible discontinuity in PES! Check PES.dat!"
223
- end if
224
- end if
207
+ ! Calculating backward second derivative formula. We are not using the newest energy
208
+ ! but the previous three. Discontinuity would not affect this formula.
209
+ second_der_back = ((en_diff(2 ) - 2 * en_diff(3 ) + en_diff(4 )) / dt** 2 )
210
+ ! We have central and backward second derivative formulas and compare them.
211
+ ! If the change is too large, we either almost hit the CI or we have discontinuity.
212
+ der_check = abs ((second_der - second_der_back) / second_der)
213
+ ! If they differ by more then 130%, we have aa unphysical change of curvature --> certain discontinuity.
214
+ if (der_check > 1.3 ) then
215
+ write (stdout, * ) " ERROR: Change of curvature --> discontinuity in PES!"
216
+ write (stdout, * ) " Probability set to 0!"
217
+ prob(ist1) = 0.0D0
218
+ ! 30% threshold was set empirically and should capture most discontinuities
219
+ ! yet it can also be a conical intersection. Thus, we just issue an warning
220
+ ! and let the user to evaluate on his own.
221
+ else if (der_check > 0.3 ) then
222
+ write (stdout, * ) " WARNING: Possible discontinuity in PES! Check PES.dat!"
223
+ end if
224
+ end if
225
225
226
226
end if
227
227
end do
228
228
229
229
! LZ warning
230
230
if (sum (prob) > 1 ) then
231
- write (stdout, * ) " WARNING: Sum of hopping probabilities > 1. Breakdown of LZ assumptions"
231
+ write (stdout, * ) " WARNING: Sum of hopping probabilities > 1. Breakdown of LZ assumptions"
232
232
end if
233
233
234
234
! Hop?
0 commit comments