Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Baeck-An couplings for SH #186

Merged
merged 26 commits into from
Dec 3, 2024
Merged

Baeck-An couplings for SH #186

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
26 commits
Select commit Hold shift + click to select a range
2374cb3
Preparing SH for Baeck-An couplings
JanosJiri Nov 7, 2024
f0fb571
Using inac=1 instead of beackan
JanosJiri Nov 7, 2024
6d3ae69
Using string for couplings input
JanosJiri Nov 7, 2024
d34d1d9
Adding the new keyword to tests
JanosJiri Nov 7, 2024
f0a3947
Adding two missing MPICH tests
JanosJiri Nov 7, 2024
d2ade90
Changing adjmom to mom_adjust keyword that takes string as input
JanosJiri Nov 8, 2024
918f683
Added energy history tracking (analogously to LZSH)
JanosJiri Nov 8, 2024
1f4ac58
Better way of energy history
JanosJiri Nov 8, 2024
172acc3
Baeck-An couplings calculated
JanosJiri Nov 11, 2024
98e235c
Clean auxilliary plotting
JanosJiri Nov 11, 2024
9e4d031
Fix for couplings = 'none' in TERAPI-SH-S0
JanosJiri Nov 11, 2024
98fac22
Writing nacme_all.dat only for inac==0.or.2
JanosJiri Nov 11, 2024
48e4fa4
Test for Baeck-An couplings
JanosJiri Nov 11, 2024
de550d1
Formatting and trying to solve unknown problems GFrotran 7
JanosJiri Nov 11, 2024
55dec86
New reference for the test, still not mending GFortran 7
JanosJiri Nov 11, 2024
e608c59
v_int missing in Baeck-An surface hop for decoherence correction
JanosJiri Nov 12, 2024
3f53c9e
Not a single commit where I would not forget to autoformat the code...
JanosJiri Nov 12, 2024
1ecb7ae
Correcting NaI potential with correct factor.
JanosJiri Nov 14, 2024
3882145
Fixing the NaI test
JanosJiri Nov 14, 2024
3f7e0dd
Creating subroutine sh_files_init()
JanosJiri Nov 18, 2024
1870900
Autoformatting..
JanosJiri Nov 18, 2024
4aa27c8
Merge remote-tracking branch 'origin/master' into baeck-an
JanosJiri Nov 18, 2024
4e3e36d
Merge branch 'master' into baeck-an
JanosJiri Nov 18, 2024
f1cc175
Changing mom_adjust to velocity_rescaling
JanosJiri Nov 30, 2024
adad854
Adding and modifying comments.
JanosJiri Dec 3, 2024
1cf13d1
open nacme_all.dat only if couplings='analytic'
JanosJiri Dec 3, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 2 additions & 0 deletions sample_inputs/input.in.sh
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -22,6 +22,8 @@ inose=0, ! NVE ensemble (thermostat turned OFF)
&sh
istate_init=2, ! initial electronic state (1 is ground state)
nstate=3, ! number of electronic states
couplings='analytic', ! non-adiabatic coupling terms 'analytic', 'baeck-an', 'none'
velocity_rescaling='nac_then_velocity' ! momentum adjustment along either 'nac_then_velocity' or 'velocity'
deltaE=2.0, ! maximum energy difference (eV) between states for which we calculate NA coupling
PopThr=0.001, ! minimum population of either state, for which we compute NA coupling
EnergyDifThr=0.50, ! maximum energy difference between two consecutive steps
Expand Down
64 changes: 44 additions & 20 deletions src/files.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -170,26 +170,7 @@ subroutine files_init(isbc, phase, ndist, nang, ndih)
end if

if (ipimd == 2) then
open (UPOP, file=chfiles(UPOP), access=chaccess, action='write')
open (UPROB, file=chfiles(UPROB), access=chaccess, action='write')
open (UPES, file=chfiles(UPES), access=chaccess, action='write')
open (UNACME, file=chfiles(UNACME), access=chaccess, action='write')
open (UDOTPROD, file=chfiles(UDOTPROD), access=chaccess, action='write')

