Unexpected behaviour during simulation time

Hello!

I have recently started using the LAMMPS package to run some MD simulations on a few simple systems, but I am running into issues during simulation time. Essentially, I am setting up a system with a few defined bonds, then running a simulation with Brownian motion using the fix Brownian command, along with the verlet run_style. I would also like to create bonds in the system, if two particles come within a certain distance from each other. I have tried doing this through both manipulating the pair_coeff as well as using fix create/bonds command, however I am getting behaviour that I do not expect from both instances. Namely, certain vertices will approach within the required distance, are of the right types for the bond, have sufficient space for additional bonds to be created, but the bond is not made. The commands seem to be working in other instances, as some bonds are being made, they just seem to be passing up opportunities in some instances and not in others (I have also checked the probability of making bonds which I believe by default is 1). This leads me to believe that I am overlooking something in the way bonds are defined in the system, however I am unfortunately not able to find out where the issue lies.

Below is a snippet of the code.

lmp.command('units real')
lmp.command('timestep 1')
lmp.command('dimension 3')
lmp.command('atom_style full')
lmp.command('pair_style gauss 10')
lmp.command('bond_style harmonic 4.0')
lmp.command('angle_style harmonic')
lmp.command('boundary p p p')

lmp.command('variable x equal 12')
lmp.command('region box block -${x} ${x} -${x} ${x} -${x} ${x}')
lmp.command('create_box 3 box bond/types 1 angle/types 1 extra/bond/per/atom 3 extra/angle/per/atom 5 extra/special/per/atom 10')

lmp.command('create_atoms 1 single -2 -2 0')
lmp.command('create_atoms 1 single 2 -2 0')
lmp.command('create_atoms 1 single 2 2 0')
lmp.command('create_atoms 1 single -2 2 0')
lmp.command('group g1 type 1')

lmp.command('create_atoms 2 single 6 -2 0')
lmp.command('create_atoms 2 single 6 2 0')
lmp.command('group g2 type 2 2')

lmp.command('create_atoms 3 single 2 6 0')
lmp.command('create_atoms 3 single -2 6 0')
lmp.command('group g3 type 2 3')

lmp.command('mass * 10.0')

lmp.command('pair_coeff * * 0 0.2 5')
lmp.command('pair_coeff 2 3 200 0.2 5')

lmp.command('bond_coeff 1 400 4.0')
lmp.command('angle_coeff 1 800 90')

lmp.command('group g8 type 1:3')

lmp.command('create_bonds many g8 g8 1 3.5 4.0')

lmp.command('create_bonds single/angle 1 1 2 3')
lmp.command('create_bonds single/angle 1 1 4 3')
lmp.command('create_bonds single/angle 1 4 1 2')

lmp.command('create_bonds single/angle 1 3 2 5')
lmp.command('create_bonds single/angle 1 3 6 5')
lmp.command('create_bonds single/angle 1 2 3 6')

lmp.command('create_bonds single/angle 1 4 8 7')
lmp.command('create_bonds single/angle 1 4 3 7')
lmp.command('create_bonds single/angle 1 3 4 8')

lmp.command('fix f1 all brownian 50 1 gamma_t 1')

lmp.command('run_style verlet')

lmp.command('dump m0 all movie 1000 ex1.avi type type size 640 480')
lmp.command('run 6000000')

Any help is much appreciated. Thank you.

You cannot create explicit bonds with a pair_coeff command and there is no fix bond/create in your quote (BTW, why are you entering command via python individually and not just a plain input file? You are not using any python features).

Hi! Thank you for your reply.

I have also tried using the fix bond/create command, that I defined in the following way:

lmp.command(‘fix f7 all bond/create 1 2 3 0.5 2’)

(In this case I also defined another harmonic bond type with slightly different parameters to the one I included in the code snippet.) However, this doesn’t change the behaviour. I have also tried simplifying the system to the point where I am only considering 3 particles connected in a line and tried forming a bond between the two side atoms and I am still unable to form the bond between them. Hence why I believe I am making a mistake in my understanding of the different commands or perhaps the parameters I set for the system.

As for you other question, I found it slightly easier to get going with the way I have set it up now, but plan on moving to a text file later on.

Thank you

You may want to compile and use the LAMMPS shell program instead for interactive testing.

For a proper testing and discussion, however, you need to provide a complete input deck in a regular text file, the matching log file, a more specific description of what you expect to get and what get instead and how you determine that.

Hi!

Thank you, I will look into the shell program you mentioned.
Unfortunately I am not able to upload the files, but I have included them here:
Input file:

atom_style full
pair_style gauss 10
bond_style harmonic 4.0
angle_style harmonic
boundary p p p

