canonical_configuration | add --seed for reproducibility

src/canonical_configuration/main.f90

@@ -11,6 +11,7 @@ program canonical_configuration
 use type_forceconstant_secondorder, only: lo_forceconstant_secondorder
 use type_jij_secondorder, only: lo_jij_secondorder
 use type_mdsim, only: lo_mdsim
+use lo_randomnumbers, only: lo_mersennetwister
 
 use options, only: lo_opts
 use semirandom, only: generate_semirandom_configurations
@@ -23,19 +24,37 @@ program canonical_configuration
 type(lo_mpi_helper) :: mw
 type(lo_mem_helper) :: mem
 integer, dimension(:), allocatable :: alloy_permutation
+type(lo_mersennetwister), save :: tw
 
 ! Get the necessary things first
 init: block
+    real(r8) :: seed
     ! Get CLI options
     call opts%parse()
     call mw%init()
     call mem%init()
     if (.not. mw%talk) opts%verbosity = -100
 
+    ! initialize random numbers
+    if (opts%seed < 0) then
+        ! fully random
+        seed = walltime()
+        if (mw%talk) print *, '... walltime() used to initialize random state'
+    else
+        ! use inverse of seed because integer part is ignored in lo_randomnumbers.f90
+        if (mw%talk) print *, '... RANDOM SEED: ', opts%seed
+        if (opts%seed == 0) then
+            seed = 0.0_r8
+        else
+            seed = 1.0_r8/float(opts%seed)
+        end if
+    end if
+    call tw%init(iseed=mw%r, rseed=seed)
+
     ! Just make sure no clever person starts running this in parallel.
-    ! if ( mw%n .gt. 1 ) then
-    !     call lo_stop_gracefully(['Do not run this in parallel.'],lo_exitcode_param,__FILE__,__LINE__)
-    ! endif
+    if (mw%n .gt. 1) then
+        call lo_stop_gracefully(['Do not run this in parallel.'], lo_exitcode_param, __FILE__, __LINE__)
+    end if
 
     ! Read structures
     if (mw%talk) write (*, *) '... reading infiles'
@@ -115,7 +134,11 @@ program canonical_configuration
     if (mw%talk) write (*, '(1X,A,I5,A)') '--> ', opts%nconf, ' structures'
 end if
 
-write (*, "(1X,A)") '  no.   +/-      ek(K)          ep(K)          <ek/ep>         T(K)          <T>(K)      <msd>(A)'
+if (opts%modes) then
+    write (*, "(1X,A)") '  no.   +/-      ek(K)          ep(K)          <ek/ep>         T(K)          <T>(K)      <msd>(A)         frequency(THz)'
+else
+    write (*, "(1X,A)") '  no.   +/-      ek(K)          ep(K)          <ek/ep>         T(K)          <T>(K)      <msd>(A)'
+end if
 dumpconf: block
     type(lo_mdsim) :: sim
     type(lo_crystalstructure) :: p
@@ -188,9 +211,9 @@ program canonical_configuration
             else
                 invert = .false.
             end if
-            call fcss%initialize_cell(p, uc, fc, opts%temperature, opts%zpm, .false., opts%mindist, mw, imode=imode, invert=invert)
+            call fcss%initialize_cell(p, uc, fc, opts%temperature, opts%zpm, .false., opts%mindist, mw, imode=imode, invert=invert, tw=tw)
         else
-            call fcss%initialize_cell(p, uc, fc, opts%temperature, opts%zpm, .false., opts%mindist, mw)
+            call fcss%initialize_cell(p, uc, fc, opts%temperature, opts%zpm, .false., opts%mindist, mw, tw=tw)
         end if
 
         ! dump to file
@@ -272,8 +295,15 @@ program canonical_configuration
         else
             invert_char = '+'
         end if
