Unexpected behavior while simulating the Temperature vs Density curve of Nitrogen gas

LAMMPS LAMMPS Beginners
,
1

Hi,
I’m trying to simulate the T vs d curve of the Nitrogen gas using a 5 point model of the Nitrogen molecule. This model is made of 2 Nitrogen atoms, and three “ghost sites” ( one in the center and two at the ends). The N2-N2 pair potential is represented in terms of the interaction of these sites. The site-site interaction strengths are taken from this work.

I created a simple model with 1000 such N2 molecules and did rigid/NPT calculations at different temperatures to get the T vs density curve. I also simulated a 3 Point model of the Nitrogen molecule using the files given here (I extracted the tabulated form of the LJ terms given there and used the pairpotential type as tabular instead of using the inbuilt LJ function. ). I then compared the results of these two simulations to the standard data .
T_vs_density.

As we see here the 3 point model is giving excellent agreement to the standard values while the five point model is showing unexpected behavior below a temperature of 250 K. I’m using the same parameters and similar simulation conditions for both the 3 point and 5 point model. However, the 5 point model is giving totally wrong results, despite the fact that the 5 point model gives more accurate potential energy of the N2-N2 interaction.

Here is a comparison of the N2-N2 pair-potential of the two models. We can see that they agree quite well.
Comparison_5p3P.

I tried so many times and can not find out what is wrong with the 5 point model, LAMMPS is not giving any errors or warnings, however the results seems to be inaccurate. For example this is the step-pressure plot of the rigid/NPT simulation of the 5 point model at T = 100K and Pressure = 0.98 atm.
Press_NPT. We can see the pressure oscillating too much. also, when we observe the trajectory file, we can see that after some steps the atoms are not moving around but they are moving back and forth with respect to a point. All this suggest that there is something wrong, can anyone help me identify the issue? Here are the input files for both the 3 point and 5 point models with the respective structure files and the tabulated potential files.

Is there any way to extract the molecule-molecule potential from LAMMPS? I wonder if there is any issue with my tabulated potential for the 5 point model. If I could somehow extract the potential between two rigid molecules, I could verify that the sum of the site-site interactions are leading to accurate N2-N2 potential.

Input files for the 3 point model: SimFiles3P.zip (175.7 KB)

Input files for the 5 point model:
Simfiles_5P.zip (269.5 KB)

I was working on this single issue since almost 10 days and any help and suggestions are highly appreciated!

Update: I was previously using a timestep of 1fs, I repeated the calculations, but the caklculation failed after 150,000 steps and started ouputing NaN.

LAMMPS Version: 7 Feb 2024 - Update 1
OS: Linux “Ubuntu 24.10” 6.11.0-19-generic x86_64

2

You may want to check if a smaller integration time step helps with the 5-point model.

Having said that, molecular dynamics does not have any general guarantee that a “more detailed” interaction parameterisation will do better. In quantum mechanics there are variational principles guaranteeing that a more complete basis set yielding a lower ground state energy always gives an approximation that is closer to the “actual wave function”.

No such guarantee exists in molecular dynamics – rather, MD works because in the regime where molecules exist and matter, most of the funky quantum stuff either is negligible or is close enough to mean field that cheap spherical potentials are good enough. More terms aren’t always better, and more terms sometimes reintroduce interactions that make cancellation of error worse, not better.

As a simple but commonplace example – SPC/E water does not have LJ sites at the hydrogens, while some species of TIP3P do. But SPC/E is as good as TIP3P for most applications while being cheaper, with the exception that protein force fields parameterised against TIP3P will be strictly less accurate if combined with SPC/E. Also, four- and five-point water models exist but have never had sufficient accuracy gains over three-point water to see extensive use.

3

Hi Srtee,
Thank you for your response.

I tried this and repeated the simulation with a timestep of 0.25 fs (previous 1 fs ). In this simulation, the calculation collapsed after around 150,000 time-step and started ouputing NaN as output for all measurables.

