Do we need a BAOAB-splitting Langevin integrator?

[ghost]

The current Langevin integrator implemented in GPUMD uses the OBABO-splitting reported by Bussi et al. in [1]. But according to Leimkuhler et al. [2], the Langevin integrator with a BAOAB-splitting could provide much better stability and configuration-sampling ability when using a large time step. So this splitting methodology could be useful in some cases, especially when the users are interested in the conformational changes of complex systems. If we need this new integrator, I would like to open a PR and work around it.

[1] Bussi, G., & Parrinello, M. (2007). Accurate sampling using Langevin dynamics. Physical Review E, 75(5), 056707.
[2] Fass, J., Sivak, D. A., Crooks, G. E., Beauchamp, K. A., Leimkuhler, B., & Chodera, J. D. (2018). Quantifying configuration-sampling error in Langevin simulations of complex molecular systems. Entropy, 20(5), 318.

[brucefan1983]

Good suggestion, go ahead if you are interested in this.

[ghost]

Thanks a lot! I have written the code to implement the new integrator, you can find them in my GPUMD fork. The new integrator is named ENSMBLE_BAO and could be invoked from the input file with the keyword nvt_bao. Due to the different splitting scheme, I could not use the velocity_verlet function from the base class. So I divided this function into different parts and put them in the new source file. I think this is something that could be improved.

And here are my test results: I tested the new integrator on a 20×20×20 silicon lattice system (with 428 atoms in it), under a temperature of 3500 K. I run multiple simulations using different timesteps. After these simulations, I calculated the RDF curves from the trajectories. There is the run.in file:

potential     ../../../potentials/tersoff/Si_Fan_2019.txt 0
velocity      300

ensemble      npt_ber 300  3500 100 0 0 0 53.4059 53.4059 53.4059 2000
time_step     1
dump_thermo   100
run           200000

ensemble      npt_ber 3500 3500 100 0 0 0 53.4059 53.4059 53.4059 2000
time_step     1
dump_thermo   100
run           50000

ensemble      nvt_bao 3500 3500 400
time_step     0.5
dump_thermo   100
dump_position 200
run           1000000

First, I checked that the new integrator could generate the same RDF curves as the RDF curves generated by the OBABO integrator implemented in GPUMD (I know they look poor since the simulation time is short) : Then, I compared the RDF curves generated by the two integrators under different timesteps. It seems that with large timesteps, the new integrator could generate RDF curves with lower errors, which stands for better stability: But I am too busy to test the heat functionality, and since the existing tests are really simple, the tested code is not guaranteed to be bug-free. So I think only if I or someone run more tests on the new code, the code could meet the standard of merging.

To reproduce my tests, you could clone and compile my fork, then use the following xyz.in: xyz.in.zip

[ghost]

I run the example thermal_transport_nemd_and_hnemd/ballistic with the new integrator while keeping other parameters unchanged. It seems that the heat_bao method works fine, but the temperature control is poor. This is due to the splitting scheme, so maybe we should not perform NEMD with this integrator. Here is the jupyter notebook:

[brucefan1983]

Very impressive! I will have a cloer look during the weekend.

[brucefan1983]

Interesting results. The temperature control in NEMD is indeed less ideal and you can leave it out for the time being. As you are adding something new, there is no risk of breaking any existing things and I encourage you continue to make a PR. Then we can discuss the implementation.

[ghost]

I would like to! But since I have token the gpu_correct_momentum function from your ENSEMBLE_LAN code, which may contain bugs as we discussed in https://github.com/brucefan1983/GPUMD/issues/170 , should I open the PR now or wait until the bugs are fixed (then I could use the correct gpu_correct_momentum function directly)?

[brucefan1983]

Ok, let me fix that bug first. Thanks.

[ghost]

I have created a new PR https://github.com/brucefan1983/GPUMD/pull/173 with the new gpu_correct_momentum function. Could you please review the code?