Hi everyone in the Material Science Discourse community,
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.25 - 1.28 g/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 (568.3 KB)
run.equil.npt (2.0 KB)
Thank you in advance for your time and expertise!
Hi @KamdemArnaud,
This question would more likely find an answer in the Science Discussion section of the forum.
But concerning your simulation procedure, here are some insights:
- Your system is very small. 855 atoms is very little at the molecular scale, expect huge fluctuations in pressure and density.
- You have a highly anisotropic system since you are elongated along your molecules axis and did not replicate them much in the plane. As such using an isotropic barostat will perform poorly. Are you sure your pressure converged to equilibrium in all direction (pxx, pyy and pzz)?
- It is unnecessary to perform energy minimization at the end of the finite temperature simulation unless you are interested in a minimum energy gerometry (which is not very interesting at finite T). However, it is important to get statistics of your equilibrium properties like density, pressure tensor and box size lengths to see if you reached equilibrium. I see no such things in your script.
I haven’t opened your files (I am on my phone right now), but I may be able to offer useful information (additional to the ones of Germain) for you to check if your force field parameters and set up are all good (and eliminate this as a possible source of problem).
There is a charmm-gui website, namely https://charmm-gui.org/, which has many tools on it that allow creating input files for classical simulations made using the CHARMM force fields. Specifically in your case, go on “input generator” and “force field converter”. This tool should allow you to create input files specifically for lammps (it will output a 1 molecule topology LAMMPS data file + set of commands (input setup) including the force field parameters). It supports versions 36 and 36m. If you check the initial page of the tool, you can find the publication underlying it. Apparently it concerns a work in which people were evaluating what should be the most suitable set of commands that should be used in many different molecular simulation softwares in order to be able to reproduce the “reference” CHARMM force field results for lipid systems. They mention that there is no reason for the setup to change for other systems (ie not lipid ones). Apparently this was motivated by different researchers finding different results when using different molecular simulation softwares despite the force field being the same (maybe people were using slightly different sets of command lines or something). Anyways, it should function as a reference for you to compare with the one you have.
PS: I think you need to input files for CHARMM (the molecular simulation software I mean). You should be able to get these from the tool found on “input generator” => Ligand & residue builder (or something on these lines).
- Concerning the size of the system i used a small system because i read somewhere in this forum that it is sometimes good to scale down the system in which we already face much difficulties simulating. This is in order to better understand it. Nevertherless, i also equilibrated larger systems of size comparable to the ones in the literature, https://doi.org/10.1155/2023/7890912 and still obtained low densities around 1.32g/cm3.
- I also changed the fix to accomodate anisotropic barostat:
fix fxnpt all npt temp 298.0 298.0 100.0 x 2.0 2.0 500.0 y 2.0 2.0 500.0 z 1.0 1.0 500.0 drag 2.0
but the pressure seems to not converge to equilibrium.
Here is input file :
boundary p p p
units real
neigh_modify delay 2 every 1
atom_style full
bond_style harmonic
angle_style charmm
dihedral_style charmm
pair_style lj/charmm/coul/long 8 12
pair_modify mix arithmetic
kspace_style pppm 1e-6
read_data crystal.data
special_bonds charmm
timestep 0.5
I have a 3-step equilibraion process (energy minimisation,NVT,NPT).
dump dumpeq1 all custom 50 traj_eq1_min.lammpstrj id mol type x y z ix iy iz
thermo_style custom step temp pe etotal lx ly lz pxx pyy pzz vol density
thermo 50
thermo_modify norm yes
– Equilibration: part 1: initial minimization –
minimize 1.0e-5 1.0e-7 100000 400000
undump dumpeq1
write_data system_after_min.data
– equilibration part 2: NVT Ensemble —
timestep 0.5
Give each atom a random initial velocity consistent with a system at 300K.
velocity all create 298.0 12345
dump dumpeq2 all custom 200 traj_eq2_nvt.lammpstrj id mol type x y z ix iy iz
Now continue the simulation using fix nvt (Nose-Hoover).
fix fxnvt all nvt temp 298.0 298.0 100.0
timestep 0.5
run 40000
undump dumpeq2
write_data system_after_eq2_nvt.data
unfix fxnvt
– equilibration part 3: NPT Ensemble —
dump dumpeq3 all custom 200 traj_eq2_npt.lammpstrj id mol type x y z ix iy iz
fix fxnpt all npt temp 298.0 298.0 100.0 x 2.0 2.0 500.0 y 2.0 2.0 500.0 z 1.0 1.0 500.0 drag 2.0
timestep 0.5
run 100000
undump dumpeq3
write_data system_after_npt.data
unfix fxnpt
print “WORK DONE”
Here is my data file :
crystal.data (1.4 MB)
The density after equilibrating is sliightly larger,1.32g/cm3 but still too low compared to the expected value of 1.6g/cm3
Well apart from basic methodology remarks on your simulation, there is little we can do since the investigation relies on what you have actually done compared to the methodology and data available from the paper. As it states, they are available upon reasonable request so you might consider contacting the authors of the article. But that said, without providing any parameter, and with bad reference to the files this kind of paper is what I would then consider as not reproducible.
Like @ceciliaalvares, mentioned, it is reasonable to suspect there are differences in the forcefield you used compared to the one used in the paper and reducing the sources of error by using tools that can directly produce LAMMPS input files can be a good call to understand what went wrong.
1 Like
I am closing this topic here now since it has veered quite a bit off-topic. It continues here
