Charge equilibration + NVE

Dear all,

I'm currently working with COMB3 and the qeq/comb charge equilibration scheme. If I'm not wrong, the free energy isn't conserved during the charge equilibration process, which leads to incompatibility with simulations in the microcanonical ensemble. Technically, your atoms would progressively freeze, energy leaking out of the system - that is, if the charge gets correctly equilibrated. My question is, are there any workarounds? I can't think of any way to re-inject the estimated energy loss after the qeq that would make any physical sense.

Thanks,
Arthur

Dear all,

I'm currently working with COMB3 and the qeq/comb charge equilibration
scheme. If I'm not wrong, the free energy isn't conserved during the
charge equilibration process, which leads to incompatibility with

True, due to the damping in the damped dynamics scheme that qeq/comb uses.
But the effect is very small, almost negligible.

simulations in the microcanonical ensemble. Technically, your atoms
would progressively freeze, energy leaking out of the system - that is,

Only technically, since the effect is so small that it should not be
noticeable.

if the charge gets correctly equilibrated. My question is, are there any
workarounds? I can't think of any way to re-inject the estimated energy

In fix_qeq_comb,cpp find the variable "heatpq" and make it zero, then you
will have no damping.

Ray

Dear all,

I'm currently working with COMB3 and the qeq/comb charge equilibration
scheme. If I'm not wrong, the free energy isn't conserved during the
charge equilibration process, which leads to incompatibility with

True, due to the damping in the damped dynamics scheme that qeq/comb uses.
But the effect is very small, almost negligible.

this makes me think of born-oppenheimer dynamics with DFT codes. the
origin of the energy loss there is due to the fact that the previously
converged wavefunction is use as initial guess for reoptimizing it to
self-consistency. unless you converge very tightly, you will have a
small residual force pointing towards the previous position. the
strategy to work around this is to use an extrapolation to extract the
initial guess from previous steps. rather then always guessing
"backwards", with a properly tuned extrapolation, you'd get a balanced
distribution, errors will cancel and energy conservation is much
improved, especially with relaxed convergence and thus simulations are
faster.

axel.

Thanks for your responses. Unfortunately, it looks like the effect of the energy not being conserved is not as negligible as thought. I did some testing, using version (11 Jul 2014-ICMS).

The structure is a very simple 2x2x2 supercell of CuO. The simulation consists of 40ps in the NPT ensemble at 1atm/300K, and a 30ps-long run in the micro-canonical ensemble. I've used different values of Nevery for the fix qeq/comb, and also did one run at fixed charges (~+-0.9, using the average values computed with the qeq scheme after an NPT run).

The result is that the kinetic energy leaks out of the system, resulting in a progressive immobilization of the atoms (see ke=f(t).pdf attachment). I've also attached the input structure, script and ff parameters (actually the comb3 parameters from lammps, but only the Cu-O lines) if you want to reproduce the simulation.

By the way, thanks for the workarounds!

Arthur

ke=f(t).pdf|attachment (24.6 KB)

comb3parameters.dat (3.57 KB)

input_f=10.dat|attachment (3.81 KB)

structure.dat (3.67 KB)

The system size may be too small. You can turn off fix qeq/comb damping and try again.

Ray

I’ve compiled a lammps dev build (29 July) without damping (heatpq = 0.0) and gave it a try (the same small supercell, 40ps NPT/30ps NVE, Nevery=10, initial charge values = 0 for each atom). During the run the charge values never converged, and as a consequence I suppose, the system never really relaxed.

I’ve also tried to run the previous simulation (that is with damping, NPT and NVE with charge equilibration and Nevery=10) with a much bigger sample (8x8x8 supercell). The variation of the total energy is indeed negligible. I don’t understand the influence of the size of the sample, though. How do you explain it? Is there some kind of cascade effect when the sample is small enough, due to the order of magnitude of the energy loss being “close enough” to the one of the total energy of the system?

In the next days I’ll try without damping, on a big system.

Thanks,
Arthur

Hi Arthur, I thought it might be useful to others if I communicate through this forum, instead of having private communications that we had had.
You wrote: “I’ve also tried to run the previous simulation (that is with damping, NPT and NVE with charge equilibration and Nevery=10) with a much bigger sample (8x8x8 supercell). The variation of the total energy is indeed negligible.” Could you pls provide the structure file, and the lammps input file. I am really surprised and curious that I get so different result for rutile! Like you, I also used the comb potential that came with lammps distribution.
Cheers!
– Abrar

I’ve compiled a lammps dev build (29 July) without damping (heatpq = 0.0) and gave it a try (the same small supercell, 40ps NPT/30ps NVE, Nevery=10, initial charge values = 0 for each atom). During the run the charge values never converged, and as a consequence I suppose, the system never really relaxed.

I’ve also tried to run the previous simulation (that is with damping, NPT and NVE with charge equilibration and Nevery=10) with a much bigger sample (8x8x8 supercell). The variation of the total energy is indeed negligible. I don’t understand the influence of the size of the sample, though. How do you explain it? Is there some kind of cascade effect when the sample is small enough, due to the order of magnitude of the energy loss being “close enough” to the one of the total energy of the system?

In the next days I’ll try without damping, on a big system.

Thanks,
Arthur