What does _state.cpu_local_snapshot.particles.net_force store with the rigid constraint?

[XinyanYang95]

Hi All,

I'm running MD simulations with rigid bodies. Each rigid body contains several constituent particles. I want to fix the positions of some rigid particles, so I only added center particles of moving rigid bodies to the integrator. The simulation runs well. Fixed rigid particles are at their initial locations and other rigid particles are moving as expected.

Then I want to collect forces on those fixed rigid particles. I checked _state.cpu_local_snapshot.particles.net_force, there are forces on constituent particles and center particles. However, summing all forces on constituent particles does not equal to forces on center particles. So I'm confused what are the forces on those rigid bodies. Should I sum up forces on constituent particles as well as the center particle to get the total force on a rigid body?

Best,
Xinyan

[XinyanYang95]

BTW, I'm running GPU simulations with -np 1 (single CPU).

[XinyanYang95]

Does the net_force include forces from rigid constraint? If I sum up net_forces of constituent particles and center particles, the forces from rigid constraint will be canceled out?

[joaander]

net_force should be 0 on all constituent particles. constrain.Rigid adds the forces from all constituents to the net force on the central particle. It also zeros the net force on all constituents:

hoomd-blue/hoomd/md/ForceComposite.cc

Lines 748 to 838 in abfaac9

	// sum up forces and torques from constituent particles
	for (unsigned int constituent_index = 0; constituent_index < h_molecule_length.data[ibody];
	++constituent_index)
	{
	unsigned int idxj = h_molecule_list.data[molecule_indexer(constituent_index, ibody)];
	assert(idxj < m_pdata->getN() + m_pdata->getNGhosts());
	assert(idxj == central_idx || constituent_index > 0);
	if (idxj == central_idx)
	continue;
	// force and torque on particle
	Scalar4 net_force = h_net_force.data[idxj];
	Scalar4 net_torque = h_net_torque.data[idxj];
	vec3<Scalar> f(net_force);
	// 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);
	// only add forces for local central particles
	if (central_idx < m_pdata->getN())
	{
	// if the central particle is local, the molecule should be complete
	if (h_molecule_length.data[ibody] != h_body_len.data[type] + 1)
	{
	std::ostringstream error_msg;
	error_msg << "Composite particle with body tag " << central_tag
	<< " is incomplete.";
	throw std::runtime_error(error_msg.str());
	}
	// sum up center of mass force
	h_force.data[central_idx].x += f.x;
	h_force.data[central_idx].y += f.y;
	h_force.data[central_idx].z += f.z;
	// sum up energy
	h_force.data[central_idx].w += net_force.w;
	// fetch relative position from rigid body definition
	vec3<Scalar> dr(h_body_pos.data[m_body_idx(type, constituent_index - 1)]);
	// rotate into space frame
	vec3<Scalar> dr_space = rotate(orientation, dr);
	// torque = r x f
	vec3<Scalar> delta_torque(cross(dr_space, f));
	h_torque.data[central_idx].x += delta_torque.x;
	h_torque.data[central_idx].y += delta_torque.y;
	h_torque.data[central_idx].z += delta_torque.z;
	/* from previous rigid body implementation: Access Torque elements from a single
	particle. Right now I will am assuming that the particle and rigid body reference
	frames are the same. Probably have to rotate first.
	*/
	h_torque.data[central_idx].x += net_torque.x;
	h_torque.data[central_idx].y += net_torque.y;
	h_torque.data[central_idx].z += net_torque.z;
	if (compute_virial)
	{
	// sum up virial
	Scalar virialxx = h_net_virial.data[0 * net_virial_pitch + idxj];
	Scalar virialxy = h_net_virial.data[1 * net_virial_pitch + idxj];
	Scalar virialxz = h_net_virial.data[2 * net_virial_pitch + idxj];
	Scalar virialyy = h_net_virial.data[3 * net_virial_pitch + idxj];
	Scalar virialyz = h_net_virial.data[4 * net_virial_pitch + idxj];
	Scalar virialzz = h_net_virial.data[5 * net_virial_pitch + idxj];
	// subtract intra-body virial prt
	h_virial.data[0 * m_virial_pitch + central_idx] += virialxx - f.x * dr_space.x;
	h_virial.data[1 * m_virial_pitch + central_idx] += virialxy - f.x * dr_space.y;
	h_virial.data[2 * m_virial_pitch + central_idx] += virialxz - f.x * dr_space.z;
	h_virial.data[3 * m_virial_pitch + central_idx] += virialyy - f.y * dr_space.y;
	h_virial.data[4 * m_virial_pitch + central_idx] += virialyz - f.y * dr_space.z;
	h_virial.data[5 * m_virial_pitch + central_idx] += virialzz - f.z * dr_space.z;
	}
	}
	// zero net virial
	h_net_virial.data[0 * net_virial_pitch + idxj] = 0.0;
	h_net_virial.data[1 * net_virial_pitch + idxj] = 0.0;
	h_net_virial.data[2 * net_virial_pitch + idxj] = 0.0;
	h_net_virial.data[3 * net_virial_pitch + idxj] = 0.0;
	h_net_virial.data[4 * net_virial_pitch + idxj] = 0.0;
	h_net_virial.data[5 * net_virial_pitch + idxj] = 0.0;
	}
	}
	}