Contributing a fix for ASPC polarizability treatment

Dear LAMMPS developers, and users,

sorry for the long mail.

I would like to contribute a fix to the LAMMPS code. The fix can treat
polarizable molecules using the time-reversible always-stable predictor
corrector algorithm by Jiri Kolafa (https://dx.doi.org/10.1002/jcc.10385).

I have finished the implementation for Drude-polarizability, which I have
verified against Jiri Kolafa’s original implementation (thanks to Jiri), and
used it quite extensively over the past months.
I would like to extend the fix to inducible point-dipoles and have started some
work in that direction (i.e., test-implementation for non-periodic systems).

Because I would like to see the fix being integrated into LAMMPS, I’d like to
discuss possible design issues.

For both Drude- and inducible point-dipoles, the electric field at the point of
the Drude particle/dipole is needed. For Drude that can be simply achieved by
calculating the electric force and dividing by the charge of the Drude
particle. This is what I am currently doing. Things get more complicated for
dipolar interactions, and one needs to explicitly calculate the field.
Also, it’s (currently) tricky if the particles have non-zero vdW interactions
(which Drude particles by definition don’t have), because electric and vdW
forces are calculated together as pair-forces.

So, my question is: Would it make sense to implement member functions to the
pair and kspace classes that explicitly calculate the electric field at the
positions of the particles?
If so, would you rather create another array for storing the values of the
electric field, or “re-use” the force array (might be a bit unpredictable in
the long run, but introduces a memory overhead if you don’t calculate them(?)).

The implementation of a corresponding, let’s call her “compute_efield”,
member function would also be helpful for implementing constant potential
methods (that is something else people in our group are looking into).

Another question (cry for help?) would be if someone familiar with the dipolar
Ewald summation code would be willing to help creating a separate compute_efield
function for that!? I can probably do it, but it might take me much longer than
someone familiar with that part of the code.

All my preliminary code is located here:

https://github.com/fuulish/ASPC

The master branch contains the Drude dipole implementation, and there is a
branch dipole which contains the test-implementation of inducible point-dipoles.

Looking forward to your response and thanks in advance,

Frank

Dear LAMMPS developers, and users,

sorry for the long mail.

I would like to contribute a fix to the LAMMPS code. The fix can treat
polarizable molecules using the time-reversible always-stable predictor
corrector algorithm by Jiri Kolafa (https://dx.doi.org/10.1002/jcc.10385).

I have finished the implementation for Drude-polarizability, which I have
verified against Jiri Kolafa's original implementation (thanks to Jiri),
and
used it quite extensively over the past months.
I would like to extend the fix to inducible point-dipoles and have started
some
work in that direction (i.e., test-implementation for non-periodic
systems).

Because I would like to see the fix being integrated into LAMMPS, I'd like
to
discuss possible design issues.

For both Drude- and inducible point-dipoles, the electric field at the
point of
the Drude particle/dipole is needed. For Drude that can be simply achieved
by
calculating the electric force and dividing by the charge of the Drude
particle. This is what I am currently doing. Things get more complicated
for
dipolar interactions, and one needs to explicitly calculate the field.
Also, it's (currently) tricky if the particles have non-zero vdW
interactions
(which Drude particles by definition don't have), because electric and vdW
forces are calculated together as pair-forces.

So, my question is: Would it make sense to implement member functions to
the
pair and kspace classes that explicitly calculate the electric field at the
positions of the particles?

​hmm... would it be possible to include this computation into the call to
Pair::ev_tally()?
for that we already have a callback mechanism in place, so that no changes
to many pair styles would be required.​

If so, would you rather create another array for storing the values of the
electric field, or "re-use" the force array (might be a bit unpredictable
in
the long run, but introduces a memory overhead if you don't calculate
them(?)).

​i would not worry about cost of extra storage that scales ​with the number
of particles.
it can be easily implemented on-demand, even without code changes by using
fix property/atom.

The implementation of a corresponding, let's call her "compute_efield",
member function would also be helpful for implementing constant potential
methods (that is something else people in our group are looking into).

​the best way to get feedback for this, would be to do a minimal prototype
implementation and submit it as a pull request on github and then request a
review from steve.​

​axel.​

Hi Axel,

response inline below…

In many cases one could. It looks like many pair styles already have separate variables for coulomb and vdW forces anyhow. If no point dipoles are involved, one could include the individual force parts in a call to Pair::ev_tally() and tally to corresponding coulomb and vdW force arrays, and/or electric field.
What bothers me here, is that the interface to ev_tally would need to be changed to something it was not designed to do and where no explicit control mechanisms are in place. I.e., how do I decide from e- or vflag whether the forces should be tallied, or the field be calculated?
One would then probably want to introduce another flag for force/field tallying, which would also mean changing the interface to Pair::compute() (i.e., changes to all pair styles).
How about a separate, optional tallying mechanism (call it ff_tally - forces and field tallying) for the Pair class? Then one would probably need to change the interface to Pair::compute(), but would keep the ev_tally mechanisms for what they were designed for.
All in all, changing the interface seems like a thing to avoid, but I don’t see a way around it.
Also, one needs separate calculations for charge-charge, charge-dipole, and dipole-dipole interactions, in one central method, and I’d rather leave those specifics to the pair style.

Not all pair styles would need to provide a compute_efield method. Not all provide a compute_inner/middle/outer either. So there might be some incompatibility, but it’s also not so difficult to extract the computation of the electric field from the pair forces. And, in this case no interfaces need changing.

good point. thanks!

Then I will do that. Thanks,

Frank

Hi Axel,

​[...]

​hmm... would it be possible to include this computation into the call to
Pair::ev_tally()?
for that we already have a callback mechanism in place, so that no
changes to many pair styles would be required.​

In many cases one could. It looks like many pair styles already have
separate variables for coulomb and vdW forces anyhow. If no point dipoles
are involved, one could include the individual force parts in a call to
Pair::ev_tally() and tally to corresponding coulomb and vdW force arrays,
and/or electric field.

​my main point was to *avoid* any API changes. from within the the tally
you would have your callback from a compute called, which in turn can do
the particular pairwise calculation​, that you want to do.

What bothers me here, is that the interface to ev_tally would need to be
changed to something it was not designed to do and where no explicit
control mechanisms are in place. I.e., how do I decide from e- or vflag
whether the forces should be tallied, or the field be calculated?

​you have to write a compute, that provides a
Compute::pair_tally_callback() function, this compute will also set a flag,
that will trigger the call to Pair::ev_tally(), typically requesting the
total energy is sufficient. please see the examples in src/USER-TALLY

your fix aspc would then internally create the compute that would collect
the data you want.

One would then probably want to introduce another flag for force/field
tallying, which would also mean changing the interface to Pair::compute()
(i.e., changes to all pair styles).

​nope.​

How about a separate, optional tallying mechanism (call it ff_tally -
forces and field tallying) for the Pair class? Then one would probably need
to change the interface to Pair::compute(), but would keep the ev_tally
mechanisms for what they were designed for.

Pair::ev_tally() designed​ to tally pairwise information, so adding an
extension to it, is the most elegant way to allow computing and tallying
additional properties that require additional computations or additional
properties without adding code or overhead to the normal case of not
needing/wanting it. the only significant limitation currently is, that it
is only available in Pair::ev_tally() and thus not for properties in
certain other pair styles and for properties based on bonded interactions.

we put that mechanism into Pair::ev_tally() exactly to avoid adding another
(performance limiting) API

All in all, changing the interface seems like a thing to avoid, but I
don't see a way around it.
Also, one needs separate calculations for charge-charge, charge-dipole,
and dipole-dipole interactions, in one central method, and I'd rather leave
those specifics to the pair style.

Not all pair styles would need to provide a compute_efield method. Not all
provide a compute_inner/middle/outer either. So there might be some
incompatibility, but it's also not so difficult to extract the computation
of the electric field from the pair forces. And, in this case no interfaces
need changing.

​ultimately, it is not me that you have to convince, but steve.​ he is
usually hesitant to add APIs (and rightfully so), but not against it, if
there are good arguments for it.

​axel.​

That is a good point. Sorry, I missed the bottom/callback part of Pair::ev_tally(). If I put the electric field calculation into the callback, then I need to redo the pairwise calculations (almost) exactly again, because only fpair (and not fcoul/fvdw) gets passed to Pair::ev_tally().

In the cases that I am interested in, I need to get the electric field before the final force calculation. Doing a calculation of the pairwise forces, then redoing the electric field calculation through the callback in Pair::ev_tally() seems like putting an additional load on the already heavy polarizable force-field calculation.

Wouldn’t it then make sense to have a separate compute that calculates the electric field, but not included in the Pair::ev_tally() callback, that hence does not need to do the pairwise force calculation as well? I could call that from within my fix as well, couldn’t I?

OK, gotcha.

I would like to avoid API changes as well. That was one motivation for actually asking you guys.

Best,

Frank