RDF of results by HOOMD-blue simulation.

[Elsaid Younes]

Dear all,

I am calculating the radial distribution function of a particle of a radius of 75, but the results always like this "as the attached picture"

Best,

Elsaid

[Bradley Dice]

Hi Elsaid,

Can you provide some frames from your trajectory (ideally as a GSD file) and the Python code you're using to analyze it? I can't tell what is wrong just from the picture.

Best,

Bradley

[Elsaid Younes]

Hi,

GSD file:

import numpy as np
import gsd.hoomd
import hoomd
import hoomd.md


hoomd.context.initialize("--mode=cpu")


s = gsd.hoomd.Snapshot()
s.particles.N = 2961
s.particles.typeid = [0]*2940 + [1]*21
s.particles.types = ['A','B']
s.particles.diameter = [4.0]*2940 + [150.0]*21
s.configuration.box = [906.0,906.0,906.0,0,0,0]
s.particles.charge = [+1]*2940 + [-140]*21
s.particles.position = result

with gsd.hoomd.open(name='file.gsd', mode='wb') as f:
    f.append(s)
   

Simulation code:

import hoomd
import freud
import gsd.hoomd
import numpy as np
import matplotlib.pyplot as plt
import hoomd
from hoomd import md

hoomd.context.initialize('--mode=cpu')

system = hoomd.init.read_gsd('file.gsd')

nl = hoomd.md.nlist.tree()

# NVT integration
groupB = hoomd.group.type(type ='B', name= 'mic', update= False);

hoomd.md.integrate.mode_standard (dt = 0.005)
nvt = hoomd.md.integrate.nvt (group = groupB, kT = 1.2, tau = 1.0)
nvt.randomize_velocities (seed = 1)


pppm = hoomd.md.charge.pppm(groupB, nl)
pppm.set_params(Nx=64, Ny=64, Nz=64, order=6, rcut=400.0, alpha=  0.007066536138444733)

# run the simulation
hoomd.run (500)

rdf = freud.density.RDF(bins=200, r_max=100)
w6 = freud.order.Steinhardt(l=6, wl=True)
def calc_rdf(timestep):
    hoomd.util.quiet_status()
    snap = system.take_snapshot()
    hoomd.util.unquiet_status()
    rdf.compute(system=snap, reset=False)
def calc_W6(timestep):
    hoomd.util.quiet_status()
    snap = system.take_snapshot()
    hoomd.util.unquiet_status()
    w6.compute(system=snap, neighbors={'num_neighbors': 21})
    return np.mean(w6.particle_order)
# Equilibrate the system a bit before accumulating the RDF.
hoomd.run(500)
hoomd.analyze.callback(calc_rdf, period=100)
logger = hoomd.analyze.log(filename='output.log',
quantities=['w6'],
period=100,
header_prefix='#',
overwrite=True)
logger.register_callback('w6', calc_W6)
hoomd.run(500)
# Store the computed RDF in a file
np.savetxt('rdf.csv', np.vstack((rdf.bin_centers, rdf.rdf)).T,
delimiter=',', header='r, g(r)')
rdf_data = np.genfromtxt('rdf.csv', delimiter=',')
plt.plot(rdf_data[:, 0], rdf_data[:, 1])
plt.title('Radial Distribution Function')
plt.xlabel('$r$')
plt.ylabel('$g(r)$')
plt.show()



w6_data = np.genfromtxt('output.log')
plt.plot(w6_data[:, 0], w6_data[:, 1])
plt.title('$w_6$ Order Parameter')
plt.xlabel('$t$')
plt.ylabel('$w_6(t)$')
plt.show()

Best,

Elsaid

--
You received this message because you are subscribed to a topic in the Google Groups "freud-users" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/freud-users/ZJuJSBtTYsw/unsubscribe.
To unsubscribe from this group and all its topics, send an email to freud-users...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/freud-users/36efee7d-87a5-40ea-810d-b1394106a4cen%40googlegroups.com.

[vramasub]

Hi Elsaid,

I can't run your script exactly because I don't know what the variable `result` is supposed to be, but if you are confused because your RDF has nonzero values for R < 75 the reason is because you are not excluding the ions (the 'A' particle type) from the snapshot that you're passing to the RDF. If you only want the RDF for the 'B' particles, you need to filter the snapshot to only include those particles before passing it to freud in your `calc_rdf` function. See HOOMD's documentation for how to get subsets of positions etc from a snapshot.

Be aware that active HOOMD development is now on the 3.x beta line, so you may want to switch over in the near future. On an unrelated note, be aware that the length scale for HOOMD's PPPM the alpha parameter (the inverse Debye length) and the particle diameter isn't really relevant there. It looks like you already set this appropriately since 1/0.007066536138444733 ~ 141 ~150, but worth mentioning just in case.