if (idebug > 1) then
open (UBKL, file=chfiles(UBKL), access=chaccess, action='write')
open (UWFCOEF, file=chfiles(UWFCOEF), access=chaccess, action='write', recl=250)
if (phase == 1) then
open (UPHASE, file=chfiles(UPHASE), access=chaccess, action='write')
end if
end if

if (pot == '_tera_') then
open (UCHARGES, file=chfiles(UCHARGES), access=chaccess, action='write')
open (UDOTPRODCI, file=chfiles(UDOTPRODCI), access=chaccess, action='write')
open (UDIP, file=chfiles(UDIP), access=chaccess, action='write')
open (UTDIP, file=chfiles(UTDIP), access=chaccess, action='write')
end if
call sh_files_init(phase, chaccess)
end if

if (ipimd /= 2 .and. pot == '_tera_') then
Expand Down Expand Up @@ -230,6 +211,49 @@ subroutine files_init(isbc, phase, ndist, nang, ndih)

end subroutine files_init

subroutine sh_files_init(phase, chaccess)
use mod_general, only: idebug, pot
integer, intent(in) :: phase
character(len=10), intent(in) :: chaccess

open (UPOP, file=chfiles(UPOP), access=chaccess, action='write')
open (UPROB, file=chfiles(UPROB), access=chaccess, action='write')
open (UPES, file=chfiles(UPES), access=chaccess, action='write')
open (UDOTPROD, file=chfiles(UDOTPROD), access=chaccess, action='write')

if (idebug > 1) then
open (UBKL, file=chfiles(UBKL), access=chaccess, action='write')
open (UWFCOEF, file=chfiles(UWFCOEF), access=chaccess, action='write', recl=250)
if (phase == 1) then
open (UPHASE, file=chfiles(UPHASE), access=chaccess, action='write')
end if
end if

if (pot == '_tera_') then
open (UCHARGES, file=chfiles(UCHARGES), access=chaccess, action='write')
open (UDOTPRODCI, file=chfiles(UDOTPRODCI), access=chaccess, action='write')
open (UDIP, file=chfiles(UDIP), access=chaccess, action='write')
open (UTDIP, file=chfiles(UTDIP), access=chaccess, action='write')
end if

end subroutine sh_files_init

subroutine nacmefile_init()
use mod_general, only: irest
character(len=10) :: chaccess

! Here we ensure, that previous files are deleted
if (irest == 0) then
chaccess = 'SEQUENTIAL'
else
chaccess = 'APPEND'
end if

! open file
open (UNACME, file=chfiles(UNACME), access=chaccess, action='write')

end subroutine nacmefile_init

subroutine print_file_headers()
use mod_general, only: ipimd, natom
use mod_system, only: names
Expand Down
1 change: 0 additions & 1 deletion src/potentials_sh.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -125,7 +125,6 @@ subroutine force_nai(x, y, z, fx, fy, fz, eclas)
E2 = (VA + VX) / 2.0D0 + dsqrt((VA - VX)**2.0D0 + 4.0D0 * VXA**2.0D0) / 2.0D0
! nonadiabatic coupling vector in the reduced system
d12 = -(VXA * (dVA - dVX) + (-VA + VX) * dVXA) / (VA**2.0D0 - 2.0D0 * VA * VX + VX**2.0D0 + 4.0D0 * VXA**2.0D0)
d12 = d12 / ANG ! one more conversion necessary
JanosJiri marked this conversation as resolved.
Show resolved Hide resolved
! derivatives of energies
dE1 = (dVA + dVX) / 2.0D0 - (2.0D0 * (VA - VX) * (dVA - dVX) + 8.0D0 * VXA * dVXA) / &
(4.0D0 * dsqrt((VA - VX)**2.0D0 + 4.0D0 * VXA**2.0D0))
Expand Down
193 changes: 165 additions & 28 deletions src/surfacehop.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -33,20 +33,23 @@
integer :: substep = 100

