Re: [hoomd-users] Problem with repulsive LJ implementation

[Joshua Anderson]

Amir,

Use r_cut=2**(1/6).

------
Joshua A. Anderson, Ph.D.
Research Area Specialist, Chemical Engineering, University of Michigan

On Oct 5, 2021, at 12:06 PM, Amir Abbasi <che...@gmail.com> wrote:



Dear hoomd-blue users,

I encountered a problem implementing the repulsive LJ potential (which is also called WCA potential or shifted lennard-jones potential) using hoomd-blue v2.9.7 on a GPU machine and as you may already know this potential has a cut_off distance between particles equal to 2^(1/6) or in other words cut_off radius equal to 2^(1/6)/2.

Now, I am trying to do the tutorial for lj particle simulation using langevin thermostat implementing the shifted LJ potential as follows:

import hoomd
import hoomd.md
hoomd.context.initialize("");
import numpy
from matplotlib import pyplot

hoomd.init.create_lattice(unitcell=hoomd.lattice.sc(a=2.50), n=5);
nl = hoomd.md.nlist.cell();
slj = hoomd.md.pair.slj(nlist=nl,r_cut = 2**(1/6)/2)
slj.pair_coeff.set('A','A', epsilon=1.0, sigma=1.0)
hoomd.md.integrate.mode_standard(dt=0.000005);
all = hoomd.group.all();
hoomd.md.integrate.langevin(group=all, kT=0.20, seed=42);

notice(2): Group "all" created containing 125 particles

notice(2): Notice: slj set d_max=1.0

notice(2): integrate.langevin/bd is using specified gamma values

hoomd.analyze.log(filename="log-output.log",
                  quantities=['potential_energy', 'temperature'],
                  period=100,
                  overwrite=True);
hoomd.dump.gsd("trajectory_tut_hard.gsd", period=2e3, group=all, overwrite=True);
hoomd.run(2e6);
%matplotlib inline
data = numpy.genfromtxt(fname='log-output.log', skip_header=True);
pyplot.figure(figsize=(4,2.2), dpi=140);
pyplot.plot(data[:,0], data[:,1]);
pyplot.xlabel('time step');
pyplot.ylabel('potential_energy');

notice(2): -- Neighborlist exclusion statistics -- :

notice(2): Particles with 0 exclusions : 125

notice(2): Neighbors included by diameter : yes

notice(2): Neighbors excluded when in the same body: no

** starting run **

Time 00:00:14 | Step 355019 / 2000000 | TPS 35501.8 | ETA 00:00:46

Time 00:00:24 | Step 703449 / 2000000 | TPS 34843 | ETA 00:00:37

Time 00:00:34 | Step 1005058 / 2000000 | TPS 30160.8 | ETA 00:00:32

Time 00:00:44 | Step 1352421 / 2000000 | TPS 34736.3 | ETA 00:00:18

Time 00:00:54 | Step 1698837 / 2000000 | TPS 34641.6 | ETA 00:00:08

Time 00:01:02 | Step 2000000 / 2000000 | TPS 34597.5 | ETA 00:00:00

Average TPS: 34068.5 --------- -- Neighborlist stats: 0 normal updates / 6667 forced updates / 0 dangerous updates n_neigh_min: 0 / n_neigh_max: 1 / n_neigh_avg: 0.096 shortest rebuild period: 100 -- Cell list stats: Dimension: 13, 13, 13 n_min : 0 / n_max: 2 / n_avg: 0.0568958 ** run complete **

potential energy plot shows two sharp positive peaks. When I am visualizing the trajectory (using ovito) I see that particles are overlapping. Could you please let me know what I am missing in my code?

Cheers,

Amir

--
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/071291fb-450f-4e4c-899a-552e04318109n%40googlegroups.com.

[Amir Abbasi]

Thanks Joshua, I have tried that before. Actually it was my first attempt. It turned out that the problem is with visualization using ovito.

Does ovito interpret the distances in gsd file as hoomd-blue? I will try it with fresnel as well to see if this is the source for the problem.

Cheers,

Amir

[Joshua Anderson]

Amir,

Please refer to the OVITO or fresnel documentation for information on setting the display radius of particles. Neither OVITO nor fresnel will read your job script to automatically determine the proper radius based on your pair potential parameters.

If you desire your HOOMD-blue simulation to have particles interact with a hard core diameter of 1/2 instead of 1, then set sigma and r_cut appropriately:
lj.params[('A', 'A')] = dict(epsilon=1, sigma=0.5)
lj.r_cut[('A', 'A')] = 2**(1/6) * lj.params[('A', 'A')]['sigma']

------
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/ee5c4938-a3a2-404b-b363-6dce3ca1c515n%40googlegroups.com.