Flying ice cube in `npt_scr` for low temperatures if no initial velocity

[elindgren]

The npt_scr ensemble can lead to a flying ice cube effect of the system at low temperatures (<100K) if one forgets to set the initial velocity in the run.in file. Flying ice cube here refers to the system having overall non-zero linear momentum, and thus the entire system translates (GIF 1). If one initializes the velocities using the velocity keyword the problem goes away (GIF2). The effect can also be avoided by using the correct_velocities keyword.

@brucefan1983 Is this expected behavior from an MD system with zero initial velocities? Does it have to do with energy transfer between modes being slow at these low temperatures? If not, how is the thermostat implemented? Are all velocities corrected in the same direction?

GIF 1: Flying ice cube

GIF 2: No flying ice cube

Here is the run.in file for generating GIF1.

potential ../C_2022_NEP3.txt
time_step 1.0
ensemble npt_scr 10 10 100 0 0 0 0 0 0 5 5 5 5 5 5 200 
dump_thermo 100
dump_position 100
run 10000

Here is the run.in file for generating GIF2.

potential ../C_2022_NEP3.txt
time_step 1.0
velocity 10
ensemble npt_scr 10 10 100 0 0 0 0 0 0 5 5 5 5 5 5 200 
dump_thermo 100
dump_position 100
run 10000

[elindgren]

The files I used to generate the example above.

flying_ice_cube_mre.zip

[elindgren]

ping @erhart1

[elindgren]

Me and Fredrik Eriksson investigated this some more, and our guess is that it has do with the scaling factor in npt_scr becoming large when the initial velocities are strictly zero, leading to a numerical issue that remains for subsequent timesteps.

We printed the kinetic energy and scaling factor from withing ensemble_npt_scr.cu

double ek[1];
  thermo.copy_to_host(ek, 1);
  int ndeg = 3 * (number_of_atoms - N_fixed);
  ek[0] *= ndeg * K_B * 0.5;
  double sigma = ndeg * K_B * temperature * 0.5;
  double factor = resamplekin(ek[0], sigma, ndeg, temperature_coupling, rng);
  factor = sqrt(factor / ek[0]);
  std::cout << "Ek: " << ek[0] << " Factor: " << factor << "\n";

Here is the output: ... Run 3 steps. Ek: 4.51601e-12 Factor: 151220 1 steps completed. Ek: 0.0985378 Factor: 1.42342 2 steps completed. Ek: 0.177958 Factor: 1.25371 3 steps completed.

Here is the run.in file

```bash
potential ../C_2022_NEP3.txt
time_step 1.0
ensemble npt_scr 10 10 100 0 0 0 0 0 0 5 5 5 5 5 5 200 
run 3

[brucefan1983]

The results are expected since when the "velocity" keyword is missing, the initial velocities are zero and the initial temperature (total kinetic energy) is zero. Most thermostats work by rescaling the velocities based on a factor which involves the instant kinetic energy in the denominator. Therefore, they cannot work starting from zero temperature. Langevin thermostat might be an expectation. So even if you want to simulate low temperatures and even aim to heat up the system starting from close to zero temperature, you need to start from a finite one, such as 1 K.

Zheyong

On Wed, Jun 28, 2023 at 11:41 PM Eric Lindgren @.***> wrote:

Me and Fredrik Eriksson investigated this some more, and our guess is that it has do with the scaling factor in npt_scr becoming large when the initial velocities are strictly zero, leading to a numerical issue that remains for subsequent timesteps.

We printed the kinetic energy and scaling factor from withing ensemble_npt_scr.cu https://github.com/brucefan1983/GPUMD/blob/f585028c93942ced16a215a27dace776389af0fe/src/integrate/ensemble_npt_scr.cu#L265

double ek[1]; thermo.copy_to_host(ek, 1); int ndeg = 3 (number_of_atoms - N_fixed); ek[0] = ndeg K_B 0.5; double sigma = ndeg K_B temperature * 0.5; double factor = resamplekin(ek[0], sigma, ndeg, temperature_coupling, rng); factor = sqrt(factor / ek[0]); std::cout << "Ek: " << ek[0] << " Factor: " << factor << "\n";

Here is the output: ... Run 3 steps. Ek: 4.51601e-12 Factor: 151220 1 steps completed. Ek: 0.0985378 Factor: 1.42342 2 steps completed. Ek: 0.177958 Factor: 1.25371 3 steps completed.

Here is the run.in file

potential ../C_2022_NEP3.txt
time_step 1.0
ensemble npt_scr 10 10 100 0 0 0 0 0 0 5 5 5 5 5 5 200
run 3

—
Reply to this email directly, view it on GitHub
<https://github.com/brucefan1983/GPUMD/issues/467#issuecomment-1611670757>,
or unsubscribe
<https://github.com/notifications/unsubscribe-auth/AF546OPSSAZETHZP4OIWPRDXNRGCLANCNFSM6AAAAAAZXCIIDQ>
.
You are receiving this because you were mentioned.Message ID:
***@***.***>

[elindgren]

The results are expected since when the "velocity" keyword is missing, the initial velocities are zero and the initial temperature (total kinetic energy) is zero. Most thermostats work by rescaling the velocities based on a factor which involves the instant kinetic energy in the denominator. Therefore, they cannot work starting from zero temperature. Langevin thermostat might be an expectation. So even if you want to simulate low temperatures and even aim to heat up the system starting from close to zero temperature, you need to start from a finite one, such as 1 K. Zheyong

Makes sense! I even tried starting from a very small temperature, 1e-16, and that was enough to not see the flying ice cube effect. I don't need to start from such low temperatures though, I had just forgotten the keyword :zipper_mouth_face: I just wanted to make sure what I saw is reasonable.

I'll go ahead and close the issue then.