long-range part of lj/long missing?

Martin_Fitzner · April 10, 2017, 5:28pm

Dear LAMMPS users,

I seem to be getting unreasonable results when using the pair style lj/long/coul/long with the 31Mar2017 version of lammps (and older versions).

In particular I was simulating the wetting behavior of water droplets on rigid substrates. When I use a lj/cut potential for the substrate-water interaction I see that there is a strong dependence of the wetting-transition temperature, which converges with Rcut. However when I compute the transition temperature with lj/long/coul/long the result doesn’t agree with the largest hard cut-off simulation but rather with the simulation that has the same hard cut-off Rcut as the real space cut-off Rreal in the lj/long one. It seems as if the long-range part of the forces is missing (although I notice a computational slow-down when using lj/long).

To confirm this I looked at a simpler system: a single atom above a slab of rigid atoms, interacting with lj/cut or lj/long/coul/long respectively. I computed binding curves of this atom for different choices of Rcut/Rreal. In the picture (link: ) you can see that the curves for lj/cut and lj/long/coul/long are identical (slightly shifted in x-direction because otherwise they exactly overlap) if Rcut=Rreal which again suggests that the long-range part is missing. Am I right in the assumption that the binding-curve for lj/long/coul/long should in principle be independent of the real space cut-off used?

Here is the part of my input concerning the interactions:

#Eps=0.68, Sigma=2.5027, nCut varies

pair_style lj/long/coul/long long off (v_nCut*v_Sigma) pair_coeff 1 2 {Eps} {Sigma} # C-mW pair_coeff 2 2 0.0 {Sigma} # C-C
pair_coeff 1 1 0.0 ${Sigma} # mW-mW
kspace_style pppm/disp 1.0e-4
kspace_modify force/disp/real 0.0001
kspace_modify force/disp/kspace 0.002
kspace_modify mix/disp none
neigh_modify delay 0 every 1 check yes

which results in the following output about PPPMdisp:

PPPMDisp initialization …

Using 2 structure factors
Dispersion G vector (1/distance) = 1.41955e-15
Dispersion grid = 2 2 2
Dispersion stencil order = 5
Dispersion estimated absolute RMS force accuracy = 0.000120517
Dispersion estimated absolute real space RMS force accuracy = 0.000120517
Dispersion estimated absolute kspace RMS force accuracy = 0
Dispersion estimated relative force accuracy = 3.62934e-07
using double precision FFTs
3d grid and FFT values/proc dispersion = 216 2

I read the corresponding parts of the manual and also tried to choose the kspace settings by hand, but to no avail, i.e. there seems to be no long-range part no matter what I do. Is there something wrong with the way I am using these commands or could there be a bug?

Best wishes,

Martin

akohlmey · April 11, 2017, 1:04am

martin,

i cannot reproduce your claim. with a simple input modified from the melt example, i see the expected behavior.
i.e. for pppm/disp the distribution of energy between evdwl and elong changes as the real space cutoff changes, but the total remains similar and also the pressure is the same within the limit of accuracy, while it changes for lj/cut and approaches the pppm/disp value.

please try the attached input and see, if you get the same result to confirm that your executable is compiled correctly.

axel.

$ grep -A1 ^Step log.lammps
Step TotEng E_vdwl E_long Press
0 0 0 0 0

in.lj-long-bulk (1.02 KB)

Martin_Fitzner1 · April 11, 2017, 3:29pm

Dear Axel,

Upon running the script you sent I got identical output to yours.

However when looking at evdwl and elong in my test system, elong seems to be constant and numerically zero at any height of the atom above the substrate. Do you reckon something might be wrong with the kspace parameters? Maybe the PPPMDisp output indicates something is wrong?

PPPMDisp initialization …
Using 2 structure factors
Dispersion G vector (1/distance) = 2.8391e-15
Dispersion grid = 2 2 2
Dispersion stencil order = 5
Dispersion estimated absolute RMS force accuracy = 6.89468e-05
Dispersion estimated absolute real space RMS force accuracy = 6.89468e-05
Dispersion estimated absolute kspace RMS force accuracy = 0
Disperion estimated relative force accuracy = 2.07631e-07
using double precision FFTs
3d grid and FFT values/proc dispersion = 343 8

I attached a minimal example of my test including the structure file.

Martin.

$ grep -A2 Step log.lammps
Step TotEng E_vdwl E_long
0 616357.11 616357.11 -1.1818496e-48
1 616357.11 616357.11 -1.1818496e-48

in.test.long (1.52 KB)

test.data (888 KB)

