Mdstep refactor (#108)

* abin.F90: print_footer(), decrease indentation
* Use abstract interface
* normalize pot_ref

src/abin.F90

@@ -22,31 +22,30 @@ program abin
    use mod_arrays
    use mod_files, only: stdout
    use mod_general, only: sim_time, pot, pot_ref, iremd, ipimd, &
-      & md, nwrite, nstep, ncalc, it, inormalmodes, istage, irest
+      & nwrite, nstep, ncalc, it, inormalmodes, istage, irest
    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
    use mod_kinetic, only: temperature
-   use mod_utils, only: abinerror, archive_file, get_formatted_date_and_time
+   use mod_utils, only: abinerror, archive_file
    use mod_transform
-   use mod_mdstep
+   use mod_mdstep, only: mdstep
    use mod_minimize, only: minimize
    use mod_analysis, only: analysis, restout
    use mod_interfaces
    use mod_en_restraint
-   use mod_plumed
    use mod_terampi_sh, only: move_new2old_terash
    use mod_mpi, only: get_mpi_rank, mpi_barrier_wrapper
-   use mod_remd
+#ifdef USE_MPI
+   use mod_remd, only: remd_swap, nswap
+#endif
    implicit none
-   ! TODO: These should probably be defined and stored in some module, not here
    real(DP) :: dt = 20.0D0, eclas = 0.0D0, equant = 0.0D0
    logical :: file_exists
-   integer, dimension(8) :: time_start, time_end
+   integer, dimension(8) :: time_start
    integer :: my_rank
-   real(DP) :: total_cpu_time
 
-   call date_and_time(VALUES=time_start)
+   call date_and_time(values=time_start)
 
    ! INPUT AND INITIALIZATION SECTION
    call init(dt)
@@ -73,33 +72,34 @@ program abin
    write (stdout, '(A)') 'Job started at: '//trim(get_formatted_date_and_time(time_start))
    write (stdout, *) ''
 
-   ! Transform coordinates and velocities Path Integral MD
+   ! Transform coordinates and velocities for Path Integral MD
    ! (staging or normal modes)
    if (istage == 1 .or. inormalmodes > 0) then
       call initialize_pi_transforms(x, y, z, vx, vy, vz)
    end if
 
-   ! Note that 'amt' equals 'am' for non-PI simulations
+   ! Note that 'amt' equals physical atomic masses in non-PI simulations
    px = amt * vx
    py = amt * vy
    pz = amt * vz
 
    if (ipimd == 3) then
-
       call minimize(x, y, z, fxc, fyc, fzc, eclas)
-
-   else
+      call print_footer(time_start)
+      call finish(0)
+      stop 0
+   end if
 
       write (stdout, *)
-      write (stdout, *) '#      Step     Time [fs]'
+      write (stdout, '("#",10X,A,11X,A)') 'Step','Time [fs]'
 
       ! ---------------- PROPAGATION-------------
 
       ! Without this Barrier, ranks > 0 do not write geom.dat in force_clas
       ! I don't know why the hell not.
       call mpi_barrier_wrapper()
 
-      ! getting initial forces and energies
+      ! Get initial forces and energies
       call force_clas(fxc, fyc, fzc, x, y, z, eclas, pot)
       if (ipimd == 1) then
          call force_quantum(fxq, fyq, fzq, x, y, z, amg, equant)
@@ -110,13 +110,13 @@ program abin
          call lz_rewind(en_array_lz)
       end if
 
-      ! if we use reference potential with RESPA
+      ! Get reference forces and energies for ab initio  MTS RESPA
       if (pot_ref /= '_none_') then
          call force_clas(fxc, fyc, fzc, x, y, z, eclas, pot_ref)
          call force_clas(fxc_diff, fyc_diff, fzc_diff, x, y, z, eclas, pot)
       end if
 
-      ! setting initial values for surface hopping
+      ! Set initial values for surface hopping
       if (ipimd == 2) then
          if (irest /= 1) then
             call get_nacm()
@@ -163,21 +163,8 @@ program abin
 
          end if
 
-         ! CALL the integrator, propagate through one time step
-         select case (md)
-         case (1)
-            call respastep(x, y, z, px, py, pz, amt, amg, dt, equant, eclas, fxc, fyc, fzc, fxq, fyq, fzq)
-         case (2)
-            call verletstep(x, y, z, px, py, pz, amt, dt, eclas, fxc, fyc, fzc)
-            ! include entire Ehrenfest step, in first step, we start from pure initial state so at first step we dont
-            ! need NAMCE and take forces just as a grad E
-         case (3)
-            call respashake(x, y, z, px, py, pz, amt, amg, dt, equant, eclas, fxc, fyc, fzc, fxq, fyq, fzq)
-         case (4)
-            call doublerespastep(x, y, z, px, py, pz, amt, amg, dt, equant, eclas, fxc, fyc, fzc, fxq, fyq, fzq)
-         end select
-
-         sim_time = sim_time + dt
+         ! PROPAGATE through one time step
+         call mdstep(x, y, z, px, py, pz, amt, dt, eclas, fxc, fyc, fzc)
 
          vx = px / amt
          vy = py / amt
@@ -219,7 +206,7 @@ program abin
          end if
 #endif
 
-         ! --- Ttrajectory analysis ---
+         ! --- Trajectory analysis ---
          ! In order to analyze the output, we have to perform the back transformation
          ! Transformed (cartesian) coordinates are stored in trans matrices.
 
@@ -236,30 +223,30 @@ program abin
             call QtoX(x, y, z, transx, transy, transz)
             call FQtoFX(fxc, fyc, fzc, transfxc, transfyc, transfzc)
             call analysis(transx, transy, transz, transxv, transyv, transzv,  &
-                 &       transfxc, transfyc, transfzc, eclas, equant)
+                 &       transfxc, transfyc, transfzc, eclas)
 
          else if (inormalmodes > 0) then
 
             call UtoX(x, y, z, transx, transy, transz)
             call UtoX(vx, vy, vz, transxv, transyv, transzv)
             call UtoX(fxc, fyc, fzc, transfxc, transfyc, transfzc)
             call analysis(transx, transy, transz, transxv, transyv, transzv,  &
-                 &        transfxc, transfyc, transfzc, eclas, equant)
+                 &        transfxc, transfyc, transfzc, eclas)
          else
 
