Polymer simulation breaks

Hi,

I am simulating an LJ polymer system, a single chain with 500 monomers in a confined box.
Boundary conditions are fixed. I fix the temperature at 1.0, yet the system starts from a very low
temperature (0.0002) and the temperature keeps going up steadily, it does not stay fix at 1, it just keeps going up.
I wonder what I should do to keep the temperature fix. I am using NVT. The following is the sample of
my input file.

The other issue is that sometimes the simulation breaks in the middle, beofore reaching the
desired number of steps (I use very small time steps for this reason) ,
and I get the error message of (example): “ERROR on proc 0: Bond atoms 294 295 missing on proc 0 at step 102658”

I’ll appreciate any help in these regards.
Than you.

Input file:

units lj
boundary f f f

atom_style molecular
bond_style harmonic

angle_style none
pair_style lj/cut 1.12246
pair_modify shift yes

read_data dna_n500_rect_p0.02

bond_coeff 1 1.5 1.0

pair_coeff * * 1.0 1.0

log log.lj500

timestep 0.000014

thermo 100
thermo_style custom step temp etotal ke pe epair press vol
thermo_modify flush yes
thermo_modify norm yes

fix 1 all nvt temp 1.0 1.0 0.5

dump 1 all atom 200000 lj500.dump

run 200000

hi jane,

Hi,

I am simulating an LJ polymer system, a single chain with 500 monomers in a
confined box.

i don't see any confinement.

Boundary conditions are fixed. I fix the temperature at 1.0, yet the system
starts from a very low
temperature (0.0002) and the temperature keeps going up steadily, it does
not stay fix at 1, it just keeps going up.

the kinetic energy (i.e. temperature) will be at 0.0 unless
you initialize it. using the nvt integrator will gradually ramp
it up until it reaches the target.

from then on, energy conservation depends on how
accurate the integration of your equations of motion
with finite differences is. this depends on the highest
characteristic frequencies in your system. the time
step has to be chosen small enough to represent
that well.

I wonder what I should do to keep the temperature fix. I am using NVT. The
following is the sample of
my input file.

you should first study what is needed to have reasonable
energy conservation. i would suggest to initialize
velocities using the velocity command and then start
with fix nve instead of fix nvt and see how well your
system conserves the total(!) energy as a function of
the length of time step. then you can move to fix nvt
with those time step setting. the "strength" of the
thermalization of fix nvt is constrolled by the time constant
(0.5 in your case), reducing it should accelerate the
thermal exchange, but also may make your calculation
less realistic as this interferes with the trajectories.
depending on what your system is representing, you
may also be better off using fix nve + fix langevin.

The other issue is that sometimes the simulation breaks in the middle,
beofore reaching the
desired number of steps (I use very small time steps for this reason) ,
and I get the error message of (example): "ERROR on proc 0: Bond atoms 294
295 missing on proc 0 at step 102658"

since you have fixed i.e. non-periodic boundaries,
atoms may get "lost". if this happens for an atom
that has a bond, you will get this error.

your system looks a bit odd. you have repulsive
only non-bonded interactions, so there is nothing
that will keep them together, yet you don't employ
periodic boundaries, so you should always end up
with "lost" atoms. it only is a matter of time.

cheers,
    axel.

Thanks for your email and explanations. I did my simulation with fix wall and used nve style, the total energy is conserved with the time step that I use, but after some time the same error happens and the simulation breaks again, If I am using reflecting
walls, why are the particles going out of the box again? Thank you.

Here is a part of the log file and then the input file:

