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.

Sign up for GitHub

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

Mdstep refactor #108

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
f882836
WIP: mdstep refactor
danielhollas Apr 2, 2022
846dd52
abin.F90: print_footer(), decrease indentation
danielhollas Apr 3, 2022
a170cac
Use abstract interface
danielhollas Apr 3, 2022
8161401
Fix MPI build
danielhollas Apr 3, 2022
3135eb8
normalize pot_ref
danielhollas Apr 4, 2022
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
104 changes: 49 additions & 55 deletions src/abin.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand All @@ -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)
Expand All @@ -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()
Expand Down Expand Up @@ -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
Copy link
Contributor Author

@danielhollas danielhollas Apr 3, 2022 •
edited
Loading

Choose a reason for hiding this comment

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

Here's the nice part enabled by this refactor.

call mdstep(x, y, z, px, py, pz, amt, dt, eclas, fxc, fyc, fzc)

vx = px / amt
vy = py / amt
Expand Down Expand Up @@ -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.

Expand All @@ -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
Expand All @@ -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)'
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()
Expand All @@ -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
14 changes: 5 additions & 9 deletions src/analysis.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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, &
Expand All @@ -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.)
Expand Down Expand Up @@ -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
Expand All @@ -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
Expand All @@ -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
Expand Down
1 change: 1 addition & 0 deletions src/ekin.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
2 changes: 0 additions & 2 deletions src/en_restraint.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Copy link
Contributor

Choose a reason for hiding this comment

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

Btw do you have the answer why this works? Maybe the variable is initialized once the subroutine is called and at that time the natom is known?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Exactly. These are so called automatic arrays that get allocated on the stack when the function is called.
Instead of relying on the global natom variable, a slightly better way would be to use the size() intrinsic to get the size of the input arrays, like this:

real(DP) :: a(size(x, 1), size(x, 2))

Copy link
Contributor Author

Choose a reason for hiding this comment

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

The disadvantage in comparision to allocatable arrays that are allocated on the heap is that stack size is usually limited. So if we were simulating big system, we might get a segfault / stack overflow. So in general manually allocating arrays of size (natom, nwalk) is probably a bit better for the future.

On Unix, the stack size can be increased by

ulimit -s unlimited

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Actually, reading that SO link I posted above, looks like whether automatic arrays are placed on the stack or on the heap depends on the compiler and compiler options. So if this ever became an issue in the future, tweaking compiler options would solve it.

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
Expand Down
3 changes: 1 addition & 2 deletions src/force_terash.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)

Expand Down
74 changes: 0 additions & 74 deletions src/forces.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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:
Copy link
Contributor Author

Choose a reason for hiding this comment

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

Moved to mdstep.F90

! 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
7 changes: 0 additions & 7 deletions src/fortran_interfaces.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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(:, :)
Expand Down
Loading