total energy is not conserved in NVE + qeq/comb

LAMMPS LAMMPS Mailing List Mirror
1

Hi all, this has been reported earlier by Arthur France-Lanord. The issue is total energy does not remain conserved in NVE runs with COMB3 potential while charge equilibration is on (NVE+qeq/comb). The total energy decreases progressively, and so does the temperature.
Pls find attached two plots of total energy and temperature (vs time) from a NVE run with rutile (TiO2). I used the comb3 potential that is provided in lammps distribution (stable version of 1 Feb 2014). Pls also find attached the files required to run this on lammps (in.*, structure, potential: ffield.comb3, taken from lammps distribution).
The rutile structure file has been inspected carefully, and is hopefully correct. Per Tao Liang’s suggestion, the system size has been so chosen that each dimension is equal or greater than twice of the cutoff length for Coulomb interaction set in comb3 (which is 11 A, the 35th element in any entry in ffield.comb3, and read as params[nparams].coulcut in pair_comb3.cpp). The lammps input file is a modified version of one I received from Ray Shan earlier.
Is there anything in the mathematics of charge equilibration (Rappe and Goodard 1991, Rick and Stuart 1994) that bars the total energy from being conserved in NVE in comb3? I was curious to know what is the scenario in case of qeq/reax (which uses Nakano 1997). I need to build lammps with relevant packages to check that. Meanwhile, any comment on all these will be much appreciated.
Best regards.
– Abrar

rutile_total_energy_NVE_qeq.gif

rutile_temperature_NVE_qeq.gif

in.rutile_NVE (814 Bytes)

ffield.comb3 (324 KB)

data.rutile_5x5x8 (90.7 KB)

2

Hi all, this has been reported earlier by Arthur France-Lanord. The issue is
total energy does not remain conserved in NVE runs with COMB3 potential
while charge equilibration is on (NVE+qeq/comb). The total energy decreases
progressively, and so does the temperature.

[...]

Is there anything in the mathematics of charge equilibration (Rappe and
Goodard 1991, Rick and Stuart 1994) that bars the total energy from being
conserved in NVE in comb3? I was curious to know what is the scenario in
case of qeq/reax (which uses Nakano 1997). I need to build lammps with
relevant packages to check that. Meanwhile, any comment on all these will be
much appreciated.

please note that i have no practical experience with COMB and COMB3,
but your report reminds me a lot of the problems that people in the ab
initio MD community have with doing born-oppenheimer dynamics (as
opposed to car-parrinello MD). the reason for the issue is that one
uses the wavefunction from the previous step as initial guess for the
convergence of the wavefunction. while this reduces the number of scf
cycles to reach a typical convergence, it also adds a small bias
resulting in a force in the exact opposite direction of the motion. to
minimize this force, one has to converge the wavefunction extremely
tightly, and thus making this kind of MD very costly.

that is, until people developed extrapolation scheme like ASPC (
http://dx.doi.org/10.1002/jcc.10385 ) where the initial guess is
extrapolated from the last few time steps. in the ideal case the guess
leads to about a gaussian distribution of the initial guess around the
final, converged result.
this leads to better energy conservation and faster simulations, since
the initial guess is generally even better than reusing the previous
data. the only price to pay, is additional storage for multiple
wavefunctions. the same principle applies to models with
polarizability and should be similarly transferable to charge
equilibration.

consequently, without extrapolation, you have two options: 1) use a
tighter convergence on the charge equilibration, that should lead to a
more accurate trajectory, less energy loss and most computational
effort, 2) reduce the frequency of charge equilibration. that will
make it less accurate, but also less time consuming and reduces the
unphysical biasing force.

axel.

p.s.: the implementation of ASPC in LAMMPS should be fairly straight
forward (best as a generic fix that could be used by multiple other
fixes). i currently don't have the time to do it myself, but if
anybody feels comfortable hacking LAMMPS and needs a project, i am
happy to give a few pointers.

3

> Hi all, this has been reported earlier by Arthur France-Lanord. The
issue is
> total energy does not remain conserved in NVE runs with COMB3 potential
> while charge equilibration is on (NVE+qeq/comb). The total energy
decreases
> progressively, and so does the temperature.

It is a well-known issue with the extended Lagragian QEq method (see Rick,
Stuart and Berne) that total energy is not conserved. It has to do with
the formulation, particularly the damping to keep the electronic degree of
freedom at sufficiently low temperature (a BO-MD concept), and has nothing
to do with COMB nor COMB3.

[...]

> Is there anything in the mathematics of charge equilibration (Rappe and
> Goodard 1991, Rick and Stuart 1994) that bars the total energy from being
> conserved in NVE in comb3?

Again, not COMB3, but all damped dynamics (aka extended Lagragian) methods
- see the Rick paper for more details. If you solve charges using the
matrix inversion method (Rappe & Goddard) then you don't have damping and
energy is better conserved.

I was curious to know what is the scenario in
> case of qeq/reax (which uses Nakano 1997). I need to build lammps with
> relevant packages to check that. Meanwhile, any comment on all these
will be
> much appreciated.

You don't really need ReaxFF (pair_style reax or reax/c) to see this
effect. You can observe the total energy drift with qeq/dynamic (extended
Lagragian method) in /examples/qeq. In the attached pdf files, you can see
with qeq/dynamic total energy drift is worse than the cases with no qeq and
qeq/point, and qeq/reax.