molecular dynamics does not have any general guarantee that a “more detailed” interaction parameterisation will do better

Thank you, but I don’t want a more accurate results, I just want a result that is at-least as good as the 3 point model. I’m doing molecular dynamics to study the impact of atomic collisions and the formation of Van der Wals complexes. To accurately model this for the Nitrogen gas system, it is necessary to account for the anisotropy of the potential energy surface caused due to the presence of electron lonepairs.

Also, there is an element of fundamental interest here. Since the simulation is failing so catastrophically by giving NaN error. Understanding the reason behind this is necessary to effectively use LAMMPS.

4

From the data you present, it looks like the 5-point model is not suitable for high-density situations. It is difficult to pinpoint a specific reason for that. Usually, it is just the balance of attractive and repulsive interactions, e.g. between LJ and Coulomb. You have to check its publication for the range of thermodynamic conditions it was parameterized for. Classical potentials are not necessarily transferable to the full range of conditions.

1 Like
5

I’ve checked the paper you referenced for the five site model and it appears to use a Buckingham-type potential, with exponential repulsion and reciprocal sixth-power attraction.

As separation goes to zero an exponential term stays bounded but a reciprocal term does not – therefore such potentials have an unphysical, divergent attraction whenever particles get too close together.

You will have to “patch in” a reciprocal repulsion, which needs to be parameterised. See the references in Buckingham potential - Wikipedia for an example.

1 Like
6

Hi,
Thank you, I’m not using the function given in the paper directly. But I’m using a Tabulated potential based on this paper. I made sure that the cutoff distance is well above this point of attraction. There is no attraction at close range in my tabulated potential.

7

Hi,
Thank you, I’m not using the function given in the paper directly. But I’m using a Tabulated potential based on this paper. I made sure that the cutoff distance is well above this point of attraction. There is no attraction at close range in my tabulated potential.

8

Hi, in the five point model I use, there is no charge and I’m using a tabulated potential form and not using LJ type potential. This 5 point and 3 point model both dies down to zero after a few angstroms. So what is inherently different in the 5 point model that makes it not stable for high density simulations? There should be a reason right? Is there any way to print some quantities that will help me identify this?

9

You are asking this in the wrong place. This is not my area of research and this is not a LAMMPS problem but a problem of the underlying science, i.e. the models you are comparing. That is something you need to discuss with people that know and care about your research. You will not likely find them in a public forum. For example, I don’t want to become your adviser and figure out your problems for you. I understand your desired to resolve this, but you have to understand that this requires effort and time and interest and that is not likely to find here.

10

Well – I have no idea what is happening either. You will have to do basic troubleshooting, such as checking that NVE and NVT integrators run okay. The fact that you have NaN pressures could indicate a problem in your tabulation. Additionally, trying to impose a linear geometry on so many sites in a row could be strongly susceptible to numerical instability.

11

Hi,
I have compared the models extensively myself and I’m not here to achieve that.

My problem is that when using a Tabulated potential form, LAMMPS simulation is breaking down, without giving any clear warning or error messages. I’m trying to understand why is LAMMPS behaving this way. Comparing the two models is not my intention.
This is certainly a LAMMPS issue as the potential I’m using have been extensively analyzed and it is accurate in determining the potential enregy surface of N2-N2 interaction. However, when I use the same potential inside LAMMPS, the simulation is giving ‘NaN’ ouputs,

I’m trying to see what is casing this behavior:

  • does Lammps break down if ghots atoms are placed too close to each other?
  • does it break down if atoms are given too low masses?
  • Is there anything else one should be careful while using ghost atoms?

