Hello everyone,
I am having some problem to equilibrate a cellulose molecule using Reax/c potential by LAMMPS 14 May 2016 version. I will talk about the simplest case of cellulose having 2 units of glucose rings consisting of 42 atoms but my problem is consistent with higher length scales. For convenience I have attached an image of my system below (2ru.png).
I saw this paper (http://www.pnas.org/content/112/29/8971.abstract) where the authors have used timestep of 0.5 fs but with Reax. Since Reax is getting outdated , I stuck to reax/c and used timesteps 0.05 and 0.5 fs.
First I got segmentation fault which I corrected by checking the mailing list and introducing “safezone 16” and “mincap 1000” to my pair_style reax/c. Should these values be more??? How do I judge what the value will be for a specific system?
My input file looks like this :(also attached below)
INPUT FILE
dimension 3
boundary p p p
units real
atom_style full
newton on
read_data data2rus.cellulose
pair_style reax/c NULL checkqeq no safezone 16 mincap 1000
pair_coeff * * ffield.reax C H O
neighbor 2.0 bin
neigh_modify delay 50 every 10 check yes
timestep 0.05
min_style hftn
min_modify line quadratic
minimize 1.0e-20 1.0e-15 100000000 1000000000
group bead1 molecule 1
group bead2 molecule 2
velocity all create 300.0 4928459 dist gaussian mom yes rot yes loop local units box
dump 1 bead1 custom 20000 releasedall_equilibrate_bead1.lammpstrj id type x y z
dump 2 bead2 custom 20000 releasedall_equilibrate_bead2.lammpstrj id type x y z
dump 3 all custom 20000 releasedall_equilibrate.lammpstrj id type x y z
fix 1 all npt temp 300.0 300.0 0.1 iso 0.0 0.0 1000
run 5000000
unfix 1
fix 2 all nvt temp 300.0 300.0 50
run 5000000
unfix 2
fix 3 all nve
run 1000000
I am constantly getting two errors :
-
While running the simulation for timestep 0.05, the simulation is running but the structure is bending and twitching a lot as I observed from the trajectories. There were no lost atoms though. I am attaching that trajectory file here titled “releasedall_equilibrate.lammpstrj”
-
I also tried changing the minimize style from cg to hftn but still I don’t think I am getting correct dynamics. While running the simulation for timestep 0.5 fs, I am getting an output like that is attached below. I am also attaching my data file. Can you please help me out as to why this is happening?
Thank you so much for your time and help.
LOG FILE
LAMMPS (14 May 2016)
Reading data file …
orthogonal box = (1.5 -6 -14.5) to (14.5 4.5 0.5)
2 by 1 by 2 MPI processor grid
reading atoms …
42 atoms
Finding 1-2 1-3 1-4 neighbors …
Special bond factors lj: 0 0 0
Special bond factors coul: 0 0 0
0 = max # of 1-2 neighbors
0 = max # of 1-3 neighbors
0 = max # of 1-4 neighbors
1 = max # of special neighbors
Reading potential file ffield.reax with DATE: 2010-02-19
WARNING: Resetting reneighboring criteria during minimization (…/min.cpp:168)
Neighbor list info …
1 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 = 3 2 3
Setting up hftn style minimization …
Unit style: real
Memory usage per processor = 14481.6 Mbytes
Step Temp Press Volume PotEng KinEng TotEng
0 0 37881.052 2047.5 -4253.6126 0 -4253.6126
50 0 4739.5241 2047.5 -4418.512 0 -4418.512
100 0 -1555.7856 2047.5 -4465.0979 0 -4465.0979
145 0 -1473.3792 2047.5 -4465.0982 0 -4465.0982
Loop time of 131.148 on 4 procs for 144 steps with 42 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 =
-4253.6125843 -4465.0981979 -4465.0981979
Force two-norm initial, final = 637.488 6.74525
Force max component initial, final = 196.553 3.38862
Final line search alpha, max atom move = 0 0
Iterations, force evaluations = 144 2358
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
in.equilibrate (986 Bytes)
data2rus.cellulose (2.5 KB)
releasedall_equilibrate.lammpstrj (192 KB)