Hello,
I have made a model which has 3 layers. The top layer (T) is a CNTCOOH (acid treated carbon nanotube) structure, the middle layer (M) is one chain of cellulose molecule and the bottom layer (B) is again another CNTCOOH. So T and B are the same. I am trying to pull out the layer T and B out of M and check the forces on M.
For that I have made 3 groups which are

T and B together is called cnt in my input script

M layer is called nfc.
I am using Reax/c potential (ffield.reax.mattsson) using Lammps 14May2016 version.
The error that I am observing is while applying a velocity (0.06 A/fs) to the T and B layers, the CNTCOOH structures is sliding for some time but then the structure is getting crumbled. The regular hexagonal arrangement in the lattice is getting distorted. Also, on applying a lower velocity (0.01 A/fs), using the same input script otherwise, the structure is almost not moving and the visualization in VMD is similar to what will happen to a cylinder when pressed with a force on its two circular ends. I do not understand what is happening and what mistake did I make in my script.
I cheked that one previous literature (http://www.pnas.org/content/112/29/8971) used a much lesser velocity (0.001A/fs) to pull out cellulose from a bundle.
To add to this, I am observing similar phenomenon when using ffield.reax.cho and ffield.reax.rdx. If I remove the M layer, then also I observe similar crumbling to occur to the CNTCOOH.
I am pasting my input script and log file below and also attaching my data file.
Can someone please help me find out the error? Thank you so much for your advice.
Input Script
dimension 3
boundary s s s
units real
atom_style full
newton on
read_data data.CNT_NFC_CNT_full
###SPECIFYING THE GROUPS
region 1 block INF INF INF 5 INF INF units box
region 2 block INF INF 5.5 13 INF INF units box
region 3 block INF INF 13.2 24 INF INF units box
group cnt_upper region 1
group nfc region 2
group cnt_lower region 3
group cnt union cnt_upper cnt_lower
###SPECIFYING THE POTENTIAL
pair_style reax/c NULL checkqeq yes safezone 16 mincap 1000
pair_coeff * * ffield.reax C H O
fix 1 all qeq/reax 1 0.0 10.0 1.0e6 reax/c
neighbor 2.0 bin
neigh_modify delay 50 every 10 check yes
dump 1 all custom 500 releasedall_equilibrate.lammpstrj id type x y z
compute reax all pair reax/c
variable eb equal c_reax[1]
variable ea equal c_reax[2]
variable elp equal c_reax[3]
#variable emol equal c_reax[4]
variable ev equal c_reax[5]
variable epen equal c_reax[6]
variable ecoa equal c_reax[7]
variable ehb equal c_reax[8]
variable et equal c_reax[9]
variable eco equal c_reax[10]
variable ew equal c_reax[11]
thermo_style custom step cpu temp pe evdwl v_eb v_ea v_elp v_ev v_epen v_ecoa v_ehb v_et v_eco v_ew
thermo 200
thermo_modify flush yes
timestep 0.5
min_style cg
min_modify line quadratic
minimize 1.0e20 1.0e15 100000000 1000000000
fix 11 nfc store/force
fix 22 nfc setforce 0.0 0.0 0.0
fix 33 cnt setforce 0.0 0.0 0.0
dump 11 nfc custom 1000 force_disp.myforce id x f_11[1]
initial velocities
velocity cnt create 5 482748 units box
fix 2 all nve
fix 3 cnt temp/rescale 100 5 5 0.02 1.0
#Set velocity in X direction
velocity cnt set 0.06 0.0 0.0 sum yes units box
run 100000
Log file
LAMMPS (14 May 2016)
dimension 3
boundary s s s
units real
atom_style full
newton on
read_data data.CNT_NFC_CNT_full
orthogonal box = (30.779 40.75 25.7885) to (109.221 59.25 25.7885)
1 by 1 by 1 MPI processor grid
reading atoms …
1226 atoms
0 = max # of 12 neighbors
0 = max # of 13 neighbors
0 = max # of 14 neighbors
1 = max # of special neighbors
###SPECIFYING THE GROUPS
region 1 block INF INF INF 5 INF INF units box
region 2 block INF INF 5.5 13 INF INF units box
region 3 block INF INF 13.2 24 INF INF units box
group cnt_upper region 1
508 atoms in group cnt_upper
group nfc region 2
210 atoms in group nfc
group cnt_lower region 3
508 atoms in group cnt_lower
group cnt union cnt_upper cnt_lower
1016 atoms in group cnt
###SPECIFYING THE POTENTIAL
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: 20100219
fix 1 all qeq/reax 1 0.0 10.0 1.0e6 reax/c
neighbor 2.0 bin
neigh_modify delay 50 every 10 check yes
dump 1 all custom 500 releasedall_equilibrate.lammpstrj id type x y z
compute reax all pair reax/c
variable eb equal c_reax[1]
variable ea equal c_reax[2]
variable elp equal c_reax[3]
#variable emol equal c_reax[4]
variable ev equal c_reax[5]
variable epen equal c_reax[6]
variable ecoa equal c_reax[7]
variable ehb equal c_reax[8]
variable et equal c_reax[9]
variable eco equal c_reax[10]
variable ew equal c_reax[11]
#variable ep equal c_reax[12]
#variable efi equal c_reax[13]
#variable eqeq equal c_reax[14]
thermo_style custom step cpu temp pe evdwl v_eb v_ea v_elp v_ev v_epen v_ecoa v_ehb v_et v_eco v_ew
thermo 200
thermo_modify flush yes
timestep 0.5
min_style cg
min_modify line quadratic
minimize 1.0e20 1.0e15 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 = 10 5 2
Memory usage per processor = 317.628 Mbytes
Step CPU Temp PotEng E_vdwl eb ea elp ev epen ecoa ehb et eco ew
0 0 0 169845.37 166333.6 285253.5 11562.739 862.52293 17432.86 57.069022 1.2048587e17 1140.6161 6539.1027 5394.0373 89000.259
200 101.0088 0 210656.25 207575.92 278610.1 1919.8982 552.03905 6835.8967 559.72153 8.9475027e19 866.23725 5296.8837 6078.4291 67773.646
330 192.14328 0 211129.58 208053.64 278580.53 1890.5539 547.34688 6608.8476 554.4239 9.239688e19 872.9901 5146.3114 6095.8767 67638.227
Loop time of 192.143 on 1 procs for 330 steps with 1226 atoms
100.0% CPU use with 1 MPI tasks x no OpenMP threads
Minimization stats:
Stopping criterion = energy tolerance
Energy initial, nexttolast, final =
169845.372726 211129.575245 211129.575281
Force twonorm initial, final = 13427.3 202.995
Force max component initial, final = 4110.57 45.6169
Final line search alpha, max atom move = 6.85051e05 0.00312499
Iterations, force evaluations = 330 1616
MPI task timing breakdown:
Section  min time  avg time  max time %varavg %total
data.CNT_NFC_CNT_full (66.4 KB)