Hard sphere particle simulation queries

[arj...@ncsu.edu]

Hi guys

I am a novice in HOOMD Blue. I am trying to simulate hard sphere suspension(particles have brownian motion) system with a simple code (See attached). I am also attaching two images which indicates the initial condition and final state of the system.Could someone help  through the following queries?

1) How do I implement a pure repulsive hard sphere potential in HOOMD. In 'Final state.png' image I am observing unphysical overlaps which is not what I require.

2) How do I control the radius of a particle in HOOMD to implement in the following code.

3) How do I control the number of particles in the simulation box

4) Can I put particles at random positions for my initial state unlike seen in 'Initial Condition .png'

Please bear with me if some of the questions might sound pretty trivial to you.

Any help would be of invaluable help to me.

Regards
Alan

[Matthew Spellings]

Hello Alan,

  I think the main difficulty you're having is due to not defining any interactions between particles in your system. 

1. For truly hard particles, you'll need to use an HPMC integrator rather than md.integrate.brownian. If you want to still use md.integrate.brownian with differentiable forces, you can use something like md.pair.lj.
   
2. You can construct the unit cell you're passing into init.unitcell and explicitly give each particle a different diameter. Alternatively, you can use the snapshot API (link) to set per-particle diameters.
3. The number of particles is controlled by the number of unit cells you initialize with, in this case the 'n=10' argument of create_lattice
4. You can use any random initialization you'd like to code into your python script; the snapshot API is helpful for this as well.

Matthew Spellings

[Michael Howard]

To add onto Matthew's great answers, to get nearly-hard spheres with md.pair.lj, set r_cut=2**(1./6.)*sigma, which is the minimum of the potential with particle diameter sigma, so that only the repulsive part of the potential is active. This is what people usually call the Weeks-Chandler-Andersen (WCA) potential, since it comes out of their perturbation theory. By setting sigma, you are setting the diameters of the particles through their interactions. You may still get some overlaps, but you can decrease these by increasing epsilon. Otherwise, you will need a potential that diverges at contact.

Regards,

Mike

[arj...@ncsu.edu]

Thank you guys!!!

I have found deprecated module which is exactly what I am looking for initialization of random particle positions. But now there is error in my cell list that I cannot figure out.

The code is the following

import hoomd
import hoomd.md
import hoomd.deprecated
#import hoomd.hpmc
#import ex_render
hoomd.context.initialize('');
#hoomd.init.create_lattice(unitcell=hoomd.lattice.sq(a=2.0), n=10);
#initialise random particles
hoomd.deprecated.init.create_random(N=1000, phi_p=0.5,name='A', min_dist=0.7)
all = hoomd.group.all();
#dump data for the visulaisation
d = hoomd.dump.gsd("trajectory.gsd", period=10, group=all, overwrite=True);
# setup pair potential inbetween particles
nl = hoomd.md.nlist.cell();
lj = hoomd.md.pair.lj(r_cut=2.0, nlist=nl)
lj.pair_coeff.set('A','A', epsilon=0.25, sigma=2, alpha=0, r_cut=2.0);
# setup forces, integrators and run simulation
hoomd.md.integrate.mode_standard(dt=0.001);
hoomd.md.integrate.brownian(group=all, kT=1,seed=250);
hoomd.run(2e4);


The ERROR Is

** starting run **
**ERROR**: Particle with unique tag 856 is no longer in the simulation box.

**ERROR**: Cartesian coordinates:
**ERROR**: x: -302.26 y: -48.2951 z: -147.622
**ERROR**: Fractional coordinates:
**ERROR**: f.x: -29.2649 f.y: -4.25583 f.z: -14.037
**ERROR**: Local box lo: (-5.07746, -5.07746, -5.07746)
**ERROR**:           hi: (5.07746, 5.07746, 5.07746)
Traceback (most recent call last):
  File "test2.py", line 18, in <module>
    hoomd.run(2e4);
  File "/home/alan/miniconda3/lib/python3.6/site-packages/hoomd/__init__.py", line 194, in run
    context.current.system.run(int(tsteps), callback_period, callback, limit_hours, int(limit_multiple));
RuntimeError: Error computing cell list



Can you provide some insights into what I am doing wrong and how to sort it out?

Thanks
Alan

[Michael Howard]

There are several reasons:

1. Your sigma is 2.0, but you set min_dist = 0.7, so the particles are overlapping in the initial configuration and then experiencing large forces, causing the simulation to blow up as the error message is telling you. You could see this if you looked at a snapshot of the initial configuration before starting the run. You need to increase min_dist. (In general, if you are new to hoomd, I would recommend following Matthew's suggestion to use snapshots and python to setup your simulation rather than the deprecated commands, which are legacy code and will be removed at some point.)

2. You are using dt=0.001 with the BD integrator and the default gamma value. This might be too large of a timestep. You probably need to increase gamma (with set_gamma) so that you truly have overdamped dynamics and also probably decrease dt.

3. You set r_cut=2.0 instead of the r_cut=2.0**(1./6.)*sigma that I recommended. Your potential is then discontinuous, and so you will have energy & momentum drift, which can cause problems.

As an aside, your epsilon is also very small relative to kT, and you will get significant particle overlaps at equilibrium. You need to *increase* epsilon if you want the potential to be more hard sphere like. But, increasing epsilon can require a smaller dt.

Regards,

Mike

[arj...@ncsu.edu]

Thanks Mike.

I am trying to make my own gsd file for initial configuration. I have been searching the information stored by a gsd file. It does not seem to store particle sizes or radius. Can I store particle sizes along with particle positions in a gsd file?

Regards
Alan

[Michael Howard]

Yes, it stores the diameter, with the default value being 1:

http://gsd.readthedocs.io/en/latest/schema-hoomd.html#chunk-particles/diameter

However, on its own this information doesn't do anything. You still need to use an appropriate pair potential in hoomd to implement the interactions between the particles as you integrate the equations of motion. gsd does not store this information.

Regards,

Mike