LAMMPS (11 Aug 2010)
Scanning data file …
1 = max bonds/atom
Reading data file …
orthogonal box = (-10.77 -10.77 -10.77) to (10.77 10.77 10.77)
1 by 1 by 1 processor grid
200 atoms
199 bonds
Finding 1-2 1-3 1-4 neighbors …
2 = max # of 1-2 neighbors
2 = max # of 1-3 neighbors
4 = max # of 1-4 neighbors
6 = max # of special neighbors
Setting up run …
Memory usage per processor = 2.11097 Mbytes
Step Temp TotEng KinEng PotEng E_pair Press Volume
0 0.2 0.6485 0.2985 0.35 0.35 0.06001632 9993.9483
100 0.20073464 0.6485 0.29959645 0.34890355 0.34890237 0.059899487 9993.9483
200 0.20188633 0.6485 0.30131535 0.34718465 0.34717995 0.059713693 9993.9483
300 0.2034477 0.6485 0.3036457 0.3448543 0.34484374 0.059459919 9993.9483
400 0.20540926 0.6485 0.30657332 0.34192668 0.34190791 0.059139435 9993.9483
500 0.20775951 0.6485 0.31008107 0.33841893 0.33838962 0.058753785 9993.9483
600 0.21048504 0.6485 0.31414893 0.33435107 0.33430888 0.058304774 9993.9483
.
.
.
279100 0.31543995 0.6485 0.47079413 0.17770587 0.00085911196 0.00061173048 9993.9483
279200 0.31546967 0.6485 0.47083848 0.17766152 0.00082795155 0.00060259702 9993.9483
279300 0.31549915 0.6485 0.47088248 0.17761751 0.00079713413 0.00059340655 9993.9483
279400 0.31552838 0.6485 0.4709261 0.1775739 0.00076670282 0.0005841654 9993.9483
279500 0.31555731 0.6485 0.47096928 0.17753072 0.00073669957 0.00057487978 9993.9483
279600 0.31558591 0.6485 0.47101197 0.17748802 0.00070716502 0.00056555581 9993.9483
279700 0.31561417 0.6485 0.47105414 0.17744585 0.00067813845 0.00055619946 9993.9483
279800 0.31564204 0.6485 0.47109574 0.17740426 0.00064965769 0.00054681659 9993.9483
279900 0.3156695 0.6485 0.47113673 0.17736327 0.00062175907 0.00053741291 9993.9483
280000 0.31569653 0.6485 0.47117706 0.17732293 0.00059447734 0.00052799399 9993.9483
ERROR on proc 0: Bond atoms 199 200 missing on proc 0 at step 280062

INPUT FILE:
units lj

boundary f f f

atom_style molecular
bond_style harmonic

angle_style none
pair_style lj/cut 1.12246
pair_modify shift yes

read_data dna_n200_rect_p0.02

fix walls all wall/reflect xlo -10.77 xhi 10.77 units box
fix walls all wall/reflect ylo -10.77 yhi 10.77 units box
fix walls all wall/reflect zlo -10.77 zhi 10.77 units box

bond_coeff 1 1.5 1.0

pair_coeff * * 1.0 1.0

log in.lj.log-200

timestep 0.000014

thermo 100
thermo_style custom step temp etotal ke pe epair press vol
thermo_modify flush yes
thermo_modify norm yes

velocity all create 0.2 49284

fix 1 all nve

dump 1 all atom 500000 in.lj3.dump

run 500000

Thanks for your email and explanations. I did my simulation with fix wall
and used nve style, the total energy is conserved with the time step that I

that is good news.

use, but after some time the same error happens and the simulation breaks
again, If I am using reflecting
walls, why are the particles going out of the box again? Thank you.

there are multiple reasons why you may get the "bond atom missing"
error. an atom leaving the box was just the most likely explanation
based on the information about your simulation that you provided.

a second possibility is that atoms move too fast for the neighborlist
rebuild. you can tweak the neighbor list settings with the neigh_modify
command. the documentation of this command gives some hints at
how to tweak the settings. you want to avoid "dangerous builds" but
at the same time avoid too frequent neighbor list builds. the "skin"
parameter also has an impact, as it decides how many additional
particles are considered beyond the cutoff. since you use a very
short cutoff, you may want to increase the skin a bit. or reduce
delay or every. or ???

axel.

Thanks a lot for all your comments. I did modify the neighbor list and got the simulation to run longer but it still breaks
after some time. But the other thing is, I am trying to implement the WCA potential (purely repulsive part of LJ) b/w the monomers in the polymer system in order to have a self avoiding walk polymer and so that’s why I had the LJ cut off
at 1.12246 (the following commands)
pair_style lj/cut 1.12246
pair_modify shift yes

but if I use the following:

pair_style lj/cut 2.5
pair_modify shift yes
pair_coeff * * 1.0 1.0 1.122

The simulation will run fine and does not break. I am not sure if the above setting is doing what
I want (the atoms having purely repulsive interactions). I’ll appreciate your help regarding this.

Thank you.

Thanks a lot for all your comments. I did modify the neighbor list and got
the simulation to run longer but it still breaks
after some time. But the other thing is, I am trying to implement the WCA
potential (purely repulsive part of LJ) b/w the monomers in the polymer
system in order to have a self avoiding walk polymer and so that's why I had
the LJ cut off
at 1.12246 (the following commands)

