Problem in equilibrating Cellulose molecule--- molecule bending too much--- Reax/c potential used --- Bond angle and Dihedral angle too low

Hello everyone,

I am having some problem to equilibrate a cellulose molecule using Reax/c potential from mattsson parameter file by LAMMPS 14 May 2016 version. I have seeked help and confirmed from the authors of this paper (http://www.pnas.org/content/112/29/8971) that they have used the mattsson parameter file (even though its not mentioned in the paper and they actually used REAX instead of REAX/C I am using) and so I am using that. But I have also tried my problem with the “ffield.reax.cho” and the results obtained are similar to what I will describe here now using “ffield.reax.mattsson”. I will try to tell my problem in detailed steps.

As many of you may be aware, cellulose is made of long units of glucose rings as shown in the link posted (https://en.wikipedia.org/wiki/Cellulose#/media/File:Cellulose-Ibeta-from-xtal-2002-3D-balls.png). 5 complete units of glucose (and a bit of the 6th) are shown in that link. So I am trying to select each such unit as a group and trying to find out the bond length, bond angle and dihedral angle between them. I am dumping the atomic trajectories and finding the center of mass of each unit and then evaluating the Bond lengths and the angles, a process which is similar to the one described in this paper (http://pubs.acs.org/doi/abs/10.1021/jp5105938)

So to find out each bond length (BL), bond angle (BA) and dihedral angle (DA), I am using cellulose structures having 2, 3 and 4 repeat units of glucose respectively. I am plotting each of the parameters with respect to different frames in the simulations and whenever I am seeing that the parameters are reaching a more or less constant saturation value (which is staying put for a time scale around 30 ps or so), I am surmising that this much be the equilibrium BL, BA or DA.

Using this, my BL for 2 unit cellulose structure is coming around 5.1 Angstorm which compares well with previous literature {(http://pubs.acs.org/doi/abs/10.1021/jp5105938) and (http://link.springer.com/article/10.1007/s10570-014-0481-2))}.

But BA is coming around 80 degrees (instead of around 170 degrees which previous literatures suggest) and DA is coming around 60 degress (instead of around 180 degrees which previous literatures suggest). This means that there my cellulose molecules having 3 and 4 units of glucose is bending a lot. But no atoms are getting lost and the BA or DA values, however less they are, still reach a constant saturation at the values I mentioned with respect to different frames in the simulation. But physically, I think so much bending is unrealistic. So can you please suggest what I should do.

I am posting a sample input file and log file here. Both are from the case of cellulose having 3 glucose units from where I compute the BA.

Thank you so much for the help. I appreciate any feedback which helps me improve my simulation here.

INPUT FILE

dimension 3
boundary p p p
units real
atom_style full
newton on
read_data data3rus.cellulose

pair_style reax/c NULL checkqeq yes safezone 16 mincap 1000
pair_coeff * * ffield.reax C H O
fix 2 all qeq/reax 1 0.0 10.0 1.0e-6 reax/c

neighbor 2.0 bin
neigh_modify delay 0 every 10 check yes

thermo 500
thermo_style custom step temp press vol pe ke etotal
timestep 0.5

min_style hftn
min_modify line quadratic
minimize 1.0e-20 1.0e-15 100000000 1000000000

group bead1 molecule 1
group bead2 molecule 2
group bead3 molecule 3

velocity all create 300.0 4928459 dist gaussian mom yes rot yes loop local units
box

dump 1 bead1 custom 1000 released_bent_equilibrate_bead1.lammpstrj id type x y z
dump 2 bead2 custom 1000 released_bent_equilibrate_bead2.lammpstrj id type x y z
dump 3 bead3 custom 1000 released_bent_equilibrate_bead3.lammpstrj id type x y z
dump 4 all custom 1000 releasedall_equilibrate.lammpstrj id type x y z

fix 3 all nve
run 70000
unfix 3

fix 4 all nvt temp 300.0 300.0 50
run 30000
unfix 4

write_restart equil.restart1

LOG FILE

LAMMPS (14 May 2016)
dimension 3
boundary p p p
units real

atom_style full
newton on

read_data data3rus.cellulose
orthogonal box = (-7.5 -0.5 -3.5) to (11.5 9.5 5.5)
2 by 2 by 1 MPI processor grid
reading atoms …
63 atoms
0 = max # of 1-2 neighbors
0 = max # of 1-3 neighbors
0 = max # of 1-4 neighbors
1 = max # of special neighbors

pair_style reax/c NULL checkqeq yes safezone 16 mincap 1000
pair_coeff * * ffield.reax C H O
Reading potential file ffield.reax with DATE: 2010-02-19
fix 2 all qeq/reax 1 0.0 10.0 1.0e-6 reax/c

neighbor 2.0 bin
neigh_modify delay 0 every 10 check yes

thermo 500
thermo_style custom step temp press vol pe ke etotal
timestep 0.5

min_style hftn
min_modify line quadratic
minimize 1.0e-20 1.0e-15 100000000 1000000000
WARNING: Resetting reneighboring criteria during minimization (…/min.cpp:168)
Neighbor list info …
2 neighbor list requests
update every 1 steps, delay 0 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 12
ghost atom cutoff = 12
binsize = 6 -> bins = 4 2 2
Memory usage per processor = 244.309 Mbytes
Step Temp Press Volume PotEng KinEng TotEng
0 0 55804.077 1710 -6818.6685 0 -6818.6685
300 0 -1362.2427 1710 -7060.1794 0 -7060.1794
Loop time of 55.6193 on 4 procs for 299 steps with 63 atoms

99.8% CPU use with 4 MPI tasks x no OpenMP threads

Minimization stats:
Stopping criterion = trust region too small
Energy initial, next-to-last, final =
-6818.66852702 -7060.17943307 -7060.17943402
Force two-norm initial, final = 573.207 13.4574
Force max component initial, final = 117.105 4.31466
Final line search alpha, max atom move = 0 0
Iterations, force evaluations = 299 3082

MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total

Hello everyone,

I am having some problem to equilibrate a cellulose molecule using Reax/c
potential from mattsson parameter file by LAMMPS 14 May 2016 version. I have
seeked help and confirmed from the authors of this paper
(http://www.pnas.org/content/112/29/8971) that they have used the mattsson
parameter file (even though its not mentioned in the paper and they actually
used REAX instead of REAX/C I am using) and so I am using that. But I have
also tried my problem with the "ffield.reax.cho" and the results obtained
are similar to what I will describe here now using "ffield.reax.mattsson". I
will try to tell my problem in detailed steps.

As many of you may be aware, cellulose is made of long units of glucose
rings as shown in the link posted
(Cellulose - Wikipedia).
5 complete units of glucose (and a bit of the 6th) are shown in that link.
So I am trying to select each such unit as a group and trying to find out
the bond length, bond angle and dihedral angle between them. I am dumping
the atomic trajectories and finding the center of mass of each unit and then
evaluating the Bond lengths and the angles, a process which is similar to
the one described in this paper
(http://pubs.acs.org/doi/abs/10.1021/jp5105938)

So to find out each bond length (BL), bond angle (BA) and dihedral angle
(DA), I am using cellulose structures having 2, 3 and 4 repeat units of
glucose respectively. I am plotting each of the parameters with respect to
different frames in the simulations and whenever I am seeing that the
parameters are reaching a more or less constant saturation value (which is
staying put for a time scale around 30 ps or so), I am surmising that this
much be the equilibrium BL, BA or DA.

Using this, my BL for 2 unit cellulose structure is coming around 5.1
Angstorm which compares well with previous literature
{(http://pubs.acs.org/doi/abs/10.1021/jp5105938) and
(Coarse-grained simulation of cellulose Iβ with application to long fibrils | Cellulose)}.

But BA is coming around 80 degrees (instead of around 170 degrees which
previous literatures suggest) and DA is coming around 60 degress (instead of
around 180 degrees which previous literatures suggest). This means that
there my cellulose molecules having 3 and 4 units of glucose is bending a
lot. But no atoms are getting lost and the BA or DA values, however less
they are, still reach a constant saturation at the values I mentioned with
respect to different frames in the simulation. But physically, I think so
much bending is unrealistic. So can you please suggest what I should do.

I am posting a sample input file and log file here. Both are from the case
of cellulose having 3 glucose units from where I compute the BA.

this is rather useless without the data file. often problems of the
kind you describe are due to issues with the starting geometry and
mistakes in the topology.

the paper you are comparing to describes a densely packed crystal. is
your structure similar dense (via PBC)?
or just the sequence of connected glucose molecules in vacuum?

axel.