Fix Rigid Superdiffusive Behavior

LAMMPS LAMMPS Mailing List Mirror
1

Hello,

I’m simulating a rigid mixture of CO2 and Acetonitrile (primarily ACN, with a little bit of CO2). I’ve started by running phase equilibria calculations in another program (GOMC, using Gibbs Ensemble Monte Carlo), and then I have been starting from the compositions calculated from those calculations and running lammps MD simulations to calculate the diffusion coefficient. Eventually, I want to include ions in the calculation, but I’m currently working on the simple binary system. The densities that I am calculating from lammps match well with my GEMC calculated densities.

I’m getting a strange behavior in the diffusion coefficient, where for some portions of the trajectory the calculated mean-squared displacements have a linear regime (slope of 1 on a log-log plot) and for other portions it has a non-linear regime (slope of like 1.3 on a log log plot). There doesn’t seem to be a clear indicator to me of what the cause might be. I’ve experienced this issue on both the 7Aug2019 version of lammps, as well as the 3Mar2020 version. This behavior, when it happens, typically occurs later on in the trajectory after quite a long period of time (usually sometime over a nanosecond into the simulation). I’ve looked through the mailing list and it seems the general advice for fix rigid in these scenarios is use a smaller timestep; however, as I mention below, that didn’t fix my problem. Both of my molecules are linear, 3 site models, so fix SHAKE isn’t an option in this case.

I’ve included a plot of the MSD versus time calculated for 2 ns chunks from my simulations, as you can see, the grouping on the bottom is linear as expected, the grouping on the top is very much not.

Things I’ve tried:

  • Doubling/halving the system size
  • Halving the simulation timestep from 1.0 fs to 0.5 fs
  • Tried running with both rigid/nvt and rigid/nvt/small variants
  • Running with a separate rigid fix for each molecule type, instead of a single rigid fix for both
  • Doubling the output frequency of the dump/outputting in wrapped/unwrapped coordinates (to eliminate coding issues with the msd calculation)
  • My MSD code works for non-rigid systems in reproducing the diffusion coefficient
    I’ve attached a sample set of input files a binary simulation test set that I have been using, with 999 ACN molecules and 1 CO2 molecule. The test set that I am providing has no ions in it, but does my whole procedure of 500 ps of NVT equilibration, 1 ns of npt equilibration after which it calculates the average box size for 1 ns, and then a 500 ps nvt equilibration again at the new box size before a 5 ns production run (I include it all in case there is an issue with using change-box with fix rigid or some other problem I might be overlooking). I would appreciate if someone had a suggestion for more tests that I could pursue to trying and figure out the cause of the problem, or if this is something anyone has seen before?

I’ve been challenged by this for a few weeks, so I really would appreciate any insight into this that anyone might have! Thanks so much in advance! I can also supply my MSD code if that would be helpful, but I don’t expect anyone to want to look through it (and at this point I’m leaning towards the problem being with something I’m doing in my simulations rather than in my code).

Best,
Zeke

msd.png

data.lmps.gz (78.2 KB)

lmps.include.coeffs.gz (342 Bytes)

in.acncx.gz (758 Bytes)

2

when looking for the optimal time step, you should run without a thermostat to be able to monitor energy conservation. depending on the conditions even 0.5fs may be too large. the time integration of the rotational degrees of freedom is much more sensitive than the time integration of the translational DOFs, so time steps that are possible with fix shake are not possible with fix rigid.

another thing to look out for would be whether you have a (different) center of mass drift. that can be removed with the velocity command at or near the end of the equilibration or in pathological cases using fix momentum. whether you get a COM drift and how much also depends on the time step, but can also be influenced by the choice of origin and (initial) composition of your system.

axel.

3

Hi Dr. Kohlmeyer,

First I want to thank you for your previous reply - I did find that energy was not being conserved when I ran an NVE run with the 1.0 fs timestep, so I appreciate you pointing this out.

