Hello everyone! I have a system composed of molecules that consist of 3 beads arranged linearly, where the harmonic angular potential has a equilibrium angle of 180 degrees, and the value of k varies between 6 and 200. After an equilibration protocol, I noticed that when I calculate the angle distribution for the molecules, the value is shifted towards smaller angles, close to 175 degrees, for all values of the k parameter in the angular potential, even for k=200. It seems that the topology with an angle of 180 degrees is penalized. For intramolecular interactions, only the bond potential and the angular potential are considered (special_bond 0 0 0), and for intermolecular interactions, an LJ potential is used.

I can’t find an explanation in this simple model for why most of the molecules do not exhibit angles of 180 degrees for very large values of k. Can anyone think of an explanation?

Maintaining an exact angle of 180 degrees is impossible with the definition of the harmonic angle potential, because the force is ill-defined when trying to project it on the individual atoms.

The problem is that at some point you need to multiply the angle force force with \frac{1}{\sin(\theta)} when projecting it on individual atoms, and that term will overflow for angles close to 0 or 180 degrees. Thus the LAMMPS implementation caps the value of this prefactor. If you look at the source code, you will see a constant SMALL that determines the smallest allowed value (for \sin(\theta)). You can experiment with making that constant smaller, but ultimately, you will not succeed to get a perfect distribution, since the computed forces will overflow.

Hi Axel,
I understand that it should oscillate around its equilibrium value in the harmonic potential, but what seemed strange to me is that its value is practically zero at theta_0. However, now I comprehend where this may come from with the explanation you mentioned about the force calculation and the variable ‘Small,’ which is particularly affected at the extremes (0 and 180 degrees). So, I will experiment on the code variable SMALL.
Thank you for your insight and clarification on this matter!

I’ve been contemplating the angular forces with the internal variable in the lammps code, and I have a question that came to mind: does it make sense to compare the dynamic behavior of the system when varying the A-A-A angle? In other words, do you think the variables ‘s’ and ‘SMALL’ strongly influence the dynamics, particularly as we approach angles close to 180 degrees, compared to the case of 90 degrees?
My intention is to filter out results that are purely derived from programming artifacts and focus on understanding the underlying physical behavior.

Instead of continuously recompiling, you can study the angle distribution of a single molecule (using a Langevin thermostat for ergodicity), at various angle strengths and equilibrium values.

In addition to the code variable SMALL there is also the issue that the bond angle distribution is “censored” in that it’s not easy (or practical) to define an angle > 180 degrees. Thus the Boltzmann distribution “doubles back” on itself and you should expect an expected angle less than your preset for a close enough preset to 180 degrees and a large enough temperature.

I have conducted the analysis of varying different angles and the intensity of the potential for various values of theta_0 and k. However, this analysis was performed under NVP/NPT ensembles and for all molecules in the system. What would be gained by performing the analysis for a single molecule?

Regarding the limiting value of theta_0=180, it is expected that the distribution will exhibit an ‘anomalous’ shift at that value, as you mentioned. However, for example, when theta_0=160 or 170 and k>100, the distributions are very narrow and there are no angles near the 180-degree limit. In this case, my question is: is the difference in the dynamic behavior that I observe for the model with molecules having theta_0=160 or 170, compared to a system with theta_0=90 degrees, an issue related to lammps or a phenomenology to be explored with this model?
There are chemical systems that exhibit this behavior, where significant dynamic changes occur throughout the system due to a topological change in a molecule. So, can I relate what I’m finding with the angles to this phenomenon? Or is there something I might be missing that could be attributed to a lammps artifact?

Thank you for providing your insights and comments!

A single molecule has a configurational energy that you can calculate by hand, U(\theta), and by extension a thermal probability distribution you can calculate by hand Pr(\theta) \propto \exp[-\beta U(\theta)]. Then you can check for yourself if the distribution you observe in LAMMPS respects your known theoretical result.

The equivalent Boltzmann distribution for the angular distribution over a dense collection of molecules is much harder to calculate (not impossible – it’s what the BBGKY hierarchy was invented for – but certainly beyond my theoretical capabilities), which is ultimately why it is hard for you to ascertain whether the problem lies in LAMMPS or in your theoretical prediction (and I would struggle to ascertain that too).

Certainly if you simulate a large enough box to be in the molecular gas phase, the system partition function can be decoupled into copies of the isolated molecular partition function and you can average over molecules for faster statistics.