Add 'ignore_state' option to Launda-Zener hopping (#185)

* Tweaks to test.sh
* Move stuff into `mod_lz` for more privacy

src/abin.F90

@@ -25,7 +25,7 @@ program abin
       & nwrite, nstep, it, inormalmodes, istage, irest, STOP_SIMULATION
    use mod_init, only: init
    use mod_sh, only: surfacehop, sh_init, get_nacm, move_vars
-   use mod_lz, only: lz_hop, en_array_lz, lz_rewind, nsinglet_lz, ntriplet_lz
+   use mod_lz, only: lz_hop, en_array_lz, lz_rewind
    use mod_kinetic, only: temperature
    use mod_utils, only: del_file, archive_file
    use mod_transform, only: initialize_pi_transforms, &
@@ -68,12 +68,6 @@ program abin
    write (stdout, '(A)') 'Job started at: '//trim(get_formatted_date_and_time(time_start))
    write (stdout, *) ''
 
-   ! LZ warning for too many states
-   if (ipimd == 5 .and. (nsinglet_lz > 2 .or. ntriplet_lz > 2)) then
-      write (*, *) 'WARNING: LZ was derived for a two-state problem. More states might cause unphysical behavior.'
-      write (stdout, *) ''
-   end if
-
    ! Transform coordinates and velocities for Path Integral MD
    ! (staging or normal modes)
    if (istage == 1 .or. inormalmodes > 0) then

src/force_abin.F90

@@ -68,31 +68,6 @@ subroutine write_sh_data(nstate, tocalc)
       close (u)
    end subroutine write_sh_data
 
-   subroutine write_lz_data(nstate, nsinglet, ntriplet, tocalc)
-      integer, intent(in) :: nstate
-      integer, intent(in) :: nsinglet, ntriplet
-      integer, dimension(:), intent(in) :: tocalc
-      integer :: ist
-      integer :: u
-
-      if (nstate /= nsinglet + ntriplet) then
-         call fatal_error(__FILE__, __LINE__, &
-            & 'LZ: nstate /= nsinglet + ntriplet')
-      end if
-
-      open (newunit=u, file='state.dat')
-      write (u, '(I0)') nstate
-
-      ! Print for which state we need gradient
-      ! First we have singlets, then triplets
-      do ist = 1, nstate
-         if (tocalc(ist) == 1) write (u, '(I0)') ist
-      end do
-      write (u, '(I0)') nsinglet
-      write (u, '(I0)') ntriplet
-      close (u)
-   end subroutine write_lz_data
-
    function get_shellscript(potential) result(shellscript)
       use mod_utils, only: toupper
       character(len=*), intent(in) :: potential
@@ -231,7 +206,7 @@ subroutine force_abin(x, y, z, fx, fy, fz, eclas, chpot, walkmax)
       use mod_system, only: names
       use mod_sh_integ, only: nstate
       use mod_sh, only: tocalc, en_array, istate
-      use mod_lz, only: nstate_lz, tocalc_lz, en_array_lz, istate_lz, nsinglet_lz, ntriplet_lz
+      use mod_lz, only: nstate_lz, tocalc_lz, en_array_lz, istate_lz, write_lz_data
       use mod_qmmm, only: natqm
       use mod_utils, only: toupper, append_rank
       implicit none
@@ -276,7 +251,7 @@ subroutine force_abin(x, y, z, fx, fy, fz, eclas, chpot, walkmax)
 
             ! Landau-Zener
             if (ipimd == 5) then
-               call write_lz_data(nstate_lz, nsinglet_lz, ntriplet_lz, tocalc_lz)
+               call write_lz_data()
             end if
 
             ! Call the external program

src/init.F90

@@ -43,7 +43,7 @@ subroutine init(dt)
       use mod_potentials_sh
       use mod_sh_integ, only: phase, nstate
       use mod_sh, only: read_sh_input, print_sh_input
-      use mod_lz
+      use mod_lz, only: read_lz_input, print_lz_input, lz_init
       use mod_qmmm, only: natqm, natmm
       use mod_force_mm, only: initialize_mm
       use mod_force_h2o, only: initialize_h2o_pot, h2opot
@@ -120,7 +120,6 @@ subroutine init(dt)
       ! - remd:    parameters for Replica Exchange MD
       ! - sh:      parameters for Surface Hopping, moved to mod_sh module
       ! - system:  system-specific parameters for model potentials, masses, SHAKE constraints...
-      ! - lz:      parameters for Landau-Zener excited state dynamics.
       ! - qmmm:    parameters for internal QMMM (not really tested).
       !
       ! All namelists need to be in a single input file, and the code
@@ -142,8 +141,6 @@ subroutine init(dt)
          k, r0, kx, ky, kz, r0_morse, d0_morse, k_morse, D0_dw, lambda_dw, k_dw, r0_dw, &
          Nshake, ishake1, ishake2, shake_tol, potential_file
 
-      namelist /lz/ initstate_lz, nstate_lz, nsinglet_lz, ntriplet_lz, deltaE_lz, energydifthr_lz
-
       namelist /qmmm/ natqm, natmm, q, LJ_rmin, LJ_eps, mm_types
 
       chcoords = 'mini.dat'
@@ -370,9 +367,8 @@ subroutine init(dt)
       end if
 
       if (ipimd == 5) then
-         read (uinput, lz)
-         rewind (uinput)
-         call lz_init(pot) !Init arrays for possible restart
+         call read_lz_input(uinput)
+         call lz_init(pot)
       end if
 
       if (iremd == 1) then
@@ -630,10 +626,6 @@ subroutine check_inputsanity()
             write (*, *) 'Invalid ipimd value'
             error = 1
          end if
-         if (ipimd == 5 .and. pot == '_tera_' .and. ntriplet_lz > 0) then
-            write (*, *) 'ERROR: Landau-Zener with Singlet-Triplet transitions not implemented with TeraChem over MPI.'
-            error = 1
-         end if
          if (ipimd == 5 .and. nwalk /= 1) then
             write (*, *) 'ERROR: LZ not implemented with multiple walkers.'
             error = 1
@@ -794,7 +786,7 @@ subroutine print_simulation_parameters()
             write (stdout, *)
          end if
          if (ipimd == 5) then
-            write (stdout, nml=lz, delim='APOSTROPHE')
+            call print_lz_input()
             write (stdout, *)
          end if
          if (iqmmm > 0 .or. pot == '_mm_') then

src/landau_zener.F90

@@ -13,22 +13,27 @@ module mod_lz
    use mod_const, only: DP
    use mod_error, only: fatal_error
    use mod_files, only: stdout, stderr
-   ! TODO: Break coupling between LZ and SH modules
-   ! The following are needed for TERA-MPI interface
-   use mod_sh, only: set_current_state, inac
-   use mod_sh_integ, only: nstate
+   use mod_sh, only: set_current_state
    implicit none
    private
-   public :: lz_init, lz_hop, lz_rewind, lz_restin, lz_restout, lz_finalize !Routines
-   public :: initstate_lz, nstate_lz, nsinglet_lz, ntriplet_lz, deltaE_lz, energydifthr_lz !User defined variables, [deltaE_lz] = eV
-   public :: en_array_lz, tocalc_lz, istate_lz !Routine variables
-   !Caveat: Every time we call force_clas en_array_lz is updated
-
-   real(DP) :: deltaE_lz = 1.0D0 !Energy difference, up to which we consider possibility of LZ hops
-   real(DP) :: energydifthr_lz = 0.5D0 !Maximum energy difference (eV) between two consecutive steps, used in check_energy_lz()
+   ! Public subroutines
+   public :: lz_init, print_lz_input, read_lz_input, write_lz_data
+   public :: lz_hop, lz_rewind, lz_restin, lz_restout, lz_finalize
+   ! Other variables that must be public (but perhaps shouldn't be!)
+   ! Caveat: Every time we call force_clas en_array_lz is updated
+   public :: nstate_lz
+   public :: en_array_lz, tocalc_lz, istate_lz
+
+   ! Energy difference (eV), up to which we consider possibility of LZ hops
+   real(DP) :: deltaE_lz = 1.0D0
+   ! Maximum energy difference (eV) between two consecutive steps, used in check_energy_lz()
+   real(DP) :: energydifthr_lz = 0.5D0
    ! Initial electronic state
    integer :: initstate_lz = 1
    integer :: nstate_lz, nsinglet_lz = 0, ntriplet_lz = 0
+   ! Do not consider hopping into this state
+   ! Same meaning as in mod_sh
+   integer, public :: ignore_state = 0
 
    ! Module variables
    integer :: istate_lz
@@ -38,10 +43,24 @@ module mod_lz
    real(DP), allocatable, dimension(:, :) :: px_temp, py_temp, pz_temp
    real(DP), allocatable, dimension(:, :) :: x_prev, y_prev, z_prev
    real(DP), allocatable, dimension(:, :) :: vx_prev, vy_prev, vz_prev
+
+   namelist /lz/ initstate_lz, nstate_lz, nsinglet_lz, ntriplet_lz, ignore_state, deltaE_lz, energydifthr_lz
    save
 
 contains
 
+   subroutine read_lz_input(param_unit)
+      integer, intent(in) :: param_unit
+
+      rewind (param_unit)
+      read (param_unit, lz)
+      rewind (param_unit)
+   end subroutine read_lz_input
+
+   subroutine print_lz_input()
+      write (stdout, nml=lz, delim='APOSTROPHE')
+   end subroutine print_lz_input
+
    subroutine check_lz_parameters()
       use mod_utils, only: int_positive, int_nonnegative, real_nonnegative, real_positive
 
@@ -60,8 +79,28 @@ subroutine check_lz_parameters()
             & '(LZ): Sum of singlet and triplet states must give total number of states.')
       end if
 
+      if (ignore_state == initstate_lz) then
+         call fatal_error(__FILE__, __LINE__, &
+            & 'ignore_state == initstate_lz, cannot start simulation on ignored state')
+      end if
+
+      if (ignore_state > nstate_lz) then
+         call fatal_error(__FILE__, __LINE__, &
+            & '(LZ): ignore_state > nstate_lz')
+      end if
+
+      if (ignore_state /= 0) then
+         write (stdout, '(A,I0,A)') 'Ignoring state ', ignore_state, ' for hopping'
+      end if
+
+      if (nsinglet_lz > 2 .or. ntriplet_lz > 2) then
+         write (*, *) 'WARNING: LZ was derived for a two-state problem. More states might cause unphysical behavior.'
+         write (stdout, *) ''
+      end if
+
       call int_positive(initstate_lz, 'initstate_lz')
       call int_positive(nstate_lz, 'nstate_lz')
+      call int_nonnegative(ignore_state, 'ignore_state')
       call int_nonnegative(nsinglet_lz, 'nsinglet_lz')
       call int_nonnegative(ntriplet_lz, 'ntriplet_lz')
       call real_positive(deltaE_lz, 'deltaE_lz')
@@ -76,7 +115,7 @@ subroutine lz_init(pot)
 
       call check_lz_parameters()
 
-      !Initial state
+      ! Initial state
       istate_lz = initstate_lz
 
       allocate (tocalc_lz(nstate_lz))
@@ -105,9 +144,17 @@ subroutine lz_init(pot)
 
    end subroutine lz_init
 
+   ! TODO: Break coupling between LZ and SH modules
+   ! The following function is needed for TERA-MPI interface
    subroutine lz_init_terash()
       use mod_general, only: natom
-      use mod_sh, only: en_array, tocalc, nacx, nacy, nacz
+      use mod_sh, only: inac, en_array, tocalc, nacx, nacy, nacz
+      use mod_sh_integ, only: nstate
+
+      if (ntriplet_lz > 0) then
+         call fatal_error(__FILE__, __LINE__, &
+            & 'Landau-Zener with Singlet-Triplet transitions not implemented with TeraChem interface')
+      end if
 
       nstate = nstate_lz !Needed in init_terash
       inac = 2 !Turns off couplings calculation
@@ -125,6 +172,23 @@ subroutine lz_init_terash()
       call set_current_state(istate_lz)
    end subroutine lz_init_terash
 
+   ! This routine is used in the file interface defined in force_abin.F90
+   subroutine write_lz_data()
+      integer :: ist, u
+
+      open (newunit=u, file='state.dat')
+      write (u, '(I0)') nstate_lz
+
+      ! Print for which state we need gradient
+      ! First we have singlets, then triplets
+      do ist = 1, nstate_lz
+         if (tocalc_lz(ist) == 1) write (u, '(I0)') ist
+      end do
+      write (u, '(I0)') nsinglet_lz
+      write (u, '(I0)') ntriplet_lz
+      close (u)
+   end subroutine write_lz_data
+
    !LZ singlets hop
    subroutine lz_hop(x, y, z, vx, vy, vz, fxc, fyc, fzc, amt, dt, eclas, chpot)
       use mod_const, only: ANG, AUTOFS, PI, AUTOEV, AUTOCM
@@ -177,7 +241,8 @@ subroutine lz_hop(x, y, z, vx, vy, vz, fxc, fyc, fzc, amt, dt, eclas, chpot)
 
       do ist1 = ibeg, iend
          if (ist1 == ist) cycle
-         ! only closest states are considered for hopping
+         if (ist1 == ignore_state) cycle
+         ! only closest upper and lower states are considered for hopping
          if (ist1 > (ist + 1) .or. ist1 < (ist - 1)) cycle
 
          do ihist = 1, 4

src/surfacehop.F90

@@ -89,7 +89,10 @@ module mod_sh
    integer, allocatable :: tocalc(:, :)
    ! Current electronic state
    integer, public, protected :: istate
-   ! for ethylene 2-state SA3 dynamics
+   ! Do not consider hopping into this state
+   ! Useful e.g. for ethylene 2-state SA3 dynamics
+   ! This is a bit of an ugly hack, it would be more general to have an array
+   ! 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, &
@@ -118,7 +121,7 @@ subroutine sh_init(x, y, z, vx, vy, vz)
       ! TODO: Write a print_sh_header()
       ! with all major parameters included and explained
       if (ignore_state /= 0) then
-         write (stdout, '(A,I0)') 'Ignoring state ', ignore_state
+         write (stdout, '(A,I0,A)') 'Ignoring state ', ignore_state, ' for hopping'
       end if
 
       ! Determining the initial state
@@ -253,6 +256,10 @@ subroutine check_sh_parameters()
          write (stderr, '(A)') 'I cannot start on an ignored state'
          error = .true.
       end if
+      if (ignore_state > nstate) then
+         write (stderr, '(A)') 'Ignored state > number of computed states.'
+         error = .true.
+      end if
 
       call int_switch(nohop, 'nohop')
       call int_switch(adjmom, 'adjmom')

tests/test.sh

@@ -227,27 +227,31 @@ do
       continue
    fi
 
+   # Don't exit on errors in the rest of the script, since abin invocation can fail
+   # and this makes the code simpler.
+   set +e
    # For special cases such as REMD, we need a more complicated test setup.
    # If a file 'test.sh' is present in the test directory we will use it.
-   if [[ -f "test.sh" ]];then
+   # Otherwise, execute abin directly.
+   if [[ -f "test.sh" ]]; then
 
       # Redirection to dev/null apparently needed for CP2K tests.
       # Otherwise, STDIN is screwed up. I have no idea why.
       # http://stackoverflow.com/questions/1304600/read-error-0-resource-temporarily-unavailable
       # TODO: Figure out a different solution
       #./test.sh $ABINEXE 2> /dev/null
-      ./test.sh $ABINEXE || true
+      ./test.sh $ABINEXE
 
    else
       if [[ -f "velocities.in" ]];then
-         $ABINEXE -v "velocities.in" > $ABINOUT 2>&1 || true
+         $ABINEXE -v "velocities.in" > $ABINOUT 2>&1
       else
-         $ABINEXE > $ABINOUT 2>&1 || true
+         $ABINEXE > $ABINOUT 2>&1
       fi
 
-      #for testing restart
-      if [[ -e input.in2 ]];then
-         $ABINEXE -i input.in2 >> $ABINOUT 2>&1 || true
+      # for testing restart, only execute if previous run did not fail
+      if [[ -f input.in2 && $? -eq 0 ]]; then
+         $ABINEXE -i input.in2 >> $ABINOUT 2>&1
       fi
    fi
 
@@ -257,16 +261,11 @@ do
 
    else
 
-      # Since we're running in the -e mode,
-      # we need to "hide" this possibly failing command
-      # https://stackoverflow.com/a/11231970/3682277
-      current_error=0
-      diff_files || current_error=$?
-      if [[ $current_error -ne 0 ]];then
+      if diff_files; then
+        echo -e "\033[0;32mPASSED\033[0m"
+      else
         global_error=1
         echo -e "$dir \033[0;31mFAILED\033[0m"
-      else
-        echo -e "\033[0;32mPASSED\033[0m"
       fi
    fi