I was curious about the origin of the sensitivity in the rotational degrees of freedom that you mentioned so I started running some tests to try and understand it. My advisor and I (cc’ed) were especially confused after reading Kamberaj et. al which describes how the rotational degrees of freedom are handled in the rigid algorithm used by lammps, with a symplectic integrator neither of us would expect it to be so timestep dependent. We started to suspect that it was possibly an issue with how linear molecules are handled in fix rigid, and we would appreciate hearing your perspective on this. Again, we are using the 3Mar2020 version of lammps, and the input decks mostly remain the same as in my previous email with only the few modifications I have detailed below.

To summarize what we did:

  1. I ran a test case with using fix/rigid/small with spc/e water and a 1.0 fs timestep, the energy conservation for the nve run was perfectly fine, and I didn’t see any of the super diffusive behavior I was seeing with my rigid linear molecules when I calculated the mean-squared displacement. I ran this test case multiple times, with multiple random number seeds for the velocity creation and get similar results each time.

  2. We double checked our input deck that we sent to you in our last message, and focused on the linearity of the molecules. In our initial input, the molecules were slightly off linear (packmol, it seems, doesn’t preserve perfect linearity when it generates configurations if you allow it to rotate the molecules to put it in the box (179.99X degree angle), we instead generated a new initial configuration that had all the molecules initially oriented horizontally, but now with perfect linearity (180.0 degree angle)). We then calculated these angles through the course of our simulation run with fix rigid/small using both the nvt and nve variants. We found that in this case the angle was not 180 degrees, but was something like 179.99X degrees after even 1 step of the simulation. Importantly, this was even the case when we ran with a 0.1 fs timestep.

  3. We then calculated the inertia tensor using the following command,

compute chk all chunk/atom molecule

compute moi all inertia/chunk chk

fix wmoi all ave/time 1 1 1 c_moi[*] file tmp.out mode vector

and then diagonalized the matrix using python. We then plotted the average smallest principal moment of inertia as a function of timestep and found that with a 1.0, 0.1, and 0.01 fs timestep the smallest moment does not remain zero as we would expect, but instead starts from zero (as it would have to since our initial configuration was perfectly linear) and then grows to a small constant value, consistent with the non-linear angle that I was seeing above. I’ve included an example plot of this behavior.

  1. It was only later that I discovered reading the documentation how to directly output the principal moments directly from the rigid fix, using the following command,

compute rmi all rigid/local 1 inertiax inertiay inertiaz
dump mom all local 1 tmp.dump c_rmi[1] c_rmi[2] c_rmi[3]

but in this case, even when doing both 3) and the above command in the same run, I get in this case the expected result where one of the principal moments is exactly zero and the other two are constant. This however, doesn’t match with the results I get when I diagonalize the matrix (which I’ve checked both by hand and by code) or the fact that the molecules in my outputted trajectory don’t remain totally linear. Are these being computed at different points in the integration?

Anyway, we are a bit confused and were hoping that you might be able to shed some light on this, especially on the conflicting output between my points in 2) and 3) and my point 4). Is it possible that there is a bug somewhere in the rigid algorithm applying it to 3 atom linear molecules?

Thank you so much,
Zeke

out_1.png

4

It is unlikely (although not impossible) that there is a bug in the
rigid body integrator. please also note, that there are multiple
"flavors" that share code to a different degree.

whether your rigid objects are detected as linear depends the
precision of the geometry you use and on the constant EPSILON in the
rigid_const.h header file. if need be, you may change that constant,
or add some debug output to the setup_bodies_xxx() functions in the
rigid classes of the LAMMPS code.

when writing output, your precision may be limited (e.g. when using
default floating point output precision, which writes out just 6
digits). this may affect your manual diagonalization and may be a
problem with packmol, e.g. if you are using PDB files which use a
notoriously low precision for storing positions. you may want to look
for alternatives for the latter, e.g. using the create_atoms command
in LAMMPS itself in combination with a molecule file.

axel.