pair_style lj/cut 1.12246
pair_modify shift yes

but if I use the following:

pair_style lj/cut 2.5
pair_modify shift yes
pair_coeff * * 1.0 1.0 1.122

The simulation will run fine and does not break. I am not sure if the above
setting is doing what
I want (the atoms having purely repulsive interactions). I'll appreciate
your help regarding this.

if you place the cutoff at the pair coeff definition or
use a global cutoff has no difference in how the
forces and potential energies are computed. you
should get the same energies, if you use exactly
the same values. you can try both settings and
compare.

however, having a global cutoff of 2.5 has about
the same effect as increasing the "skin" distance
for the neighbor lists to about 1.3. usually the
neighbor lists are built based on the largest of
all defined cutoffs plus the skin distance.

cheers,
     axel.

Did you run a soft potential before starting the LJ interactions? That
way, you can ensure there are no overlaps between beads.

Also, I would suggest using a soft wall (wall/lj96) rather than wall/reflect.

-Monojoy

Thanks for the suggestions. I changed to lj/wall93 and used a soft potential twice. Yet the simulation still breaks,
Could you please take a look at my input file and see if there is anything you find unreasonable. I appreciate your feedback.

units lj
boundary f f f

atom_style molecular

neighbor 1.3 bin

angle_style none

read_data dna_n200_rect_p0.02

fix walls all wall/lj93 xlo -10.8 1.0 1.0 2.5 xhi 10.8 1.0 1.0 2.5 units box
fix walls all wall/lj93 ylo -10.8 1.0 1.0 2.5 yhi 10.8 1.0 1.0 2.5 units box
fix walls all wall/lj93 zlo -10.8 1.0 1.0 2.5 zhi 10.8 1.0 1.0 2.5 units box

pair_style soft 1.12246
pair_coeff * * 30.0 1.12246

bond_style harmonic
bond_coeff 1 1.5 1.0

pair_style soft 1.12246
pair_coeff * * 30.0 1.12246

pair_style lj/cut 2.5
pair_modify shift yes
pair_coeff * * 1.0 1.0 1.12

log lg.n200

timestep 0.000014

thermo 100
thermo_style custom step temp etotal ke pe epair press vol
thermo_modify flush yes
thermo_modify norm yes

velocity all create 0.05 49284

fix 1 all nve
#fix 2 all langevin 1.0 1.1 100.0 48279

dump 1 all atom 1000000 lj.n200.dump

run 1000000

Thanks for the suggestions. I changed to lj/wall93 and used a soft potential
twice. Yet the simulation still breaks,
Could you please take a look at my input file and see if there is anything
you find unreasonable. I appreciate your feedback.

units lj
boundary f f f

atom_style molecular

neighbor 1.3 bin

angle_style none

read_data dna_n200_rect_p0.02

fix walls all wall/lj93 xlo -10.8 1.0 1.0 2.5 xhi 10.8 1.0 1.0 2.5 units
box
fix walls all wall/lj93 ylo -10.8 1.0 1.0 2.5 yhi 10.8 1.0 1.0 2.5 units
box
fix walls all wall/lj93 zlo -10.8 1.0 1.0 2.5 zhi 10.8 1.0 1.0 2.5 units box

pair_style soft 1.12246
pair_coeff * * 30.0 1.12246

bond_style harmonic
bond_coeff 1 1.5 1.0

pair_style soft 1.12246
pair_coeff * * 30.0 1.12246

this is pointless. if you specity
pair_style a second time it will
override the first setting. the same
is true for pair_coeff.

pair_style lj/cut 2.5
pair_modify shift yes
pair_coeff * * 1.0 1.0 1.12

and again. these lines will remove
everything that was set for pair_style soft

axel.

Thanks for all the suggestions. I managed to get the simulation run, but for a lower temperature,
when I set the temperature to 0.3 (using NVT), the simulation does not break. But I wonder if I could run this
system with a higher temperature and have the system equilibrate: get constant temperature. I very much appreciate
all the help.

units lj
boundary f f f

atom_style molecular

neighbor 1.5 bin

angle_style none

read_data dna_n200_rect_p0.02

