Forces and rigid bodies

[trvsst]

I have a system consisting of many rigid bodies with some polymers attached to it. Let us call one of this units (rigid + flexible polymers attached to it) a nanoparticle.

As explained in the documentation
https://hoomd-blue.readthedocs.io/en/latest/module-md-constrain.html#hoomd.md.constrain.Rigid

The forces, energies and virials of the rigid parts are transferred to the center tag of the rigid body. Those are my questions:

• It is therefore true that the sum of the center tag and all the particles polymers (non-rigid) attached to it would give me the total force over the nanoparticle exerted by the other nanoparticles in the system. In general, that force is not zero, because the nanoparticle is not in equilibrium.

However, If I decide not to integrate the rigid center, because I want to study my nanoparticles at fixed positions, then what I find is that the force at the rigid center is zero, and the sum of all forces on the rigid particles is the negative of the one of the non-rigid particles. Basically the center contains the sum of all forces, whether rigid or not. That means that HOOMD has zeroed the force on the center so that it does stay in place (zero acceleration) but it somehow has modified the force on the rigid particle.

• How can I get the original forces exerted on the nanoparticle?
• Which local particle forces has HOOMD modified and which ones has not?

[joaander]

1. As the constrain.Rigid documentation states, the force on a given rigid body center is equal to the sum of the net forces on all the constituent particles. The md.Integrator documentation states that the net force any any given particle is the sum of the forces from all objects in the forces list. You can review the docs for the specific Force class instances you are using, but in general they all compute forces on particle i due to all other particles in the system. Using your definitions, this would include both intra- and inter-nanoparticle forces.

2. You can write a custom code outside HOOMD-blue that computes the inter-nanoparticle forces.

3. constrain.Rigid zeroes the net force and torque on the constituent particles:

   hoomd-blue/hoomd/md/ForceComposite.cc

   Lines 764 to 767 in 37f6348

   	// zero net energy on constituent particles to avoid double counting
   	// also zero net force and torque for consistency
   	h_net_force.data[idxj] = make_scalar4(0.0, 0.0, 0.0, 0.0);
   	h_net_torque.data[idxj] = make_scalar4(0.0, 0.0, 0.0, 0.0);

md.Integrator always computes all forces on all particles, regardless of the filter you apply to your integration method(s). If you find the force on a given particle is 0, that means you have provided a system state that results in a 0 force.

[trvsst]

ok, so I prepared these two scripts,

The first run creates a simulation with 8 nanoparticles, consisting of a rigid rod of 20 beads anf 4 flexible ones attached at the end. The simulation runs smoothly and after execution you can see the movie stored in rigid_movie.gsd. The simulation saves local energies and forces for the lennard jones particles. I run with large lennard jones cut-off.

https://www.dropbox.com/scl/fi/gqz9rxyvc6yn7smn3msif/rigid_body_rot.py?rlkey=zmgv8raycrz4x3wmpy1cf54a8&dl=0

The results are analyzed here:
https://www.dropbox.com/scl/fi/gqz9rxyvc6yn7smn3msif/rigid_body_rot.py?rlkey=zmgv8raycrz4x3wmpy1cf54a8&dl=0

The output is as follows, for each of the 8 nanoparticles I output the lj force on the last configuration for: the center, the force on the rigid particles and the force on the flexible ones. Here is what I do not understand:

• The force on the center is always zero, but the forces on the rigid components are not. The rigid body is clearly moving, so I do not understand why this is the case if the force on the center is zero. From the documentation "Rigid transfers forces, energies, and torques from constituent particles to the central particle and adds them to those from the interaction on the central particle itself."
• What is the meaning of the forces stored for the rigid body components
• Why the sum of the rigid have the same magnitude and opposite sign to the flexible ones (in some cases they are exactly the same, in some others this is only approximate).

[trvsst]

on a second thought, I think that I understand what happens, those are the answers to my questions, please confirm it it is true:

