I wanted to model the RMIM cation with a virtual site introduced by Acevedo et al. (Doherty, B.; Zhong, X.; Acevedo, O., Virtual Site OPLS Force Field for Imidazolium-Based Ionic Liquids. The Journal of Physical Chemistry B 2018, 122(11), 2962-2974.), as shown below. In this model, a dummy atom (ε = σ = 0) carrying a negative charge was defined to offload the negative charge to inside the plane of the ring.
To do so, the N1-X bond and the C2-N2-X angle should kept fixed (rigid) during the simulations. I tested applying a large K value for bond and angle constants, but it didn’t work.
Are there any methods or commands in LAMMPS for this application?
I wonder if the geometry might be suitable for you to apply the TIP4P styles, having N3-C2-N1 and X along its bisector. Otherwise, you can try fix shake or fix rigid/*/small to rigidify either the whole molecule or suitable subparts.
You may want to include more information in a subsequent post about the rest of the force field. For example:
what non-Coulomb forces are used, including taper / smoothing / cutoff and combination rules
what angles and bond lengths exactly are fixed, and which ones are flexible
whether the surrounding atoms like N1 are still charged, or neutral
While I’m sure this information is available in the paper, it will be a challenge for volunteers to read through and collate all the information.
Q: I’m afraid fix/rigid or fix/rigid/small are not able to keep a subpart of a molecule rigid. If you think the opposite would you elaborate more?
As mentioned in the original post and noted in the caption of Figure (d and theta are reported there), the LJ parameters for dummy X atoms are ignored. X carries a charge of -0.088 e. All other atoms are treated by an OPLS-style forcefield (of course, the combination rule would be geometric, etc.).
All other bonds and angles except the targeted N1-X and C2-N2-X are flexible (fix shake are applied to H-connected bonds)
All other atoms, including the surrounding N1 atom are charged (e.g., N1 carries a charge of +0.176 in this model)
Thank you for your volunteer contribution. Any other comments from you or others regarding this specific case (rigidifying a bond and angle within a molecule) would be appreciated.
Hello Axel,
It appears I was mistaken! Thank you for your note. I am re-reading the “fix rigid” documentation. Any specific hint to add to the procedure?
BTW, I have already a group of rigid molecules (HF of TIPS model) in my simulation box. Does adding another rigid subpart make the dynamics simulations a challenging task? If yes, do you have any suggestions for that?
This is too generic a question to provide a meaningful answer. What is challenging to one person, can be easy to another. The main general issue with using rigid fixes is that they require collective communications and thus parallel scaling to a very large number of processors is limited. This is somewhat alleviated by the rigid/small fixes at the expense of some additional limitations in how they can be used. This is all discussed in the documentation, so you can see (and decide!) for yourself. For many people those issues are not relevant because of the size of the system and kind of computers they are running on.
