Re: [hoomd-users] rigid patchy particles with bond swaps

[Joshua Anderson]

Moritz,

HOOMD always assumes periodic boundary conditions. Particle out of the box errors indicate that the numerical solution to the equations of motion are unstable. This can be due to large forces between particles and/or too large a step size dt.

In your initialization script, you place particles with:
`points.append((random.uniform(*xrange), random.uniform(*yrange), random.uniform(*zrange))) for i in range(N)`
I would guess that this places particles very close together so that the LJ forces in your simulation are very large in the first timestep of the simulation.
------
Joshua A. Anderson, Ph.D.
Research Area Specialist, Chemical Engineering, University of Michigan

> On May 3, 2022, at 9:09 AM, Moritz Hopfenmüller <moritz.hop...@gmail.com> wrote:
>
> Hi,
>
> I try to build a networkstructure out of monomers( with 2 patches) and crosslinkers ( 4 patches). So first step I build a uniform random distribution of pointparticles with 5 percent crosslinker content (see uniform.py).
> Next step I define monomers and crosslinkers as rigid bodies.
> I apply a repulsive leonnard jones potential and a three body potential (for bondswaps, see rigids.py)
> Without the three body potential the particles move out of the simulation box resulting in an error, I'm new to hoomd-blue, so I assume the periodic boundary condition has to be set explicitly like in lammps?
> With the three body potential and n=100: **ERROR**: Particle with unique tag 15830 has NaN for its position.
> With the three body potential and n=10, again:**ERROR**: Particle with unique tag 2196 is no longer in the simulation box.
> the repuslive potential should only work between the center particles of different monomers, crosslinkers and between the center particles of monomers and crosslinkers, the attractive patchy potential and the bondswap potential should only work between patches (of different rigid molecules).
> Thanks for any help!
>
> --
> You received this message because you are subscribed to the Google Groups "hoomd-users" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to hoomd-users...@googlegroups.com.
> To view this discussion on the web visit https://groups.google.com/d/msgid/hoomd-users/d4db2940-bedb-4930-8ab2-9b6c22de86c4n%40googlegroups.com.
> <rigids.py><uniform.py>

[Isha Malhotra]

Hello Moritz, 

Maybe you have solved this issue by now, I also faced the same problem. 

Try removing "lj.params[('patch', ('mono', 'patch'))] = dict(epsilon=0, sigma=0)". Just put r_cut=0 for non interacting pairs and remove epsilon=0 and sigma=0 from your code. It will work fine then.

I don't have any idea why we are not able to put epsilon=0 and sigma=0, but this solved it for me.

Best, 

Isha Malhotra

[Joshua Anderson]

This advice is not correct. You must define pair parameters for all pair type pairs in the matrix, even when they are zero. To demonstrate, thr script below fails with `hoomd.error.IncompleteSpecificationError: For <class 'hoomd.md.pair.pair.LJ'> in TypeParameter params for key ('A', 'B') value for key epsilon is required`. Uncomment the two commented lines and the script runs without error.

As a convenience, you may set the `default` type parameters to 0 so that you do not need to explicitly define all zeros. However, the behavior of a script that does this will be equivalent to one that explicitly sets the 0's. 

If you found behavior that doesn't match the documentation or my description here, please submit a minimal script that demonstrates the behavior - such as a modified version of the script I post here.

---- script

import hoomd

snapshot = hoomd.Snapshot()
snapshot.configuration.box = (10,10,10,0,0,0)
snapshot.particles.N = 3
snapshot.particles.position[:] = [(-1, 0, 0), (0, 0, 0), (1, 0, 0)]
snapshot.particles.types = ['A', 'B']
snapshot.particles.typeid[:] = [0, 0, 1]

sim = hoomd.Simulation(device=hoomd.device.CPU())
sim.create_state_from_snapshot(snapshot)

cell = hoomd.md.nlist.Cell(buffer=0.4)
lj = hoomd.md.pair.LJ(nlist=cell)
lj.params[('A', 'A')] = dict(epsilon=1, sigma=1)
# lj.params[('A', 'B')] = dict(epsilon=0, sigma=0)
# lj.params[('B', 'B')] = dict(epsilon=0, sigma=0)
lj.r_cut[('A', 'A')] = 2.5
lj.r_cut[('A', 'B')] = 0
lj.r_cut[('B', 'B')] = 0

integrator = hoomd.md.Integrator(dt=0.005, forces=[lj])
sim.operations.integrator = integrator
sim.run(0)

------

Joshua A. Anderson, Ph.D.

Research Area Specialist, Chemical Engineering, University of Michigan

To view this discussion on the web visit https://groups.google.com/d/msgid/hoomd-users/6e62536d-3de2-4b68-9e26-cb75159399abn%40googlegroups.com.

[Isha Malhotra]

If we put r_cut to be zero and give any random value to epsilon and sigma, would the pair still remain non-interacting? 

