About the Shifted Lennard-Jones pair force

[John Greenfield]

Hi, I am a new user of Hoomd. I was cheking the documentation for the
Shifted Lennard-Jones pair force and it seems that the offset (delta)
is fixed in Hommd:

delta = (di+dj)/2 + 1

(should not be sigma instead of 1?)

In LAMMPS it is possible to use positive and negative values for
delta.

Is there any way in which delta can be varied in the current compiled
version of Hoomd or is it necessary to modify any particular
subroutine?

Thanks

John Greenfield

[Carolyn Phillips]

Hi John

For the Shifted LJ, the idea is that the potential is "shifted" radially to the surface of the particle. If sigma were changed, than the shape of the potential would also change. So no, it should not be sigma. This is the same formula as the lj/expand formula in lammps

You have copied the delta formula incorrectly.

delta = (di +dj)/2 -1

Note, that if di and dj are less than 1, this is negative. So for Hoomd it is also possible to have both positive and negative values of dleta.

In Lammps, particles are implicitly given diameters by type and it is up to you to calculate the correct delta for every particle pair based on the formula above.
In Hoomd, every particle in the system can have a different diameter. As long as these diameters are specified during the initialization of the system, the SLJ pair potential automatically correctly calculatess every delta correctly between every particle pair.

Cheers,
-Carolyn

> --
> You received this message because you are subscribed to the Google Groups "hoomd-users" group.
> To post to this group, send email to hoomd...@googlegroups.com.
> To unsubscribe from this group, send email to hoomd-users...@googlegroups.com.
> For more options, visit this group at http://groups.google.com/group/hoomd-users?hl=en.
>
>
>

[John Greenfield]

Dear Carolyn,

Thank you for your prompt answer. Effectively, I copied incorrectly
delta, as you wrote it is defined by

delta = (di +dj)/2 -1

Notice that delta cannot be never negative. Its minimum value is zero,
because the diameters are scaled with respect to the smallest diameter
to which corresponds the value of 1 in lj units.

I agree with you that the idea of the Shifted LJ is to shift the LJ
potential to the surface of the particle. This is specially useful
when particles of different sizes are considered. For example, let us
suppose that we have a mixture of repulsive particles of type A of
diameter 1 and particles of type B of diameter 7. Let us define then a
vertical shifted slj potential between them defined as (you can type
this in gnuplot):

vsslj(x)= 4*e*( (s/(x-d))**12 - (s/(x-d))**6) + 1
e=1
s=1
d=4-1
p [4:5] [0:1] vsslj(x)

were the "contact" distance between one particle of species A and
another one of type B would be (1+7)/2=4. Let us consider also a
cutoff at the minimum in order to mimic repulsive particles. If you
change the previous expression by
s=0.1
d=4-s
p [4:5] [0:1] vsslj(x)

then you have a "harder" interaction potential. If were possible to
adjust delta as in LAMMPS, it would be possible to modulate the
softness or hardness of particles of different size in HOOMD.

Best

John Greenfield

[Carolyn Phillips]

Hi John,

This potential functions identically to what is in LAMMPS.

Yes. It can be negative. It can be negative because the diameters are not scaled with respect to the smallest diameter. That would be true if we combined our LJ potential by using combination sigmaA sigmaB rules and then made smaller of the two the units of the system. That is not what we do. Sigma is kept fixed at 1. Sigma is the units of distance of the system but is no longer "the diameter of the smallest particle". If you define the particles to have diameter = 0, and use shifted lennard jones, they will pass through each other like ghosts.

-Carolyn

[Carolyn Phillips]

Or... clarification, if the particles are given a diameter of 0, and you are using a shifted LJ, they will feel almost only attraction, and, as the temperature of the system is lowered, their centers will sit about 0.1units sigma apart.... not quite ghosts.

-Carolyn

[Joshua A. Anderson]

On Nov 18, 2010, at 6:07 PM, Carolyn Phillips wrote:

> Hi John,
>
> This potential functions identically to what is in LAMMPS.
>
> Yes. It can be negative. It can be negative because the diameters
> are not scaled with respect to the smallest diameter. That would be
> true if we combined our LJ potential by using combination sigmaA
> sigmaB rules and then made smaller of the two the units of the
> system. That is not what we do. Sigma is kept fixed at 1. Sigma
> is the units of distance of the system but is no longer "the
> diameter of the smallest particle".