fix walls all wall/lj93 xlo -10.8 1.0 1.0 2.5 xhi 10.8 1.0 1.0 2.5 units box
fix walls all wall/lj93 ylo -10.8 1.0 1.0 2.5 yhi 10.8 1.0 1.0 2.5 units box
fix walls all wall/lj93 zlo -10.8 1.0 1.0 2.5 zhi 10.8 1.0 1.0 2.5 units box

pair_style soft 1.12246
pair_coeff * * 30.0 1.12246

bond_style fene
bond_coeff 1 30.0 1.5 1.0 1.0

run 50000

pair_style lj/cut 2.5#1.12246
pair_modify shift yes
pair_coeff * * 1.0 1.0 1.122

log log.n200

timestep 0.000014

thermo 100
thermo_style custom step temp etotal ke pe epair press vol
thermo_modify flush yes
thermo_modify norm yes

velocity all create 0.05 4928

fix 1 all nve
fix 2 all nvt temp 0.3 0.3 0.61

dump 1 all atom 1000000 lj.n200.dump

run 1000000

Thanks for all the suggestions. I managed to get the simulation run,
but for a lower temperature,
when I set the temperature to 0.3 (using NVT), the simulation does not
break. But I wonder if I could run this
system with a higher temperature and have the system equilibrate: get
constant temperature. I very much appreciate
all the help.

units lj
boundary f f f

atom_style molecular

neighbor 1.5 bin

angle_style none

read_data dna_n200_rect_p0.02

fix walls all wall/lj93 xlo -10.8 1.0 1.0 2.5 xhi 10.8 1.0 1.0 2.5
units box
fix walls all wall/lj93 ylo -10.8 1.0 1.0 2.5 yhi 10.8 1.0 1.0 2.5
units box
fix walls all wall/lj93 zlo -10.8 1.0 1.0 2.5 zhi 10.8 1.0 1.0 2.5
units box

pair_style soft 1.12246
pair_coeff * * 30.0 1.12246

bond_style fene
bond_coeff 1 30.0 1.5 1.0 1.0

run 50000

this is useless, and LAMMPS should have printed
a warning. you don't have any fix defined that does
time integration for you, i.e. your atoms won't
move and you just compute the same forces over and
over again; 50,000 times.

pair_style lj/cut 2.5#1.12246
pair_modify shift yes
pair_coeff * * 1.0 1.0 1.122

log log.n200

timestep 0.000014

thermo 100
thermo_style custom step temp etotal ke pe epair press vol
thermo_modify flush yes
thermo_modify norm yes

velocity all create 0.05 4928

fix 1 all nve
fix 2 all nvt temp 0.3 0.3 0.61

and this is about as bad. now you integrate
the atoms twice!

again, LAMMPS should have printed a warning
about this. if you keep doing things like this,
you won't ever get a meaningful simulation.

remember, MD can be *very* tricky in that you
can run a simulation for a long time that will
be complete garbage, ... and you can have a
system that has proper settings, but just needs
to be coaxed into the proper equilibrium state.

if you want to learn this by trial an error,
then get used to spending a lot of time and
people getting quickly less than helpful, since
a lot of these mistakes can be avoid through
common sense and spending some time reflecting
about errors and reading suitable text books
on the practical side of MD and statistical
mechanics.

cheers,
    axel.

well, this is not my field but have learned some by now., yet b/c of some comparison in my project
am very willing to get this running if I can. I modified the input file, if possible please let me know
if it still has got major problems. Thank you.

units lj

boundary f f f

atom_style molecular

neighbor 1.5 bin

angle_style none

read_data dna_n200_rect_p0.02

fix walls all wall/lj93 xlo -10.8 1.0 1.0 2.5 xhi 10.8 1.0 1.0 2.5 units box
fix walls all wall/lj93 ylo -10.8 1.0 1.0 2.5 yhi 10.8 1.0 1.0 2.5 units box
fix walls all wall/lj93 zlo -10.8 1.0 1.0 2.5 zhi 10.8 1.0 1.0 2.5 units box

pair_style soft 1.12246
pair_coeff * * 30.0 1.12246

bond_coeff 1 1.5 1.0

bond_style fene
bond_coeff 1 30.0 1.5 1.0 1.0

fix 1 all nve
run 50000

pair_style lj/cut 2.5
pair_modify shift yes

fix 2 all nve
run 20000
dump 1 all atom 20000 in.lj1.dump

