From 1176b51495ecfd4c8674d67c1f8ee8f13afabfe2 Mon Sep 17 00:00:00 2001 From: Florian Knoop Date: Thu, 17 Oct 2024 18:38:52 +0200 Subject: [PATCH 1/3] canonical_configuration | add --seed for reproducibility --- src/canonical_configuration/main.f90 | 23 +++++++++++++++---- src/canonical_configuration/options.f90 | 6 +++++ src/libolle/Makefile | 2 ++ .../type_forceconstant_secondorder.f90 | 4 +++- .../type_forceconstant_secondorder_aux.f90 | 14 ++++++----- 5 files changed, 37 insertions(+), 12 deletions(-) diff --git a/src/canonical_configuration/main.f90 b/src/canonical_configuration/main.f90 index c5aaca52..0527efcf 100644 --- a/src/canonical_configuration/main.f90 +++ b/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,31 @@ 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 + seed = walltime() + else + ! use inverse of seed because integer part is ignored in lo_randomnumbers.f90 + seed = 1.0_r8/float(opts%seed) + end if + if (mw%talk) print *, '... RANDOM SEED : ', opts%seed + 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' @@ -188,9 +201,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 diff --git a/src/canonical_configuration/options.f90 b/src/canonical_configuration/options.f90 index db68e36c..7b14db18 100644 --- a/src/canonical_configuration/options.f90 +++ b/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 diff --git a/src/libolle/Makefile b/src/libolle/Makefile index bc62be77..b8ea8bdf 100644 --- a/src/libolle/Makefile +++ b/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\ diff --git a/src/libolle/type_forceconstant_secondorder.f90 b/src/libolle/type_forceconstant_secondorder.f90 index 4ed7c70a..ddcfb273 100644 --- a/src/libolle/type_forceconstant_secondorder.f90 +++ b/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 diff --git a/src/libolle/type_forceconstant_secondorder_aux.f90 b/src/libolle/type_forceconstant_secondorder_aux.f90 index d56d22ff..056ce5fb 100644 --- a/src/libolle/type_forceconstant_secondorder_aux.f90 +++ b/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 From 9843732533439b7bec86128647bd0677d2fa28c7 Mon Sep 17 00:00:00 2001 From: Florian Knoop Date: Tue, 22 Oct 2024 13:54:30 +0200 Subject: [PATCH 2/3] canonical_configuration | only print seed when used --- src/canonical_configuration/main.f90 | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/src/canonical_configuration/main.f90 b/src/canonical_configuration/main.f90 index 0527efcf..9a54a9a7 100644 --- a/src/canonical_configuration/main.f90 +++ b/src/canonical_configuration/main.f90 @@ -37,12 +37,18 @@ program canonical_configuration ! 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 - seed = 1.0_r8/float(opts%seed) + 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 - if (mw%talk) print *, '... RANDOM SEED : ', opts%seed call tw%init(iseed=mw%r, rseed=seed) ! Just make sure no clever person starts running this in parallel. From ccec69a9be1089a8337f8cc5a7f4191777b4e489 Mon Sep 17 00:00:00 2001 From: Florian Knoop Date: Tue, 22 Oct 2024 13:54:45 +0200 Subject: [PATCH 3/3] canonical_configuration | print mode frequencies when displacing each mode --- src/canonical_configuration/main.f90 | 17 ++++++++++++++--- 1 file changed, 14 insertions(+), 3 deletions(-) diff --git a/src/canonical_configuration/main.f90 b/src/canonical_configuration/main.f90 index 9a54a9a7..1af0b7d1 100644 --- a/src/canonical_configuration/main.f90 +++ b/src/canonical_configuration/main.f90 @@ -134,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) T(K) (K) (A)' +if (opts%modes) then + write (*, "(1X,A)") ' no. +/- ek(K) ep(K) T(K) (K) (A) frequency(THz)' +else + write (*, "(1X,A)") ' no. +/- ek(K) ep(K) T(K) (K) (A)' +end if dumpconf: block type(lo_mdsim) :: sim type(lo_crystalstructure) :: p @@ -291,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