-            call analysis(x, y, z, vx, vy, vz, fxc, fyc, fzc, eclas, equant)
+            call analysis(x, y, z, vx, vy, vz, fxc, fyc, fzc, eclas)
 
          end if
 
          if (modulo(it, nwrite) == 0) then
-            write (stdout, '(I20,F15.2)') it, sim_time * AUtoFS
+            write (stdout, '(I15,F15.2)') it, sim_time * AUtoFS
             call flush (OUTPUT_UNIT)
          end if
 
          ! Time step loop
       end do
 
-      ! DUMP restart file at the end of a run
+      ! 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
       it = it - 1
@@ -276,24 +263,9 @@ program abin
          call restout(x, y, z, vx, vy, vz, it)
       end if
 
-      ! minimization endif
-   end if
-
-   write (stdout, *) ''
-   write (stdout, '(A)') 'Job finished successfully!'
-
+   call print_footer(time_start)
    call finish(0)
 
-   call cpu_time(total_cpu_time)
-   write (stdout, '(A)') 'Total cpu time [s] (does not include ab initio calculations)'
-   write (stdout, *) total_cpu_time
-   write (stdout, '(A)') 'Total cpu time [hours] (does not include ab initio calculations)'
-   write (stdout, *) total_cpu_time / 3600.
-
-   write (stdout, '(A)') 'Job started at:  '//trim(get_formatted_date_and_time(time_start))
-   call date_and_time(VALUES=time_end)
-   write (stdout, '(A)') 'Job finished at: '//trim(get_formatted_date_and_time(time_end))
-
 contains
 
    subroutine clean_temp_files()
@@ -302,4 +274,26 @@ subroutine clean_temp_files()
       call system('rm -f ERROR engrad*.dat.* nacm.dat hessian.dat.* geom.dat.*')
    end subroutine clean_temp_files
 
+   function get_formatted_date_and_time(time_data) result(formatted_string)
+      character(len=25) :: formatted_string
+      integer, dimension(8), intent(in) :: time_data
+      formatted_string = ''
+      ! time_data must be get from date_and_time() intrinsic
+      ! e.g. 1:48:39   3.11.2020
+      write (formatted_string, "(I2.2,A1,I2.2,A1,I2.2,A2,I2,A1,I2,A1,I4)") time_data(5), ':', &
+         time_data(6), ':', time_data(7), '  ', time_data(3), '.', time_data(2), '.', &
+         time_data(1)
+   end function get_formatted_date_and_time
+
+   subroutine print_footer(time_start)
+      integer, dimension(8), intent(in) :: time_start
+      integer, dimension(8) :: time_end
+
+      call date_and_time(values=time_end)
+      write (stdout, *) ''
+      write (stdout, '(A)') 'Job finished successfully!'
+      write (stdout, '(A)') 'Job started at:  '//trim(get_formatted_date_and_time(time_start))
+      write (stdout, '(A)') 'Job finished at: '//trim(get_formatted_date_and_time(time_end))
+   end subroutine print_footer
+
 end program abin

