Hello everyone,
I want to rotate a polymer molecule using Reax/c potential by LAMMPS 14 May 2016 version. I equilibrate the molecule first by using this input script:
INPUT SCRIPT
dimension 3
boundary p p p
units real
atom_style full
newton on
read_data data4rus.poly
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
group bead4 molecule 4
velocity all create 300.0 4928459 dist gaussian mom yes rot yes loop local units box
fix 3 all nve
run 60000
unfix 3
fix 4 all nvt temp 300.0 300.0 50
run 30000
unfix 4
write_restart equil.restart1
THEN I ROTATE SPECIFIC GROUPS OF THE POLYMER BY READING FROM THE RESTART FILE OF PREVIOUS SIMULATION
read_restart equil.restart1
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
thermo 500
thermo_style custom step temp press vol pe ke etotal
timestep 0.5
group leftbead union bead1 bead2
group rightbead union bead3 bead4
dump 6 all custom 1000 releasedall_twist.lammpstrj id type x y z
fix 7 leftbead move rotate 1.2347501113 2.3272528705 0.7340133595 -2.7806280192 -0.3742854366 3.8499468824 2
run 10000
unfix 7
fix 8 rightbead move rotate 7.5062698084 2.0859288896 -0.4973933233 -1.7017111925 3.7343404625 1.5023427287 2
run 10000
unfix 8
THE PROBLEM IS WHEN I DO THIS, I AM GETTING AN EXTREMELY HIGH TEMPERATURE RIGHT WHEN move rotate starts to work as can be seen in the log file below :
LOG FILE
LAMMPS (14 May 2016)
Reading restart file …
restart file = 14 May 2016, LAMMPS = 14 May 2016
WARNING: Restart file used different # of processors (…/read_restart.cpp:717)
orthogonal box = (-8.5 -2.5 -5.5) to (20.5 12.5 3.5)
2 by 2 by 1 MPI processor grid
84 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
42 atoms in group leftbead
42 atoms in group rightbead
Neighbor list info …
2 neighbor list requests
update every 10 steps, delay 50 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 12
ghost atom cutoff = 12
binsize = 6, bins = 5 3 2
Setting up Verlet run …
Unit style : real
Current step: 90465
Time step : 0.5
Memory usage per processor = 178.43 Mbytes
Step Temp Press Volume PotEng KinEng TotEng
90465 314.56126 4071.2654 3915 -9473.7256 77.824731 -9395.9008
90500 57893281 1.6725898e+08 3915 -7001.7738 14323217 14316215
91000 57893281 1.6725898e+08 3915 -7001.7738 14323217 14316215
91500 57893281 1.6725898e+08 3915 -7001.7738 14323217 14316215
92000 57893281 1.6725898e+08 3915 -7001.7738 14323217 14316215
92500 57893281 1.6725898e+08 3915 -7001.7738 14323217 14316215
93000 57893281 1.6725898e+08 3915 -7001.7738 14323217 14316215
93500 57893281 1.6725898e+08 3915 -7001.7738 14323217 14316215
94000 57893281 1.6725898e+08 3915 -7001.7738 14323217 14316215
94500 57893281 1.6725898e+08 3915 -7001.7738 14323217 14316215
Can someone please point out where I am going wrong? Thank you so much for your help.
Regards,