Again, I’m not looking to compare the effectiveness of the two models, I’m specifically trying to understand what is going wrong when I use the 5 point model. I mentioned the case of 3 point model to let know that a similar approach has known to be effective in LAMMPS.

  • Is there any way to extract the reason for this? For example can I dump the force between two rigid bodies as a function of distance?
  • Or would it help if I redo the simulation with NPT or NVT ensemble instead of the rigid/not rigid/nvt ensemble?
  • AS I indicated, the molecules are not moving around in a free diffusion way, they just move back and forth at a fixed point. Is this bacause of using the ‘fix’ term?
12

Hi, I do not understand what you mean by this. Do you mean designing the N2 molecule as a linear chain is causing this issue? But I using the rigid/npt integrator, so it should not matter right?

13

Sorry, but I disagree. LAMMPS is executing the data you are feeding it and will just do what you tell it to do. Whether something is numerically stable or not is not a LAMMPS issue unless you can prove beyond a doubt that this is caused by some code in LAMMPS not doing what is documented.

But how do we know that you are applying it correctly and that your simulation parameters are suitable and so on? We know that you are using functionality in LAMMPS that has been used a lot over the years and it is to a large percentage covered by our internal testing, so it is highly unlikely (but not impossible) that there is a bug in this part of the code. That makes it more likely a problem of your input. LAMMPS cannot prevent you from providing input leading to unphysical results, it can only check for syntax issues. That makes all issues stemming from the specific details of a model not LAMMPS issues, but user issues.

Ghost atoms are determined from the geometry. If ghosts are too close together, the same would apply to their original non-ghost positions. You cannot avoid ghost atoms, they are part of the domain decomposition algorithm of LAMMPS and if there were such fundamental problems with them, it would have been noticed in the 30+ years that LAMMPS exists.

Low masses of individual atoms matter and determine the choice of time step. For rigid bodies they don’t, since the total mass and the moments of inertia of the rigid body matter. For near-linear molecules (i.e. where not all atoms are exactly on a line within the hardcoded epsilon in fix rigid), the time integration of the rotations can be tricky and require a very short timestep.

You have to create a special input for this placing two molecules in a given orientation and then scanning the potential versus distance using fix move.

Impossible to say from the outside, but your comments indicate a model with “virtual” charge sites and those require using fix rigid or else those virtual sites will “fly away”.

No idea what you mean by this. Since you keep discussing this all “in the abstract” there is not much to comment. And frankly, I personally don’t even care about seeing the details of your inputs since I firmly believe that this is something you have to figure out and not us. Validating models can be frustrating and confusing. Reasons for failures are many.

14

Hi @Planck056,

Please note that the 3 points model you are comparing to in different than your reference. When you say:

You are assuming a conversion to table form of a long-range interaction model is straightforward and working because your 3 points table form gives you agreement with some experimental data on a given range. This is not the general case and even more so when using the same procedure to convert a classical model on the one hand and a quantum one on the other. So stating:

You are actually not very convincing because, when we look at the details, the two source model are rather different, and you are incorrect in assuming this model “has been tested”. You are comparing pears and oranges, while referring to bananas.

Moreover, the ghost sites are supposed to represent charge “location” converting delocalized pseudo-charges to table loses all the interest of moving the charges due to the (non-negligible!) long range effects. If people could simply use table instead of computationally expensive long-range solvers, they would have been doing so since a long time.

It is no surprise if at some point you observe non-physical molecule aggregation (or so do we guess from your density curve) with such workflow. Before a LAMMPS issue, it is more likely that your parametrization procedure is incorrect in the first place and gave you “correct results” for the 3 points model by chance. It needs more thorough thinking and design. You should refer to the literature on how to conceive such models.

In addition to that are the issue @akohlmey and @srtee mentioned regarding the rigid integrator.[1] So this makes it all a complex issue that maybe needs to be boiled down to a more simple case.

  1. I would add a warning with the use of the term “ghost” atoms which is a technical term in the LAMMPS community. ↩︎

15