akohlmey · April 12, 2017, 3:17pm

Dear Axel,

Upon running the script you sent I got identical output to yours.

However when looking at evdwl and elong in my test system, elong seems to
be constant and numerically zero at any height of the atom above the
substrate. Do you reckon something might be wrong with the kspace
parameters? Maybe the PPPMDisp output indicates something is wrong?

yes, i cannot even run your input with the current version of LAMMPS (11
April 2017):

$ ~/compile/lammps-icms/src/lmp_omp -in in.test.long
LAMMPS (11 Apr 2017)
  using 1 OpenMP thread(s) per MPI task
Reading data file ...
  orthogonal box = (-0.731996 -0.731996 -43.175) to (175.602 175.602 66.825)
  1 by 1 by 1 MPI processor grid
  reading atoms ...
  21601 atoms
1 atoms in group mW
21600 atoms in group C
Setting atom values ...
  1 settings made for x
Setting atom values ...
  1 settings made for y
WARNING: No fixes defined, atoms won't move (../verlet.cpp:55)
PPPMDisp initialization ...
  Using 2 structure factors
ERROR: Cannot compute initial g_ewald_disp (../pppm_disp.cpp:3609)
Last command: run 1 pre no post no

this means, that your setup is so, that the algorithms cannot converge to
meaningful parameters.
my guess is, that you are using an older version of LAMMPS, that is
missing the necessary sanity checks and thus your generated kspace settings
are likely bogus.

axel.

Martin_Fitzner · April 20, 2017, 4:43pm

Ok, thanks for trying out the script. Interestingly, when I use the latest LAMMPS version from github (11 April 2017) I still don’t get an error and the output is identical to the older versions.

I have now tried to find suitable parameters (in particular the ewald parameter) by screening many different choices but did not succeed. It seems that I can choose a parameter that gives the correct long-range part but fails at short range. This leads me to the question: Is it possible that there is no set of parameters that make lj/long work with my system?

I wonder which particular “property” of my system could cause this problem? Again, I’m dealing with a surface slab (so in principle x,y periodic, not necessarily in z) with something on top of that.

Martin.

akohlmey · April 21, 2017, 9:44pm

yes, i cannot even run your input with the current version of LAMMPS (11
April 2017):
Ok, thanks for trying out the script. Interestingly, when I use the latest
LAMMPS version from github (11 April 2017) I still don't get an error and
the output is identical to the older versions.

what kind of compiler and optimization settings are you using and what
kind of machine are you running on?
i've taken a closer look and ran some debugging tools and found two
locations where a (theoretically harmless) division by zero can happen and
updated the code to not trigger it.

however, bother ewald/disp and pppm/disp are written in somewhat unique
programming style that make it much more difficult to understand and debug
than most of LAMMPS. it also makes it much more sensitive to compiler bugs.

this means, that your setup is so, that the algorithms cannot converge to
meaningful parameters.
my guess is, that you are using an older version of LAMMPS, that is
missing the necessary sanity checks and thus your generated kspace settings
are likely bogus.

I have now tried to find suitable parameters (in particular the ewald
parameter) by screening many different choices but did not succeed. It
seems that I can choose a parameter that gives the correct long-range part
but fails at short range. This leads me to the question: Is it possible
that there is no set of parameters that make lj/long work with my system?

I wonder which particular "property" of my system could cause this
problem? Again, I'm dealing with a surface slab (so in principle x,y
periodic, not necessarily in z) with something on top of that.

i don't know. the lj/long and */disp code hasn't seen much maintenance
recently and also - as i already mentioned - is not exactly
maintainer-friendly, so it is quite possible, that something is not fully
in sync with the overall kspace code or some heuristics are no longer valid
due to changes in either pair.cpp or kspace.cpp. outside of (essentially)
cosmetic changes, there is little than i can do at this moment beyond what
i've done so far.

if there was somebody with reasonable programming skills, a meticulous
mindset and some time, i would recommend to re-implement, all of the pair
lj/long and related pair styles as well as ewald/disp and pppm/disp to be
in a more maintainable code style that also more closely follows the flow
of control of the current ewald and pppm implementation. these sounds more
drastic than it actually is, since one would mostly have to disentangle the
various shortcuts and unusual constructs with simpler and less compact
versions and then make certain, that known-to-be-good results are still
reproduced. at that point, it would be easier to identify any possible
logic errors that are currently mostly hidden by the style of the existing
implementation.

this part of the code is a great example for the dangers of premature
optimization, favoring complex/compact code constructs over simpler/obvious
ones.

axel.