Hope that's helpful.

Best,

Vyas

[Elsaid Younes]

Hi Vyas,

Thanks for your reply.

In the GSD:

result wich  is in positions:

a = np.array([451.0,451.0,451.0])
z = np.array([376.0,376.0,376.0])

data1 = np.linspace(-a, a, 2940)
np.random.shuffle(data1[:,0])
np.random.shuffle(data1[:,1])

data2 = np.linspace(-z, z, 21)
np.random.shuffle(data2[:,0])
np.random.shuffle(data2[:,1])

result = [*data1,*data2]

For the Debye parameter I converted it into A^-1 like the length units in my system. I attached what I have done to calculate the Debye parameter.0

Best,

Elsaid

To view this discussion on the web visit https://groups.google.com/d/msgid/freud-users/c4cdc6a7-a59e-46c6-b705-807fb5468731n%40googlegroups.com.

[vramasub]

Hi Elsaid,

I'm still not sure why you think the RDF is wrong. If the problem is the length scale, please see my previous answer: you need to limit your RDF to only the particles with radius 75, not the smaller particles with radius 2. You can do that by selecting the 'B' particles from the HOOMD snapshot.

Best,

Vyas

[Elsaid Younes]

Hello,

The mistake was ';' after 'hoomd.group...' line, when I deleted it, the RDF was calculated.

But many parameters there, could change the RDF, but I do not know  the reason.

These parameters are:

the pppm parameters; 'Nx,Ny,Nz,order ', bins and the period.

Can you tell me the meaning of them and the ideal values for my system?

Also can I know the meaning of 'w6'? How can I skip the 'w6' calculation in order to speed up the process?

Best regards,

Elsaid

To view this discussion on the web visit https://groups.google.com/d/msgid/freud-users/abc64b5f-6e79-4383-8324-fc48cf26422en%40googlegroups.com.

[vramasub]

Hi,

I'm not sure I follow, semicolons at the end of lines in Python should have no effect, but I am glad that you have a working script now.

The bins parameter governs the distances at which the RDF is sampled. If you use too many bins the calculation will eventually slow down, but you can usually afford to make it reasonably large. If you choose too small a value the spacing between bins may be bigger than the features in the RDF (the peaks and valleys), so I would suggest choosing a number of bins such that the bin size is small enough to capture variations on the order of your particle diameters since that's the smallest length scale at which your RDF could have features. The period is how frequently data is added to the RDF during your simulation. You should see how the RDF plot looks, if the period is too large the RDF will look very jagged because it's not well-sampled.

The PPPM parameters are unfortunately harder to choose. Typical approaches involve using a model system (e.g. a pair of ions), setting a desired accuracy in the electrostatic calculations, then selecting a grid size that that will achieve that accuracy. Choosing larger numbers of grid points will slow the Ewald summation part, so keep that in mind. I would suggest looking at some literature for suggestions on approaches for selecting the parameters.

w6 is a Steinhardt order parameter as originally defined in this paper. It is a method for characterizing the crystalline ordering in a system. You can remove w6 by deleting the following from your script:

# This line creates the w6 order parameter calculator

w6 = freud.order.Steinhardt(l=6, wl=True)

# The code below defines a callback function that uses the calculator.

def calc_W6(timestep):
    hoomd.util.quiet_status()
    snap = system.take_snapshot()
    hoomd.util.unquiet_status()
    w6.compute(system=snap, neighbors={'num_neighbors': 21})
    return np.mean(w6.particle_order)

# The lines below tell HOOMD to use the function above to calculate w6 and store in the output.log file.

logger = hoomd.analyze.log(filename='output.log',
quantities=['w6'],
period=100,
header_prefix='#',
overwrite=True)
logger.register_callback('w6', calc_W6)

Best,

Vyas

[Elsaid Younes]

Dear all,

How can I ask freud to calculate r.d.f. only for certain type (if I have types )of particles in GSD file)?

Best,

Elsaid

--

You received this message because you are subscribed to a topic in the Google Groups "freud-users" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/freud-users/ZJuJSBtTYsw/unsubscribe.
To unsubscribe from this group and all its topics, send an email to freud-users...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/freud-users/513617ab-ca87-4815-b996-bde2cc23c586n%40googlegroups.com.

[Tommy Waltmann]

Hi Elsaid,

All you have to do is pick the indices of the particles of the type you want to use. Using your gsd file from earlier in the thread as an example, if you wanted to compute the RDF only for your 'A' particles,

