Hello everyone,
I am trying to equilibrate a cellulose molecule (just 2 repeat units) using REAX/c (mattsson) potential in Lammps 14May-2016 version. A similar paper doing this is given here (https://link.springer.com/article/10.1007/s10570-014-0481-2).
My question is with the intermolecular van der waals forces and hydrogen bonding captured in this paper by a standard lj12/6 potential. The paper quotes " The intermolecular potential must describe van der Waals interactions and hydrogen bonding. Because the coarse-grained beads are neutral, electrostatic interactions are not considered. We represent the van der Waals interactions with a 6-12 Lennard-Jones potential"
I also want to capture the vdw and hydrogen bond energy by a standard lj potential. When I am carrying my simulations,the value of epsilon is much much higher. My value of epsilon is coming around 60kcal/mole (~2.4 ev) using the hydrogen bond energy contribution of REAX force field. The value in the existing paper came around 0.1 kcal/mole. This is orders of magnitude higher and made me think of the effectiveness of my results. The sigma value is comparable though (5.04 Angstorm). We can output a separate hydrogen bond energy using Reax, which the existing paper wasn’t using, but will this value be so different? Also the vdw forces are much higher (~1900kcal/mole). Can someone please help me understand where I am going wrong?
I am attaching the input files and log files for your understanding.
Thank you so much for your help.
Regards,
Input file
dimension 3
boundary s s s
units real
atom_style full
newton on
read_data data.2beads
###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
###SPECIFYING THE GROUPS
region 1 block -5.0 0.15 -5.0 5.0 -5.0 5.0 units box
region 2 block 0.3 5.0 -5.0 5.0 -5.0 5.0 units box
group bead id <= 42
group bead1 region 1
group bead2 region 2
displace_atoms bead2 move 0.5 0 0 units box
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 500
thermo_modify flush yes
timestep 0.5
dump 1 bead1 custom 500 releasedall_equilibrate_bead1.lammpstrj id type x y z
dump 2 bead2 custom 500 releasedall_equilibrate_bead2.lammpstrj id type x y z
dump 4 bead custom 500 releasedall_equilibrate_bead.lammpstrj id type x y z
min_style cg
min_modify line quadratic
minimize 1.0e-20 1.0e-15 100000000 1000000000
velocity all create 300.0 4928459 dist gaussian mom yes rot yes loop local units box
fix 2 all nvt temp 300.0 300.0 50
run 50000
unfix 2
write_restart equil.restart1
Log file
LAMMPS (14 May 2016)
dimension 3
boundary s s s
units real
atom_style full
newton on
read_data data.2beads
orthogonal box = (-5.5 -5.5 -4.5) to (5.5 3 4.5)
1 by 1 by 1 MPI processor grid
reading atoms …
42 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 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
###SPECIFYING THE GROUPS
region 1 block -5.0 0.15 -5.0 5.0 -5.0 5.0 units box
region 2 block 0.3 5.0 -5.0 5.0 -5.0 5.0 units box
group bead id <= 42
42 atoms in group bead
group bead1 region 1
21 atoms in group bead1
group bead2 region 2
21 atoms in group bead2
displace_atoms bead2 move 0.5 0 0 units box
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 500
thermo_modify flush yes
timestep 0.5
dump 1 bead1 custom 500 releasedall_equilibrate_bead1.lammpstrj id type x y z
dump 2 bead2 custom 500 releasedall_equilibrate_bead2.lammpstrj id type x y z
dump 4 bead custom 500 releasedall_equilibrate_bead.lammpstrj id type x y z
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 = 2 2 1
Memory usage per processor = 18.3661 Mbytes
Step CPU Temp PotEng E_vdwl eb ea elp ev epen ecoa ehb et eco ew
0 0 0 -4454.1246 -4213.8495 -6490.7237 -51.936002 0.35154844 82.751806 0.0019810209 0 -23.720825 143.24071 -14.868565 2141.0535
200 1.458657 0 -4666.7758 -4431.8118 -6452.1606 -108.23699 7.7570673 59.108315 0.004845264 0 -26.835471 139.09144 -11.93034 1961.39
Loop time of 1.45871 on 1 procs for 200 steps with 42 atoms
100.0% CPU use with 1 MPI tasks x no OpenMP threads
Minimization stats:
Stopping criterion = linesearch alpha is zero
Energy initial, next-to-last, final =
-4454.12456652 -4666.77584578 -4666.7758412
Force two-norm initial, final = 637.489 39.6633
Force max component initial, final = 200.621 19.0886
Final line search alpha, max atom move = 1.34355e-11 2.56464e-10
Iterations, force evaluations = 200 1161
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
data.2beads (2.48 KB)
in.2beads (1.61 KB)
ffield.reax.mattsson (14.8 KB)