@bcdandurand this is a promising idea. I will look at the implementation and give it a try.

@bcdandurand Thanks for contributing the work, this is potentially very helpful for large scale simulations with long-range electrostatics on a large number of MPI processes. Have you observed any speedup with this implementation in your test runs?

I can build with CMake and test with a modified version of bench/in.rhodo and can run with the GPU package with pair/only yes. I have a couple of (early) comments:

1. The implementation logic is clear and concise. Communicating only the grid-based charge density (R->K) and potential (K->R) is an improvement to the existing VerletSplit to siginficantly reduce the amount of data transfer between R and K processes.
2. The change where the R-K communication is done requires changes to the logic in the run() function compared to VerletSplit. The addition of rk_flag and 3 new virtual functions to the core KSpace class makes the VerletSplitRK::run() looks cleaner than having to check if force->kspace style name matches pppm/rk in the init() function and to cast force->kspace to (PPPM_RK*).
3. PPPM_RK currently comments out compute(), but there are quite a few classes that invoke force->kspace->compute(). Do you think this function should be enabled? Also, is it possible to support compute_group_group()?
4. Please add sections for verlet/split/rk and pppm/rk in the corresponding doc pages (run_style.rst and kspace_style.rst).
5. I think you could remove the commented out sections with #if 0 .. #endif in verlet_split_rk.cpp.

@ndtrung81 Thank you very much for reviewing this contribution.
First, for the question about observed speedup, yes, improvement is observed where there would be communication bottleneck between R- and K-processes. At small node counts and for R to K ratio small, there would not be much bottleneck and thus little improvement. Otherwise, for larger node count or large R to K ratio, improvement is noticeable. For example, if I run VISCOSITY using non-default settings

replicate       4 4 4
kspace_style   pppm/rk 1.0e-4
run_style      verlet/split/rk

then the simulation takes about 47 seconds, versus about 77 seconds with the non-default (otherwise baseline) settings

replicate       4 4 4
kspace_style   pppm 1.0e-4
run_style      verlet/split

In both runs, I use 96 R processes and 4 K processes with a sbatch script similar to

#!/bin/bash

#SBATCH --ntasks=100                            # Number of tasks to run (for MPI tasks this is number of cores)
#SBATCH --nodes=1                               # Maximum number of nodes on which to launch tasks
#SBATCH --exclusive
#SBATCH --mem=0

mpirun -n 100 ~/lammpsfork/lammps/build/lmp -partition 96 4 -screen vsrk_r96k4 -log none -in in.cos.1000SPCE

@ndtrung81
PPPM::compute: Currently, the PPPM::compute would be invoked from an instance of PPPM_RK. Whether this is ideal default behavior is up for debate. It probably is not sensible, since the same computation would be done redundantly by both the communicator of R processes and the communicator of K processes. I'll reply again after I've thought of what the use of pppm_rk requires outside of the use of verlet/split/rk.

Commented out code: Otherwise, the currently commented out implementation of PPPM_RK::compute might be refactored as a comment, to indicate the contents of a decomposed alternative to PPPM::compute.
Other commented code might likewise be refactored as comments if deemed useful, otherwise delete.

compute_group_group: I'm not as familiar with this (though I'm aware of its existence), I can look into this next.

Please advise as to the preferred way to provide updated changes in addition to pushing to the bcdandurand fork.

doc pages: I'll look into this and enquire as needed.

I'm looking at incorporating an analogous decomposition for PPPM::compute_group_group. For examples that would use this, I see
examples/PACKAGES/tally/in.pe
examples/PACKAGES/tally/in.force

I'm looking at incorporating an analogous decomposition for PPPM::compute_group_group. For examples that would use this, I see examples/PACKAGES/tally/in.pe examples/PACKAGES/tally/in.force

Testing in.pe under the script modifications:

kspace_style    pppm/rk 1.0e-4
replicate       4 4 4
run             500

with sbatch command, e.g.,
mpirun -n 60 ~/lammpsfork/lammps/build/lmp -partition 30 30 -screen verlet_r30k30_pppm_rk -log none -in in.pe