• The statement in the documentation refers to the net forces, but not to the lj forces. In other words, the lj forces are not transferred to the rigid center.
• The forces stored in the constituent rigid body are the lj forces exerted by the other particles.
• The sum is almost zero, because of Newton third law. Despite a generous cut-off most of the interactions are self-nanoparticle interactions.

[joaander]

Check rg.forces to find the forces on the central particles that constrain.Rigid computes.

[trvsst]

I post this code as I think it may be helpful for others. The code maybe made a little nicer, and perhaps, if others request it I will work further on it.

This is a straight-forward analysis that reads forces and energies(although energies are not used) for a system of n_rigid nanoparticles each consisting of n_apr particles forming a rigid body with n_bpr flexible particles attached at one end. Thus, each nanoparticle includes 1+n_apr+n_bpr particles, the first being the rigid center.

The actual simulation that generated the file rigid_local.gsd can be obtained for the code provided in the next comment

The analysis checks

• The forces on the center obtained from constrain.Rigid are the same as the sum of the forces, lj + bond, on all constituent particles.
• The total force on the nanoparticle is the sum of lj+bond for the rigid and flexible particles is equal to lj for the rigid and flexible because the bond are all internal forces and therefore action and reaction cancels them out. Only non-bonded interactions from other nanoparticles can generate a net force.

import gsd.hoomd

import numpy as np

n_rigid = 8
n_apr = 20
n_bpr = 4
n_period = n_apr + n_bpr + 1

fl = 'rigid_local.gsd'
lst_data = gsd.hoomd.open(fl, 'r')

c_tags = [ ind*n_period for ind in range(n_rigid)]
r_tags = [ind+1 for ind in range(n_apr)]
f_tags = [ind+1+n_apr for ind in range(n_bpr)]

quant_hoomd =list(lst_data[0].log.keys())
print('This program checks the relation of forces')
print('Quantities that are stored')
print(quant_hoomd)
print('\n')
print('results per nanoparticle')
print('\n')

n_confs = 10
force_c = np.zeros([3, n_rigid, 3])
force_r = np.zeros([3, n_rigid, 3])
force_f = np.zeros([3, n_rigid, 3])
for ind_q, val in enumerate(quant_hoomd):
    lst_array = n_confs*[0]
    for ind_l in range(n_confs):
        lst_array[ind_l]=lst_data[ind_l].log[val]
    mat=np.array(lst_array)
    if val  == 'particles/md/pair/LJ/forces':
        for ind_p in range(n_rigid):
            ctr = c_tags[ind_p]
            ind_r_ini = ctr + r_tags[0]
            ind_r_end = ctr + r_tags[-1]+1
            fr = np.sum(mat[-1][ind_r_ini:ind_r_end], axis=0)
            ind_f_ini = ctr + f_tags[0]
            ind_f_end = ctr + f_tags[-1]+1
            fl = np.sum(mat[-1][ind_f_ini:ind_f_end], axis=0)
            force_c[0, ind_p][:] = mat[-1, ctr,:]
            force_r[0, ind_p][:] = fr[:]
            force_f[0, ind_p][:] = fl[:]
   if val == 'particles/md/constrain/Rigid/forces':
        for ind_p in range(n_rigid):
            ctr = c_tags[ind_p]
            ind_r_ini = ctr + r_tags[0]
            ind_r_end = ctr + r_tags[-1]+1
            fr = np.sum(mat[-1][ind_r_ini:ind_r_end], axis=0)
            ind_f_ini = ctr + f_tags[0]
            ind_f_end = ctr + f_tags[-1]+1
            fl = np.sum(mat[-1][ind_f_ini:ind_f_end], axis=0)
            force_c[1, ind_p][:] = mat[-1, ctr,:]
            force_r[1, ind_p][:] = fr[:]
            force_f[1, ind_p][:] = fl[:]
    if val == 'particles/md/bond/Harmonic/forces':
        for ind_p in range(n_rigid):
            ctr = c_tags[ind_p]
            ind_r_ini = ctr + r_tags[0]
            ind_r_end = ctr + r_tags[-1]+1
            fr = np.sum(mat[-1][ind_r_ini:ind_r_end], axis=0)
            ind_f_ini = ctr + f_tags[0]
            ind_f_end = ctr + f_tags[-1]+1
            fl = np.sum(mat[-1][ind_f_ini:ind_f_end], axis=0)
            force_c[2, ind_p][:] = mat[-1, ctr,:]
            force_r[2, ind_p][:] = fr[:]
            force_f[2, ind_p][:] = fl[:]

