Removed triclinic restrictions on slab ewald/pppm.#411
Conversation
|
Unfortunately, it isn't this simple. This PR fails the most basic sanity check. An orthogonal box should give the same energy and force as a triclinic box when all the tilt factors are equal to 0.0, since the triclinic box isn't really tilted. However, if I add this line below to an input script using this PR I get significantly different forces and energies than the orthogonal case: change_box all triclinic xy final 0.0 xz final 0.0 yz final 0.0 A second good check for triclinic systems is that when there are periodic boundary conditions in change_box all triclinic xy final ${Ly} xz final ${Lz} yz final 0.0
A third good check for kspace slab systems is that this input script and data file (attached below) give the analytic solution of forces in magnitude equal to |
|
For the Kspace slab analytic solution, just place two oppositely charged particles in a slab system. They form two infinite sheets of charge due to periodic boundary conditions in the x and y directions. If the charges are far apart enough in the z-direction, the force in the z-direction for the orthogonal case converges to |
|
Sorry, in the above post that should be tilt factor = 0.0, not 1.0. I've edited the post above. If I run the script I attached with a pure orthogonal box, I get an energy = 6.5933251. If I run it with triclinic and all zero tilt factors, I get an energy = 2.6182985. The energies should match in these two cases but don't. |
|
This is certainly odd behavior. I can confirm it even with your first
example. The force is radically different even by setting the box using:
run 0
print "f = ${fz}"
change_box all triclinic xy final 0.0001 xz final 0.0 yz final 0.0
run 0
print "f = ${fz}"
It's definitely related to the extension in the z-direction, since when
I set
"kspace_modify slab 1.0".
the forces before and after switching to triclinic are nearly identical:
f = 3.14494261727358
f = 3.1449426178097
…On 3/11/17 9:19 PM, Stan Moore wrote:
Sorry, in the above post that should be tilt factor = 0.0, not 1.0.
I've edited the post above. If I run the script I attached with a pure
orthogonal box, I get an energy = 6.5933251. If I run it with
triclinic and all zero tilt factors, I get an energy = 2.6182985. The
energies should match in these two cases but don't.
—
You are receiving this because you authored the thread.
Reply to this email directly, view it on GitHub
<#411 (comment)>, or
mute the thread
<https://github.com/notifications/unsubscribe-auth/ACI7d7yrbNikO1oeu6WLQc9UMS4D1i7bks5rk1Y8gaJpZM4MaEec>.
|
|
i looked into this a little bit today. turns out, the slab correction does work correctly with the provided modification, but that is not the origin of the difference. the issue is, that activating the slab correction changes the volume of the cell and that means, that all computations done before applying the slab correction, need to account for that correctly, and that is where things get a bit complicated, since for triclinic cells, all of those computations are done in lamda coordinates. |
The force-style YAML tests (kspace-*_tri_slab.yaml) exercise the pair+kspace force in a triclinic slab box, but they never invoke "compute group/group", whose own slab correction (KSpace::slabcorr_groups) was separately enabled for triclinic cells when EW3DC was generalized to non-orthogonal boxes. This adds unittest/commands/test_kspace_group_group_slab.cpp, a GoogleTest that builds the two-charged-sheet reproducer from lammps#411 (charges +/-1 in a box periodic in x,y and non-periodic in z, Lx=1 Ly=2) and checks the group A <-> group B interaction for ewald and pppm: - the orthogonal reference matches the analytic infinite-sheet result, z force = 2*pi*q^2/(Lx*Ly) = pi and energy = 2*pi*q^2*d/(Lx*Ly) = 5*pi; - the triclinic xy-tilted box reproduces the orthogonal result (Ewald to ~1e-13, PPPM within mesh tolerance), confirming the triclinic group/group slab path. Verified passing on 1, 2, and 4 MPI ranks. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Since the slab correction applies specifically to the z-direction, it is not influenced by tilt. I have confirmed the correction by calculating it manually for a system with 2 charges (+1 at 0,0,1.5 and -1 at 0.0,-1.5). The files are attached:
ewtest.zip
Four cases are tested -- an orthogonal and a triclinic box with and without the slab correction.
Although changing the box tilt alters the energies and forces, the correction term is identical for both, as it should be.