variable x equal 12
region box block -${x} ${x} -${x} ${x} -${x} ${x}
create_box 3 box bond/types 2 angle/types 1 extra/bond/per/atom 3 extra/angle/per/atom 5 extra/special/per/atom 10

create_atoms 1 single -2 -2 0
create_atoms 1 single 2 -2 0
create_atoms 1 single 2 2 0
create_atoms 1 single -2 2 0
group g1 type 1

create_atoms 2 single 6 -2 0
create_atoms 2 single 6 2 0
group g2 type 2 2

create_atoms 3 single 2 6 0
create_atoms 3 single -2 6 0
group g3 type 2 3

mass * 10.0

pair_coeff * * 0 0.2 5

bond_coeff 1 400 4.0
bond_coeff 2 400 4.0
angle_coeff 1 800 90

group g8 type 1:3

create_bonds many g8 g8 1 3.5 4.0

create_bonds single/angle 1 1 2 3
create_bonds single/angle 1 1 4 3
create_bonds single/angle 1 4 1 2

create_bonds single/angle 1 3 2 5
create_bonds single/angle 1 3 6 5
create_bonds single/angle 1 2 3 6

create_bonds single/angle 1 4 8 7
create_bonds single/angle 1 4 3 7
create_bonds single/angle 1 3 4 8

fix f1 all brownian 50 1 gamma_t 1

fix f2 all bond/create 1 2 3 0.5 2

run_style verlet

dump m0 all movie 1000 ex1.avi type type size 640 480
run 6000000

Corresponding log:

LAMMPS (7 Jan 2022)
OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (src/comm.cpp:98)
  using 1 OpenMP thread(s) per MPI task
units real
timestep 1
dimension 3
atom_style full
pair_style gauss 10
bond_style harmonic 4.0
angle_style harmonic
boundary p p p
variable x equal 12
region box block -${x} ${x} -${x} ${x} -${x} ${x}
region box block -12 ${x} -${x} ${x} -${x} ${x}
region box block -12 12 -${x} ${x} -${x} ${x}
region box block -12 12 -12 ${x} -${x} ${x}
region box block -12 12 -12 12 -${x} ${x}
region box block -12 12 -12 12 -12 ${x}
region box block -12 12 -12 12 -12 12
create_box 3 box bond/types 2 angle/types 1 extra/bond/per/atom 3 extra/angle/per/atom 5 extra/special/per/atom 10
Created orthogonal box = (-12 -12 -12) to (12 12 12)
  1 by 1 by 1 MPI processor grid
create_atoms 1 single -2 -2 0
Created 1 atoms
  using lattice units in orthogonal box = (-12 -12 -12) to (12 12 12)
  create_atoms CPU = 0.001 seconds
create_atoms 1 single 2 -2 0
Created 1 atoms
  using lattice units in orthogonal box = (-12 -12 -12) to (12 12 12)
  create_atoms CPU = 0.001 seconds
create_atoms 1 single 2 2 0
Created 1 atoms
  using lattice units in orthogonal box = (-12 -12 -12) to (12 12 12)
  create_atoms CPU = 0.001 seconds
create_atoms 1 single -2 2 0
Created 1 atoms
  using lattice units in orthogonal box = (-12 -12 -12) to (12 12 12)
  create_atoms CPU = 0.001 seconds
group g1 type 1
4 atoms in group g1
create_atoms 2 single 6 -2 0
Created 1 atoms
  using lattice units in orthogonal box = (-12 -12 -12) to (12 12 12)
  create_atoms CPU = 0.001 seconds
create_atoms 2 single 6 2 0
Created 1 atoms
  using lattice units in orthogonal box = (-12 -12 -12) to (12 12 12)
  create_atoms CPU = 0.001 seconds
group g2 type 2 2
2 atoms in group g2
create_atoms 3 single 2 6 0
Created 1 atoms
  using lattice units in orthogonal box = (-12 -12 -12) to (12 12 12)
  create_atoms CPU = 0.001 seconds
create_atoms 3 single -2 6 0
Created 1 atoms
  using lattice units in orthogonal box = (-12 -12 -12) to (12 12 12)
  create_atoms CPU = 0.001 seconds
group g3 type 2 3
4 atoms in group g3
mass * 10.0
pair_coeff * * 0 0.2 5
bond_coeff 1 400 4.0
bond_coeff 2 400 4.0
angle_coeff 1 800 90
group g8 type 1:3
8 atoms in group g8
create_bonds many g8 g8 1 3.5 4.0
  generated 0 of 3 mixed pair_coeff terms from geometric mixing rule
