Charge equilibration convergence for COMB potential when initial charge not set

Lammps version: lammps-11Mar13

I have been testing the COMB potential to describe adsorption of oxygen atoms on a Ti surface. Ultimately I would like to describe the deposition of oxygen onto the Ti surface using MD.

For the bare Ti surface I set up a 2D slab and was able to optimise the structure using COMB and charge equilibration. The final charges on Ti ions are close to zero (as expected).

However, when I perform a calculation on the same slab with an oxygen atom placed above the surface the optimisation fails. Looking into the charge equilibration iterations I see that it is not converging, e.g.

Charge equilibration on step 1
   iteration: 0, enegtot 4.05076, enegmax 1373.58, fq deviation: 7.60988
   iteration: 1, enegtot -3365.83, enegmax 1.21214e+06, fq deviation: 6715.48
   iteration: 2, enegtot 9.47103e+11, enegmax 1.60328e+15, fq deviation: 8.88244e+12
   iteration: 3, enegtot -nan, enegmax nan, fq deviation: nan
   iteration: 4, enegtot -nan, enegmax nan, fq deviation: nan
   iteration: 5, enegtot -nan, enegmax nan, fq deviation: nan

Even if I put the O atom far above the surface of Ti (14 Angstroms) I find I still need to set the initial charge in order to get the charge equilibration to converge. However, I discovered that if I set the initial charge of the oxygen atom to -1 and put some positive charge on the remaining Ti ions (so that overall neutral) I can get it to converge, e.g.

Charge equilibration on step 1
   iteration: 0, enegtot 0.437564, enegmax 3.49434, fq deviation: 0.186271
   iteration: 1, enegtot 0.443535, enegmax 3.21395, fq deviation: 0.163696
   iteration: 2, enegtot 0.454393, enegmax 2.71604, fq deviation: 0.125754
   iteration: 3, enegtot 0.468788, enegmax 2.08327, fq deviation: 0.0987725
   iteration: 4, enegtot 0.485291, enegmax 1.40373, fq deviation: 0.0940399
   iteration: 5, enegtot 0.502544, enegmax 1.09408, fq deviation: 0.101073
   iteration: 6, enegtot 0.519354, enegmax 1.01511, fq deviation: 0.11542
   iteration: 7, enegtot 0.534752, enegmax 0.963302, fq deviation: 0.115947

Should this behaviour be expected? It poses a problem for my original aim as I would like to introduce neutral O atoms and simulate their interaction with the surface using MD. Any suggestions would be appreciated.

Hi Keith,

The default charge equilibration step size may be too big for your
system. You can change a variable in the source code to reduce the
step size so that charge equilibrations do not go too wild - in the
same manner as reducing timesteps for MD runs. Around line 202 of
fix_qeq_comb.cpp, please reduce the "dtq = 0.01" to 0.002.