[lammps-users] Ewald Parameters

LAMMPS LAMMPS Mailing List Mirror
1

Dear LAMMPS,

I am trying to understand the Ewald implementation in LAMMPS. It seems the only thing you can specify is a precision.

I think there are generally 3 parameters that go into the method. Using the notation in the DLPOLY 2.17 manual, page 61

http://www.google.com/url?sa=t&source=web&cd=1&ved=0CBYQFjAA&url=ftp%3A%2F%2Fftp.dl.ac.uk%2Fccp5%2FDL_POLY%2FDL_POLY_CLASSIC%2FDOCUMENTS%2FUSRMAN2.17.pdf&rct=j&q=dlpoly%20manual&ei=RLpRTe62C8GclgfAkvCsCg&usg=AFQjCNHnJN6N8nCB7w8_Ih6PLMiGb0JW6A&sig2=G0dpR6oIZzeX3kmJKBU8qg&cad=rja

these are kappa, alpha and the number of lattice vectors. Is LAMMPS using some scheme to determine these parameters from the precision? If so, what is this scheme?

I see in ewald.cpp

g_ewald = (1.35 - 0.15*log(precision))/rcut
kmax = something

Can someone help me to understand ewald.cpp and which variables correspond to those in the DLPOLY manual (a description is not given in the LAMMPS manual)?

Is there a way for me to know what these values are that are being used? I do see that g_ewald, and “vectors: actual 1d max” are printed. Is the 1d value kmax?

The reason I want to know this is that I’d like to relate the values LAMMPS is using back to kappa, alpha, kmax as used in the DLPOLY manual as I’m trying to check results to another code. With this other code I have some differences, about 0.3 for rcut = 14 Angstroms and 0.02 for rcut = box/2. These are generally good, but I would feel more confident in both codes if I could match them as close as possibly (within rounding error and precision of the Boltzmann constant, etc. that I’ve used to convert units).

Sorry for my naive question.

Thank you in advance.

Best.

Josh

2

There are actually 2 settings in LAMMPS that effect
the Ewald (or PPPM) discretization: the precision
and the cutoff. Internally the code uses both of those
to set Gewald and the Kspace resolution (# of vectors,
grid spacing). You can override that with the
kspace_modify gewald command if you want to
set that param directly. I would just experiment
with those 3 settings on a simple system and see
what values LAMMPS prints out after the setup
phase, and you should quickly understand how the
various params are related. For PPPM, you can also
specify the grid size directly via kspace_modify mesh
but that isn't an option for Ewald.

Steve

3

Thanks Steve!

I do understand about the options for PPPM and appreciate how you can control things with those options for PPPM. However, I’d like to first get Ewald going as I want to compare to Ewald in another code. I can play with the cutoff and this would allow me to change the # of vectors and whatever g_ewald corresponds to. However the cutoff for the electrostatics would also change the short range Coulomb term would it not?

One way I can do it, is just have the other code use the same Ewald parameters as what LAMMPS is using (with the same cutoffs, etc.).

So, then my question becomes, how do the lammps parameters g_ewald and the vectors that are outputted correspond to kappa and kmax (as written in any standard solid state physics book as the manual states).

Does kappa in the Ewald equations (as written in any standard solid state physics book) = g_ewald in ewald.cpp?

Then does kmax (as written in any standard solid state physics book) = the 1d vector as printed in the vector line in the output?

My real question is that I don’t understand what g_ewald means in LAMMPS and what the vectors outputted mean (i.e., what they correspond to in the Ewald equations).

