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.

Thanks,

Keith