reaction field term and Kirkwood approximation

Dear lammps users,

1) Just wanted to ask if the following reaction field term in Coulomb interactions
is implemented already:

(Erf-1)/(2*Erf+1) Rij^2/Rc^3

I didn't notice it in the pair forces, but i might have overlooked something.

2) If it is not, can you please tell me in what .cpp file the
'dielectric ...' command is read from the input script?
I guess the right way to add this would be to add another variable in force.h
(similar to dielectric)
and write a pair style based on
pair_lj_cut_coul_cut.cpp

3) If someone is using it, do you know any papers where the motivation
for that term is described? Probably that's a common knowledge because
I have several papers in which this term is used and there is no reference to any
of Kirkwood's papers or explanation why it is there.

Thank you in advance for your help.

Sincerely yours,
Denis.

Dear lammps users,

1) Just wanted to ask if the following reaction field term in Coulomb
interactions
is implemented already:

(Erf-1)/(2*Erf+1) Rij^2/Rc^3

I didn't notice it in the pair forces, but i might have overlooked
something.

no. each pair style implements these separately.
i had a reaction field coulomb implemented in a
haphazard way a long time ago. i didn't find it
worth the effort. the key problem of cutting of
long-range interactions for an inhomogeneous
system is not solved by any of these approaches.

i strongly recommend to first implement a variant
of pair style coul/cut and then use pair style
hybrid/overlay to add reactionfield coulomb
to any nonbonded pair style. that will reduce
effort and increase flexibility.

but then again, you may want to have a look at
the coul/wolf pair styles, which seem to be a
better approach to the same idea.

2) If it is not, can you please tell me in what .cpp file the
'dielectric ...' command is read from the input script?

input.cpp

grep is your friend!

I guess the right way to add this would be to add another variable in
force.h

there is nothing to add.

(similar to dielectric)
and write a pair style based on
pair_lj_cut_coul_cut.cpp

see above.

3) If someone is using it, do you know any papers where the motivation
for that term is described? Probably that's a common knowledge because
I have several papers in which this term is used and there is no
reference to any
of Kirkwood's papers or explanation why it is there.

reaction field is a fairly old concept and at the time it was
competing with plain cutoff coulomb on system. the simplest
way to describe it is that you model the coulomb screening
from the environment (i.e. the periodic cells) as a continuum
dielectric. the result is a smooth decay of the coulomb interactions.

however, unlike in the wolf sum, you don't pay attention to
the residual multipole after the screening (programs like
CHARMM and GROMOS and derivatives of them, countered
that, but applying cutoffs and PBC not to atoms, but to charge
groups, that were usually chosen to be neutral).

you probably want to look up stuff from the mid-to-late 1980's,
when people started "bio" simulations in explicit water
and packages like charmm, gromos and so on were
born together with the force fields and the whole shebang.

once computers got more powerful, this was all superseded
by ewald summation and then PME/SPME/PPPM.

hope this helps,
      axel.