"lj.params[('B', 'B')] = dict(epsilon=0, sigma=0)

lj.r_cut[('B', 'B')] = 0"

"lj.params[('B', 'B')] = dict(epsilon=1, sigma=1)

lj.r_cut[('B', 'B')] = 0"

Are the above two scripts equivalent? 

Thanks. 

To view this discussion on the web visit https://groups.google.com/d/msgid/hoomd-users/72017F5A-F747-4650-860D-7F6AD405FE82%40umich.edu.

[Joshua Anderson]

Isha,

The answer depends on how you define equivalent. The text of the script will not be identical, nor will the hash. Different parameters will be set in memory when run. Are you asking if they will give identical output for pair forces? The answer to that is: it depends on the other forces present in the system and the execution device. On the GPU, HOOMD evaluates force sums in an undefined order, so forces may differ slightly between executions.

However, with respect to the logical evaluation of a single pair interaction, when r_cut is 0, HOOMD ignores the parameters set and evaluates 0 energy and force as defined in the documentation (see the r >= r_cut cases): https://hoomd-blue.readthedocs.io/en/v3.2.0/module-md-pair.html

Is there something unclear about the documentation that leads you to think the behavior is otherwise?

------

Joshua A. Anderson, Ph.D.

Research Area Specialist, Chemical Engineering, University of Michigan

To view this discussion on the web visit https://groups.google.com/d/msgid/hoomd-users/CA%2BdjDeVnkDSEjLwROoT42KUOZT%2BM4xVXGcfOj8bpO0nhbFSZ5w%40mail.gmail.com.

[Isha Malhotra]

Joshua, 

In my script, if I use 

"lj.params[('B', 'B')] = dict(epsilon=0, sigma=0)

lj.r_cut[('B', 'B')] = 0"

I get NAN for particle positions, but if I use:

"lj.params[('B', 'B')] = dict(epsilon=1, sigma=1)

lj.r_cut[('B', 'B')] = 0"

then I don't get any such error and my simulation runs fine. I am not sure why I get this error with epsilon=0 and sigma=0. 

This is the exact error that I get:

"**ERROR**: Particle with unique tag 10594 has NaN for its position."

Thanking you, 

Isha Malhotra

You received this message because you are subscribed to a topic in the Google Groups "hoomd-users" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/hoomd-users/vXt3lKE7zZk/unsubscribe.
To unsubscribe from this group and all its topics, send an email to hoomd-users...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/hoomd-users/F7255A34-DD96-471C-B10A-BC5289E206E3%40umich.edu.

[Joshua Anderson]

Thanks for sharing the error message. I'm going to need a complete and minimal script that reproduces this behavior to look into this further. Also tell me which version of HOOMD-blue you are using.

I am not able to reproduce it with my test. The script below purposely places two non-interacting particles at the same position and sets r_cut=0 between them. The simulation runs normally and produces the expected output (0 force on particle 2 in all cases):

---- script.py

import hoomd

snapshot = hoomd.Snapshot()

snapshot.configuration.box = (6,6,6,0,0,0)
snapshot.particles.N = 3
snapshot.particles.position[:] = [(-1, 0, 0), (0, 0, 0), (0, 0, 0)]

snapshot.particles.types = ['A', 'B']
snapshot.particles.typeid[:] = [0, 0, 1]

sim = hoomd.Simulation(device=hoomd.device.CPU(), seed=10)

sim.create_state_from_snapshot(snapshot)

cell = hoomd.md.nlist.Cell(buffer=0.4)
lj = hoomd.md.pair.LJ(nlist=cell)
lj.params[('A', 'A')] = dict(epsilon=1, sigma=1)

lj.params[('A', 'B')] = dict(epsilon=0, sigma=0)

lj.params[('B', 'B')] = dict(epsilon=0, sigma=0)
lj.r_cut[('A', 'A')] = 2.5
lj.r_cut[('A', 'B')] = 0
lj.r_cut[('B', 'B')] = 0

integrator = hoomd.md.Integrator(dt=0.005, forces=[lj])
sim.operations.integrator = integrator

langevin = hoomd.md.methods.Langevin(kT=0.1, filter=hoomd.filter.All())
integrator.methods.append(langevin)

sim.run(0)

print(lj.forces)

sim.run(1000)

print(lj.forces)

---- output

[[-24.   0.   0.]
 [ 24.   0.   0.]
 [  0.   0.   0.]]
[[-0.26216708  1.6042809  -1.56972645]
 [ 0.26216708 -1.6042809   1.56972645]
 [ 0.          0.          0.        ]]

------

Joshua A. Anderson, Ph.D.

Research Area Specialist, Chemical Engineering, University of Michigan

To view this discussion on the web visit https://groups.google.com/d/msgid/hoomd-users/CA%2BdjDeUTM0Xpfv7ORL78ozbq%2BY6p62E_VOqEHDD7cA5x9_YUnQ%40mail.gmail.com.