R.D.F from electrostatic simulation

17 views
Skip to first unread message

Elsaid Younes

unread,
Sep 15, 2021, 8:02:30 AMSep 15
to hoomd...@googlegroups.com
Dear all,

Here is the r.d.f. diagram I got from electrostatic simulations by HOOMD-blue:

image.png
The diameter of the particle is 150 so I expected the first peak would be after 150. 
I show my GSD file then script:

GSD file

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



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

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

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

list1 = data1.tolist()

# print(len(list1))

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

list2 = data2.tolist()
# print(len(list2))
result1 = [*data1,*data2]
result = [*list1, *list2]
# print(len(result))
# print(any(result))


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 = [11.838686524544189]*2940 + [-1657.4161134361864]*21
s.particles.position = result1

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

script:
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
all = hoomd.group.all();
hoomd.md.integrate.mode_standard (dt = 0.005)
nvt = hoomd.md.integrate.nvt (group = all, kT = 4.11433*(10**-24), tau = 1.0)
nvt.randomize_velocities (seed = 1)

pppm = hoomd.md.charge.pppm(all, nl)
pppm.set_params(Nx=150, Ny=150, Nz=250, order=6, rcut=450.0, alpha=2.91419808943119)

# run the simulation
hoomd.run (100)

rdf = freud.density.RDF(bins=150.0, r_max=450.0)
def calc_rdf(timestep):
     
    hoomd.util.quiet_status()
    snap = system.take_snapshot()
   
    points = snap.particles.position[-21:,]
    rdf.compute(system=(snap.box, points), reset=False)
 
# Equilibrate the system a bit before accumulating the RDF.
hoomd.run(100)
hoomd.analyze.callback(calc_rdf, period=10)
hoomd.run(100)
# Store the computed RDF in a file
np.savetxt('rdf_411_100s_bins=150_period=10_N(x,y,z)=256_order=2.csv', np.vstack((rdf.bin_centers, rdf.rdf)).T,
delimiter=',', header='r, g(r)')
rdf_data = np.genfromtxt('rdf_273_100s_bins=150_period=10_N(x,y,z)=256_order=2.csv', delimiter=',')
plt.plot(rdf_data[:, 0], rdf_data[:, 1])
plt.title('Radial Distribution Function')
plt.xlabel('$r$')
plt.ylabel('$g(r)$')
plt.show()0


Best,
Elsaid

Michael Howard

unread,
Sep 15, 2021, 2:25:44 PMSep 15
to hoomd-users
There are a few problems with this script:

1. Your initialization script does not ensure there are no overlaps of the two types. You pick both sets of positions at random from overlapping lattices, so it is likely you have some.

2. Your temperature is ludicrously small, so your particles basically do not move. Given that you observe overlap in the RDF based on the nominal diameter, this is probably due to #1. I have already given you advice about how to choose a system of units where parameters are of order unity (https://groups.google.com/g/hoomd-users/c/bo97Kj4qRzQ), so I am not going to rehash that here and you can figure out how to fix this yourself.

3. Even if you fix #1 and #2, your simulation script does not have a pair potential to exclude volume (only the electrostatic potential). You could use a WCA potential to mimic hard-sphere interactions, which could be done using the LJ potential (https://hoomd-blue.readthedocs.io/en/latest/module-md-pair.html#hoomd.md.pair.LJ). However, if you add this potential properly and you have overlap in your initial configuration, the simulation will probably blow up.

I think multiple people have recommended this to you before, but you should try doing a simulation of hard spheres before you add the electrostatic potential to make sure you understand what you are doing.

Regards,
Mike
Reply all
Reply to author
Forward
0 new messages