Bond style hydrid using harmonic and table

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. I define additional bonds using “fix bond/create” (Reference - http://lammps.sandia.gov/threads/msg59117.html) and “special_bonds extra” and use “bond_style hybrid” to reproduce the harmonic and l/j interactions between beads.

However when I am running NVT, I see that the temperature and the Ebond energy is not rising from 0. The visualized trajectory shows a chain not moving at all.

I am attaching the input script and log file here along with the data file. Can someone please help and tell me where I am going wrong? Is the sytax for bond_style hybrid or special_bonds or bonds/create wrong? Or is it something else? Please do provide your feedback. Thank you so much

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

fix 1 all bond/create 1 1 1 6.5 1

##POTENTIAL SPECIFICATION

pair_style lj/cut 12.0
pair_modify shift yes
pair_coeff 1 1 0 4.42

special_bonds lj/coul 1.0 1.0 1.0 extra 1

bond_style hybrid harmonic table spline 51
bond_coeff 1 harmonic 90.283 5.6393
bond_coeff 1 table ljvalues.table INTRABEAD

angle_style harmonic
angle_coeff 1 366068.3768 163.272

dihedral_style harmonic
dihedral_coeff 1 37442.11445 1 1

min_style cg
min_modify line quadratic
minimize 1.0e-20 1.0e-15 100000000 1000000000

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

thermo_style custom step cpu temp pe ke ebond eangle edihed etotal
thermo 1000
thermo_modify flush yes
reset_timestep 0
timestep 0.25

fix 3 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 2 bin
#neigh_modify delay 50 every 10 check yes

fix 1 all bond/create 1 1 1 6.5 1

##POTENTIAL SPECIFICATION

pair_style lj/cut 12.0
pair_modify shift yes
pair_coeff 1 1 0 4.42

special_bonds lj/coul 1.0 1.0 1.0 extra 1
2 = max # of 1-2 neighbors
3 = max # of special neighbors

bond_style hybrid harmonic table spline 51
bond_coeff 1 harmonic 90.283 5.6393
bond_coeff 1 table ljvalues.table INTRABEAD
WARNING: 1 of 51 force values in table are inconsistent with -dE/dr.
Should only be flagged at inflection points (…/bond_table.cpp:380)

angle_style harmonic
angle_coeff 1 366068.3768 163.272

dihedral_style harmonic
dihedral_coeff 1 37442.11445 1 1

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 = 14
ghost atom cutoff = 14
WARNING: Proc sub-domain size < neighbor skin, could lead to lost atoms (…/domain.cpp:925)
Memory usage per processor = 6.9135 Mbytes
Step Temp E_pair E_mol TotEng Press Volume
0 0 0 99830.089 99830.089 -2.7258398e+10 6.77032e-05
21 0 0 99822.296 99822.296 -61.021985 5.89032e-05
Loop time of 7.39098e-05 on 1 procs for 21 steps with 4 atoms

47.4% CPU use with 1 MPI tasks x no OpenMP threads

Minimization stats:
Stopping criterion = energy tolerance
Energy initial, next-to-last, final =
99830.0893103 99822.2955809 99822.2955809
Force two-norm initial, final = 6.74872 3.80275e-08
Force max component initial, final = 4.77206 2.78264e-08
Final line search alpha, max atom move = 1 2.78264e-08
Iterations, force evaluations = 21 32

MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total

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. I define
additional bonds using "fix bond/create" (Reference -
LAMMPS Molecular Dynamics Simulator) and "special_bonds extra"
and use "bond_style hybrid" to reproduce the harmonic and l/j interactions
between beads.

​there is no reason to use either fix bond/create or bond_style hybrid in
your simulation. the e-mail you quote is irrelevant to your system.

you already have a suitable bond topology, so there is no need to create
*additional* bonds.
also, with your special_bonds setting of 1.0 1.0 1.0, you do​ include the
1-2 LJ interactions as specified by your pair style, so there is no need to
specify them again with the bond style.
furthermore, bond style hybrid does an "either or" setup and not "one plus
the other", thus your tabulated potential will wipe out the harmonic
interaction, which is definitely not what you want for a spring-bead
polymer.

However when I am running NVT, I see that the temperature and the Ebond
energy is not rising from 0. The visualized trajectory shows a chain not
moving at all.

​that is because your minimization is not working as expected due to the
symmetry of your system. you need to carefully evaluate your bond topology,
it may not be what you expect. you can see what happens by replacing fix
nvt by fix nve. after a few steps, the symmetry is broken due to the
assign​ment of random velocities and then all the stored potential energy
is quickly released, which would overload the nose-hoover thermostat.

the remedy for this is simple. insert a smaller randomization of positions
*before* the minimization, e.g. with:

displace_atoms all random 0.01 0.01 0.01 443322

​this will then break the symmetry and minimize can relax your system to a
local potential energy minimum and after that, the nose-hoover thermostat
can handle your system well.

axel.​