2D Polymer: Area Expansion under NPT Equlibration

Dear Dr. Axel Kohlmeyer and LAMMPS Users,

Hello All. I am currently equilibrating a 2D polymer network subject to NPT ensemble. I would like to obtain a stress free state of this 2D network prior to conducting a stress strain test.

I first used the soft cosine potential that govern the interaction between non-connected particles. The connecting particles that make up the chains are governed by the FENE potential.

I first ran a NVE ensemble and then switched to NPT to relax the 2D system of polymer network and find its equilibrated stress free state.

Using LJ units, I conducted the NPT ensemble at T = 1.0 and p = 0.0. (This is the NPT ensemble that I use to equilibrate my 3D polymer network)

However, the volume (2D area) of the 2D polymer network keeps inflating although the pressure reaches p = 0.0 and T = 1.0. (Interestingly, when I use 3D polymer network, I do not have this problem as the polymer equilibrates to its melt density using Kremer-Grest parameters under NPT: p = 0.0 and T = 1.0)

I am not understanding why the 2D volume is inflating when, in fact, the pressure (e.g pxx and pyy) reaches steadily at 0.0 thereby nothing to push the boundaries of the 2D region outward. Moreover, I see that the free energy of the system remains constant as I observe the total potential energy remains steady.

I searched the previous mailing list, and I see that bad initial geometry (overlapping beads) or not running NVE ensemble may cause this volume expansion. However, I am using soft cosine potential which prevents any sudden repulsion due to overlapping particles.

I would greatly appreciate if I could request for an advice that determines the cause for the continuous volume expansion.

I have attached a small simulation that replicates this issue of this 2D volume expansion under NPT.

Thank you so much for your time.

Sincerely,

Masato Koizumi

2D_NPT_Equilibration.in (871 Bytes)

2D_Polymer.dat (533 Bytes)

I don’t have an easy way to run your simulation today. (I’m replying to this on my phone.)
But I have had some experience simulating a periodic 1d polymer (periodic in the x direction only) at zero pressure.
I will say that took an extremely long time to equilibrate (roughly 10^8 timesteps, I think). I speculate that this was because there was very little interaction between the polymer and it’s periodic images (in the +x and -x) direction. (My polymer was also rather long.)
I am guessing that if I had been running a simulation which was periodic in the x,y,z directions pressure equilibration would have been much faster (because there would have been more collisions between atoms in neighboring periodic images and/or more bonds crossing the periodic boundaries). (This is especially true if you are using isotropic pressure equilibration. Otherwise, if you are using “fix npt” or “fix nph” with “aniso”, then this might no longer be true.) Either way, I suspect that equilibration would also be faster at higher pressures (because the density would be higher and there would be more frequent collisions).
I don’t know if my experience is relevant to yours. Your system might be expanding for any number of reasons.

good luck.

andrew

In your example, you have defined a model with purely repulsive particles, and the additional constraints P = 0 and T > 0. The only stable solution is V->Infinity i.e. the polymer evaporates. If you want a finite volume under these conditions, you will have to provide some sort of cohesive energy term. The same model in 3D might behave differently simple because entanglement effects greatly slow down the evaporation rate.

Hi everybody! Just to expand on what Aidan said, the equilibrium density in ANY system with purely repulsive “intermolecular” interactions and (P = 0, T > 0) is zero. This has nothing to do with the fact that your system is 2D and polymeric. Afaik it would hold in any D for any type of particle.

I suggest running the soft-pushoff phase of the run at constant volume, then switching to the “true” pair interactions (e.g. Lennard-Jones), then equilibrating at constant volume a bit more, and only then switching to NPT with P = 0.

Best,
Rob

Hi everybody! Just to expand on what Aidan said, the equilibrium density in ANY system with purely repulsive “intermolecular” interactions and (P = 0, T > 0) is zero. This has nothing to do with the fact that your system is 2D and polymeric. Afaik it would hold in any D for any type of particle.

I suggest running the soft-pushoff phase of the run at constant volume, then switching to the “true” pair interactions (e.g. Lennard-Jones), then equilibrating at constant volume a bit more, and only then switching to NPT with P = 0.

Best,
Rob

One tiny clarification. If you are simulating infinitely long polymers (ie. polymers that cross the periodic boundary), then the simulation box will not expand forever (because that would break the polymer). But for long polymers it can still take a long time for the box size to equilibrate. (example)

This is not the case for Masato’s polymer. I finally had a look at Masato’s files. He only has 10 particles in the entire simulation. (One polymer which is a chain of 10 consecutively bonded particles.) Masato: Disregarding the other issues, I would build a much larger system before attempting to draw meaningful conclusions from an NPT simulation.

Hello to Rob and Aiden!!

Andrew
Masato, for what it’s worth, this is a good book:
https://www.elsevier.com/books/understanding-molecular-simulation/frenkel/978-0-12-267351-1

Dear All,

Thank you all for your advices! I really appreciate your time for the suggestion!

Thank you Dr. Andrew Jewett for your time in looking into my files and providing with feedback!

I am retrying my simulation and I will see how it goes.

Sincerely,

Masato Koizumi