Wondering about LAMMPS issue with npt + npt/rigid

LAMMPS LAMMPS Beginners
1

Hey there,

I recently found myself doing some classical MD simulations of a metal-organic framework (host) loaded with a given amount of N2 molecules (guest). My goal was to simulate the dynamics in the NPT ensemble and evaluate the resulting structure of my system once the equilibrium state had been reached. My initial microstate contained information on position of the framework (based on its crystal structure at ambient conditions) and a given fixed amount of N2 molecules (which were in somewhat “random positions”). I wanted to simulate a point of the experimental isotherm, meaning that my goal was to equilibrate the system containing that given, fixed amount of N2 molecules at the pressure which is observed to be the corresponding equilibrium one at that loading.

In short terms, the potential I was using for my N2 molecules required them to be rigid. At the same time, I wanted the atoms of my framework to move relative to one another. Since I wanted to simulate the dynamics in the NPT ensemble, the first idea that had come to mind was to define a fix npt/rigid to set up the dynamics of the guest and a fix npt to set up the dynamics of my host (example on how this would look like below). However, in the manual there is a mention about how we should not set up a ‘hybrid’ simulation where a part of the system evolves according to npt/rigid and another part with npt, the reason being “because there can only be one fix which monitors the global pressure and changes the simulation box dimensions.”

fix 16 guest rigid/npt/small molecule temp 77 77 100.0 aniso 0.0538 0.0538 1000.0
fix 1 host npt temp 77 77 100.0 aniso 0.0538 0.0538 1000.0
run 4000000
unfix 1
unfix 16

I am aware of the ways I could do my simulation (using for example nvt/rigid for the dynamics of my molecules), but I kind of just wanted to check if I understood the issue with using npt/rigid for one part of the system and npt for the other. As I am not defining some weird partial pressure computation for any of the fixes by means of a fix_modify, the pressure compute inherent to both fixes above would cherish the contributions of all the atoms composing the system. However, the standard of a fix npt/rigid is to decrease the number of degrees of freedom so that each rigid body is accounted as one entity in the computation of pressure. At the same time, the same would not done for the fix npt in the scope of its calculation of pressure (meaning it would account each atom of my rigid molecules of the guest as one atom), right? This is the (only) issue, right? Because be it in the equations shown in Shinoda et al paper (mentioned in the fix npt page of LAMMPS manual) or the basic Nose-Hoover approach to treat the dynamics of the system in the NPT ensemble, I would suppose that as long as I am using the same Text, Pext, Tdamp, Pdamp, I would have the value of my “ficticious variables” as well as their time evolution be the same in the context of both my fixes. Same would therefore hold for the predicted change in volume at each timestep.

PS: By “ficticious variables” I mean those invented to mimic the interaction with the surroundings in the context of simulating the system in an NPT ensemble. For example, looking at the equations for the basic Nose-Hoover approach shown below, it would be zĂȘta and ĂȘta.

image
image1251×475 14.8 KB