for ind_p in range(n_rigid):
    fc = force_c[0, ind_p]+force_c[2, ind_p]
    fr = force_r[0, ind_p]+force_r[2, ind_p]
    ff = force_f[0, ind_p]+force_f[2, ind_p]
    print('nanoparticle number', ind_p)
    print('force center rigid is the same as bond+lj on constituens','force center', force_c[1, ind_p], 'bond+lj force constituents', fr)
    print('net force is the same as lj only','net force: rigid+flexible bond+lj', fr+ff, 'lj only', force_r[0, ind_p]+force_f[0, ind_p])

[trvsst]

Here is the code it uses hoomd 4.0 or higher and gsd 3.0 or higher.

The file rigid_from_snap.gsd can be visualized in ovito and helps understand the system

The file rigid_local.gsd stores local energies and forces and can be analyzed with the script in the previous message

import numpy as np
import hoomd
import gsd.hoomd

#this system defined a rigid body consisting of n_apr constituents and attaching n_bpr flexible particles at one end
# in this example
# n_apr=20
# n_bpr=4


#kBT 
kT_temp=1.0
# timestep
dt_param=0.001
# time steps
num_steps = 10000
# initial separation
n_separation=2
# lj cut_off
lj_cut = 5.0

# define some numbers
ang = np.pi/3
nm = (1.0/np.sqrt(3), 1.0/np.sqrt(3), 1/np.sqrt(3))
quat = (np.cos(0.5*ang), nm[0]*np.sin(0.5*ang), nm[1]*np.sin(0.5*ang), nm[2]*np.sin(0.5*ang))


cpu = hoomd.device.CPU()
sim = hoomd.Simulation(device=cpu, seed=32)

cell = hoomd.md.nlist.Cell(buffer=0.4)
# lj interactions
lj = hoomd.md.pair.LJ(nlist=cell)

lj.params[('A','A')] = dict(epsilon=1,sigma=1)
lj.r_cut[('A','A')] =  lj_cut

lj.params[('B','B')] = dict(epsilon=2,sigma=1)
lj.r_cut[('B','B')] = lj_cut

lj.params[('A','B')] = dict(epsilon=2, sigma=1)
lj.r_cut[('A','B')] = lj_cut

lj.params[('cA','cA')] = dict(epsilon=0,sigma=0)
lj.r_cut[('cA','cA')] = 0.0

lj.params[('cA','A')] = dict(epsilon=0,sigma=0)
lj.r_cut[('cA','A')] = 0.0

lj.params[('cA','B')] = dict(epsilon=0,sigma=0)
lj.r_cut[('cA','B')] = 0.0

# bond terms
bond_term = hoomd.md.bond.Harmonic()
bond_term.params['AB']=dict(k=5.0, r0=1)
bond_term.params['BB']=dict(k=3.0, r0=1)

# number of rigid bodies
n_r = 8
# number of centers
n_c = n_r
# number of A particles per rigid body
n_apr = 20
# number of A particles 
n_a = n_r*n_apr
# number of flexible particles B per rigid body
n_bpr = 4
# total number of B particles
n_b = n_r*n_bpr
# total number of particles
n_p = n_b + n_a + n_c

position = n_p*[0]

# allocate the rigid bodies
# rigid body tags are n*(n_a+n_b+1) where n = 0, 1, 2, n_r-1
n_period = n_apr + n_bpr +1
for indx in range(n_r):
    ind_r = indx*n_period
    ind_n = (indx+1)*n_period
    # center tag
    position[ind_r]= (3*n_separation*indx, 0, 0)
    # each rifid body is placed on the x-axis
    # particles a,b tag are positioned along the z-axis
    position[(ind_r+1):ind_n]=[(3*n_separation*indx, 0, p-0.5*(n_apr+1)) for p in range(0,n_apr+n_bpr)]

