gran/hooke/history and harmonic bonds

Hi there,

I am willing to simulate a bunch of identical spheres with gran/hooke/history potential and harmonic bonds.

On the documentation I found

It is permissible to have itype = jtype. Rmin must be <= the pairwise cutoff distance between itype and jtype atoms, as defined by the pair_style command.

I put itype=jtype=1 [the only type of atoms I have]; as far as I understand, since the cutoff in this case is equal to 1 diameter, I cannot bond spheres with a length-at-rest greater than 1 diameter.

Am I wrong?

Actually I would like to bind spheres even though they are not interacting (in other words, their relative distance is greater than 1 diameter). Is there perhaps a way to override this?

It seems to me that the only way to have what I am willing to is to introduct a longer-ranged potential.

Please find attached the example if you would like to have a look.

Relevant rows (for the bond-part; diameter=1 -> Rmin=0.9):

create_box 4 reg bond/types 1 extra/bond/per/atom 6
bond_style harmonic
bond_coeff * 0.0 0.0
fix 11 all bond/create 1000 1 1 0.9 1

Many thanks

Happy holidays to the ones who are celebrating them!

Mario

in-test (1.63 KB)

Hi there,

I am willing to simulate a bunch of identical spheres with
gran/hooke/history potential and harmonic bonds.

On the documentation I found

It is permissible to have itype = jtype. Rmin must be <= the pairwise
cutoff distance between itype and jtype atoms, as defined by the pair_style
command.

I put itype=jtype=1 [the only type of atoms I have]; as far as I
understand, since the cutoff in this case is equal to 1 diameter, I cannot
bond spheres with a length-at-rest greater than 1 diameter.

Am I wrong?

Actually I would like to bind spheres even though they are not interacting
(in other words, their relative distance is greater than 1 diameter). Is
there perhaps a way to override this?

It seems to me that the only way to have what I am willing to is to
introduct a longer-ranged potential.

Please find attached the example if you would like to have a look.

​there are several issues with your approach.

1) as you already noticed, the bond creation distance cannot be larger than
the neighbor list cutoff of the pair style, since the fix uses the pairwise
neighborlist to speed up the search for possible candidates. this can
indeed be remedied by introducing a pairwise (dummy) interaction with a
longer cutoff, e.g. through:

pair_style hybrid/overlay gran/hooke/history \{Kn\} {Kt} \{Gn\} {Gt}
${mu} 1 lj/cut 2.5
pair_coeff * * gran/hooke/history
pair_coeff * * lj/cut 0.0 1.0

2) there is no special_bonds command, which means that the default settings
apply. those will result in pairwise interactions being removed for atoms
that are bonded. this is clearly incompatible with your setup where there
is a force constant for the bonds of 0.0. for each bond created, you would
have a sudden change in the interactions. either you have to adjust the
force constant for the bonds to be similar to the non-bonded interaction at
the time of creating the bond, or you have to use a different special_bonds
setting, i.e. lj/coul 1.0 1.0 1.0

3) similarly, you are not reserving space for those special neighbors and
thus will get a crash as soon as bonds are created. depending on the amount
of bonding, you may need between 20 and 100 special neighbors per atom.

4) creating the bonds will the atoms are moving and you are compressing the
system to the state you are interested in, will not likely result in a
meaningful behavior. atoms can move around without being affected by the
created bonds and will thus "collect" more and more bonds over time. it
also won't work in parallel, as bonds may become too long for the
communication cutoff and you'll "lose" bonds over time, as atoms diffuse
around.

​5) it is difficult to see any physically meaningful simulation to come out
of what​ you seem to be doing. in general, it'll be much better to separate
the step of system preparation from the actual simulation. if i understand
your intentions correctly, you should do the following steps:
  a) run the simulation without fix bond/create to generate the
"compressed" system and write out a data file, so you don't have to repeat
that step again and again.
  b) do a short run without time integration and a changed potential, e.g.
lj/cut, that has a long enough interaction cutoff, so you can assign the
bonds. then restore the desired pair style and settings and, again, write
out a data file.
  c) set up your simulation with the bonds including proper parameter
settings and re-equilibrate. store the result as a data file suitable for
production runs.
  d) do your production simulation.
  trying to do all those steps within the same input is going to be more
tedious, error-prone and inefficient (e.g. running with (many) more
neighbors required for the bond creation will significantly slow down the
non-bonded interaction calculation, as you then have many neighbors that
are irrelevant, but you need to compute their distance before you can
know it.)

​axel.​

I hereby do thank you for the inspiring and very detailed answer, Axel.

Warm regards and happy new year,

Mario