fix 3 all langevin 1.0 1.0 100.0 48279
run 10000
dump 2 all atom 50000 in.lj2.dump

fix 4 all nvt temp 1.0 1.0 0.61
run 300000
dump 3 all atom 300000 in.lj3.dump

well, this is not my field but have learned some by now., yet b/c of

well, whether this is your field or not is irrelevant.
if you want to do MD simulations, you have to learn to
do them properly, or you will be wasting you time (and
that of others). computer simulation is brutally stupid.
it does *exactly* what you tell it to do, and that may
not be what you *want* it to do. one has to understand
how the "mind" of the software is working and then how
the method works in order to get good and usable results.
the chance to get out garbage is surprisingly high,
considering how simple most of the models are.

some comparison in my project
am very willing to get this running if I can. I modified the input
file, if possible please let me know
if it still has got major problems. Thank you.

yes, it still has. you probably should study
existing inputs and examples a bit more and
stick your nose into a textbook on MD a bit
deeper. there are still both lammps specific
problems, i.e. you have not fully realized
what individual command in lammps are doing
and how you "assemble" an MD input and there
are choices of parameters that ignore good
practice guidelines that are outlined in most
of the standard text books. sorry to be a
bit harsh on this, but after pointing out
the same type of problem repeatedly, you
are exhausting the credit that i am willing
to give a person with little experience
very quickly.

units lj

boundary f f f

atom_style molecular

neighbor 1.5 bin

angle_style none

read_data dna_n200_rect_p0.02

fix walls all wall/lj93 xlo -10.8 1.0 1.0 2.5 xhi 10.8 1.0 1.0 2.5
units box
fix walls all wall/lj93 ylo -10.8 1.0 1.0 2.5 yhi 10.8 1.0 1.0 2.5
units box
fix walls all wall/lj93 zlo -10.8 1.0 1.0 2.5 zhi 10.8 1.0 1.0 2.5
units box

pair_style soft 1.12246
pair_coeff * * 30.0 1.12246

bond_coeff 1 1.5 1.0

you cannot set bond parameters before you define the bond style.

bond_style fene
bond_coeff 1 30.0 1.5 1.0 1.0

fix 1 all nve

you should initialize the kinetic energy
and have some thermostatting here, e.g.
with fix langevin to get a good
initialization and equipartitioning of
the kinetic energy.

run 50000

it is near impossible to predict for how
many steps this part of the calculation
needs to be done, if don't know your system.
it takes some experience and monitoring
to know what is a good choice. this holds
true for all run statements. so rather than
running multiple steps in one go, you may be
better off to build and understand your
system in stages and use restart files to
not have to re-do everything from the beginning
all the time. BTW: you haven't defined a time
step, so lammps will pick the default. not sure
if that is what you want.

pair_style lj/cut 2.5
pair_modify shift yes

if you define a new pair style,
you *have* to define new pair
interaction coefficients.
the old settings are lost.

each style and fix or dump command
creates an object inside of lammps
and will than hold the associated
information. if you override or delete
the object, that information will be
gone and has to be initialize again.

fix 2 all nve

you are defining a second fix
to time integrate atom positions
without removing the first through
unfix. since you had used fix nve
before, there is no need to define
it again. unless you use the "clear"
command, which will completely wipe
out the system, every "object" that
was defined stays active

run 20000
dump 1 all atom 20000 in.lj1.dump

fix 3 all langevin 1.0 1.0 100.0 48279

you are giving a *very* long time constant here
with 100.0. that will render this fix nearly
inactive. there are suggestions for what is a
good choice in the lammps manual.

run 10000
dump 2 all atom 50000 in.lj2.dump

for dumps the same rule applies than what applies
to fixes. this will create a second trajectory file,
while the fist will stay active and is going to be
written to. you can turn off fixes with "unfix"
and dumps with "undump".

fix 4 all nvt temp 1.0 1.0 0.61

... and now you define yet another time integration
fix *and* still have fix langevin active. ...and
here you use a *much* shorter time constant. in general
you would want to do it exactly the other way around:
use a somewhat shorter time constant for the thermostat
to initially bring the system close to equilibrium,
and then use a weaker thermal coupling, to get the
more realistic trajectory with no interference of
the thermostat.

run 300000
dump 3 all atom 300000 in.lj3.dump

if you define a dump object *after* the last run
command, it will have no effect. if you want to
store the final state, you need to write a restart.

axel.