Triclinic PPPM

I have some LAMMPS code for triclinic PPPM, MSM, and Ewald. If anyone wants to be a beta tester, shoot me an email and I will send you a copy.

Thanks,

Stan

Dear LAMMPS users and developers,

Last year, I gave a Master mini-project on the computation of quartz elastic moduli using LAMMPS (more specifically using Aidan Thompson’s example scripts).
At the end of the mini-project, using files included in elastic-quartz-Backup.7z and the windows version of LAMMPS dated July 1st, 2012, we got results that were more or less in agreement with experimental values:
C11 C12 C13=C23 C14=-C24=C56 C33 C44=C55 C66 C15 C16
experiment 86.2 7.0 11.8 -18.0 105.4 58.2 39.9 (0) (0)
Calculation (cutoff=8.8) 86.9 10.2 16.1 16.3 109.2 48.6 37.8 -4.7 -10.9
Calculation (cutoff=8.9) 85.9 9.25 15.3 16.8 106.8 48.1 38.3 -4.85 -11.1

I have tried to reproduce these results with the June 7th, 2013 version of LAMMPS using ewald/disp instead of ewald/n since the later does not exist any-more (files in elastic-quartz-new.7z).

The problem is that I get very incorrect values if I use standard version of the pair style potential (be it with my Linux g++ or my Intel+MKL versions of LAMMPS), whereas I do get the same results as before with both binaries if I just add -sf omp on the command line and do not change a single character in the input files!

I understand that OMP implementations are not the same as the standard implementation, so that I wonder whether there could be a bug in the standard ewald/disp or buck/coul/long triclinic implementation that would not exist in the OMP implementation ???

The script following my signature just tries to optimize the geometry of the unit cell with minimize commands.
It reproduces the problem with only one input file (but it is not as easy to see which are the good results as for the elastic moduli).
In fact, it is even worse since with this same unchanged input, I get quite different results depending on the binary that I use (serial, MKL, g++ or openmpi) both with -sf omp and without!

PS1: I find that the sensitivity to the cutoff value in the pair_style buck/coul/long line is very high: it seems that values below 8.8 or higher than 9.0 give wrong results. Would you have any idea of the origin of this sensitivity to the cutoff ?

PS2: Any idea as for why C15 and C16 are not zero as they should, due to the 32 symmetry?

PS3: I have had to use the 7-zip compress format (available from sourceforge) since our university mailer blocks .zip and .tgz attachments

elastic-quartz-Backup.7z (50.2 KB)

elastic-quartz-new.7z (21.4 KB)

michel,

i have three comments, but not a lot of help to offer.
1) ewald/disp used to be called ewald/n it was renamed in an effort to
have more consistent naming conventions (we had a couple of those
during the last year or so).
2) ewald/disp was very recently upgraded to support point dipoles,
perhaps something got broken in the process.
3) there is no /omp version of ewald/disp (yet), so that hints at a
bug that got fixed by accident when writing the /omp version, or a bug
that is present in both variants, but doesn't show in the second. i
will try to take a look.

axel.

p.s.: have you tried using pppm? this should support triclinic cells now.

Dear Axel, Thanks for your interest in my problem. I have indeed tried PPPM, but - for pppm/disp, either with or without -sf omp, I get PPPMDisp initialization … ERROR: KSpace style does not yet support triclinic geometries (…/kspace.cpp:135) - for pppm with -sf omp, I get PPPM initialization … ERROR: KSpace style does not yet support triclinic geometries (…/kspace.cpp:135) - for ppm without -sf omp, I get wrong results for the elastic constants (after a much longer time than for ewald/disp).

Dear LAMMPS users and developers,

[...]

PS1: I find that the sensitivity to the cutoff value in the pair_style
buck/coul/long line is very high: it seems that values below 8.8 or higher
than 9.0 give wrong results. Would you have any idea of the origin of this
sensitivity to the cutoff ?

please try with tighter accuracy for the ewald sum. 1.0e-4 is rather
sloppy. it is "ok" for running MD, where there is a lot of error
cancellation, but may not be good enough for your project.

i've also noticed that the estimator algorithm makes some assumptions
about the distribution of charges which may not always apply (e.g. it
doesn't work well for some of the coarse grain MD tests, that i have
run), so may want to converge kspace sum and alpha manually in the
old-fashioned way. :wink:

axel.