src/analysis.F90

@@ -19,7 +19,7 @@ module mod_analysis
 contains
 
    ! Contains all analysis stuff
-   subroutine analysis(x, y, z, vx, vy, vz, fxc, fyc, fzc, eclas, equant)
+   subroutine analysis(x, y, z, vx, vy, vz, fxc, fyc, fzc, eclas)
       use mod_analyze_ext, only: analyze_ext
       use mod_estimators, only: estimators
       use mod_general, only: it, ipimd, icv, nwrite, nwritef, nwritev, &
@@ -33,13 +33,7 @@ subroutine analysis(x, y, z, vx, vy, vz, fxc, fyc, fzc, eclas, equant)
       real(DP), intent(inout) :: x(:, :), y(:, :), z(:, :)
       real(DP), intent(in) :: fxc(:, :), fyc(:, :), fzc(:, :)
       real(DP), intent(inout) :: vx(:, :), vy(:, :), vz(:, :)
-      real(DP), intent(in) :: eclas, equant
-      real(DP) :: energy
-
-      ! eclas is the ab initio energy averaged per bead,
-      ! equant is additional harmonic energy between PI beads (from force_quantum)
-      ! TODO: Print equant or energy somewhere
-      energy = eclas + equant
+      real(DP), intent(in) :: eclas
 
       if (modulo(it, nwrite) == 0 .and. idebug > 0) then
          call remove_rotations(x, y, z, vx, vy, vz, am, .false.)
@@ -322,7 +316,7 @@ end subroutine restout
    ! It is called from subroutine init.
    subroutine restin(x, y, z, vx, vy, vz, it)
       use mod_general, only: icv, ihess, nwalk, ipimd, natom, &
-                             iremd, pot, sim_time
+                             iremd, pot, update_simtime
       use mod_nhc, only: readNHC, inose, nhc_restin
       use mod_mpi, only: get_mpi_rank
       use mod_estimators
@@ -336,6 +330,7 @@ subroutine restin(x, y, z, vx, vy, vz, it)
       real(DP), intent(out) :: x(:, :), y(:, :), z(:, :)
       real(DP), intent(out) :: vx(:, :), vy(:, :), vz(:, :)
       integer, intent(out) :: it
+      real(DP) :: sim_time
       integer :: iat, iw
       integer :: my_rank
       character(len=100) :: chtemp
@@ -357,6 +352,7 @@ subroutine restin(x, y, z, vx, vy, vz, it)
 
       open (111, file=chin, status="OLD", action="READ")
       read (111, *) it, sim_time
+      call update_simtime(sim_time)
       read (111, '(A)') chtemp
       call checkchar(chtemp, chcoords)
       do iw = 1, nwalk

src/ekin.F90

@@ -66,6 +66,7 @@ subroutine temperature(px, py, pz, amt, eclas)
       end if
    end subroutine temperature
 
+   ! WARNING: We assume physical masses here!
    real(DP) function ekin_v(vx, vy, vz)
       use mod_general, only: nwalk, natom
       use mod_system, only: am

src/en_restraint.F90

@@ -42,15 +42,13 @@ end subroutine en_rest_init
 
    subroutine energy_restraint(x, y, z, px, py, pz, eclas)
       use mod_general, only: natom, nwalk, dt0 ! dt0 is the time step
-      use mod_const, only: AMU
       use mod_system, only: am
       use mod_sh, only: en_array
       use mod_terampi_sh, only: force_terash
       real(DP), intent(inout) :: x(:, :), y(:, :), z(:, :)
       ! TODO: Eclas is not modified in this routine, but probably should be
       real(DP), intent(inout) :: eclas
       real(DP), intent(inout) :: px(:, :), py(:, :), pz(:, :)
-      ! DH: I am surprised this works, since natom is not a constant
       real(DP), dimension(natom) :: fxgs, fygs, fzgs, fxes, fyes, fzes
       real(DP) :: eclasexc, eclasground, Egrad, deltaE, lambda, lsum, deltaEnext, convercrit, deltaD
       integer :: ios, iat, iat2, iw

src/force_terash.F90

@@ -276,9 +276,8 @@ subroutine send_terash(x, y, z, tc_comm)
       call MPI_SSend(bufints, nstate * (nstate - 1) / 2 + nstate, MPI_INTEGER, 0, TC_TAG, tc_comm, ierr)
       call handle_mpi_error(ierr)
 
