Added energy history tracking (analogously to LZSH)

src/force_abin.F90

@@ -230,7 +230,7 @@ subroutine force_abin(x, y, z, fx, fy, fz, eclas, chpot, walkmax)
       use mod_general, only: ipimd, iqmmm, it
       use mod_system, only: names
       use mod_sh_integ, only: nstate
-      use mod_sh, only: tocalc, en_array, istate
+      use mod_sh, only: tocalc, en_array, istate, inac, en_hist_array
       use mod_lz, only: nstate_lz, tocalc_lz, en_array_lz, istate_lz, nsinglet_lz, ntriplet_lz
       use mod_qmmm, only: natqm
       use mod_utils, only: toupper, append_rank
@@ -297,6 +297,15 @@ subroutine force_abin(x, y, z, fx, fy, fz, eclas, chpot, walkmax)
                if (abort) cycle
                eclas = en_array(istate)
 
+               ! Store the energy history for Baeck-An couplings
+               if (inac == 1) then
+                  ! Move old energies by 1 and storing the new energy
+                  en_hist_array(:, 4) = en_hist_array(:, 3)
+                  en_hist_array(:, 3) = en_hist_array(:, 2)
+                  en_hist_array(:, 2) = en_hist_array(:, 1)
+                  en_hist_array(:, 1) = en_array(:)
+               end if
+
             else if (ipimd == 5) then
 
                ! Move old energies by 1

src/force_terash.F90

@@ -48,7 +48,7 @@ subroutine receive_terash(fx, fy, fz, eclas, tc_comm)
       use mod_qmmm, only: natqm
       use mod_io, only: print_charges, print_dipoles, print_transdipoles
       use mod_sh_integ, only: nstate
-      use mod_sh, only: check_CIVector, en_array, istate, nacx, nacy, nacz
+      use mod_sh, only: check_CIVector, en_array, istate, nacx, nacy, nacz, inac, en_hist_array
       use mod_lz, only: en_array_lz
       use mpi
       real(DP), intent(inout) :: fx(:, :), fy(:, :), fz(:, :)
@@ -80,6 +80,14 @@ subroutine receive_terash(fx, fy, fz, eclas, tc_comm)
          en_array_lz(:, 2) = en_array_lz(:, 1); 
          !Store the new one
          en_array_lz(:, 1) = en_array(:)
+      ! SH Baeck-An arrays
+      else if ((ipimd == 2) .and. (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)
+         !Store the new one
+         en_hist_array(:, 1) = en_array(:)
       end if
 
       if (idebug > 0) write (*, '(a)') 'Receiving transition dipoles from TC.'

src/potentials_sh.F90

@@ -71,7 +71,7 @@ subroutine nai_init(natom, nwalk, ipimd, nstate)
    end subroutine nai_init
 
    subroutine force_nai(x, y, z, fx, fy, fz, eclas)
-      use mod_sh, only: en_array, istate, nacx, nacy, nacz
+      use mod_sh, only: en_array, istate, nacx, nacy, nacz, inac, en_hist_array
       use mod_const, only: ANG, AUTOEV
       real(DP), intent(in) :: x(:, :), y(:, :), z(:, :)
       real(DP), intent(out) :: fx(:, :), fy(:, :), fz(:, :)
@@ -135,6 +135,15 @@ subroutine force_nai(x, y, z, fx, fy, fz, eclas)
       ! saving electronic energies for SH
       en_array(1) = E1
       en_array(2) = E2
+
+      ! store the energy history for Baeck-An couplings
+      if (inac == 1) then
+         ! Move old energies by 1 and storing the new energy
+         en_hist_array(:, 4) = en_hist_array(:, 3)
+         en_hist_array(:, 3) = en_hist_array(:, 2)
+         en_hist_array(:, 2) = en_hist_array(:, 1)
+         en_hist_array(:, 1) = en_array(:)
+      end if
 
       ! saving classical energy for Verlet
       eclas = en_array(istate)

src/surfacehop.F90

@@ -20,7 +20,7 @@ module mod_sh
 
    ! Variables holding SH state
    ! TODO: Make these protected and write Set methods
-   public :: en_array, tocalc
+   public :: en_array, en_hist_array, tocalc
    public :: nacx, nacy, nacz
 
    ! Initial electronic state
@@ -578,12 +578,19 @@ subroutine calc_nacm(pot, nac_accu)
       end if
    end subroutine calc_nacm
 
+   ! calculate Baeck-An couplings
+   subroutine calc_baeckan()
+
+
+
+   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
       real(DP), intent(in) :: vx(:, :), vy(:, :), vz(:, :)
       real(DP), intent(out) :: vx_old(:, :), vy_old(:, :), vz_old(:, :)
-      integer :: ist1, ist2, iat
+      integer :: ist1, ist2, iat, ihist
 
       do ist1 = 1, nstate
          en_array_old(ist1) = en_array(ist1)
@@ -648,7 +655,6 @@ subroutine surfacehop(x, y, z, vx, vy, vz, vx_old, vy_old, vz_old, dt, eclas)
 
       ! First, calculate NACME
       if (inac == 0) then ! Analytic ab initio couplings
-         write (stdout, '(A)') 'Analytic couplings calculated' !TODO JJ: remove later
          ! For TeraChem MPI / FMS interface, NAC are already computed!
          if (pot /= '_tera_' .and. pot /= '_nai_') then
             nacx = 0.0D0
@@ -661,7 +667,11 @@ subroutine surfacehop(x, y, z, vx, vy, vz, vx_old, vy_old, vz_old, dt, eclas)
          ! 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
+         ! TODO JJ: once I make a function to calculate Baeck-An couplings, I need to add if condition that it calculates only if
+         ! step is bigger than xx
          write (stdout, '(A)') 'Baeck-An couplings calculated' !TODO JJ: remove later
+         write (stdout, *) "when coups calc", en_hist_array(1, :)  !TODO JJ: remove later
+         write (stdout, *) "Current energy:", en_array(1) !TODO JJ: remove later
       end if
 
       ! smaller time step
@@ -676,6 +686,7 @@ subroutine surfacehop(x, y, z, vx, vy, vz, vx_old, vy_old, vz_old, dt, eclas)
          ! pop0 is later used for Tully's fewest switches
          pop0 = get_elpop(ist)
 
+         ! TODO JJ: I guess this interpolation will go only for NAC.
          ! 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, &