I tried with down to 1e-10 (especially for the single file for geometry optimization), but that did not really change things on the limited number of tests that I did. I came back to the value of 1E-4 in the files that I sent, because it was the value used by my student (who used a relatively slow laptop…). Thanks for the advice. It looks like a mini-project subject for next semester! :wink: It seems, from the archives of this list, that there have been other people that used LAMMPS with BKS potential for alpha-quartz. I call on them to tell us whether they have noticed the same phenomenon and maybe even how they solved it! :wink:

Stan and/or Rolf may want to take a look at this, as it indicates
some issue with ewald/disp. However, the only reason to use
any of the “disp” kspace styles is if you want to include long-range
LJ interactions, which your script appears not to. So you can
just try ewald or pppm, both of which can do triclinic boxes.

Steve

Hi,

I’ve just read through the problem and my suspicion is that the problem is not
in the ewald/disp part of the code:

If I understood the inquiry right, then simulation values come out correctly
when the -sf omp option is used in the command line, but incorrectly if
this option is not used. Thus the problem is likely somewhere in the parts
of the code that is replaced by an corresponding omp version. As the ewald/disp/omp
does not exist, the bug is probably somewhere else, possibly in the pair style.

Best,

Rolf

Hasen’t Paul recently added a table for the Coulomb potential to
the buck/lcoul/long pair style? I’ve seen that this table is not included
in the OMP version yet. Maybe the bug was introduced when adding
the table for Coulomb interactions.

Michel: Could you please try to check whether adding the command

pair_modify table 0

will solve the problem?

Best,

Rolf

Thanks for the answer. Since the Buckingham part of the Buck/coul/long potential includes a 1/r^6 term, it seems to me I need the disp, don’t I? Anyway, I have already tried without the disp and I got either wrong result for moderate values of the kspace_style parameter or “ERROR: Triclinic box skew is too large (…/domain.cpp:158)” for value such as 1E-10 (in spite of the fact that I now have “fix 3 all box/relax tri 0.0” thanks to your previous advice).

Hasen't Paul recently added a table for the Coulomb potential to
the buck/lcoul/long pair style? I've seen that this table is not included
in the OMP version yet. Maybe the bug was introduced when adding
the table for Coulomb interactions.

Michel: Could you please try to check whether adding the command

pair_modify table 0

will solve the problem?

yes. that makes the /omp and the regular version agree (again).

tabulation is done in single precision and thus introduces a bit of an
error (or a bit more).

axel.

Stan and/or Rolf may want to take a look at this, as it indicates
some issue with ewald/disp. However, the only reason to use
any of the "disp" kspace styles is if you want to include long-range
LJ interactions, which your script appears not to. So you can
just try ewald or pppm, both of which can do triclinic boxes.

Thanks for the answer.

Since the Buckingham part of the Buck/coul/long potential includes a 1/r^6
term, it seems to me I need the disp, don't I?

no. that is only required/used if you choose buck/long/coul/long as pair style.

axel.

I had a little trouble testing that since I could not find the syntax of this buck/long/coul/long on Fortunately, Stan kindly helped me (e.g. pair_style buck/long/coul/long long long 8.85) and I could make some tests. I have still not found at the moment a combination of pair_style cutoff and kspace_style accuracy that would give me elastic moduli close to the experimental ones (except for C15 and C16). However, I note that in the original van Beest, Kramer, van Santen (BKS) Phys. Rev. Lett., it is mentioned that: “With all force fields derived in this way, we have performed unrestrained lattice minimizations of alpha quartz. It has been found that the energy-minimized structure does not always have the alpha-quartz symmetry.” and “Surprisingly, the best description of the cluster data (lowest chi^2) is given by a set of parameters that yield the wrong low-temperature symmetry of quartz. Along with the symmetry change, the values of elastic constants change dramatically.” which may give an answer to PS1 and PS2 of my first mail… It seems to me I once read that other potentials have been parametrized for quartz after 1990. I think I will try to have a look at these! :wink:

Just a clarification (more for others than for you Michel, since
you’ve already figured it out), but the doc/pair_buck.html
page is for all the Buckingham variants with a cutoff for the exp()
and 1/r^6 terms. The doc/pair_buck_long.html page is
for the styles with the dispersion option to make 1/r^6 also be
long-range. This would be more clear I suppose if we called them all
buck/cut or buck/long. But the vanilla buck pre-dated buck/long.

Steve