So to sum it up, with damped dynamics (fix qeq/dynamic and fix qeq/comb)
you will always see total energy drift and that is just the way it is.

Ray

qeqreax.pdf (19.1 KB)

qeqbuck.pdf (19.4 KB)

4

Hi Axel and Ray,
thanks again for your committed support.
I am trying the first two suggestions given by Axel: using higher tolerance, and changing frequency of qeq. I had tried with higher tolerance before with FeO (13824 atoms, each dimension > 22 A), and found out that changes in tolerance did not affect the total energy at all! In the FeO runs, the values of tolerances were 1e-4 and 1e-5. For the qeq calculations to converge to higher tolerance, I had to increase the value of ‘loopmax’ in fix_qeq_comb.cpp (in the MD steps, the charges converged in little over 100 iterations). I’ll soon update the results from the ensuing rutile runs.
Thanks Ray for providing the two plots. Are these from lammps examples? I couldn’t find any directory /examples/qeq. Also I grep-ed with qeq/dynamic and qeq/point from top of /examples, but it did not yield any file.
In the plot ‘Buck example’, I see the total energy is increasing, as opposed to the rutile COMB3 case where it was decreasing. In ‘ReaxFF example’, the total energy is kind of oscillating.
I have a question about the plots of total energy with ‘No qeq’. For ReaxFF, it looks pretty good, however, for Buck example, I cannot be sure since the energy scale is too big. I notice an initial drop of several eV (I assume the energy is in eV in the plots). When I ran rutile NVE without qeq, there was an initial drop of 2 eV in the total energy (pls find attached a plot). To provide more specifics, there was one initial step of qeq followed by minimization with box/relaxation, and then NVE run was ensued (I followed Ray’s input file). In the attached graph, I have plotted total energy only from the NVE part. Pls find attached the relevant lammps files for this run (rutile NVE with no qeq).
Thanks again for your time.
– Abrar

rutile_energy_NVE_no_qeq.gif

in.rutile_NVE (813 Bytes)

ffield.comb3 (324 KB)

data.rutile_5x5x8 (90.7 KB)

5

Hi Axel and Ray,
thanks again for your committed support.
I am trying the first two suggestions given by Axel: using higher
tolerance, and changing frequency of qeq. I had tried with higher tolerance
before with FeO (13824 atoms, each dimension > 22 A), and found out that
changes in tolerance did not affect the total energy at all! In the FeO
runs, the values of tolerances were 1e-4 and 1e-5. For the qeq calculations
to converge to higher tolerance, I had to increase the value of 'loopmax'
in fix_qeq_comb.cpp (in the MD steps, the charges converged in little over
100 iterations). I'll soon update the results from the ensuing rutile runs.
Thanks Ray for providing the two plots. Are these from lammps examples? I
couldn't find any directory /examples/qeq.

In $YOUR_LAMMPS_DIR/examples/qeq after the 8Sep14 version.

Also I grep-ed with qeq/dynamic and qeq/point from top of /examples, but
it did not yield any file.

They are fix qeq/variants. See http://lammps.sandia.gov/doc/fix_qeq.html for
more details.

Ray

6

Ah I see. As I mentioned in my first mail, I have been using the stable version from 1 Feb 2014. Yes, I understood they are from qeq/variants, but since the directory /qeq is absent in the Feb version, I didn’t find any example where qeq/dynamic or qeq/point had been used. I have just downloaded the development version of 9 Oct, where I see the directory /qeq, and usage of the above fixes.
Anyway, now I have the results from the NVE rutile run with Nevery = 10 (upto 10 ps), and precision = 1e-4 (upto 2.5 ps) respectively. They do not show any improvement in energy conservation. For Nevery=10, it has indeed gone worse. Pls find attached the graphs, and the input files. The structure file is the same I used before, and so is the potential file (ffield.comb3). For convenience, I have attached them too.
In my previous mail, I informed with precision 1e-4 and 1e-5, I did not find any difference in total energy in case of FeO. Now, with precision 1e-3 and 1e-4, we see some difference in total energy in case of rutile. I’ve just started a rutile run with precision 1e-5, and will check how much it differs from 1e-4 run.
Best regards.
– Abrar

Hi Axel and Ray,
thanks again for your committed support.
I am trying the first two suggestions given by Axel: using higher tolerance, and changing frequency of qeq. I had tried with higher tolerance before with FeO (13824 atoms, each dimension > 22 A), and found out that changes in tolerance did not affect the total energy at all! In the FeO runs, the values of tolerances were 1e-4 and 1e-5. For the qeq calculations to converge to higher tolerance, I had to increase the value of ‘loopmax’ in fix_qeq_comb.cpp (in the MD steps, the charges converged in little over 100 iterations). I’ll soon update the results from the ensuing rutile runs.
Thanks Ray for providing the two plots. Are these from lammps examples? I couldn’t find any directory /examples/qeq.

In $YOUR_LAMMPS_DIR/examples/qeq after the 8Sep14 version.

Also I grep-ed with qeq/dynamic and qeq/point from top of /examples, but it did not yield any file.

They are fix qeq/variants. See http://lammps.sandia.gov/doc/fix_qeq.html for more details.

Ray

rutile_NVE_energy_different_qeq_frequency.gif

rutile_NVE_energy_different_precision.gif

in.rutile_NVE_hightolerance (814 Bytes)

in.rutile_NVE_Nevery10 (815 Bytes)

data.rutile_5x5x8 (90.7 KB)

ffield.comb3 (324 KB)