Mdstep refactor #108
Mdstep refactor #108
Conversation
danielhollas
added
the
refactor
| @@ -276,77 +276,3 @@ subroutine force_quantum(fx, fy, fz, x, y, z, amg, quantum_energy) | |||
| quantum_energy = equant | |||
| end subroutine force_quantum | |||
|
|
|||
| ! Based on: | |||
There was a problem hiding this comment.
Moved to mdstep.F90
| @@ -206,15 +219,6 @@ subroutine initialize_pi_masses(amg, amt) | |||
| ! write(*,*)'amg2',(amg(iat,iw)*2*sin(pi*(iw-1)/nwalk),iw=1,nwalk) | |||
| end if | |||
| end do | |||
|
|
|||
| else if (inormalmodes == 1) then | |||
There was a problem hiding this comment.
This case is covered by default, see top of the function
Codecov Report
@@ Coverage Diff @@
## master #108 +/- ##
==========================================
- Coverage 87.08% 86.96% -0.12%
==========================================
Files 41 42 +1
Lines 5805 5815 +10
==========================================
+ Hits 5055 5057 +2
- Misses 750 758 +8
|
| end select | ||
|
|
||
| sim_time = sim_time + dt | ||
| ! PROPAGATE through one time step |
There was a problem hiding this comment.
Here's the nice part enabled by this refactor.
| ! nstep_ref - internal steps with reference potential pot_ref | ||
| integer, public :: nstep_ref = 1 | ||
|
|
||
| abstract interface |
There was a problem hiding this comment.
Gettting a bit fancy here, here's a nice explanation:
https://fortran-lang.org/learn/best_practices/callbacks
There was a problem hiding this comment.
That's a cool trick. 'Hello pointers my old friend.'
There was a problem hiding this comment.
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.
danielhollas
force-pushed
the
mdstep-refactor
branch
from
9c0ee43 to
3135eb8
Compare
| allocate (ishake2(natom * 3 - 6)) | ||
| ishake1 = 0 | ||
| ishake2 = 0 | ||
| allocate (ishake1(natom * 3 - 6), source=0) |
There was a problem hiding this comment.
New trick, allocation and initialization in one statement.
| call finish(0) | ||
|
|
||
| call cpu_time(total_cpu_time) | ||
| write (stdout, '(A)') 'Total cpu time [s] (does not include ab initio calculations)' |
There was a problem hiding this comment.
I removed the CPU time printing, because it's not indicative of anything as it does not account for external program execution.
| 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 |
There was a problem hiding this comment.
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?
There was a problem hiding this comment.
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))There was a problem hiding this comment.
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 unlimitedThere was a problem hiding this comment.
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.
Trying to move more stuff out of the big init subroutine into smaller chunks.
I also wanted to make a nice interface to various integrators. I am hitting a problem here that various integrators have slightly different interface, e.g. some of them need extra force arrays. This problem in general and different ways to solve it are outlined here:
https://fortran-lang.org/learn/best_practices/type_casting
Even before this refactor this problem manifested for the MTS RESPA (
doublerespastep()) integrator, where we have additional forces fx_diff, which we do not pass in as arguments but import them from global module. This is equivalent to approach no. 3 from the link above, and I have utilized this approach for the other integrators as well so that I could make a minimal common interface (see top of mdstep module). Long terms it would be nice to define our own derive types that would hold the simulation state and that we could pass explicitly into all functions (approach no. 2). I have outline how that could look in #107.There are couple more unrelated refactors, mostly to pass data explicitly via subroutine arguments, instead of importing them from modules. This makes it more clear what depends on what.