-        write (*, "(1X,I4,5X,A,5(2X,F13.5),2X,F15.8)") iconf, invert_char, &
-            ek/lo_kb_hartree/1.5_r8, ep/lo_kb_hartree/1.5_r8, ratio, temp, avgtemp/iconf, avgmsd*lo_bohr_to_A/iconf
+        if (opts%modes) then
+            write (*, "(1X,I4,5X,A,5(2X,F13.5),2X,F15.8,2X,F10.3)") iconf, invert_char, &
+                ek/lo_kb_hartree/1.5_r8, ep/lo_kb_hartree/1.5_r8, ratio, temp, &
+                avgtemp/iconf, avgmsd*lo_bohr_to_A/iconf, fcss%omega(imode)*lo_frequency_Hartree_to_THz
+        else
+            write (*, "(1X,I4,5X,A,5(2X,F13.5),2X,F15.8)") iconf, invert_char, &
+                ek/lo_kb_hartree/1.5_r8, ep/lo_kb_hartree/1.5_r8, ratio, temp, &
+                avgtemp/iconf, avgmsd*lo_bohr_to_A/iconf
+        end if
     end do
 
     ! And, at the end, maybe dump the fake simulation

src/canonical_configuration/options.f90

@@ -26,6 +26,7 @@ module options
     real(r8) :: dielcutoff2 = -lo_huge
     real(r8) :: dielcutoff3 = -lo_huge
     logical :: modes = .false.
+    integer :: seed = -lo_hugeint
 contains
     procedure :: parse
 end type
@@ -119,6 +120,10 @@ subroutine parse(opts)
                  help='Print displacements for every individual mode.', hidden=.true., &
                  required=.false., act='store_true', def='.false.', error=lo_status)
     if (lo_status .ne. 0) stop
+    call cli%add(switch='--seed', &
+                 help='Seed the random number generator.', &
+                 hidden=.false., required=.false., act='store', def='-1', error=lo_status)
+    if (lo_status .ne. 0) stop
 
     cli_manpage
     cli_verbose
@@ -153,6 +158,7 @@ subroutine parse(opts)
     call cli%get(switch='-dc2', val=opts%dielcutoff2)
     call cli%get(switch='-dc3', val=opts%dielcutoff3)
     call cli%get(switch='--modes', val=opts%modes)
+    call cli%get(switch='--seed', val=opts%seed)
 
     ! Convert input to atomic units right away
     opts%maximum_frequency = opts%maximum_frequency*lo_frequency_THz_to_hartree

src/libolle/Makefile

@@ -528,6 +528,7 @@ $(OBJECT_PATH)type_forceconstant_firstorder_io.o:\
 
 $(OBJECT_PATH)type_forceconstant_secondorder.o:\
 	type_forceconstant_secondorder.f90\
+	type_forceconstant_secondorder_aux.f90\
     $(OBJECT_PATH)konstanter.o\
     $(OBJECT_PATH)gottochblandat.o\
     $(OBJECT_PATH)mpi_wrappers.o\
@@ -543,6 +544,7 @@ $(OBJECT_PATH)type_forceconstant_secondorder_dynamicalmatrix.o:\
     $(OBJECT_PATH)type_blas_lapack_wrappers.o
 	$(FC) $(FFLAGS) -c type_forceconstant_secondorder_dynamicalmatrix.f90 -o $@
 $(OBJECT_PATH)type_forceconstant_secondorder_aux.o:\
+	type_forceconstant_secondorder.f90\
 	type_forceconstant_secondorder_aux.f90\
     $(OBJECT_PATH)type_forceconstant_secondorder.o\
     $(OBJECT_PATH)type_distancetable.o\

src/libolle/type_forceconstant_secondorder.f90

@@ -9,6 +9,7 @@ module type_forceconstant_secondorder
 use gottochblandat, only: tochar, walltime, lo_clean_fractional_coordinates, lo_chop, lo_sqnorm
 use mpi_wrappers, only: lo_mpi_helper, lo_stop_gracefully
 use lo_memtracker, only: lo_mem_helper
+use lo_randomnumbers, only: lo_mersennetwister
 use lo_longrange_electrostatics, only: lo_ewald_parameters
 use type_crystalstructure, only: lo_crystalstructure
 use type_qpointmesh, only: lo_qpoint