snap = system.take_snapshot()

points = snap.particles.position[:2940, :]

rdf.compute(system=(snap.configuration.box, points), reset=False)

Best,

Tommy

[Elsaid Younes]

Hi Tommy,

it gives me this error:

AttributeError: 'hoomd._hoomd.SnapshotSystemData_float' object has no attribute 'configuration'

Bes,

Elsaid

You received this message because you are subscribed to the Google Groups "freud-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to freud-users...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/freud-users/13ebcfb9-4e45-4f97-adad-ae80e72f5ea8n%40googlegroups.com.

[Thomas Waltmann]

Hi Elsaid,

Getting the box from a hoomd snapshot is done via snap.box, while getting the box from a gsd snapshot is done through snap.configuration.box. I mistakenly wrote the above example for the wrong type of snapshot.

Best,

Tommy

[Joshua Anderson]

snapshot.box is the HOOMD v2 behavior. In v3, it is snapshot.configuration.box - matching GSD's API.
------
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/freud-users/CAEZTVZA9%2BrMMBR9VTrUPS%2Bwp8--Xng%2Bmn4JhRkksaUoZ_9Mg0A%40mail.gmail.com.

[Elsaid Younes]

Hi,

   snap = system.take_snapshot()

    points = snap.particles.position[:21,:]
    rdf.compute(system=(snap.box, points), reset=False)

It should measure r.d.f. for the type 'B' paticles. The diameter of the type 'B' = 150??

Best,

Elsaid 

[Tommy Waltmann]

Elsaid,

You have chosen to use the first 21 particles in your simulation to compute the RDF. Based on your previous scripts, I think you want the last 21 particles in your simulation.

-Tommy

[Elsaid Younes]

Hi,

    points = snap.particles.position[len -21]

gives me the error:

-----
-- Neighborlist stats:
0 normal updates / 1 forced updates / 0 dangerous updates
n_neigh_min: 0 / n_neigh_max: 1553 / n_neigh_avg: 761.291
shortest rebuild period: 100
** run complete **
** starting run **
Traceback (most recent call last):

  File "/home/elsaid/untitled4.py", line 56, in <module>
    hoomd.run(100)

  File "/home/elsaid/anaconda3/lib/python3.8/site-packages/hoomd/__init__.py", line 201, in run
    context.current.system.run(int(tsteps), callback_period, callback, limit_hours, int(limit_multiple));

  File "/home/elsaid/untitled4.py", line 46, in calc_rdf
    points = snap.particles.position[len -21]

TypeError: unsupported operand type(s) for -: 'builtin_function_or_method' and 'int'

Best,

Elsaid

To view this discussion on the web visit https://groups.google.com/d/msgid/freud-users/4fed608a-ebc5-4bc9-9f7e-5ec1cca1c96fn%40googlegroups.com.

[Tommy Waltmann]

Elsaid,

You have a python syntax error in your code, can you spot where it is and fix it? Please only ask questions on this mailing list related to the freud library.

Best,

Tommy

[Elsaid Younes]

Hi,

After I modified the index of the points, it runs as before. But I still have the same problem;

 (the particles seem to have a diameter less than 150)

    points = snap.particles.position[-21:,]

Best,

Elsaid

To view this discussion on the web visit https://groups.google.com/d/msgid/freud-users/4353808c-fcdf-4153-ad6e-e52d02a9abf9n%40googlegroups.com.

[Tommy Waltmann]

Elsaid,

Based on the scripts you have shown in this thread previously, you have only run your simulation for ~1000 time steps, that means your computed RDF is highly dependent on the initial configuration of your system. I don't know what your initial configuration is, so I can't say whether your particles are supposed to be within 100 distance units of each other or not.

Also, keep in mind that in your simulation, the only interactions present are electrostatic. Electrostatic interactions in hoomd do not use the diameters from the system state to compute forces.

-Tommy

[Elsaid Younes]

Hi,

That means that I have to use an alternative package other than Freud? or do you need a full-detailed script to decide?

Best,

Elsaid

To view this discussion on the web visit https://groups.google.com/d/msgid/freud-users/531ef00c-179c-4dbc-8cd8-f05347d35c01n%40googlegroups.com.

[Tommy Waltmann]

Elsaid,

Freud is a package that handles analyzing simulation output, not running simulations, so you can still use freud to analyze your simulation output. I think you want to know whether you should be using hoomd-blue or not, and I can't answer that for you. You have to think about whether the interactions supported by hoomd-blue are sufficient for you to model the system you are trying to study. Please refer to the hoomd documentation https://hoomd-blue.readthedocs.io/en/latest/