Respectfully I must disagree, Carolyn. Your implementation of slj does
completely not meet this explanation. sigma is a parameter that you do
allow to be set, and it has units of distance. As John has pointed
out, setting any sigma other than 1.0 results in a potential that is
clearly _not_ shifted to the surface of the particle.

There are two possible routes to obtain the "hard" s=0.1 potential
John gave as an example.
1) Set sigma=1.0 (i.e., remove it as a parameter for slj) and rescale
the units of distance by a factor of 10. Thus, the particles would
have diameters of 10 and 70.
2) Modify the slj potential (and the rest of the code) to compute
delta=(di+dj)/2 - s as suggested.

Frankly, neither of these is a very optimal solution. The first would
result in box sizes of tens of thousands of distance units if you
simulated a lot of these particles. That may result in serious
cancellation errors when computing forces out near the box edges.

As for the 2nd - well, slj has already wreaked havoc on the r_cut
handling in the neighbor list and pair force routines. Adding sigma
into the computation of delta is just not possible from a software
engineering perspective. The neighbor list would need a way to query
sigma per each particle type pair from pair.slj. And that is only the
beginning! What happens as more shifted pair potentials are added to
the mix? Each of them with their own preferred definition of delta?
This just isn't a scalable solution. Basically, to answer your
original question about what needs to be modified, John, is that the
delta computation shows up in quite a few places and would need to be
changed in each one in a non-trivial way.

The only scaleable solution for delta is to keep it as it is and
implement every "shifted X" potential in the future to have a natural
diameter of 1.

I only see two realistic avenues for escape from this circle.
1) What issues are you trying to solve, John? Is it that you need a
harder shifted potential as in your example? What if we added a
shifted potential that has a much higher power of n in place of the
12. That results in a very step hard potential that still has a
minimum near a natural diameter of 1.

2) As an alternate, perhaps we could expand the definition of delta
slightly. As I said earlier, it cannot possibly be linked in with the
slj parameter specification per type pair. BUT, we could define delta
to be (di+dj)/2 - S, where S is a globally chosen constant. This would
only take a moderate amount of work to implement. The biggest problem
with it is that pair.slj (or any other shifted potential) will only
work when _all_ sigmas (or natural diameter) settings for all type
pairs are equal to S.
--------
Joshua A. Anderson, Ph.D.
Chemical Engineering Department, University of Michigan

[Carolyn Phillips]

Right. Sorry John, I missed that you really were trying to do something outside of the cases our implementation of the slj potential was designed for.

[John Greenfield]

Dear Joshua and Carolyn,

I appreciate a lot the time you spent trying to understand my
question.

My interest in this issue is because I would like to study jamming and
glass transition of intruders in hard spheres. I have already
developed a serial code to perform MD of hard spheres, but takes a lot
of MD steps to measure the diffusion of the intruders (specially at
high volume fractions). I chose hard spheres because several theories
have been developed for this model. The hard sphere limit should be
recovered in the limit of S->0, when sigma=S, in the slj potential. A
good alternative should be that proposed by Joshua in 2). In that
case, the hardness would be the same for arbitrary di and dj and
modulated by S (that is, the minimum of the lj potential would be at
the same position measured from the surce for all species). One
disadvantage of this treatment is that the time steps should be very
smalls because of the very large forces near the contact, but I think
that could be compensated with the speed of GPUs. Nevertheless,
perhaps adding this feature in future HOOMD versions is far from
trivial due to the neighbor lists use.

Best,

John Greenfield

[Carolyn Phillips]

Hi John,

Definitely one problem with using "classical" MD to study hard particles is that the "harder" the particles, the smaller the time step must be. Josh summed up some approaches nicely. Personally, I find that one benefit of shifting only is that I don't have to worry about my time step size. Whatever time step worked for diameter 1, particles, works for larger particles too. Now the effective time elapse of my simulation is dropping linearly because in the same time step size, my particles move a smaller fraction of their own diameter. On the other hand, if you want to keep your particles effectively fixed in diameter, and then increase the "hardness", you will have to continually lower your time step. On the other hand, I have not pushed my particles to extreme diameters where the other types of numeric problems in a simulation code may come into play.

If you are lucky, your particles don't need to be very very hard to act hard enough. If your particles *must* be very very hard, then I think MD (or at least classical MD) is not a very good tool for the job. Event-driven MD or Monte Carlo techniques become more appropriate.

-Carolyn