@@ -190,7 +191,7 @@ module subroutine fake_forceconstant(fc, uc, ss, debye_temperature, maximum_freq
         real(r8), intent(in), optional :: maximum_frequency
         integer, intent(in), optional :: verbosity
     end subroutine
-    module subroutine initialize_cell(fcss, ss, uc, fc, temperature, quantum, exact, closest_distance, mw, nosync, imode, invert)
+    module subroutine initialize_cell(fcss, ss, uc, fc, temperature, quantum, exact, closest_distance, mw, nosync, imode, invert, tw)
         class(lo_forceconstant_secondorder), intent(inout) :: fcss
         type(lo_crystalstructure), intent(inout) :: ss
         type(lo_crystalstructure), intent(in) :: uc
@@ -203,6 +204,7 @@ module subroutine initialize_cell(fcss, ss, uc, fc, temperature, quantum, exact,
         logical, intent(in), optional :: nosync
         integer, intent(in), optional :: imode
         logical, intent(in), optional :: invert
+        type(lo_mersennetwister), intent(inout), optional :: tw
     end subroutine
     module subroutine setsumtozero(fc)
         class(lo_forceconstant_secondorder), intent(inout) :: fc

src/libolle/type_forceconstant_secondorder_aux.f90

@@ -1,6 +1,5 @@
 submodule(type_forceconstant_secondorder) type_forceconstant_secondorder_aux
 use konstanter, only: lo_exitcode_blaslapack
-use lo_randomnumbers, only: lo_mersennetwister
 use gottochblandat, only: lo_planck, lo_harmonic_oscillator_free_energy, lo_negsqrt
 use type_distancetable, only: lo_distancetable
 use type_blas_lapack_wrappers, only: lo_dsyevr, lo_dsyevd
@@ -363,7 +362,7 @@ subroutine set_values(alpha, fc)
 end subroutine
 
 !> use the harmonic model to initialize a cell
-module subroutine initialize_cell(fcss, ss, uc, fc, temperature, quantum, exact, closest_distance, mw, nosync, imode, invert)
+module subroutine initialize_cell(fcss, ss, uc, fc, temperature, quantum, exact, closest_distance, mw, nosync, imode, invert, tw)
     !> force constant for this (super) cell
     class(lo_forceconstant_secondorder), intent(inout) :: fcss
     !> supercell to be thermally populated
@@ -388,10 +387,12 @@ module subroutine initialize_cell(fcss, ss, uc, fc, temperature, quantum, exact,
     integer, intent(in), optional :: imode
     !> use negative mode amplitude so that u_i = \sum_s -Q_s X_si)
     logical, intent(in), optional :: invert
+    !> maybe there is already a random number generator
+    type(lo_mersennetwister), intent(inout), optional :: tw
     real(r8) :: inv_prefactor
 
     ! Not sure about save attribute here.
-    type(lo_mersennetwister), save :: tw
+    type(lo_mersennetwister), save :: tw_local
     integer :: solrnk
     logical :: sync
 
@@ -403,10 +404,11 @@ module subroutine initialize_cell(fcss, ss, uc, fc, temperature, quantum, exact,
     end if
 
     init: block
+        if (present(tw)) tw_local = tw
         ! Seed rng if needed
-        if (tw%initialized .eqv. .false.) then
+        if (tw_local%initialized .eqv. .false.) then
             ! make sure to seed it with the rank to make it different on all MPI ranks.
-            call tw%init(iseed=mw%r, rseed=walltime())
+            call tw_local%init(iseed=mw%r, rseed=walltime())
         end if
 
         ! Decide on syncinc thingies
@@ -481,7 +483,7 @@ module subroutine initialize_cell(fcss, ss, uc, fc, temperature, quantum, exact,
                 dstloop: do iter = 1, maxiter
                     ! Generate displacements
                     modeloop: do i = 1, ss%na*3
-                        call tw%rnd_boxmuller_pair(1.0_r8, 0.0_r8, x1, x2)
+                        call tw_local%rnd_boxmuller_pair(1.0_r8, 0.0_r8, x1, x2)
                         l = 0
                         do a1 = 1, ss%na
                         do j = 1, 3