Source: Melchionna, S., Ciccotti, G., Holian, B. L. 1993, Molec. Phys., 78, 533. ( https://doi.org/10.1080/00268979300100371)

Thanks for the help,
Best,
Cecilia

1 Like
2

No. The situation is quite simple. There is only one “correct” pressure for your system and that includes contributions from all atoms, regardless of whether they are part of a rigid body or not. Similarly, the kinetic component of the pressure (i.e. the temperature) must be that for the entire system. Consequently, to adjust the system to react to the difference between target and actual pressure you need to adjust the box size, but there is only one box and any change to the box will affect all atoms. Thus, there must be only one fix that adjusts the box size in each dimension (i.e. you may have fix npt adjust x and y and fix deform adjust z, but you cannot have fix npt and fix npt/rigid adjust x, y, or z at the same time). At the same time, you must not have multiple fixes do time integration (i.e. update velocities and positions) of the same atoms. Thus you can have only one barostat in your system. Full stop.

1 Like
3

Okay, I do understand what you said (hopefully). But if the two fixes are computing the pressure in the precise same way (since both of them concern all atoms because I am not specifying some fix_modify to any of them) and the parameters of the barostat are precisely the same, both of these fixes would converge to a same set of equations adjusting the dimensions of my box, right? Therefore I would not have some “unphysical” consequence of volume change suffered by any given atom, no? (unless one of the fixes is reducing the degrees of freedom of the rigid particles and the other isnt)

Also, I am not having multiple fixes do time integration on the same atoms (as you can see above, one fix would be acting in the molecules and the other in the framework)

1 Like
4

Yes, you would. In fact, it is likely that your simulation will crash. You are looking at this with the rationale of a human, you need to keep in mind that computer programs do not reason. They just blindly execute what is programmed. When you have two fixes changing the box dimension , each will maintain its own fictitious degrees of freedom of Nose-Hoover chains, each will scale the box and then the result will be inconsistent with each fix’es internal state, which will result in incorrect box adjustments and so forth.

The rule is simple and there is no arguing about it. Only one fix may modify the box.

What is the problem in having fix rigid/npt in combination with fix nvt?

5

There is absolute no problem with me having fix rigid/npt in combination with fix nvt (or vice versa), except that I have already run several simulations, which by now amount for a huge computing time and also analysis-time, for my PhD using fix rigid/npt to evolve my guest and npt to evolve my framework :smiling_face_with_tear: :smiling_face_with_tear:
Therefore if indeed there was something serious with npt/rigid for a part of the system and npt for the other, I would be obligated by my own sense of ethics to redo everything, which will take me ages (and also my supervisor will love it).

I actually followed my reasoning (the one I just told you) when I was setting up the methodology to address the dynamics of my system and I never even dreamed that there would be any whatsover problem with the barostatting since all the parameters of both fixes were the same and the calculation of pressure inherent to each of them would converge to the same thing in the sense that both of them would concern all atoms. I stumbled into the NPT + NPT/rigid problem 3 days ago for another thing I was checking


On the weekend, I selected some random simulations to redo using NPT for the host and nvt/rigid for the molecules to see if I would see some difference. Actually the simulations do run quite nicely (they dont crash at all) and the resulting structure was exactly the same in this new approach as in my initial approach (just like the volume of the box).

So, this was good. Still, I thought of checking with you to hear what you would say about it.

If the volume calculation of one of this fixes was overwriting the one of the first, it would really save the day, so I was hoping you would tell me that :grimacing:
If not, I will suck it up and just redo my calculations to be sure it’s all good

6

I am very surprised you received the exact same results. Since either of the fixes would scale the box, basically your box and pressure fluctuations should be much larger with two npt fixes. You may just have gotten very lucky. In other cases where people made the same mistake in their input and asked about this here, it often lead to the simulations crashing with some error about non-numeric box dimensions or similar.

What surprises me even more is that LAMMPS should have printed some quite visible warning about this (assuming you are using a reasonably recent version of LAMMPS).

I cannot see your results, don’t know your research, and thus I will not give you a thumbs up or down on whether you can salvage the time you already spent on this.

If there were not a few legitimate use cases of two box changing fixes at the same time (e.g. with fix npt + fix deform as I mentioned earlier), LAMMPS would not just print a warning but stop with an error.

7

Yeah, it runs quite well and there is no warning, but it’s okay. I will run all the simulations again with the npt + nvt/rigid to be sure that structure and volume are the same for every single case, since we have no way to be sistematically sure. Night is a child

Thanks anyways for the help!!

PS: I am sending a google drive link with a sample of input + log.lammps file just in case you want to see by yourself that it doesnt crash neither that there is a warning.
https://drive.google.com/file/d/1VqpDXw5uwV_wHx9bUjlOPncHR6wwWRKh/view?usp=sharing

I am running it on the 7Aug2019 version. Maybe I could have used a newer one. Oh well


8

I would recommend using rigid/npt plus nvt. That ensures that the rigid bodies are properly moved when the box is scaled.

9

I will try your suggested approach (rigid/npt plus nvt) and mine (with nvt/rigid plus npt) and then re-edit this comment pointing out differences I observed in structure and volume within the two approaches once I am finished (I mean just for the sake of possible social utility for someone else setting a similar methodology and reading this post).

I think in my case though, if I had then to choose to act on momentums of atoms belonging solely to one of the two groups group (i.e., those of the host or of the guest) for the sake of regulating the pressure of the system, I would go for the host since in my case I wanted to see if my classical model would picture this “phase transformation” of the host that is observed experimentally when it is at high pressure and high loading.
(I do see the point of the justification for using your approach though. And I do understand that in the end of the day, in the two approaches, the total pressure of the system (which cherishes contribution of all atoms), would be the one that I set up in the barostat anyways as well as that the atoms of the host will ultimately feel the effects of me acting on the momentum’s of the guest’s atoms as they interact).

10

For the sake of completeness. There is a discussion toward the end of the fix rigid command — LAMMPS documentation page about this scenario.

11

Thanks Axel, I believe I saw that page. In this thread I was just trying to understand more specifically what LAMMPS does when I do a npt for a set of atoms of my system and npt/rigid for another.

Anyways, I want to finish some other stuff before I go back and redo every single simulation I had done in the past with the misguided npt + npt/rigid approach, but for the sake of not coming back to the thread 40 years later and allowing for some closure, I did the analysis for 5 different classical models I was investigating for my system (composed of host + guest). For each of these models, I evaluated total volume of my simulation box, total potential energy and structure (by means of RDFs) for the three set-ups: (i) npt/rigid (for the guest) + nvt (for the host), (ii) nvt/rigid (for the guest) + npt (for the host) (iii) npt/rigid (for the guest) + npt (for the host).

Below I am leaving the results of average values of these properties (ofc, after my system has well equilibrated - I evaluated the evolution of properties with time) and the respective standard deviations. When evaluating ‘alikeness’ of RDF curves for the three possibilities, I checked it for every single existing pair of atom types X-Y in my system and, as it turns out, all of them are ultra equivalent within the three approaches (set-ups previousy mentioned) within each and all the five models. For what is worth it, I included first neighbours in the RDF check-up whenever the given pair happens to be bonded. Since this amounts to many figures, I am posting a random one that is the one I found the most mismatching (qualitatively speaking, judged by my eyes) - as you can see the curves are quite alike as the difference is solely in peak amplitude and ultra tiny.
PS: Ep in kcal/mol, V in angs^3

So ultimately it seems that I had my luck with the 2019 version of LAMMPS upon having performed npt/rigid for a part of my system and npt for the other. I hope the luck stands when I assess all the other classical models I investigated later on!

image
image1067×441 35.1 KB

Thanks again,
Best