I highly appreciate your time and willingness to help, I guess my intention is not clear. I’m not trying to compare the two models or find any bugs in LAMMPS. I’m sure I’m doing something wrong somehere and I’m trying to find it. So far I have tried everything and I can not find this.

  • But how do we know that you are applying it correctly and that your simulation parameters are suitable and so on
    Even I do not know, that is what I’m trying to find. I’m trying to find if there is anything wrong in the way I implemented my potential that is resulting in this numerical issue. I’m not trying to claim that Lammps has any bugs. I’m trying to understand which part of my input file us causing the issue. So far I have compared these two models and i can not find what is specifically wrong that is only affecting the 5 point model. I’m looking here for advices that would help me identify the issue in my input files.

  • Ghost atoms are determined from the geometry. If ghosts are too close together, the same would apply to their original non-ghost positions
    I used the word ghost atoms by mistake. I meant ‘fake atoms’, that is molecular sites with no mass but only act as source of force. Does LAAMPS assume an atomic radius even for such fake atoms and does it run into issue if such fake atoms are too close? Again I’m not claiming that there is any issue with LAMMPS.

  • For near-linear molecules , the time integration of the rotations can be tricky and require a very short timestep
    I tried reducing the timestep to 0.25 fs, but in this case the simulation started giving NaN. When the timestep was 1 fs the calculations were running and some values where given, even though they were inaccurate

-You have to create a special input for this placing two molecules in a given orientation and then scanning the potential versus distance using fix move.
Thank you, this is certainly helpful and I will do this.

16

Hi, thank you for your valuable feedback. When I said the models has been tested, I meant they are reproducing the properties of the pair potentials: such as the depth of the energy well, position of this energy well, etc. Actually, I did not create the 3 point model myself. It has been created by this group . I just used that 3 point model as a base, and expanded it to the 5 point model. Unfortunately, it was not working and I’m trying to understand what I’m doing wrong.

You should refer to the literature on how to conceive such models.
The five point model is in good agreement with the literature values of the N2-N2 potential and the position of the well.

So this makes it all a complex issue that maybe needs to be boiled down to a more simple case.

I have tried doing the simulation with fewer molecules, but then the simulation is working but it becomes difficult to tell anything as there will be too much oscillations in the measurable if the number of atoms are too few anyway.

I’m sure there is something wrong in the parametrization, but I do not understand what is wrong. The only way I can test this is by plotting the N2-N2 potential energy curve from my tabulated values and this is working fine. LAMMPS is not giving any errors, so how can I understand what is wrong?

17

No idea what you mean by this
here is a video showing the motion of the atoms in the 5 point model ( we are only visualizing the center point of the N2 molecule here)

As we can see the atoms are not ‘moving’ around as we would see in a typical simulation. They seem fixed at those points and just vibrating back and forth. Also some atoms are forming a conglomeration.

18

I understand your frustration, but - again - this is something you have to work out. It is not a LAMMPS issue but a user issue and you cannot expect people here dig in as deep as needed to sort this out, even more so since you are not even providing means to easily reproduce it.
You have to do what everybody else has to do in in science: eliminate degrees of freedom until you have the minimum system to reproduce the issue (and don’t assume it is just one issue, it could also be a combination of multiples). Then you need to apply “the scientific method™”, i.e. formulate a hypothesis about what could be causing an issue, and then construct an input that would allow to confirm that hypothesis (or not). Rinse and repeat until you have figured it out.

You cannot have tried “everything”. There are far too many permutations and had you done all of those, you would know the cause.

19

I know about the TRAPPE model and used it myself. But you miss my point, you are using a tabulated version of the model that you modified. Not the TRAPPE model. The same goes for your 5 points model.

Have you tried using the same analytical functions, charges and long-range solvers to actually use the published model and not a tabulated version?

20

Hi, thank you for the advice. I would prefer to use the actual functional form, I have read that using the Python form will let us do this. I did not take this route as the manual said it will be much slower than the tabulated form. I will try this and will let you know. Thank you for the suggestion.