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.

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

Mdstep refactor #108

Merged
merged 5 commits into from
Apr 7, 2022
Merged
Changes from 1 commit
Commits
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
Prev Previous commit
Next Next commit
abin.F90: print_footer(), decrease indentation
  • Loading branch information
danielhollas committed Apr 3, 2022

Verified

This commit was signed with the committer’s verified signature.
paring-chan 파링
commit 846dd5285e3c33a606430c49a58625d8de4b2209
74 changes: 40 additions & 34 deletions src/abin.F90
Original file line number Diff line number Diff line change
@@ -27,27 +27,27 @@ program abin
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_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
use mod_mdstep, only: mdstep, md
#ifdef USE_MPI
use mod_remd, only: remd_swap
#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)
@@ -86,21 +86,22 @@ program abin
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)
@@ -111,13 +112,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()
@@ -176,8 +177,6 @@ program abin
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

vx = px / amt
vy = py / amt
vz = pz / amt
@@ -218,7 +217,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.

@@ -251,14 +250,14 @@ program abin
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
@@ -275,24 +274,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)'
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 removed the CPU time printing, because it's not indicative of anything as it does not account for external program execution.

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()
@@ -301,4 +285,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
3 changes: 1 addition & 2 deletions src/force_terash.F90
Original file line number Diff line number Diff line change
@@ -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)

19 changes: 10 additions & 9 deletions src/mdstep.F90
Original file line number Diff line number Diff line change
@@ -5,6 +5,7 @@
module mod_mdstep
use mod_const, only: DP
use mod_utils, only: abinerror
use mod_general, only: update_simtime
use mod_transform
implicit none
private
@@ -36,21 +37,17 @@ module mod_mdstep
end type

abstract interface
Copy link
Contributor Author

Choose a reason for hiding this comment

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

Gettting a bit fancy here, here's a nice explanation:
https://fortran-lang.org/learn/best_practices/callbacks

Copy link
Contributor

@suchanj suchanj Apr 7, 2022

Choose a reason for hiding this comment

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

That's a cool trick. 'Hello pointers my old friend.'

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Yeah. The good thing here is that the logic for selecting the integrator is localized in one place, whereas before it was both in init.F90 and abin.F90.

subroutine integrator(x, y, z, px, py, pz, amt, dt, f, e)
import :: DP, energies, forces
subroutine integrator(x, y, z, px, py, pz, amt, dt, E_pot, fx, fy, fz)
import :: DP
real(DP), dimension(:, :), intent(inout) :: x, y, z, px, py, pz
real(DP), dimension(:, :), intent(in) :: amt
type(forces), intent(inout) :: f
type(energies), intent(inout) :: e
real(DP), intent(in) :: dt
real(DP), intent(inout) :: E_pot
real(DP), dimension(:, :), intent(inout) :: fx, fy, fz
end subroutine integrator
! subroutine verle(x, y, z, px, py, pz, amt, dt, eclas, fxc, fyc, fzc)
! subroutine respa(x, y, z, px, py, pz, amt, amg, dt, equant, eclas, fxc, fyc, fzc, fxq, fyq, fzq)
! subroutine shake(x, y, z, px, py, pz, amt, amg, dt, equant, eclas, fxc, fyc, fzc, fxq, fyq, fzq)
! subroutine dresp(x, y, z, px, py, pz, amt, amg, dt, equant, eclas, fxc, fyc, fzc, fxq, fyq, fzq)
end interface

procedure(integrator), public :: mdstep
procedure(integrator), pointer, public :: mdstep

contains

@@ -276,6 +273,7 @@ subroutine verletstep(x, y, z, px, py, pz, amt, dt, eclas, fxc, fyc, fzc)

if (inose > 0) call thermostat(px, py, pz, amt, dt / 2)

call update_simtime(dt)
end subroutine verletstep

! RESPA ALGORITHM 10.12.2012
@@ -322,6 +320,7 @@ subroutine respastep(x, y, z, px, py, pz, amt, amg, dt, equant, eclas, &

call thermostat(px, py, pz, amt, dt / (2 * nabin))

call update_simtime(dt)
end subroutine respastep

! RESPA ALGORITHM WITH RATTLE 10.12.2012
@@ -406,6 +405,7 @@ subroutine respashake(x, y, z, px, py, pz, amt, amg, dt, equant, eclas, &

if (inose == 1) call shiftNHC_yosh(px, py, pz, amt, dt / (2 * nabin))

call update_simtime(dt)
end subroutine respashake

! Double RESPA algorithm using reference low-cost potential pot_ref with smaller time step
@@ -467,6 +467,7 @@ subroutine doublerespastep(x, y, z, px, py, pz, amt, amg, dt, equant, eclas, &

if (inose > 0) call thermostat(px, py, pz, amt, dtsm / 2)

call update_simtime(dt)
end subroutine doublerespastep

end module mod_mdstep
16 changes: 11 additions & 5 deletions src/modules.F90
Original file line number Diff line number Diff line change
@@ -1,8 +1,7 @@
!-File with core modules created by Daniel Hollas,9.2.2012
! File with global simulation parameters

! We are using modules to initialize some variables and
! for passing global variables to different subroutines.
!------------------------------------------------------------------------------------

! mod_array_size contains various array limits.
! Modify here if you need larger arrays.
@@ -62,14 +61,21 @@ module mod_general
integer :: iknow = 0
! Linux Process ID, populated automatically for the current ABIN process
integer :: pid
! Future variables for adaptive timestep in SH
real(DP) :: dt0, sim_time = 0.0D0
! Initial time step (for future adaptime timestep functionality)
real(DP) :: dt0
! Total simulation time
real(DP), protected :: sim_time = 0.0D0
! Energy restrain MD by Jiri Suchan
integer :: en_restraint = 0
save
contains
subroutine update_simtime(dt)
use mod_const, only: DP
real(DP) :: dt
sim_time = sim_time + dt
end subroutine
end module

! Some information about simulated system, especially for distributions and shake
! TODO: Move this to a separate file, and think hard what should be inside this module.
module mod_system
use mod_const, only: DP