I'll soon have an update with the addition of an implementation of PPPM_RK::compute_group_group and ancillary structures. While it separates the execution between R- and K-processes, and produces the same output for tally/in.pe as with the baseline use of kspace_style pppm 1.0e-4, it does so in a manner that does not appear to be so useful in terms of speedup. If it is to be of any use, there needs to be a modified version of the ComputeGroupGroup class that is analogous to the VerletSplitRK variant of VerletSplit. Developing the variant of ComputeGroupGroup, if it's worth it, should probably be part of a later pull request. I need to do a bit more cleanup with the addition of PPPM_RK::compute_group_group before I push to the bcdandurand forked repo. (I won't be able to do so for a couple days.)

I can then turn to the doc pages after the above is finished.

@bcdandurand Before working on additional features, you should probably first pay attention to clean integration of your code into the LAMMPS code base. You pull request branch currently has a merge conflict with upstream and fails to compile on all our test platforms.

@akohlmey @ndtrung81 et al.
I'm working on adding sections for verlet/split/rk and pppm/rk in the corresponding doc pages (run_style.rst and kspace_style.rst). This is in tandem with deciding what to do about PPPM_RK::compute_group_group. The issue I'm bringing up in this comment is also analogous with an issue with the implementation of PPPM_RK::compute, which is already integrated with this fork.

The issue is due to the assumption in PPPM_RK that the initial MPI communicator is partitioned (see PPPM_RK::setupRKBlock). While it's possible to implement PPPM_RK::compute and PPPM_RK::compute_group_group within the partitioned framework, the invocation of these two functions, within, e.g., Verlet::run and ComputeGroupGroup::kspace_contribution, respectively, which do not assume a partitioning of the MPI communicator, would be detrimental to runtime performance. This is due to the default behavior where identical execution occurs on R-processes and K-processes. If there are 100 processes, of which 96 are R-processes and 4 are K-processes, the speedup will be bottlenecked as if only 4 processes are being used.

Currently, there is only beneficial use of MPI parallelism if kspace_style pppm/rk is paired with run_style verlet/split/rk. We would need an analogous group/group/rk compute style implemented in, say, a ComputeGroupGroupRK class for any benefit from a call to PPPM_RK::compute_group_group to be realized.

My initial thought to prevent users from using an ineffective combination of styles (from a parallel computing perspective) is to have an rk_flag in the Verlet class (or higher up in Integrate) and other relevant classes (like Compute), which would be visible to instances of KSpace. An error can be thrown if there is not consistency between the rk_flag values of these classes when their respective styles are used.

The issue described in this comment is also relevant to updating the manual, because users would need to have consistency in the rk option for run_style, kspace_style, etc.

I hope this is clear, I'll clarify where helpful.

@bcdandurand If you don't have time to complete the work on compute group/group/rk at the moment, we can leave it to a later PR. If possible, please wrap up the documentation for this PR and get it ready to the final round of review and to merge.

@bcdandurand If you don't have time to complete the work on compute group/group/rk at the moment, we can leave it to a later PR. If possible, please wrap up the documentation for this PR and get it ready to the final round of review and to merge.

@ndtrung81 Understood. The cluster I normally work on is under maintenance, I expect to have access to it again on Friday. I'll then address the documentation as soon as I can.

Thanks @bcdandurand for submitting these enhancements. And for the thorough explanation in the PR description of what you've done. Also thanks to @ndtrung81 for looking it over carefully and testing it out.

I have just a couple hi-level suggestions.

(1) In the doc pages for run_style and/or kspace_style for the new variant styles could you provide some numbers for sample speed-ups. This is to motivate users to try out the new options. Ideally, the numbers would compare vanilla verlet vs verlet/split (previous) vs verlet/split/rk. Could be explained in a few sentences, or add a simple table with 3 columns of numbers.

(2) If you like, you can add a "citation" in the source code for your pending paper which will trigger a prompt for the user to cite it, whenever the verlet/split/rk option is used.