Dear all.
I am trying to test the qeq/point fix in LAMMPS by calculating the charges on two atoms forming a dimer. The electronegativities and self-Coulomb repulsions of both elements are defined by the following qeq file:
1 3\.04 5\.5 0 0 0
2 7\.05 13\.96 0 0 0
The script looks as follows (no relaxation, just static distribution of charges):
units metal
boundary p p p
atom\_style charge
region box block \-10 10 \-10 10 \-10 10 units lattice
create\_box 2 box
create\_atoms 1 single \-0\.88 0 0 units lattice
create\_atoms 2 single 0\.88 0 0 units lattice
mass \* 1\.0
group at1 type 1
group at2 type 2
set type 1 charge 0\.1
set type 2 charge \-0\.1
pair\_style coul/long 12\.0
pair\_coeff \* \*
kspace\_style pppm 1e\-4
kspace\_modify compute no
fix 1 all qeq/point 1 12\.0 1e\-4 500 qeq\_point
compute qq1 at1 property/atom q
compute q1 at1 reduce ave c\_qq1
compute qq2 at2 property/atom q
compute q2 at2 reduce ave c\_qq2
variable qtot equal c\_q1\+c\_q2
thermo\_style custom step c\_q1 c\_q2 v\_qtot
run 0
The calculated charges on the two atoms in the dimer areqA = -qB = 0.219.By hand calculation, following the paper of Rappe and Goddard (1991) and requiring charge neutrality, gives:
qA = -qB = (chiB0-chiA0)/(JAA0-2*JAB+JBB0),
where chiA0 = 3.04 eV, chiB0 = 7.05 eV, JAA0 = 5.5 eV/e^2, JBB0 = 13.96 eV/e^2 (see the qeq file above). If I take JAB=14.4/RAB, where RAB = 1.76 A, and substitute all into the equation above, I get qA = 1.29, which is quite different from the value obtained from LAMMPS.
I checked the file fix_qeq_point.cpp, which is probably missing a factor of 14.4 on the line
H.val[m_fill] = 0.5/r;
although JAA0 and JBB0 in the qeq file do contain this factor. This line should be:
H.val[m_fill] = 0.5*14.4/r;
After this, I get the correct charge of 1.29. I tried also qeq/shielded to avoid the Coulomb catastrophe at small RAB, where this factor is correctly in (at the bottom of the file fix_qeq_shielded.cpp). However, it is stored in the variable EV_TO_KCAL_PER_MOL, which is slightly misleading because 1 eV = 23 kcal/mol (http://www.colby.edu/chemistry/PChem/Hartree.html). This variable should be probably renamed to something like ONE_BY_4PI_EPS0, which is 14.4.
Maybe I am wrong or there was a reason for not including the 14.4 factor in fix_qeq_point.cpp and having it in fix_qeq_shielded.cpp. However, this may be quite confusing for many users.
Best,
Roman