-      ! temporary hack
-      bufdoubles(1) = sim_time ! * AUtoFS !* dt
       ! Send Time
+      bufdoubles(1) = sim_time
       call MPI_Send(bufdoubles, 1, MPI_DOUBLE_PRECISION, 0, TC_TAG, tc_comm, ierr)
       call handle_mpi_error(ierr)
 

src/forces.F90

@@ -275,77 +275,3 @@ subroutine force_quantum(fx, fy, fz, x, y, z, amg, quantum_energy)
    quantum_energy = equant
 end subroutine force_quantum
 
-! Based on:
-! Efficient stochastic thermostatting of path integral molecular dynamics
-! J. Chem. Phys. 133, 124104 ?2010?
-subroutine propagate_nm(x, y, z, px, py, pz, m, dt)
-   use mod_const, only: DP, PI
-   use mod_array_size
-   use mod_general, only: nwalk, natom
-   use mod_nhc, only: temp
-   implicit none
-   real(DP), intent(inout) :: x(:, :), y(:, :), z(:, :)
-   real(DP), intent(inout) :: px(:, :), py(:, :), pz(:, :)
-   real(DP), intent(in) :: m(:, :), dt
-   real(DP) :: omega(size(x, 2)), omega_n, om, omt, c, s
-   real(DP) :: x_old(size(x, 1), size(x, 2))
-   real(DP) :: y_old(size(x, 1), size(x, 2))
-   real(DP) :: z_old(size(x, 1), size(x, 2))
-   real(DP) :: px_old(size(px, 1), size(px, 2))
-   real(DP) :: py_old(size(px, 1), size(px, 2))
-   real(DP) :: pz_old(size(px, 1), size(px, 2))
-   integer :: iat, iw, k
-
-   ! expecting m = am = physical masses
-
-   x_old = x
-   y_old = y
-   z_old = z
-   px_old = px
-   py_old = py
-   pz_old = pz
-
-   omega_n = TEMP * NWALK
-
-   do iw = 2, nwalk
-      k = iw - 1
-      omega(iw) = 2 * omega_n * sin(k * PI / nwalk)
-   end do
-
-   ! First, propagate centroid
-   iw = 1
-   do iat = 1, natom
-      X(iat, iw) = X(iat, iw) + dt * PX(iat, iw) / M(iat, iw)
-      Y(iat, iw) = Y(iat, iw) + dt * PY(iat, iw) / M(iat, iw)
-      Z(iat, iw) = Z(iat, iw) + dt * PZ(iat, iw) / M(iat, iw)
-   end do
-
-   ! eq 23 from J. Chem. Phys. 133, 124104 2010
-   ! exact propagation of a free ring polymer in normal mode coordinates
-   do iw = 2, nwalk
-      om = omega(iw)
-      omt = omega(iw) * dt
-      c = cos(omt)
-      s = sin(omt)
-      do iat = 1, natom
-         ! Propagate positions
-         x(iat, iw) = x_old(iat, iw) * c &
-                      + px_old(iat, iw) * s / m(iat, iw) / om
-         y(iat, iw) = y_old(iat, iw) * c &
-                      + py_old(iat, iw) * s / m(iat, iw) / om
-         z(iat, iw) = z_old(iat, iw) * c &
-                      + pz_old(iat, iw) * s / m(iat, iw) / om
-
-         ! propagate momenta
-         px(iat, iw) = px_old(iat, iw) * c &
-                       - x_old(iat, iw) * s * m(iat, iw) * om
-         py(iat, iw) = py_old(iat, iw) * c &
-                       - y_old(iat, iw) * s * m(iat, iw) * om
-         pz(iat, iw) = pz_old(iat, iw) * c &
-                       - z_old(iat, iw) * s * m(iat, iw) * om
-
-      end do
-
-   end do
-
-end subroutine propagate_nm

src/fortran_interfaces.F90

@@ -49,13 +49,6 @@ subroutine force_wrapper(x, y, z, fx, fy, fz, eclas, chpot, walkmax)
          character(len=*), intent(in) :: chpot
       end subroutine force_wrapper
 
-      subroutine propagate_nm(x, y, z, px, py, pz, amg, dt)
-         import :: DP
-         real(DP), intent(inout) :: x(:, :), y(:, :), z(:, :)
-         real(DP), intent(inout) :: px(:, :), py(:, :), pz(:, :)
-         real(DP), intent(in) :: amg(:, :), dt
-      end subroutine propagate_nm
-
       subroutine oniom(x, y, z, fx, fy, fz, eclas, iw)
          import :: DP
          real(DP), intent(in) :: x(:, :), y(:, :), z(:, :)