Make sure the last step is written when using dE_S0S1 (#184)

* Remove ncalc parameter
* Use new global flag `STOP_SIMULATION`

docs/abin-v1.0-doc.tex

@@ -397,7 +397,6 @@ \subsection{\&general}
 \textit{nrest} & \textbf{1}   & Print restart file every \textit{nrest} step.   \\
 & &  Larger value is recomended, since restart files can be quite big. \\
 \textit{nwritev} & \textbf{0} & Print velocities every \textit{nwritev} step.   \\
-\colorbox{black!20}{\textit{ncalc}} & \textbf{1}   & Calculate observables (energies, temperature) every \textit{ncalc} step. 	\\
 
 \end{tabular}
 \newpage

src/abin.F90

@@ -22,7 +22,7 @@ program abin
    use mod_arrays
    use mod_files, only: stdout
    use mod_general, only: sim_time, pot, pot_ref, iremd, ipimd, &
-      & nwrite, nstep, ncalc, it, inormalmodes, istage, irest
+      & 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
@@ -140,9 +140,13 @@ program abin
       ! if it is present (we delete it below before we stop the program).
       call mpi_barrier_wrapper()
 
+      if (STOP_SIMULATION) exit
+
       inquire (file="EXIT", exist=file_exists)
       if (file_exists) then
          write (stdout, *) 'Found file EXIT. Writing restart file and exiting.'
+         ! TODO: Not sure why we are writing the restart file here, since it's being
+         ! written after the main loop as well.
          if (istage == 1) then
             call QtoX(vx, vy, vz, transxv, transyv, transzv)
             call QtoX(x, y, z, transx, transy, transz)
@@ -210,11 +214,6 @@ program abin
       ! In order to analyze the output, we have to perform the back transformation
       ! Transformed (cartesian) coordinates are stored in trans matrices.
 
-      ! Enter this section only every ncalc step
-      if (modulo(it, ncalc) /= 0) then
-         cycle
-      end if
-
       call temperature(px, py, pz, amt, eclas)
 
       if (istage == 1) then
@@ -247,10 +246,12 @@ program abin
    end do
 
    ! Write restart file at the end of a run
-   ! Because NCALC might be >1, we have to perform transformation to get the most
-   ! recent coordinates and velocities
+   ! TODO: This it variable manipulation is very brittle :-(
    it = it - 1
 
+   ! Because NCALC might be >1, we have to perform transformation to get the most
+   ! recent coordinates and velocities
+   ! TODO: ncalc parameter has been removed so this is probably not necessary anymore.
    if (istage == 1) then
       call QtoX(vx, vy, vz, transxv, transyv, transzv)
       call QtoX(x, y, z, transx, transy, transz)

src/ekin.F90

@@ -9,7 +9,7 @@ module mod_kinetic
 contains
    subroutine temperature(px, py, pz, amt, eclas)
       use mod_const, only: AUtoFS, AUtoK
-      use mod_general, only: it, ncalc, sim_time, natom, nwalk, nwrite, inormalmodes, ipimd
+      use mod_general, only: it, sim_time, natom, nwalk, nwrite, inormalmodes, ipimd
       use mod_system, only: dime, f, conatom
       use mod_files, only: UTEMPER, UENERGY
       use mod_nhc, only: inose, get_nhcham
@@ -21,9 +21,7 @@ subroutine temperature(px, py, pz, amt, eclas)
       real(DP), intent(in) :: amt(:, :)
       real(DP), intent(in) :: eclas
       real(DP) :: est_temp, temp1, ekin_mom
-      real(DP) :: it2
 
-      it2 = it / ncalc
       ekin_mom = 0.0D0
       do iw = 1, nwalk
          do iat = 1, natom
@@ -45,7 +43,7 @@ subroutine temperature(px, py, pz, amt, eclas)
 
          ! Temperature and conserved quantities of thermostats
          write (UTEMPER, '(F15.2,F10.2)', advance="no") sim_time * AUtoFS, est_temp * autok
-         if (ipimd /= 2) write (UTEMPER, '(F10.2)', advance="no") est_temp_cumul * autok / it2
+         if (ipimd /= 2) write (UTEMPER, '(F10.2)', advance="no") est_temp_cumul * autok / it
 
          if (inose == 1) then
             write (UTEMPER, '(E20.10)', advance="no") get_nhcham(natom, nwalk) + ekin_mom + eclas
@@ -59,7 +57,7 @@ subroutine temperature(px, py, pz, amt, eclas)
             ! Printing to file energies.dat
             write (UENERGY, '(F15.2,3E20.10)', advance='no') sim_time * AUtoFS, &
                                                            & eclas, ekin_mom, eclas + ekin_mom
-            if (ipimd == 0) write (UENERGY, '(E20.10)', advance="no") entot_cumul / it2
+            if (ipimd == 0) write (UENERGY, '(E20.10)', advance="no") entot_cumul / it
             write (UENERGY, *)
          end if
 

src/estimators.F90

@@ -13,7 +13,7 @@ module mod_estimators
    ! Expecting cartesian coordinates and forces!
    subroutine estimators(x, y, z, fxab, fyab, fzab, eclas)
       use mod_general, only: ipimd, pot, natom, nwalk, it, sim_time, &
-                           & ncalc, nwrite, inormalmodes, imini, ihess, icv
+                           & nwrite, inormalmodes, imini, ihess, icv
       use mod_nhc, only: temp, inose
       use mod_system, only: am, dime
       use mod_potentials, only: hessian_harmonic_oscillator, hessian_morse, hessian_harmonic_rotor
@@ -34,15 +34,14 @@ subroutine estimators(x, y, z, fxab, fyab, fzab, eclas)
          temp = temp / nwalk
       end if
 
-      ! We calculate all quantities only every ncalc steps
-      ! also we begin to accumulate energies only after first enmini steps to avoid
+      ! We begin to accumulate energies only after first enmini steps to avoid
       ! large initial oscilations
-      itnc = (it - enmini) / ncalc
+      itnc = it - enmini
       nf = dime * natom - nshake !degrees of freedom
 
       ! We begin to accumulate averages of heat capacities only after it > imini
       ! This auxiliary variable is for cumulative averaging if imini > 0
-      it2 = (it - imini) / ncalc
+      it2 = it - imini
 
       ! We cannot accumulate heat capacity without accumulating energy first
       if (enmini > imini .and. icv == 1) then

src/force_spline.F90

@@ -24,6 +24,8 @@ module mod_splined_grid
    real(DP) function potential_cubic_spline(x) result(y)
       real(DP), intent(in) :: x
       call fatal_error(__FILE__, __LINE__, "Spline potential not implemented")
+      ! To squash compiler warnings
+      y = x
    end function potential_cubic_spline
 
    subroutine force_splined_grid(x, fx, eclas, walkmax)

src/init.F90

@@ -128,7 +128,7 @@ subroutine init(dt)
       namelist /general/ pot, ipimd, mdtype, istage, inormalmodes, nwalk, nstep, icv, ihess, imini, nproc, iqmmm, &
          nwrite, nwritex, nwritev, nwritef, dt, irandom, nabin, irest, nrest, anal_ext, &
          isbc, rb_sbc, kb_sbc, gamm, gammthr, conatom, mpi_milisleep, narchive, xyz_units, &
-         dime, ncalc, idebug, enmini, rho, iknow, watpot, h2opot, iremd, iplumed, plumedfile, &
+         dime, idebug, enmini, rho, iknow, watpot, h2opot, iremd, iplumed, plumedfile, &
          en_restraint, en_diff, en_kk, restrain_pot, &
          pot_ref, nstep_ref, nteraservers, max_mpi_wait_time, cp2k_mpi_beads
 
@@ -620,10 +620,6 @@ subroutine check_inputsanity()
             write (*, *) 'Cannot compute heat capacity for NVE simulation.'
             error = 1
          end if
-         if (ncalc > nwrite) then
-            write (*, *) 'Ncalc greater than nwrite.Setting nwrite=ncalc'
-            nwrite = ncalc
-         end if
 
          if (ipimd == 1 .and. inose <= 0) then
             write (*, *) 'You have to use thermostat with PIMD! (inose>=0)'
@@ -742,7 +738,6 @@ subroutine check_inputsanity()
          call int_nonnegative(narchive, 'narchive')
 
          call int_positive(nwalk, 'nwalk')
-         call int_positive(ncalc, 'ncalc')
 
          if (error == 1) then
             call fatal_error(__FILE__, __LINE__, 'Invalid input parameters')

src/modules.F90

@@ -27,7 +27,6 @@ module mod_general
    integer :: nwrite = 1
    ! output for XYZ coordinates, velocities and forces (Molden format)
    integer :: nwritex = 1, nwritev = 0, nwritef = 0
-   integer :: ncalc = 1
    ! How often do we print restart file and how often do we archive it
    integer :: nrest = 1, narchive = 10000
    ! Restart switch, 1 = restart simulation
@@ -49,6 +48,11 @@ module mod_general
    real(DP), protected :: sim_time = 0.0D0
    ! Energy restrain MD by Jiri Suchan
    integer :: en_restraint = 0
+   ! Global flag to stop the simulation, checked at the end of each time step.
+   ! The simulation will finish prematurely (before nstep is reached),
+   ! but otherwise successfully. For example, when reaching the dE_S0S1 threshold
+   ! in Surface Hopping simulations.
+   logical :: STOP_SIMULATION = .false.
    save
 contains
    subroutine set_natom(num_atom)

src/surfacehop.F90

@@ -1097,7 +1097,7 @@ subroutine check_energy(vx_old, vy_old, vz_old, vx, vy, vz)
    end subroutine check_energy
 
    subroutine check_S0S1_gap(en_S0, en_S1, threshold_ev)
-      use mod_interfaces, only: finish
+      use mod_general, only: STOP_SIMULATION
       use mod_const, only: AUtoEV
       real(DP), intent(in) :: en_S0, en_S1
       real(DP), intent(in) :: threshold_ev
@@ -1108,8 +1108,12 @@ subroutine check_S0S1_gap(en_S0, en_S1, threshold_ev)
       if (dE_S0S1 < threshold_ev) then
          write (stdout, *) 'S1 - S0 gap dropped below threshold!'
          write (stdout, *) dE_S0S1, ' < ', threshold_ev
-         call finish(0)
-         stop 0
+         ! TODO: This is an unfortunate hack,
+         ! but we can't directly stop here since we need the last step
+         ! to be written to the output files.
+         write (stdout, *) 'Simulation will stop at the end of this step.'
+         ! This global flag is checked in the main loop in abin.F90
+         STOP_SIMULATION = .true.
       end if
    end subroutine check_S0S1_gap
 

tests/SH_S0S1/PES.dat.ref

@@ -1 +1,2 @@
  #    Time[fs] Potential energies
+           0.50   -0.9439722351E+02   -0.9407547374E+02   -0.9399745168E+02

tests/SH_S0S1/energies.dat.ref

@@ -1 +1,2 @@
  #        Time[fs] E-potential           E-kinetic     E-Total    E-Total-Avg
+           0.50   -0.9399745168E+02    0.6700581882E-01   -0.9393044586E+02

tests/SH_S0S1/pop.dat.ref

@@ -1 +1,2 @@
  #    Time[fs] CurrentState   Populations Sum-of-Populations
+           0.50  3   0.00000   0.00016   0.99984 0.9999999

tests/SH_S0S1/prob.dat.ref

@@ -1 +1,2 @@
  #    Time[fs] CurrentState   Probabilities
+           0.50  3   0.00000   0.00023   0.00000