Like I said, with a cutoff of 14.0 angstroms (VDW and Electrostatics) I’m getting 0.3 error between the two codes and with a cutoff of 25.0 angstroms (half the box I'm using) I'm getting 0.02 difference. However, I really just want to (1) make sure that the other code is correct and (2) want to use the same parameters so I can compare results to make sure I am correctly comparing data between the two.

I’d like to compare both codes with the same kappa, kmax, and cutoff. Then I think I should be getting exactly the same result for the total energy (within the precision of conversion factors, etc.)

Thanks again.

Best.

Josh

4

I'll let Paul take a crack at your follow-up questions.

Without digging into the Ewald code,
my recollection is that kmax is strored internally, though it may not
be output. Also, I seem to recall there is a sqrt(2) difference between
Ewald and PPPM in their definition of Gewald, due to some difference
in standard nomenclature, or the way they were initially implemented.
So that might throw you off.

Steve

5

I do understand about the options for PPPM and appreciate how you can
control things with those options for PPPM. However, I'd like to first get
Ewald going as I want to compare to Ewald in another code. I can play with
the cutoff and this would allow me to change the # of vectors and whatever
g_ewald corresponds to. However the cutoff for the electrostatics would
also change the short range Coulomb term would it not?

Yes.

So, then my question becomes, how do the lammps parameters g_ewald and the
vectors that are outputted correspond to kappa and kmax (as written in any
standard solid state physics book as the manual states).

Does kappa in the Ewald equations (as written in any standard solid state
physics book) = g_ewald in ewald.cpp?

Yes. I'd have to look at the same equations you're looking at to be sure. If you look at eqns 4, 5, and 6 of JCP 109, pg 7679, LAMMPS's g_ewald corresponds to the alpha parameter there. I've seen this Ewald parameter written as alpha or kappa typically in the literature.

Then does kmax (as written in any standard solid state physics book) = the
1d vector as printed in the vector line in the output?

Yes. Again, assuming that equations you're looking at are the same as those I'm looking at, kmax is the value you should be using. Looking at Ewald::eik_dot_r(), you'll see that LAMMPS uses kmax as the summation limit in each dimension.

My real question is that I don't understand what g_ewald means in LAMMPS and
what the vectors outputted mean (i.e., what they correspond to in the Ewald
equations).

You'll want to use g_ewald (the Ewald parameter that divides real and reciprocal space) and kmax (the reciprocal space vector length in each dimension) for comparison with other codes.

Like I said, with a cutoff of 14.0 angstroms (VDW and Electrostatics) I'm
getting 0.3 error between the two codes and with a cutoff of 25\.0 angstroms \(half the box I'm using\) I'm getting 0\.02 difference. However,
I really just want to (1) make sure that the other code is correct and (2)
want to use the same parameters so I can compare results to make sure I am
correctly comparing data between the two.

Note that there is a boundary condition energy term that might be different between the two codes. (See eqn 7 of JCP 109, pg 7679.) If your system is not net neutral, you'll have to make sure that these extra terms are the same in both codes. LAMMPS uses the following energy correction for systems that are not net neutral:

    -0.5*PI*qsum*qsum / (g_ewald*g_ewald*volume)

I'd like to compare both codes with the same kappa, kmax, and cutoff. Then
I think I should be getting exactly the same result for the total energy
(within the precision of conversion factors, etc.)

The user usually specifies the cutoff and the precision, from these numbers, LAMMPS computes (kappa = g_ewald) using this simple heuristic:

    g_ewald = (1.35 - 0.15*log(precision))/cutoff;

Alternatively, you can set g_ewald from the input deck using the kspace_modify command.

Then kmax is computed like this:

  // determine kmax
  // function of current box size, precision, G_ewald (short-range cutoff)
  int nkxmx = static_cast<int> ((g_ewald*xprd/PI) * sqrt(-log(precision)));
  int nkymx = static_cast<int> ((g_ewald*yprd/PI) * sqrt(-log(precision)));
  int nkzmx = static_cast<int> ((g_ewald*zprd_slab/PI) * sqrt(-log(precision)));
  int kmax_old = kmax;
  kmax = MAX(nkxmx,nkymx);
  kmax = MAX(kmax,nkzmx);

As you've noted, LAMMPS prints kmax to the screen for you:

    fprintf(screen," vectors: actual 1d max = %d %d %d\n",kcount,kmax,kmax3d);

(It's the second number printed in that line.)

That gives you kappa, kmax, and cutoff, so you should be able to get exact agreement between two codes if these parameters are the same.

Paul

6

Also, I seem to recall there is a sqrt(2) difference between Ewald and PPPM in their definition of > Gewald, due to some difference in standard nomenclature, or the way they were initially
implemented. So that might throw you off.

I think we changed the code a while back (many years ago) to consistetize it and got rid of this sqrt(2) difference between Ewald and PPPM. So I think that they are consistent now.

Paul

7

Dear Paul and Steve,

Thank you so much. This answers all my questions. Obviously I did not look at the code close enough, as Paul clearly points out kmax is clearly written there.

Thank you for the citation as well. I think I will be able to figure it out from that (kappa or alpha).

Thanks again.

Josh

8

Dear Paul and Steve,

Thanks again for your help. I’m sorry for bothering about this question, as now it’s obvious what these parameters are looking at the code (although maybe that’s just because I know that g_ewald = alpha (or kappa) now).

The major problem I discovered was the inclusion of the surface dipole term (which LAMMPS doesn’t have) in the other code. Without this term and using the same parameters, I get 0.0026 agreement in the total energy. For my system the surface dipole term seems to be about 0.3 of the total energy, so it’s probably ok not to include this.

Furthermore, LAMMPS default Ewald parameter scheme seems to be very good, at least for my system, which I’m sure you already knew.

One thing though: You said that you could control g_ewald with kspace_modify. Although LAMMPS accepts this input, it seems to override it with the precision, so it doesn’t change the value used. I think maybe g_ewald can only be specified with kspace_modify with PPPM (which I will probably use now that I know Ewald is good). This is actually what the manual implies too.

Thanks again.

Josh

9

Paul can comment on the surface dipole term.

One thing though: You said that you could control g_ewald with
kspace_modify. Although LAMMPS accepts this input, it seems to override it
with the precision, so it doesn't change the value used. I think maybe
g_ewald can only be specified with kspace_modify with PPPM (which I will
probably use now that I know Ewald is good). This is actually what the
manual implies too.

Paul - this looks to be the case (ewald ignore the user setting for
gewald, only PPPM uses it). This might be as simple a change as
adding "if (!gewaldflag)" in front of this line in ewald.cpp, similar
to what pppm.cpp does:
  g_ewald = (1.35 - 0.15*log(precision))/cutoff;

Do you agree?

Steve

10

Yes, I agree that making that little change to ewald.cpp (and the docs) makes sense. The user should be able to set g_ewald from the input deck.

As far as the surface dipole term goes, LAMMPS uses metallic or "tin foil" boundary conditions, so that term is zero.

Paul

11

Posted a patch 3Feb11 that includes this turn-on of the user
Gewald setting for Ewald.

Steve

12

Thanks. This is great!