Hello everyone,
I am trying to run a bead-spring model in Lammps 14May-2016 version at 300K. There is harmonic potential between bonds, angles and dihedrals. Also, between each beads there is a lj12/6 non-bonded potential. The bond, angle and dihedral coefficients match close to this paper (Table 1) (https://link.springer.com/article/10.1007/s10570-014-0481-2)
However, when I am running the simulations, I see that the energy is extremely high and the temperature is not increasing at all. Reducing the time step does not help. For ease, I simplified the system to just 4 beads in a chain and am producing the log file, data file and input script below.
When I am visualizing the trajectories, I am seeing that the beads are going far apart and then again coming back together. But the final position of the chain is perpendicular to the initial position. This definitely should not happen.
Can someone please help me out with this? Thank you so much for your help.
INPUT SCRIPT
dimension 3
boundary s s s
units real
atom_style full
newton on
read_data data.chain
neighbor 2.0 nsq
neigh_modify every 10 delay 10 check yes
##POTENTIAL SPECIFICATION
pair_style lj/cut 7.0
pair_modify shift yes
pair_coeff 1 1 9.02 4.42
bond_style harmonic
bond_coeff 1 90.283 5.6393
angle_style harmonic
angle_coeff 1 366068.3768 163.272
dihedral_style harmonic
dihedral_coeff 1 37442.11445 1 1
special_bonds lj/coul 1 1 1
dump 1 all custom 500 releasedall_equilibrate.lammpstrj id type x y z
velocity all create 300.0 53244 dist gaussian mom no rot no units box
fix 1 all nve
run 100
minimize 1.0e-4 1.0e-6 100 1000
write_restart equil.restart1
thermo_style custom step cpu temp pe ke ebond eangle edihed etotal
thermo 1000
thermo_modify flush yes
reset_timestep 0
timestep 0.25
unfix 1
fix 2 all nvt temp 300.0 300.0 10
run 50000
write_restart equil.restart2
LOG FILE
LAMMPS (14 May 2016)
dimension 3
boundary s s s
units real
atom_style full
newton on
read_data data.chain
orthogonal box = (-5.5 -5 -5) to (23.5 5 5)
1 by 1 by 1 MPI processor grid
reading atoms …
4 atoms
scanning bonds …
1 = max bonds/atom
scanning angles …
1 = max angles/atom
scanning dihedrals …
1 = max dihedrals/atom
reading bonds …
3 bonds
reading angles …
2 angles
reading dihedrals …
1 dihedrals
2 = max # of 1-2 neighbors
1 = max # of 1-3 neighbors
2 = max # of 1-4 neighbors
3 = max # of special neighbors
WARNING: Proc sub-domain size < neighbor skin, could lead to lost atoms (…/domain.cpp:925)
neighbor 2.0 nsq
neigh_modify every 10 delay 10 check yes
#neighbor 0.002 bin
#neigh_modify delay 50 every 10 check yes
##POTENTIAL SPECIFICATION
pair_style lj/cut 7.0
pair_modify shift yes
pair_coeff 1 1 9.02 4.42
bond_style harmonic
bond_coeff 1 90.283 5.6393
angle_style harmonic
angle_coeff 1 366068.3768 163.272
dihedral_style harmonic
dihedral_coeff 1 37442.11445 1 1
special_bonds lj/coul 1 1 1
2 = max # of 1-2 neighbors
2 = max # of special neighbors
dump 1 all custom 500 releasedall_equilibrate.lammpstrj id type x y z
velocity all create 300.0 53244 dist gaussian mom no rot no units box
fix 1 all nve
run 100
Neighbor list info …
1 neighbor list requests
update every 10 steps, delay 10 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 9
ghost atom cutoff = 9
WARNING: Proc sub-domain size < neighbor skin, could lead to lost atoms (…/domain.cpp:925)
Memory usage per processor = 6.71467 Mbytes
Step Temp E_pair E_mol TotEng Press Volume
0 300 -12.840924 99849.356 99839.198 -2.6169035e+10 6.77032e-05
100 1.2745609e+11 0 8.481596e+08 1.9879268e+09 10253693 1312916.1
Loop time of 0.000164986 on 1 procs for 100 steps with 4 atoms
Performance: 52368.189 ns/day, 0.000 hours/ns, 606113.295 timesteps/s
64.9% CPU use with 1 MPI tasks x no OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total