python_rigid_rod.py

import itertools
import math

import gsd.hoomd
import hoomd
import matplotlib
import numpy as np
#define position of the 5 points in the rigid rod. We can change it to whatever number of points we want+ whatever spacing we want.
## doing this for 50 points so that it is relevant to the aspect ratio of 8nm/400nm
#placing all the constituent particles on x-axis
positions_x = np.arange(-30, 30, 12.0)
print(len(positions_x))

# Create the 3x3 array with the first column and zeros in the other columns
rod_positions = np.column_stack([positions_x, np.zeros(5), np.zeros(5)])
#print(rod_positions)
#rod_positions = [[-2.4,0,0],[-1.2,0,0],[0,0,0],[1.2,0,0],[2.4,0,0]]
# each instance of rigid body is placed in system global coordinates, located by position and orienttaion.
#where in the box should this rod be located
central_position = [10, 5, 0]
central_rotation = 0.9
#just positioning the rod and the particles in the box
cos_theta = math.cos(central_rotation)
sin_theta = math.sin(central_rotation)

global_positions = []
for i in range(len(rod_positions)):
    x, y = rod_positions[i][:2]

    global_positions.append(
        [[central_position[0] + (x * cos_theta - y * sin_theta)],
         [central_position[1] + (y * cos_theta + x * sin_theta)]])

## define properties of rigid rod, what are the comprising particle types
frame = gsd.hoomd.Frame()
frame.particles.types = ['rod', 'A']

#place rod particles in simulation box. place them far apart so A's don't overlap. We'll probably need to change\
#  this according to the concentration of rods in solution we want to use.
#TODO: Figure out box dimension and number of rods to use.

N_particles = 10
spacing = 50*2
K = math.ceil(N_particles**(1 / 3))
L = K * spacing
print(L)
print(L)
Lc=400*L*10**(-9)/60.0
print(N_particles/(6*10**23*Lc**3))
#L=(N_particles/(24*10**12))**(1/3)
#print(L)
#L=1410
#K = 2*math.ceil(N_particles**(1 / 3))
#L = K * spacing*5
x = np.linspace(-L / 2 , L / 2, K, endpoint=False)
position = list(itertools.product(x, repeat=3))
#print(position)
position = np.array(position) + [spacing / 2, spacing / 2, spacing / 2]
print(position)
print(len(position))
frame.particles.N = N_particles
frame.particles.position = position[0:N_particles, :]
frame.particles.typeid = [0] * N_particles
frame.configuration.box = [L, L, L, 0, 0, 0]

# define mass of each rod. Maybe need to discuss this?
frame.particles.mass = [5] * N_particles
#moment of inertia for each road. Important for rotation.
mass = 1.0
I = np.zeros(shape=(3, 3))
for r in rod_positions:
    I += mass * (np.dot(r, r) * np.identity(3) - np.outer(r, r))
frame.particles.moment_inertia = [I[0, 0], I[1, 1], I[2, 2]] * N_particles
frame.particles.orientation = [(1, 0, 0, 0)] * N_particles
rodnames=['A']*5
orientation_constituent=[(1.0,0.0,0.0,0.0)]*5
print(rodnames)
#save the centers of rods
with gsd.hoomd.open(name='dimer_centers_rods.gsd', mode='w') as f:
    f.append(frame)
#start making the rod rigid
rigid = hoomd.md.constrain.Rigid()
#start defining the properties of this rigid rod
rigid.body['rod'] = {
    "constituent_types": rodnames,
    "positions": rod_positions,
    "orientations": orientation_constituent,
}
#define a simulation object
#read the dimer centers of rods that we saved earlier.
simulation = hoomd.Simulation(device=hoomd.device.CPU(), seed=4)
simulation.create_state_from_gsd(filename='dimer_centers_rods.gsd')
#add body to the center of rods
rigid.create_bodies(simulation.state)
hoomd.write.GSD.write(state=simulation.state, mode='wb', filename='lattice_rod.gsd')
#start the simulation from the lattice rod state we saved
cpu = hoomd.device.CPU()
simulation = hoomd.Simulation(device=cpu, seed=1)
simulation.create_state_from_gsd(filename='lattice_rod.gsd')
#rigid and constrain our rod

rigid = hoomd.md.constrain.Rigid()
#properties of the rigid rod

rigid.body['rod'] = {
    "constituent_types": rodnames,
    "positions": rod_positions,
    "orientations": orientation_constituent,
}
rigid.create_bodies(simulation.state)
#add an integrator here to the simulation box
integrator = hoomd.md.Integrator(dt=0.0005, integrate_rotational_dof=True)
simulation.operations.integrator = integrator
# make sure that the integrator assumes that the bodies are rigid
integrator.rigid = rigid
kT = 0.01
#setting the temperature. I think the filter decided which particles to apply rigid body definition to?
rigid_centers_and_free = hoomd.filter.Rigid(flags=('center',))
#what kind of integrator do we need to add. Because we are doing brownian dynamics that's what we'll use here
#nvt = hoomd.md.methods.ConstantVolume(
#    filter=rigid_centers_and_free,
#    thermostat=hoomd.md.methods.thermostats.Bussi(kT=kT))


brownian_d=hoomd.md.methods.Brownian(filter=rigid_centers_and_free,kT=kT)
integrator.methods.append(brownian_d)
#neighbor list for force calculation
cell = hoomd.md.nlist.Cell(buffer=0, exclusions=['body'])
"""""
yukawa = hoomd.md.pair.Yukawa(default_r_cut=2.0, nlist=cell)
yukawa.params[('rod', 'rod')] = dict(epsilon=1.0, kappa=1.0)
yukawa.params[('A', 'rod')] = dict(epsilon=0.0, kappa=1.0)
yukawa.params[('A', 'A')] = dict(epsilon=0.0, kappa=1.0)
""" 
lj = hoomd.md.pair.ForceShiftedLJ(nlist=cell,default_r_cut=0.8)
lj.params[('rod','rod')] = dict(epsilon=0.000, sigma=0.0)
#lj.r_cut[('A', 'A')] = 0.6*(2**(1/6))
lj.params[('A', 'A')] = dict(epsilon=0.001, sigma=0.6)
lj.params[('rod', 'A')] = dict(epsilon=0.0, sigma=0.0)
#lj.r_cut[('rod', ['rod', 'A'])] = 0

integrator.forces.append(lj)
simulation.run(1000)