Equilibrating cellulose crytal to a realistic density with charmm36 force field

Science Talk
1

Hi everyone,

I’m currently working on simulating a cellulose uniaxial test and I’m reaching out to seek some advice from the community.

  • I’ve built a cellulose data file using the CHARMM36 force field. The procedure followed the general steps outlined in this paper: (Molecular dynamics simulation of mechanical properties of carbon nanotube reinforced cellulose | Research Square).
  • To achieve this, I used the cellulose-builder toolkit for constructing the cellulose crystal geometry.
  • Subsequently, I employed charmm2lammps to apply the CHARMM36 force field with the top-all36-carb.rtp and top-all36-carb.prm files.
  • I then attempted an equilibration process involving minimization, followed by an NVT ensemble run for 40,000 timesteps and an NPT ensemble run for 100,000 timesteps.

Unfortunately, the resulting cellulose density after equilibration is around 1.32g/cm³, which falls significantly short of the expected value of 1.6 g/cm³ for cellulose.

I would greatly appreciate any insights or suggestions the community might have to help me achieve the desired density in my cellulose simulation.
Below are my input and data files.
crystal.data (1.4 MB)
run.equil.npt (1.8 KB)

Thank you in advance for your time and expertise!

2

If you have not done so, please read up on a process called “jamming” that can happen to polymers which prevents them from reaching equilibrium, especially when using MD only.

In general, it is always much more difficult to have a system relax to higher density from lower density than the other way around. Once you have reached a certain density, rather complex movements have to happen to make room for another conformational change that increases the density. Those become increasingly less likely and thus you may get stuck.

There are two general ways to get around this:

  1. use a mix of MD and MC
  2. apply a significant pressure and/or run at elevated temperature until you reach a density that is somewhat higher than your desired density, keep it at that for a bit and then remove the excess pressure and temperature and equilibrate to the desired values.

In both cases you may still not get the expected density and in that case you have to carefully review all simulation and force field parameters for errors.