respa

Dear LAMMPS users,

I am trying to simulate the NVT simulation of a simple alkane in LAMMPS using the method respa to integrate the equation of motion. If I use the method verlet to integrate the equation of motion with a step of 1 fs, there are no problems, but I would like to increase the integration step for pair interactions and kspace in order to speed up my calculations. To do this, I used the integration method respa by specifying the following command:

run_style respa 2 8 and run_style respa 2 4

increasing the integration step to 8 fs, I also checked 4 fs. The calculation is started but only the first step of integration and the values of all energy are issued, and then everything ā€œstopsā€ and no information appears in the log file, however the task does not fall. Please tell me how to correctly switch from integration by the method verlet to method respa, in order to increase the integration step for pair and long-range electrostatic interactions. To what value can I increase the outer integration step for paired and long-range electrostatic interactions?

Best regards,
Victor

Please always report which version of LAMMPS you are using and what platform you are running on.

Dear LAMMPS users,

I am trying to simulate the NVT simulation of a simple alkane in LAMMPS using the method respa to integrate the equation of motion. If I use the method verlet to integrate the equation of motion with a step of 1 fs, there are no problems, but I would like to increase the integration step for pair interactions and kspace in order to speed up my calculations. To do this, I used the integration method respa by specifying the following command:

run_style respa 2 8 and run_style respa 2 4

where are the assignments of what part of the force calculation to run at what respa level?

Axel

Dear Prof. Kohlmeyer,

I use LAMMPS version 19 Mar 2020.

For assigned force I use command

timestep 8.0 or timestep 4.0
run_style respa 2 8 bond 1 angle 1 dihedral 1 pair 2 kspace 2

Victor

Š²Ń, 5 ŠøюŠ». 2020 Š³. Š² 15:36, Axel Kohlmeyer <[email protected]>:

when trying to use r-RESPA, i would recommend the following workflow.

  • start from a reasonably equilibrated restart/data file (with run style verlet).
  • use the most conservative neighbor list setting every 1 delay 0 check yes
  • use fix nve to be able to monitor energy conservation
  • do not use fix rigid or fix shake or anything equivalent
  • get a ā€œbaselineā€ timestep, i.e. empirically determine what kind of timestep is required with verlet to get a good energy conservation. that will be the largest timestep you may have for your innermost timestep
  • switch to r-RESPA with the same timestep and run_style respa 2 1 bond 1 pair 2 this should get identical energies for at least a 1000 timesteps
  • now gradually try to explore how you can decouple the time domains, i.e. start with doubling the timestep and using run_style respa 2 2 bond 1 pair 2
  • see if you can switch to doubling again and run_style respa 2 4 bond 1 pair 2
  • if you have a very large number of dihedrals in your system (usually only important when using implicit solvent), you might try running dihedrals together with pair instead of bond or do run_style respa 3 2 2 bond 1 dihedral 2 pair 3.
  • after you have optimized these settings (make sure you have sufficient energy conservation) and you if are using long-range electrostatics, you might try running kspace every other step, so double the timestep one more time and try something like run_style respa 3 4 2 bond 1 pair 2 kspace 3
  • if a factor of 4 is too large, you may also try a factor 3, but that makes it harder to compare energies on the same timesteps.
  • if you have found a satisfactory setting for the division of the timestepping, you can try to reduce the cost of neighbor list builds. keep in mind, that the neighbor list build is on the outer timestep, so be careful with increasing every, better to see, if you can set delay to 1 or 2. you can compute the average number of (outer) timesteps between neighbor list builds from the neighbor list statistics at the end of a run. best to set the delay to one or two steps fewer than the average. keep in mind that the ā€œreturn on investmentā€ is largest at the beginning and the with ā€œevery 1 check yesā€ the extra cost is primarily the O(N) loop over the atoms to check if half the skin distance has past, not a full neighbor list rebuild. also choosing a larger skin can help (it decreases the frequency of required neighbor list rebuilds, but increases the cost for pair style computation). keep in mind, that with very aggressive r-RESPA level multipliers, you may need to increase the skin or atoms may move too far between outer timesteps, even if you rebuild on every step.
  • if the pair style supports it, you may also consider using ā€œinnerā€ and ā€œouterā€ (i have not found a case where using ā€œmiddleā€ is beneficial) to divide the non-bonded forces by a factor of 2. ā€œinnerā€ would then likely need to be run with ā€œdihedralā€ and ā€œouterā€ with ā€œkspaceā€ unless you get lucky. the choice of inner cutoff and overlap region is crucial, so this is best tested independently (i.e. first run ā€œinnerā€ at the level 1 and determine what the separation of time between those can maximally be, which will give you the optimal timestep for the ā€œouterā€ level, then increase the separation between ā€œbondā€ and ā€œinnerā€ for as long as it maintains energy conservation while reducing the separation between ā€œinnerā€ and ā€œouterā€ accordingly).

in summary, the basic principle is that each ā€œdomainā€ has its ā€œnaturalā€ optimal timestep, so you have to watch your choices to not push beyond that. but since you can only specify those timesteps implicitly through integer multipliers, you need to adapt accordingly. this is particularly true for the ā€œinnerā€, ā€œmiddleā€, ā€œouterā€ divisions. those only make sense with large cutoffs. you need a sufficiently large overlap region, but that costs performance because the atom pairs in the overlap region(s) are computed twice. but with inner/middle/outer you cannot push for a larger timestep on outer than for inner/outer only, and the optimal timestep for inner will be the same as well.

axel.

1 Like

Dear Prof. Kohlmeyer,

Thank you very much for your detailed explanation.

Best regards,
Victor

Š²Ń, 5 ŠøюŠ». 2020 Š³. Š² 19:22, Axel Kohlmeyer <akohlmey@ā€¦43ā€¦4ā€¦>: