Bond forces too strong in Reax/c

Hello,

I have made a model which has 3 layers. The top layer (T) is a CNT-COOH (acid treated carbon nanotube) structure, the middle layer (M) is one chain of cellulose molecule and the bottom layer (B) is again another CNT-COOH. 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

  1. T and B together is called cnt in my input script

  2. M layer is called nfc.

I am using Reax/c potential (ffield.reax.mattsson) using Lammps 14May-2016 version.

The error that I am observing is while applying a velocity (0.06 A/fs) to the T and B layers, the CNT-COOH 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 CNT-COOH.

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.0e-6 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.0e-20 1.0e-15 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 1-2 neighbors
0 = max # of 1-3 neighbors
0 = max # of 1-4 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: 2010-02-19
fix 1 all qeq/reax 1 0.0 10.0 1.0e-6 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.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 = 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.2048587e-17 -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.9475027e-19 -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.239688e-19 -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, next-to-last, final =
-169845.372726 -211129.575245 -211129.575281
Force two-norm initial, final = 13427.3 202.995
Force max component initial, final = 4110.57 45.6169
Final line search alpha, max atom move = 6.85051e-05 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)

Hello

I’m not sure whether your velocity command is correct or not, because I don’t use that command when I move something. But as far as I remember, the velocity command just induce one time initial velocity… How about test fix move command? (http://lammps.sandia.gov/doc/fix_move.html) This command induce continuous movement to your ‘cnt’ group.

Does the A in your velocity unit means Angstrom? If it is, 0.01 A/fs and 0.06 A/fs are very very very fast speed… they are 1000m/s and 6000m/s, respectively. This will bring huge amount of energy to your system… Are you sure your velocity condition is correct?

The main problem is that a thermostat (fix temp/rescale) was also applied to group cnt that effectively slows it movement by cooling the temperature down to 5K. But it is true that fix move is the more appropriate option as it applies a constant moving velocity to the atoms. By only assigning the initial velocity via the velocity set command, velocities on the atoms can change due to interactions with atoms from other groups.

Ray