! Controls calculations of Non-adiabatic Couplings (NAC)
! 0 - Analytical NAC
! 1 - Numerical Hammers-Schffer-Tully model (currently not implemented)
! 2 - Do not compute couplings
integer :: inac = 0
! couplings = 'analytic' (inac=0) - Analytical NAC (default)
! couplings = 'baeck-an' (inac=1) - Baeck-An couplings
! couplings = 'none' (inac=2) - Do not compute couplings
integer :: inac = 0 ! for working within the code
character(len=50) :: couplings = 'analytic' ! for reading the input file
! energy history array and time-derivate couplings (sigma_ba) necessary for Beack-An couplings
real(DP), allocatable :: en_hist_array(:, :), sigma_ba(:, :) ! sigma_ba is actually the dotproduct

! 1 - Turn OFF hopping
integer :: nohop = 0

! How to adjust velocity after hop:
! 0 - Adjust velocity along the NAC vector (default)
! 1 - Simple velocity rescale
! NOTE: Simple v-rescale is invoked as a fallback
! if there is not enough momentum along the NAC vector.
integer :: adjmom = 0
! velocity_rescaling = 'nac_then_velocity' (adjmom=0) - Adjust velocity along the NAC vector, if not possible,
! try the velocity vector (default)
! velocity_rescaling = 'velocity' (adjmom=1) - Rescale along the velocity vector
integer :: adjmom = 0 ! for working within the code
character(len=50) :: velocity_rescaling = 'nac_then_velocity' ! for reading the input file
! 1 - Reverse momentum direction after frustrated hop
integer :: revmom = 0

Expand Down Expand Up @@ -95,8 +98,8 @@
! of states that are calculated but ignored.
integer :: ignore_state = 0

namelist /sh/ istate_init, nstate, substep, deltae, integ, inac, nohop, phase, decoh_alpha, popthr, ignore_state, &
nac_accu1, nac_accu2, popsumthr, energydifthr, energydriftthr, adjmom, revmom, &
namelist /sh/ istate_init, nstate, substep, deltae, integ, couplings, nohop, phase, decoh_alpha, popthr, ignore_state, &
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks for introducing the new keywords!