# box size
l_val = 3*n_separation*n_period

# define snapshot
snapshot = gsd.hoomd.Frame()
snapshot.particles.N = n_p
snapshot.particles.types = ['cA','A', 'B']
snapshot.particles.position = position[0:n_p]
snapshot.particles.orientation = n_p*[(1.0, 0.0, 0.0, 0.0)]
snapshot.configuration.box = [l_val, l_val, l_val, 0, 0, 0]
snapshot.particles.typeid = n_r*([0]+n_apr*[1]+n_bpr*[2])
snapshot.particles.mass = n_p*[1]
lst = []
for ind1 in range(n_r):
    lst += [ind1*n_period]*(n_apr+1)+[-1]*n_bpr
snapshot.particles.body = lst
snapshot.particles.moment_inertia = [(0,0,0)]*n_p

# define bonds
snapshot.bonds.N = n_r*n_bpr
snapshot.bonds.types = ['AB','BB']
snapshot.bonds.typeid = n_r*([0]+(n_bpr-1)*[1])

lst = []
for ind1 in range(n_r):
    for ind2 in range(n_bpr):
        pnt = ind1*n_period+n_apr+ind2
        lst.append((pnt, pnt+1))
snapshot.bonds.group = lst

# moment of inertia for center particles
for ind in range(0,n_p,n_period):
    snapshot.particles.moment_inertia[ind]=(1,1,0)
#    snapshot.particles.orientation[ind] = quat

sim.create_state_from_snapshot(snapshot)
g = hoomd.md.constrain.Rigid()
rg.body['cA']={'constituent_types': n_apr*['A'],
        'positions':[(0, 0, p-0.5*(n_apr+1)) for p in range(0,n_apr)],
        'orientations':n_apr*[(1.0, 0, 0, 0)]
}

intg = hoomd.md.Integrator(dt=dt_param, integrate_rotational_dof=True)
intg.forces.append(lj)
intg.forces.append(bond_term)
intg.rigid = rg
'rigid_local.gsd'
rigid_centers_and_free = hoomd.filter.Rigid(('center', 'free'))

#busi  = hoomd.md.methods.thermostats.Bussi(kT=kT_temp)
#nvt = hoomd.md.methods.ConstantVolume(filter=rigid_centers_and_free, thermostat=busi)
mttk = hoomd.md.methods.thermostats.MTTK(kT=kT_temp, tau=100*dt_param)
#mttk.translational_dof = (1.0, 1.0) 
#mttk.rotational_dof = (1.0, 1.0)
nvt = hoomd.md.methods.ConstantVolume(filter=rigid_centers_and_free, thermostat=mttk)

intg.methods.append(nvt)
sim.state.thermalize_particle_momenta(filter=rigid_centers_and_free, kT=kT_temp)
sim.operations.integrator = intg

# save every 10 configurations
trig = hoomd.trigger.Periodic(10)

# ouput intermediate configurations
fl_w='rigid_movie.gsd'
int_wri = hoomd.write.GSD(filename=fl_w, trigger=trig, mode='wb', filter=hoomd.filter.All())
sim.operations.writers.append(int_wri)

# ouput local forces and energies
fl_l='rigid_local.gsd'
logger = hoomd.logging.Logger()
logger.add(lj, quantities=['energies', 'forces'])
logger.add(rg, quantities=['forces'])
logger.add(bond_term, quantities=['energies', 'forces'])
int_local = hoomd.write.GSD(filename=fl_l, trigger=trig, mode='wb', filter=hoomd.filter.All())
int_local.logger = logger
sim.operations.writers.append(int_local)

sim.run(num_steps)

# here we obtain the forces and energies

hoomd.write.GSD.write(state=sim.state, filename='rigid_from_snap.gsd', mode='wb')