Neighbor list info ...
  update every 1 steps, delay 10 steps, check yes
  max neighbors/atom: 2000, page size: 100000
  master list distance cutoff = 7
  ghost atom cutoff = 7
  binsize = 3.5, bins = 7 7 7
  2 neighbor lists, perpetual/occasional/extra = 1 1 0
  (1) command create_bonds, occasional
      attributes: full, newton on
      pair build: full/bin
      stencil: full/bin/3d
      bin: standard
  (2) pair gauss, perpetual
      attributes: half, newton on
      pair build: half/bin/newton
      stencil: half/bin/3d
      bin: standard
WARNING: Communication cutoff 7 is shorter than a bond length based estimate of 8. This may lead to errors. (src/comm.cpp:723)
Added 10 bonds, new total = 10
Finding 1-2 1-3 1-4 neighbors ...
  special bond factors lj:    0        0        0       
  special bond factors coul:  0        0        0       
     4 = max # of 1-2 neighbors
     6 = max # of 1-3 neighbors
    12 = max # of 1-4 neighbors
    21 = max # of special neighbors
  special bonds CPU = 0.001 seconds
create_bonds single/angle 1 1 2 3
Finding 1-2 1-3 1-4 neighbors ...
  special bond factors lj:    0        0        0       
  special bond factors coul:  0        0        0       
     4 = max # of 1-2 neighbors
     6 = max # of 1-3 neighbors
    12 = max # of 1-4 neighbors
    21 = max # of special neighbors
  special bonds CPU = 0.001 seconds
create_bonds single/angle 1 1 4 3
Finding 1-2 1-3 1-4 neighbors ...
  special bond factors lj:    0        0        0       
  special bond factors coul:  0        0        0       
     4 = max # of 1-2 neighbors
     6 = max # of 1-3 neighbors
    12 = max # of 1-4 neighbors
    21 = max # of special neighbors
  special bonds CPU = 0.001 seconds
create_bonds single/angle 1 4 1 2
Finding 1-2 1-3 1-4 neighbors ...
  special bond factors lj:    0        0        0       
  special bond factors coul:  0        0        0       
     4 = max # of 1-2 neighbors
     6 = max # of 1-3 neighbors
    12 = max # of 1-4 neighbors
    21 = max # of special neighbors
  special bonds CPU = 0.001 seconds
create_bonds single/angle 1 3 2 5
Finding 1-2 1-3 1-4 neighbors ...
  special bond factors lj:    0        0        0       
  special bond factors coul:  0        0        0       
     4 = max # of 1-2 neighbors
     6 = max # of 1-3 neighbors
    12 = max # of 1-4 neighbors
    21 = max # of special neighbors
  special bonds CPU = 0.001 seconds
create_bonds single/angle 1 3 6 5
Finding 1-2 1-3 1-4 neighbors ...
  special bond factors lj:    0        0        0       
  special bond factors coul:  0        0        0       
     4 = max # of 1-2 neighbors
     6 = max # of 1-3 neighbors
    12 = max # of 1-4 neighbors
    21 = max # of special neighbors
  special bonds CPU = 0.001 seconds
create_bonds single/angle 1 2 3 6
Finding 1-2 1-3 1-4 neighbors ...
  special bond factors lj:    0        0        0       
  special bond factors coul:  0        0        0       
     4 = max # of 1-2 neighbors
     6 = max # of 1-3 neighbors
    12 = max # of 1-4 neighbors
    21 = max # of special neighbors
  special bonds CPU = 0.005 seconds
create_bonds single/angle 1 4 8 7
Finding 1-2 1-3 1-4 neighbors ...
  special bond factors lj:    0        0        0       
  special bond factors coul:  0        0        0       
     4 = max # of 1-2 neighbors
     6 = max # of 1-3 neighbors
    12 = max # of 1-4 neighbors
    21 = max # of special neighbors
  special bonds CPU = 0.001 seconds
create_bonds single/angle 1 4 3 7
Finding 1-2 1-3 1-4 neighbors ...
  special bond factors lj:    0        0        0       
  special bond factors coul:  0        0        0       
     4 = max # of 1-2 neighbors
     6 = max # of 1-3 neighbors
    12 = max # of 1-4 neighbors
    21 = max # of special neighbors
  special bonds CPU = 0.001 seconds
create_bonds single/angle 1 3 4 8
Finding 1-2 1-3 1-4 neighbors ...
  special bond factors lj:    0        0        0       
  special bond factors coul:  0        0        0       
     4 = max # of 1-2 neighbors
     6 = max # of 1-3 neighbors
    12 = max # of 1-4 neighbors
    21 = max # of special neighbors
  special bonds CPU = 0.001 seconds