I think we should keep the old variables here for backwards compatibility (so that we don't break existing input files).

Can you also update sample_inputs/input.in.sh

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Can you also update sample_inputs/input.in.sh

Done!

About keeping inac and adjmom, if I keep them in the namelist, someone can set them but they will still be reassigned by couplings and mom_adjust, which have their default values. I'm afraid this would cause more problems. I've never seen input where these two were used so I wouldn't be worried about backwards compatibility of things done in Prague. When you changed the name for decoh_alpha, it was not backwards compatible and still people quickly adapted (again at least here in Prague). So I would stick just to the new keywords, yet I leave the executive decision up to you. :)

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What you could do is to set the couplings keyword to some default value, and if the user doesn't change it in the input, don't overwrite the old keywords. I think I already do that for ipimd and mdtype. But I don't feel too strongly.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Just coming back to this, I agree it's okay to not support the old keywords and not complicate the code. We need to also document all of this in the documentation, but that's a separate conversation. :-)

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Just a note for myself: I don't know why I put the sh namelist at the module level, since it is not public, it would be better to put it into sh_init function, to better separate variables that are only needed there (such as the new couplings keyword.

nac_accu1, nac_accu2, popsumthr, energydifthr, energydriftthr, velocity_rescaling, revmom, &
dE_S0S1_thr, correct_decoherence
save

Expand All @@ -108,6 +111,7 @@
use mod_general, only: irest, natom, pot
use mod_interfaces, only: force_clas
use mod_kinetic, only: ekin_v
use mod_files, only: nacmefile_init
real(DP), intent(inout) :: x(:, :), y(:, :), z(:, :)
real(DP), intent(in) :: vx(:, :), vy(:, :), vz(:, :)
real(DP) :: dum_fx(size(x, 1), size(x, 2))
Expand Down Expand Up @@ -153,6 +157,14 @@
en_array = 0.0D0
en_array_old = en_array

! Initialize the history array we use to calculate the Baeck-An couplings
if (inac == 1) then
allocate (en_hist_array(nstate, 4)) !last 3 energies (1: current, 2: n-1, 3: n-2, 4: n-3)
en_hist_array = 0.0D0
allocate (sigma_ba(nstate, nstate)) !this is equivalent to dotproduct, but I will need to store old and new values
sigma_ba = 0.0D0
end if

allocate (tocalc(nstate, nstate))
tocalc = 0
tocalc(istate, istate) = 1
Expand All @@ -165,6 +177,9 @@
dum_fx = 0.0D0; dum_fy = 0.0D0; dum_fz = 0.0D0
call force_clas(dum_fx, dum_fy, dum_fz, x, y, z, dum_eclas, pot)

! open nacme_all.dat if NACME is calculated
if (inac == 0) call nacmefile_init()

! When restarting, initial SH WF was already read from the restart file
if (irest == 0) then
call sh_set_initialwf(istate)
Expand Down Expand Up @@ -227,18 +242,45 @@
error = .true.
end if

if (inac > 2 .or. inac < 0) then
write (stderr, '(A)') 'Parameter "inac" must be 0, 1 or 2.'
! converting input 'couplings' into inac which is used in the code
select case (couplings)
case ('analytic')
inac = 0
write (stdout, '(A)') 'Using analytic ab initio couplings.'
JanosJiri marked this conversation as resolved.
Show resolved Hide resolved
case ('baeck-an')
inac = 1
write (stdout, '(A)') 'Using approximate Baeck-An couplings.'
case ('none')
inac = 2
write (stdout, '(A)') 'Ignoring nonadaiabatic couplings.'
case default
write (stderr, '(A)') 'Parameter "couplings" must be "analytic", "baeck-an" or "none".'

Check warning on line 257 in src/surfacehop.F90

View check run for this annotation

Codecov / codecov/patch

src/surfacehop.F90#L257 

Added line #L257 was not covered by tests
error = .true.
end if
end select

! converting input 'velocity_rescaling' into inac which is used in the code
select case (velocity_rescaling)
case ('nac_then_velocity')
adjmom = 0
write (stdout, '(A)') 'Rescaling velocity along the NAC vector after hop.'
write (stdout, '(A)') 'If there is not enough energy, try rescaling along the velocity vector.'
case ('velocity')
adjmom = 1
write (stdout, '(A)') 'Rescaling velocity along the momentum vector after hop.'
case default
write (stderr, '(A)') 'Parameter "velocity_rescaling" must be "nac_then_velocity" or "velocity".'

Check warning on line 271 in src/surfacehop.F90

View check run for this annotation

Codecov / codecov/patch

src/surfacehop.F90#L271 

Added line #L271 was not covered by tests
error = .true.
end select

if (adjmom == 0 .and. inac == 1) then
write (stderr, '(A)') 'Combination of adjmom=0 and inac=1 is not possible.'
write (stderr, '(A)') 'NAC vectors are not computed if inac=1.'
write (stderr, '(A)') 'Combination of velocity_rescaling="nac_then_velocity" and couplings="baeck-an" is not possible.'
write (stderr, '(A)') 'Velocity cannot be rescaled along NAC when using Baeck-An.'
write (stderr, '(A)') 'Change velocity_rescaling="velocity" to rescale along the velocity vector.'
error = .true.
end if

if (inac == 2 .and. nohop == 0) then
write (stdout, '(A)') 'WARNING: For simulations without couplings, inac=2, hopping probability cannot be determined.'
write (stdout, '(A)') 'WARNING: For simulations without couplings(="none") hopping probability cannot be determined.'
nohop = 1
end if

Expand Down Expand Up @@ -549,6 +591,37 @@
end if
end subroutine calc_nacm

! Calculate Baeck-An couplings
! Implementation was based on these two articles:
! original by Barbatti: https://doi.org/10.12688/openreseurope.13624.2
! another implementation by Truhlar: https://doi.org/10.1021/acs.jctc.1c01080
! In the numeric implementation, we follow Barbatti and use a higher order formula.
subroutine calc_baeckan(dt)
JanosJiri marked this conversation as resolved.
Show resolved Hide resolved
use mod_general, only: it
integer :: ist1, ist2
real(DP), intent(in) :: dt
real(DP) :: de(4), de2dt2, argument

sigma_ba = 0.0D0

! we don't have sufficient history until step 4
if (it < 4) return

do ist1 = 1, nstate
do ist2 = ist1 + 1, nstate
de = en_hist_array(ist2, :) - en_hist_array(ist1, :)
! Second derivative (de2dt2) comes from Eq. 16 from https://doi.org/10.12688/openreseurope.13624.2
de2dt2 = (2.0D0 * de(1) - 5.0D0 * de(2) + 4.0D0 * de(3) - de(4)) / dt**2
JanosJiri marked this conversation as resolved.
Show resolved Hide resolved
argument = de2dt2 / de(1)
if (argument > 0.0D0) then
sigma_ba(ist2, ist1) = dsqrt(argument) / 2.0D0
end if
sigma_ba(ist1, ist2) = -sigma_ba(ist2, ist1)
end do
end do

end subroutine calc_baeckan

! move arrays from new step to old step
subroutine move_vars(vx, vy, vz, vx_old, vy_old, vz_old)
use mod_general, only: natom
Expand All @@ -568,9 +641,21 @@
end do
end do
end if

end do

! Shift the energy history for Baeck-An couplings
if (inac == 1) then
! Move old energies by 1
en_hist_array(:, 4) = en_hist_array(:, 3)
en_hist_array(:, 3) = en_hist_array(:, 2)
en_hist_array(:, 2) = en_hist_array(:, 1)
! new energy is stored before the couplings are calculated
! I avoided doing the same as with LZSH energy history tracking because then I need to modify force_abin, force_terash and
! every potential in potentials_sh (and all possible future ones). This way it is kept private and does not depend on the
! way energies are calculated.
! See commit: https://github.com/PHOTOX/ABIN/pull/186/commits/918f6837a76ec0f41b575d3ca948253eed2f30cc
end if

vx_old = vx
vy_old = vy
vz_old = vz
Expand All @@ -594,7 +679,7 @@
real(DP), dimension(nstate) :: en_array_int, en_array_newint
real(DP), dimension(natom, nstate, nstate) :: nacx_int, nacy_int, nacz_int
real(DP), dimension(natom, nstate, nstate) :: nacx_newint, nacy_newint, nacz_newint
real(DP), dimension(nstate, nstate) :: dotproduct_int, dotproduct_newint
real(DP), dimension(nstate, nstate) :: dotproduct_int, dotproduct_newint, sigma_ba_old
! Switching probabilities
real(DP) :: t(nstate, nstate)
! Cumulative switching probabilities
Expand All @@ -618,7 +703,7 @@
t_tot = 1.0D0

! First, calculate NACME
if (inac == 0) then
if (inac == 0) then ! Analytic ab initio couplings
! For TeraChem MPI / FMS interface, NAC are already computed!
if (pot /= '_tera_' .and. pot /= '_nai_') then
nacx = 0.0D0
Expand All @@ -630,6 +715,11 @@
! TODO: Should we call this with TeraChem?
! I think TC already phases the couplings internally.
call phase_nacme(nacx_old, nacy_old, nacz_old, nacx, nacy, nacz)
else if (inac == 1) then ! Baeck-An couplings
! saving the current energy to the energy history (shifting was already done in previous step in move_vars)
en_hist_array(:, 1) = en_array(:)
sigma_ba_old = sigma_ba ! saving old sigma_ba
call calc_baeckan(dt)
end if

! smaller time step
Expand All @@ -645,15 +735,26 @@
pop0 = get_elpop(ist)

! INTERPOLATION
fr = real(itp, DP) / real(substep, DP)
call interpolate(vx, vy, vz, vx_old, vy_old, vz_old, vx_newint, vy_newint, vz_newint, &
nacx_newint, nacy_newint, nacz_newint, en_array_newint, &
dotproduct_newint, fr)
if ((inac == 0) .or. (inac == 2)) then
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It's not clear to me if we need to interpolate if inac==2? But better keep the existing behaviour for now.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't think it is necessary to interpolate but the original code was interpolating also for inac==2 so I kept it because I didn't want to mess up something.

fr = real(itp, DP) / real(substep, DP)
call interpolate(vx, vy, vz, vx_old, vy_old, vz_old, vx_newint, vy_newint, vz_newint, &
nacx_newint, nacy_newint, nacz_newint, en_array_newint, &
dotproduct_newint, fr)

fr = real(itp - 1, DP) / real(substep, DP)
call interpolate(vx, vy, vz, vx_old, vy_old, vz_old, vx_int, vy_int, vz_int, &
nacx_int, nacy_int, nacz_int, en_array_int, &
dotproduct_int, fr)
else if (inac == 1) then
fr = real(itp, DP) / real(substep, DP)
call interpolate_ba(vx, vy, vz, vx_old, vy_old, vz_old, vx_newint, vy_newint, vz_newint, &
en_array_newint, dotproduct_newint, sigma_ba, sigma_ba_old, fr)

fr = real(itp - 1, DP) / real(substep, DP)
call interpolate_ba(vx, vy, vz, vx_old, vy_old, vz_old, vx_int, vy_int, vz_int, &
en_array_int, dotproduct_int, sigma_ba, sigma_ba_old, fr)

fr = real(itp - 1, DP) / real(substep, DP)
call interpolate(vx, vy, vz, vx_old, vy_old, vz_old, vx_int, vy_int, vz_int, &
nacx_int, nacy_int, nacz_int, en_array_int, &
dotproduct_int, fr)
end if

! Integrate electronic wavefunction for one dtp time step
call sh_integrate_wf(en_array_int, en_array_newint, dotproduct_int, dotproduct_newint, dtp)
Expand Down Expand Up @@ -1012,6 +1113,42 @@

end subroutine interpolate

! interpolation of time-derivative coupling calculated via Baeck-An approximation
! this routine interpolates sigma_ba between integration steps
subroutine interpolate_ba(vx, vy, vz, vx_old, vy_old, vz_old, vx_int, vy_int, vz_int, &
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I am not stoked about yet another interpolation routine. At some point it would be good to refactor them into a more general routine that could then be called on various arrays. But not in this PR anyway.

en_array_int, dotproduct_int, sigma_ba, sigma_ba_old, fr)
use mod_general, only: natom
real(DP), intent(in) :: sigma_ba(:, :), sigma_ba_old(:, :)
real(DP), intent(in) :: vx(:, :), vy(:, :), vz(:, :) ! for velocity interpolation
real(DP), intent(in) :: vx_old(:, :), vy_old(:, :), vz_old(:, :)
real(DP), intent(out) :: dotproduct_int(:, :)
real(DP), intent(out) :: en_array_int(:)
real(DP), intent(out) :: vx_int(:, :), vy_int(:, :), vz_int(:, :) ! interpolated velocities
! How far are we interpolating?
real(DP), intent(in) :: fr
real(DP) :: frd
integer :: ist1, ist2, iw, iat !iteration counters

frd = 1.0D0 - fr

do ist1 = 1, nstate
en_array_int(ist1) = en_array(ist1) * fr + en_array_old(ist1) * frd
do ist2 = 1, nstate
! interpolating dot product
dotproduct_int(ist1, ist2) = sigma_ba(ist1, ist2) * fr + sigma_ba_old(ist1, ist2) * frd
end do
end do

! interpolating velocity which is necessary for Ekin in the decoherence correction
iw = 1
do iat = 1, natom
vx_int(iat, iw) = vx(iat, iw) * fr + vx_old(iat, iw) * frd
vy_int(iat, iw) = vy(iat, iw) * fr + vy_old(iat, iw) * frd
vz_int(iat, iw) = vz(iat, iw) * fr + vz_old(iat, iw) * frd
end do

end subroutine interpolate_ba

subroutine try_hop_simple_rescale(vx, vy, vz, instate, outstate, eclas)
use mod_general, only: pot
use mod_kinetic, only: ekin_v
Expand Down
Loading
Loading