[QE-users] Random seed for MD calculations

[Han Hsu (徐翰)]

Dear QE developers and users,

I've been performing Born-Oppenheimer molecular dynamics (BOMD)
calculations for hcp Fe at high-PT conditions using Quantum Espresso (pw.x,
calculation= 'md'). For the thermostats, I am testing stochastic velocity
rescaling (svr) and Andersen, as they are the only ones (in QE) that
generate canonical ensembles. For these two thermostats, random numbers are
generated in each time step. I find, however, that in restart calculations
(restart_mode= 'restart'), the random numbers are not generated properly.
In each restart calculation, the same sequence of random numbers are
generated, as shown in detail in the next few paragraphs. Consequently, the
calculation results (e.g pressure, temperature) exhibit a periodic
behavior, with a period same as the length of each restart calculation. My
question is: Can we avoid this problem by properly setting the random seed
in restart calculations? If so, how to do that? (For restart_mode=
'from_scratch', such a problem does not occur, as also shown in the next
few paragraphs).

Note: To print out the random numbers generated in each time step, a few
more lines (write statements) have to be added to the file
PW/src/dynamics_module.f90. Also, to obtain lots of data within a short
time, I use fcc Si (8-atom cells) to run extensive tests. For these tests,
each run consists of 12000 time steps divided into 12 rounds of
calculations: A from_scratch round (steps 1–1000) followed by 11 restart
rounds (steps 1001–2000, 2001–3000, ..., 11001–12000).

For the 'svr' thermostat, random numbers (R1, R2,..., R_Nf) from gaussian
distribution are generated in teach time step. In the from_scratch
calculation (steps 1–1000), the random numbers generated in steps 1–4 are
are as below:

     Entering Dynamics:    iteration =     1
     rr**2 = R1**2 =  2.2273
     sum_of_gaussians2 = R2**2+...+R_Nf**2 = 23.3447
     Entering Dynamics:    iteration =     2
     rr**2 = R1**2 =  0.0441
     sum_of_gaussians2 = R2**2+...+R_Nf**2 = 17.0921
     Entering Dynamics:    iteration =     3
     rr**2 = R1**2 =  0.6032
     sum_of_gaussians2 = R2**2+...+R_Nf**2 = 16.4800
     Entering Dynamics:    iteration =     4
     rr**2 = R1**2 =  1.2989
     sum_of_gaussians2 = R2**2+...+R_Nf**2 = 13.2204

Next, in the restart calculation for steps 1001–2000, the random numbers
generated in steps 1001–1004 are are as below:

     Entering Dynamics:    iteration =  1001
     rr**2 = R1**2 =  0.0008
     sum_of_gaussians2 = R2**2+...+R_Nf**2 = 18.0810
     Entering Dynamics:    iteration =  1002
     rr**2 = R1**2 =  0.8306
     sum_of_gaussians2 = R2**2+...+R_Nf**2 = 19.9568
     Entering Dynamics:    iteration =  1003
     rr**2 = R1**2 =  1.0171
     sum_of_gaussians2 = R2**2+...+R_Nf**2 = 14.5561
     Entering Dynamics:    iteration =  1004
     rr**2 = R1**2 =  0.0605
     sum_of_gaussians2 = R2**2+...+R_Nf**2 = 28.3579

Next, in the restart calculation for steps 2001–3000, the random numbers
generated in steps 2001–2004 are are as below:

     Entering Dynamics:    iteration =  2001
     rr**2 = R1**2 =  0.0008
     sum_of_gaussians2 = R2**2+...+R_Nf**2 = 18.0810
     Entering Dynamics:    iteration =  2002
     rr**2 = R1**2 =  0.8306
     sum_of_gaussians2 = R2**2+...+R_Nf**2 = 19.9568
     Entering Dynamics:    iteration =  2003
     rr**2 = R1**2 =  1.0171
     sum_of_gaussians2 = R2**2+...+R_Nf**2 = 14.5561
     Entering Dynamics:    iteration =  2004
     rr**2 = R1**2 =  0.0605
     sum_of_gaussians2 = R2**2+...+R_Nf**2 = 28.3579

Next, in the restart calculation for steps 3001–4000, the random numbers
generated in steps 3001–3004 are are as below:

     Entering Dynamics:    iteration =  3001
     rr**2 = R1**2 =  0.0008
     sum_of_gaussians2 = R2**2+...+R_Nf**2 = 18.0810
     Entering Dynamics:    iteration =  3002
     rr**2 = R1**2 =  0.8306
     sum_of_gaussians2 = R2**2+...+R_Nf**2 = 19.9568
     Entering Dynamics:    iteration =  3003
     rr**2 = R1**2 =  1.0171
     sum_of_gaussians2 = R2**2+...+R_Nf**2 = 14.5561
     Entering Dynamics:    iteration =  3004
     rr**2 = R1**2 =  0.0605
     sum_of_gaussians2 = R2**2+...+R_Nf**2 = 28.3579

Clearly, the same sequence of random numbers are generated in all restart
rounds (steps 1001–2000, 2001–3000, 3001–4000, ..., etc.), while the
from_scratch round (steps 1–1000) is an exception. I am thus wondering how
to properly generate/seed the random numbers in restart calculations?

A few more remarks:
(1) I have done this test on different facilities at different institutes.
The same pattern persists.
(2) For Andersen thermostat, the same pattern persists: In the restart
calculations, at the same respective time steps (e.g. steps 1001, 2001,
3001, ...), the same atoms are "randomly" chosen to interact with the heat
bath.

Finally, my understanding is that in Quantum Espresso, generation and seed
of random numbers are controlled by Modules/random_numbers.f90 and a
subroutine set_random_seed(). To properly set the random seed in restart
calculations, perhaps the subroutine set_random_seed() should be called
somewhere? (Probably not in dynamics_module.f90, as we do not want to
reseed the random number in every time step...)

Any suggestions and advice are appreciated. Thank you very much for the
help!!


Best regards,


Han Hsu
Department of Physics
National Central University (Taiwan)

-- 
===============================
 Han Hsu (徐翰)
 Associate Professor
 Department of Physics
 National Central University
 Taoyuan 32001, Taiwan
 Email: hanhsu at ncu.edu.tw
 Web: sites.google.com/site/hanhsuphys
 Tel: +886-3-422-7151 ext 65303
 Fax: +886-3-425-1175
===============================
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.quantum-espresso.org/pipermail/users/attachments/20230810/b7b5a86a/attachment.html>