simulation box exploding

Hello,

I am trying to simulate an ideal filler dispersion within a polymer melt, and I am having some trouble keeping the simulation box from blowing up. In order to create an ideal dispersion, I am using a packing software (packmol) to randomly place the fillers in a box, and then excluding the fillers from any time integration to keep their positions fixed and truly random (the fillers are generated within a small cubic region which is contained within a larger box filled with polymer). My filler particles are basically spherical, and are modeled by 2 outer layers and a hollow center (looks like a shell). Each filler particle contains 50 or so beads that are separated by the LJ equilibrium distance (1.12 sigma). I am intentionally avoiding making each filler a rigid body because I do not want any of their positions to change over time (nor do I want to dilate them in the event of applied pressure). In the generation phase, I apply pressure (fix npt) to only the polymer in order to allow the box to contract to roughly the cubic dimensions containing the filler particles. The polymer/polymer and polymer/filler LJ interactions are attractive. I ramp the pressure from 1.0 to 0.0 over the course of 4000 tau, and as soon as the pressure gets close to zero, the simulation box explodes.

I believe the box is blowing up because of the vacuum regions within each of my filler particles, but I am not positive. I would greatly appreciate any feedback or advice on the root of this problem and potential solutions. Thank you.

If your system always has a positive pressure, then if you ramp the external

pressure to 0.0 during NPT, the box will continue to expand. You can test this by doing

NVE or NVT (not NPT) simulations at various volumes and monitoring the pressure.

Steve

Hello,

I am trying to simulate an ideal filler dispersion within a polymer melt,
and I am having some trouble keeping the simulation box from blowing up. In
order to create an ideal dispersion, I am using a packing software (packmol)
to randomly place the fillers in a box, and then excluding the fillers from
any time integration to keep their positions fixed and truly random (the
fillers are generated within a small cubic region which is contained within
a larger box filled with polymer). My filler particles are basically
spherical, and are modeled by 2 outer layers and a hollow center (looks like
a shell). Each filler particle contains 50 or so beads that are separated by
the LJ equilibrium distance (1.12 sigma). I am intentionally avoiding making
each filler a rigid body because I do not want any of their positions to
change over time (nor do I want to dilate them in the event of applied
pressure). In the generation phase, I apply pressure (fix npt) to only the
polymer in order to allow the box to contract to roughly the cubic
dimensions containing the filler particles. The polymer/polymer and
polymer/filler LJ interactions are attractive. I ramp the pressure from 1.0
to 0.0 over the course of 4000 tau, and as soon as the pressure gets close
to zero, the simulation box explodes.

I believe the box is blowing up because of the vacuum regions within each of
my filler particles, but I am not positive. I would greatly appreciate any
feedback or advice on the root of this problem and potential solutions.

difficult to say without having a complete input deck to look at.
one thing comes to mind: do your filler particles interact with each other?
i.e. do they have a potential defined between them? if yes, have you
tried excluding those interactions from the neighbor list (or just
setting them to zero)?
that might generate spurious contributions to the stress tensor.

axel.

Yes, the beads within each filler have LJ interactions, and I think that is part of the problem–I only want inter-filler interactions, and I want to turn off intra-filler interactions. The only way to do this (that I know of) with my current setup is to change the atom type in every filler, setting epsilon = 0 for self interactions (beads in same filler), and then changing epsilon as needed for inter-filler attractions. However, I have over 400 filler particles in a given system, which will create a laundry list of atom types and pair interactions to specify. Instead, I am thinking of bonding every bead (in a given filler particle) that falls within the LJ cutoff distance (we use 2.5 sigma) in a harmonic potential with a spring constant of 0. Using “special bonds 0 1 1”, this should turn off any pair interactions within 2.5 sigma, and since the spring constant is 0, bonded interactions should not contribute to the pressure (or density).

In my bonding scheme, imagine that atom 1 is 1 sigma from atom 2, and atom 2 is 1 sigma from atom 3. With my 2.5 sigma bonding criterion, there will be a bond between 1-2, 2-3, and 1-3, effectively creating a 3-membered ring. Thus, a given pair of atoms could have a 1-2 or a 1-3 topology, which leads to my question about the special bonds command–do 1-2 interactions take precedence over 1-3 and 1-4 interactions?

Yes, the beads within each filler have LJ interactions, and I think that is
part of the problem--I only want inter-filler interactions, and I want to
turn off intra-filler interactions. The only way to do this (that I know of)
with my current setup is to change the atom type in every filler, setting
epsilon = 0 for self interactions (beads in same filler), and then changing
epsilon as needed for inter-filler attractions. However, I have over 400
filler particles in a given system, which will create a laundry list of atom
types and pair interactions to specify. Instead, I am thinking of bonding
every bead (in a given filler particle) that falls within the LJ cutoff
distance (we use 2.5 sigma) in a harmonic potential with a spring constant
of 0. Using "special bonds 0 1 1", this should turn off any pair
interactions within 2.5 sigma, and since the spring constant is 0, bonded
interactions should not contribute to the pressure (or density).

there is another option: assign molecule ids to each filler and then
use the per molecule exclusion option of neigh_modify. this should be
much easier and cleaner.

In my bonding scheme, imagine that atom 1 is 1 sigma from atom 2, and atom 2
is 1 sigma from atom 3. With my 2.5 sigma bonding criterion, there will be a
bond between 1-2, 2-3, and 1-3, effectively creating a 3-membered ring.
Thus, a given pair of atoms could have a 1-2 or a 1-3 topology, which leads
to my question about the special bonds command--do 1-2 interactions take
precedence over 1-3 and 1-4 interactions?

yes.

axel.