fix f1 all brownian 50 1 gamma_t 1
fix f2 all bond/create 1 2 3 0.5 2
run_style verlet
dump m0 all movie 1000 ex1.avi type type size 640 480
run 6000000
  generated 0 of 3 mixed pair_coeff terms from geometric mixing rule
Neighbor list info ...
  update every 1 steps, delay 10 steps, check yes
  max neighbors/atom: 2000, page size: 100000
  master list distance cutoff = 7
  ghost atom cutoff = 7
  binsize = 3.5, bins = 7 7 7
  2 neighbor lists, perpetual/occasional/extra = 1 1 0
  (1) pair gauss, perpetual
      attributes: half, newton on
      pair build: half/bin/newton
      stencil: half/bin/3d
      bin: standard
  (2) fix bond/create, occasional, copy from (1)
      attributes: half, newton on
      pair build: copy
      stencil: none
      bin: none
WARNING: Communication cutoff 7 is shorter than a bond length based estimate of 8. This may lead to errors. (src/comm.cpp:723)
Per MPI rank memory allocation (min/avg/max) = 9.709 | 9.709 | 9.709 Mbytes
Step Temp E_pair E_mol TotEng Press 
       0            0            0            0            0            0 
 6000000     1585.522            0   0.91514354    33.998122    109.77392 
Loop time of 59.3543 on 1 procs for 6000000 steps with 8 atoms

Performance: 8733.998 ns/day, 0.003 hours/ns, 101087.939 timesteps/s
53.9% CPU use with 1 MPI tasks x 1 OpenMP threads

MPI task timing breakdown:
Section |  min time  |  avg time  |  max time  |%varavg| %total
---------------------------------------------------------------
Pair    | 0.64316    | 0.64316    | 0.64316    |   0.0 |  1.08
Bond    | 4.6042     | 4.6042     | 4.6042     |   0.0 |  7.76
Neigh   | 0.013939   | 0.013939   | 0.013939   |   0.0 |  0.02
Comm    | 1.2484     | 1.2484     | 1.2484     |   0.0 |  2.10
Output  | 45.248     | 45.248     | 45.248     |   0.0 | 76.23
Modify  | 6.0432     | 6.0432     | 6.0432     |   0.0 | 10.18
Other   |            | 1.553      |            |       |  2.62

Nlocal:              8 ave           8 max           8 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost:             10 ave          10 max          10 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs:              0 ave           0 max           0 min
Histogram: 1 0 0 0 0 0 0 0 0 0

Total # of neighbors = 0
Ave neighs/atom = 0
Ave special neighs/atom = 7
Neighbor list builds = 1565
Dangerous builds = 1
Total wall time: 0:00:59

The system here is composed of 3 square faces, one base and connected to it, two faces on adjacent edges. I am trying to get the system to connect the two adjacent faces along the vertices, provided they are close enough (i.e. the faces make an angle of ~ 90degrees with the base at the same time). The faces should bind in this way. However, it appears that the vertices closest to each other, are unable to bind at all (even though they are types which should be allowed to bind). Instead in simulation time, the faces move until the far vertices are close enough and then bind. This isn’t necessarily undesired, but looking through simulation videos, I can see the close vertices approaching each other, but not binding, which is where the unexpected behaviour comes from. I am thinking that there is something else preventing them from binding, but I am not sure where I am going wrong.

Apologies if this is not the format you had asked for.

Hi Axel!

I just wanted to post an update, as I have had another look through the documentation and managed to track down the special_bonds command. I believe that the default value was causing the problem. Including the following command;

lmp.command(‘special_bonds coul 0 1 0’)

solves my issue.

Thank you for your help

While this may solve the issue of forming the bonds, you need to be aware that this changes how non-bonded interactions are computed. With the default of special_bonds lj/coul 0 0 0 all non-bonded interactions between pairs that are connected via a chain of bonds a 1-2, 1-3, or 1-4 neighbors are excluded from the non-bonded pair interactions. That is the usual expected behavior for molecular force fields as then the interactions for those pairs are entirely determined by the corresponding bonded interactions and thus straightforward to parameterize. If your system does not have dihedral interactions, they the setting would usually be special_bonds lj/coul 0 0 1 and some force fields also mix in the 1-4 interactions to a fractional degree.

Thus, by changing your special_bonds setting you are now adding non-bonded forces to your model that are potentially unwanted. You can get the same effect that you are looking for with exclusion of the interactions by using special_bonds lj/coul 0.0 1.0e-20 1.0. Additional explanations are in the manual.

That’s great, thank you for the suggestion. At the moment it seems that the system is behaving closer to what I would hope. But I will experiment a little with this command to get a better idea of it and have a close read of the documentation as well.

Thanks for all your help today