Dear LAMMPS users,
I have been playing with the new fix/qEq option and a Buckingham potential for CeO2. However, I am getting some pretty weird results which make me think that maybe there is something wrong with the new fix. At the end of this email you can find the input file and the values for the qfile.
Here are my problems. When I perform a short MD simulation, I observe the following behavior:
-
The energy, lattice parameter and other parameters depend a lot (20% change!) on whether I use qeq/dynamic or qeq/point. It is my understanding that these should yield comparable results, so I am surprised.
-
I also see that the results of my simulation change a lot if I slightly modify the qeq cut-off (say from 10.5 to 10). Again, that is surprising, because I would expect the results to be quite converged for large cut-offs
-
Related to this, I have also tried fitting a potential (including the qeq parameters) to DFT data. The behavior of this qeq potential is quite weird again. Indeed, I get a much better agreement with DFT using a simple Buckingham potential with partial charges than using qeq. Considering that the qeq model has more parameters, I would expect the opposite, i.e. a better agreement or at most the same as simple Buckingham…
Am I doing something wrong or is this an issue with the fix?
Thanks!
Dario
Dear LAMMPS users,
I have been playing with the new fix/qEq option and a Buckingham potential
for CeO2. However, I am getting some pretty weird results which make me
think that maybe there is something wrong with the new fix. At the end of
this email you can find the input file and the values for the qfile.
Here are my problems. When I perform a short MD simulation, I observe the
following behavior:
1) The energy, lattice parameter and other parameters depend a lot (20%
change!) on whether I use qeq/dynamic or qeq/point. It is my understanding
that these should yield comparable results, so I am surprised.
You might want to update your LAMMPS version since there has been some
modifications to the fix qeq/* commands. I got almost identical results
using the two:
Step Temp Volume Press TotEng Lx CPU q1 q2
Dynamic: 1000 3.5722215 13607.373 -2.9574365 -994.57052
23.873976 3.6343489 -0.76570673 0.38285336
Point: 1000 4.0631692 13611.254 15.799083 -994.33882
23.876246 3.5666311 -0.76570426 0.38285213
2) I also see that the results of my simulation change a lot if I slightly
modify the qeq cut-off (say from 10.5 to 10). Again, that is surprising,
because I would expect the results to be quite converged for large cut-offs
Your unitcell length is 5.41 angstroms, so cutoffs of 10 and 10.5 can not
be considered as "large". You should extend your cutoff to at least ~3
times the unitcell size.
3) Related to this, I have also tried fitting a potential (including the
qeq parameters) to DFT data. The behavior of this qeq potential is quite
weird again. Indeed, I get a much better agreement with DFT using a simple
Buckingham potential with partial charges than using qeq. Considering that
the qeq model has more parameters, I would expect the opposite, i.e. a
better agreement or at most the same as simple Buckingham...
Your charges flip signs with your QEq parameters. Did you notice your Ce
changed from +4 to -0.8 and O from -2 to +0.4? This indicates your QEq
parameters are bad, which could be the reason you don't get good results
after using fix qeq.
Ray
Hi Ray,
Thanks for the quick reply. There was an issue in the input file (I had mixed up Ce and O). The attached one is good. Still, I get some problems.
I am using the latest beta version (21 Jan 2015). Which one are you using? In my case the energies are quite different between the two cases (point vs dynamic):
Dynamic: 1000 85.873085 12437.058 -619.73641 -1426.6116 23.168925 68.031887
Point: 1000 346.76277 12851.969 431.31735 -1489.2256 23.423757 40.307313
Regarding the cut-off, you were right! I have run the same calculation with different cut-offs, going from 10 to 20, in steps of 1. Apparently, once the cut-off is larger than ~3 unit cells, the energy converges to exactly the same. Incidentally, is this cut-off limited to 1/2 the box length? Does it have to be the same as the Coulomb cut-off?
1000 346.76277 12851.969 431.31735 -1489.2256 23.423757 37.042404
1000 434.80234 13237.311 -1037.2001 -1554.6161 23.655561 41.725277
1000 407.99734 12432.976 2067.5925 -1745.5004 23.16639 46.863295
1000 443.16304 13175.728 1546.8717 -1490.2356 23.61882 39.016055
1000 408.59714 13221.345 38.263072 -1512.4118 23.646047 35.812034
1000 408.59714 13221.345 38.263072 -1512.4118 23.646047 35.69181
1000 408.59714 13221.345 38.263072 -1512.4118 23.646047 35.633128
1000 408.59714 13221.345 38.263072 -1512.4118 23.646047 35.503157
1000 408.59714 13221.345 38.263072 -1512.4118 23.646047 35.749284
1000 408.59714 13221.345 38.263072 -1512.4118 23.646047 35.942666
1000 408.59714 13221.345 38.263072 -1512.4118 23.646047 35.921411
Thanks!
Dario
Hi Ray,
Thanks for the quick reply. There was an issue in the input file (I had
mixed up Ce and O). The attached one is good. Still, I get some problems.
I am using the latest beta version (21 Jan 2015). Which one are you using?
In my case the energies are quite different between the two cases (point vs
dynamic):
Dynamic: 1000 85.873085 12437.058 -619.73641 -1426.6116
23.168925 68.031887
Point: 1000 346.76277 12851.969 431.31735 -1489.2256
23.423757 40.307313
I am using 21Jan15 version. As said in the doc page, qeq/point and
qeq/dynamic give comparable answers if they are both converged. I do not
know what your qeq/dynamic settings are, but they should give similar
results. Please check if the charges are converged.
Regarding the cut-off, you were right! I have run the same calculation
with different cut-offs, going from 10 to 20, in steps of 1. Apparently,
once the cut-off is larger than ~3 unit cells, the energy converges to
exactly the same. Incidentally, is this cut-off limited to 1/2 the box
length? Does it have to be the same as the Coulomb cut-off?
It is not relavant to 1/2 the box size since fix qeq/* use the LAMMPS ghost
convention. Fix qeq/* build their own neighbor lists so it can be
different from the Coulomb cutoff.
Ray
I am using the following settings:
fix 1 all qeq/dynamic 1 20.0 1.0e-4 2000 qfile
100 2.3374669 11261.93 -19342.332 -2356.7843 22.414941 9.87058
and
fix 1 all qeq/point 1 20.0 1.0e-7 2000 qfile
100 242.46371 11578.24 2066.188 -2122.12 22.622861 5.827913
As you can see the energies are quite different (11% difference). In both cases the tolerance is much smaller than suggested in the manual and the results are clearly converged (i.e. if I change the tolerance for qeq/point or qeq/dynamic, the energy does not change).
Is there anything else that I should be aware of? Some incompatibility between different options?
Thanks!
Dario
I still get the same results. I’d suggest you re-make your LAMMPS executable, make sure it’s linked to the right executable, and try again. Otherwise I can not help you if I can not reproduce what you see. All in all, qeq/point and qeq/dynamic give comparable results if no errors were made.
Ray
Hi Ray,
So, I have tried what you suggested, I recompiled LAMMPS and also tried the 9Dec14 version. I then used the examples in examples/qeq to see if I get the same numbers. I DO NOT! My energies and charges are completely different (see below)!
My LAMMPS
Step PotEng q1 q2 qtot
0 -4203.3633 0.96050687 -0.48025343 -1.7053026e-13
10 -4242.2467 0.96522386 -0.48261193 -1.7053026e-13
From log file in example/qeq
Step PotEng q1 q2 qtot
0 -15435.276 0.85155361 -0.4257768 1.9326762e-12
10 -15450.727 0.85448844 -0.42724422 1.6484591e-12
This is kind of crazy. I tested the executable with some other examples/potentials (LJ, eam) and it works fine.
I guess that begs the question. Is there anything special that you have to do when compiling LAMMPS with the QEQ package? The on-line manual has no instructions, so maybe I am missing something? Do I need another package?
Thanks!
Dario
The log files in the example directory just needed to be updated.
Please download the lasted development version of LAMMPS, try with your example script with my slight modifications (pasted below), and you should see comparable results from qeq/point and qeq/dynamic.
Thanks,
Ray