Rigid body distorting during ReaxFF SiO2-SiO2 bonding

[scottrb]

Hi all,

I’m encountering some strange behavior that I can’t seem to fix. For background, I’m looking at bonding two SiO2 surfaces together. Using version 11 Feb 2026.

From bottom up, I have the following:

1. Bottom hydroxylated SiO2 “wafer”, with the lower 1 nm of it fixed with setforce 0.0 0.0 0.0 to prevent movement during the simulation.
2. Top hydroxylated SiO2 “wafer”, with the upper 1 nm of it set as a single rigid body to apply a load to compress the two wafers together.

The “free” layers of each wafer are updated with time through fix nvt. I move the top wafer towards the bottom wafer and attempt to bond them by simulating a pressure with setforce or addforce. My current protocol is shown below:

# ----------------------- DEFINE REGIONS AND GROUPS ----------------------------
region bottomrigidlayer block EDGE EDGE EDGE EDGE EDGE 10 # rigid layer, bottom 1 nm of bottom wafer
group bottomrigidlayer region bottomrigidlayer

region toprigidlayer block EDGE EDGE EDGE EDGE 92.969602 EDGE # rigid layer, top 1 nm of top wafer
group toprigidlayer region toprigidlayer

group wafers subtract all bottomrigidlayer toprigidlayer

# ----------------------- EQUILIBRATE AT T = 300 K ---------------------------- 
reset_timestep 0
timestep 0.25

variable forceperatom equal -0.002199 # 10 MPa

compute com_toprigidlayer toprigidlayer com

fix toprigidfix toprigidlayer rigid single force 1 off off on torque 1 off off off 

velocity bottomrigidlayer set 0.0 0.0 0.0
fix freezebottom bottomrigidlayer setforce 0.0 0.0 0.0

fix toprigidcompressionforce toprigidlayer addforce 0.0 0.0 v_forceperatom

fix equilibrate wafers nvt temp 300 300 25

When I use addforce on the rigid layer, it appears to split into quadrants and tilt. Using fix rigid delays this behavior until ~70 ps into the simulation, whereas using fix rigid/nvt or rigid/nve shows the behavior within the first couple ps of the simulation. Using setforce instead keeps the rigid layer parallel to the wafers, however the rigid layer accelerates and “infinitely” compresses the two wafers together until the simulation ends. Examples shown below.

Addforce, throughout the simulation:

Setforce, throughout the simulation:

The only thing I haven’t yet tried is:

velocity toprigidlayer zero linear rigid toprigidfix
velocity toprigidlayer zero angular rigid toprigidfix

, as perhaps this could help to prevent the rotation / splitting. My thought process is that since this didn’t happen with setforce, maybe there are some forces acting on the rigid layer that I have not accounted for.

Is there some way to keep the “pressure” through addforce, without this strange splitting behavior occurring? I have thought about switching to fix move and measuring pressure instead of directly applying it, but this would complicate the methodology for future steps of my simulations.

Thanks,

Scott

[akohlmey]

Note: This response was drafted with the assistance of an AI tool (Claude) after this thread sat without any replies for a few days. While every effort has been made to reason carefully from the available information, AI-assisted answers carry a real risk of error – please treat the points below as a starting point rather than authoritative guidance.

The most likely cause for the undesired behavior are the initial velocities of the atoms in the top rigid layer.

fix rigid computes the center-of-mass velocity and the angular momentum of the body from the velocities of its atoms at the beginning of each run. Your input sets the velocities of the bottom layer to zero, but not those of toprigidlayer, so those atoms still have the thermal velocities from however the system was prepared before. Summed over the layer, they give the rigid body a small, random linear and angular momentum.

With torque 1 off off off you then remove the only thing that could change that angular momentum: the interactions with the wafer below can no longer produce a torque that counteracts the rotation, so the body keeps rotating for the entire simulation. The same applies to the x- and y-components of the center-of-mass velocity with force 1 off off on. The fix rigid documentation has a note about exactly this:

Note that if you expect a rigid body not to move or rotate by using these keywords, you must ensure its initial center-of-mass translational or angular velocity is 0.0. Otherwise the initial translational or angular momentum the body has will persist.

Because the rigid layer is heavy, the rotation is slow, which is why it takes tens of ps to become visible with plain fix rigid. And because the layer spans the periodic box in x and y, any rotation is incompatible with the periodic boundaries: atoms are wrapped back into the box at the boundaries, so the tilted layer appears to break into blocks. If some atoms in the layer crossed a periodic boundary during earlier runs (i.e. have different image flags), the rigid body is assembled from pieces shifted by a box length, and a rotation moves these pieces by very different amounts, which makes this even more pronounced.

The remedy is to start the rigid layer at rest, the same way you do for the bottom layer:

velocity toprigidlayer set 0.0 0.0 0.0

This must be executed before the run command. Since fix rigid re-initializes the bodies at the start of each run (the default reinit yes setting), it does not matter whether it comes before or after the fix rigid command. You can check that it worked by adding the angular velocity of the layer to the thermo output, which should then remain zero:

variable wx equal omega(toprigidlayer,x)
variable wy equal omega(toprigidlayer,y)
variable wz equal omega(toprigidlayer,z)
thermo_style custom step temp pe v_wx v_wy v_wz

A few more comments:

• I would use plain fix rigid and not fix rigid/nvt for the pressure plate. The thermostat of rigid/nvt acts on all translational and rotational degrees of freedom of the body regardless of the force and torque settings, so it would add or remove kinetic energy from the motion of the plate in z-direction, i.e. act like an additional, artificial force on it. Thermostatting the motion of a single heavy body is not very meaningful anyway, and the free atoms are already thermostatted by your fix nvt.

• The behavior with fix setforce is expected: setforce replaces the force on every atom of the layer, including all forces from the interactions with the wafer below. Fix rigid then sums up those forces, so the body only feels the applied load and nothing resisting it, and thus accelerates indefinitely. fix addforce is the correct choice for applying a constant load.

• I reproduced the effect with a small Lennard-Jones model of the same setup (bottom layer frozen, top layer as rigid single with force 1 off off on torque 1 off off off, addforce load, fix nvt on the remaining atoms). With thermal velocities on the rigid layer atoms, the layer keeps rotating and tilts steadily until atoms are lost. After adding velocity ... set 0.0 0.0 0.0 for the rigid layer, it stays flat for the entire run, with both fix rigid and fix rigid/nvt.

[scottrb]

Thanks for the help, Axel. All fixed!