Making Silicon a rigid body

I am trying to migrate hydrogen atoms deposited on a silicon surface by applying electric field. To make silicon rigid I used the following command:

pair_style reax/c NULL
pair_coeff * * ffield2.txt Si H
fix 15 silicon setforce 0.0 0.0 0.0
fix 16 hydrogen setforce 0.0 0.0 0.0
velocity silicon set 0 0 0
velocity hydrogen set 0 0 0

fix 33 silicon_rigid rigid single Langevin 70 70 2 17891

however, I am getting an error message:

ERROR: Temperature compute degrees of freedom < 0 (src/compute_temp.cpp:100)

My objective is to move hydrogen by keeping silicon rigid. Any insight on how to do this?


I dont know if this fits your case, but there was some previous posts on the forum about the same error in which the problem was related to setting up particles as something other than point particles in the input script: apparently the computation of temperature of LAMMPS is by default on point-particles, so this poses a problem.
For more, see ERROR: Temperature compute degrees of freedom < 0 - #2 by akohlmey. Maybe it has to do with your problem.

By the way, I think that using the “setforce” command to lead to a 0 force at each iteration in your silicon phase (whose atoms already have velocity = 0) + setting it up as a rigid body with a specific thermostat acting on it may not make a lot of sense.

Maybe in your case, you could keep the “velocity silicon set 0 0 0” and do a “fix 33 silicon_rigid rigid single” without assigning equations of motion to your silicon phase. You can even turn off the non-bonded and bonded interactions on it to save some computational time (see neigh_modify command — LAMMPS documentation and delete_bonds command — LAMMPS documentation for more).
How are you setting up the dynamics of the hydrogen atoms?

In case you keep the setforce strategy, you should also be careful because using this “setforce” to set a 0 force over the atoms can lead to some bogus pressure calculation. Although I think that you should be able to surpass that by turning off the silicon-silicon bonded and non-bonded interactions.

1 Like

Hi, thanks for your reply.

I want to deposit hydrogen on silicon surface that forms Si-H bonds. Later on, I want to apply electric field so that hydrogen atoms break bonds with one silicon atom and form bond with another silicon atom. I am applying 0.00001 V/Angstorm (in real unit) electric field with (tilmestep 0.05 and run 0.5 million). However, I am observing that although hydrogen atoms are hopping on the surface, some of the silicon atoms bonded with hydrogen atoms are also coming out of their structure which is undesired for me. I want silicon atoms to retain their bonds with other silicon atoms, but hydrogen atoms break bond with one surface silicon and forms bonds with another.

Also, is there any way to create a silicon body that has a slope on its top surface instead of a rectangular one?

Hi Nazneen,

There are really two possibilities:

  1. For your system of interest in real life, silicon really does leave the surface, and the potential you’re using models this realistically, and you should accept the result.
  2. For your system of interest in real life, silicon doesn’t leave the surface, and your simulation procedure isn’t accurate.

Either way, what you think should happen doesn’t really enter the picture – the simulations are useful, or they are not. Personally I would be very skeptical of the ReaxFF forcefield and would check other publications to see what potentials have been used in the study of similar systems before.


You mustn’t disable interactions between silicon atoms, because it will change the interactions between hydrogen and surface silicon as ReaxFF is a manybody potential.

1 Like

Hmm I see… I never worked with reactive forcefields, so I am not aware of the potential form.
However, if it is the case that simply the forces between Silicon and hydrogen depend on the relative positions of Si-Si, it could be that you can still turn off the forces between Si-Si and everything would be fine, no? I mean, the relative position between Si atoms would still exist.

In any case I suppose he can still follow the advice of simply using “velocity silicon set 0 0 0” and do a “fix 33 silicon_rigid rigid single” without turning off the silicon interactions. I think he would not get any bogus pressure in that case (I dont know cause I never tried using a rigid body with leaving interactions on). He would simply have more computational effort (which apparently seems to be necessary).

But with interactions turned off ReaxFF won’t see that the silicon atoms are close to each other. It will treat each atom as a separate entity, therefore highly reactive.

Maybe then it’s better to also nvm my advice on making the whole Si rigid cause I have no idea on what impacts this could have on the eyes of the reaxff. Better to wait for more suitable advices on how to immobilize the Si atoms for the case of when reactive force fields are being used.

Indeed, if the system is meant to simulate a surface on top of a bulk system, a layered approach is needed. 1-2 immobile layers at the bottom, a thermalized zone on top of that (to model the thermal exchange with the bulk) and a sufficiently large non-thermalized zone at the top. This has been discussed in detail many times in the past, so a search through the archives can be helpful.

@srtee makes some important points. If atoms move in unexpected way, the solution is not to suppress the motions, but to understand why this happens. After all, you want to create a computer simulation and not a computer animation (i.e. you aim for the Nobel Prize and not the Oscar). Thus either the force field parameters are unsuitable (ReaxFF parameters need to be re-trained for different kinds of systems, using other force fields may be difficult, since they usually do not have charged particles and thus cannot react to an external field) or the applied electric field is too strong (a common issue in atomic scale simulation, since physically realistic field strengths often have too little impact for the accessible timescales and thus people have a tendency to boost field strengths, but there are limits to that).

You can create or delete atoms based on regions and you can build complex regions through unions or intersections of multiple primitive regions, so very complex shapes can be created.