LAMMPS WWW Site - LAMMPS Documentation - LAMMPS Mailing List Archives
Re: [lammps-users] Semi-Flexible polymers in 2D

# Re: [lammps-users] Semi-Flexible polymers in 2D

 From: Stefan Paquay Date: Tue, 27 Jun 2017 10:12:13 -0400

Oops, did I forget that? I changed the script a bit but the following deck will reproduce the issue. Atoms 49 and 48 are the first to pick up huge momentum due to them being too close.

On Mon, Jun 26, 2017 at 11:46 PM, Trung Nguyen wrote:
Hi Stefan,

is it possible to share an example (small-sized) chains_flat.data for a complete input deck that reproduces the issue?

-Trung

On Mon, Jun 26, 2017 at 5:56 PM, Stefan Paquay wrote:
So, the bottom line seems to be: If you're going to use angle potentials, make sure the bond lengths are sufficiently far from zero, because the angle becomes ill-defined if the bonds are too short.

On Mon, Jun 26, 2017 at 6:48 PM, Stefan Paquay wrote:
Hi Andrew,

It's not as easy as that. Reducing the angular potential leads to less straight polymers, which in turn leads to a smaller likelihood of the problem manifesting itself. If the polymer is almost straight, it is easiest to get a bond length close to 0, because (in 2D) there is only variation along a single axis. If the angle potential is small, the bonds explore a two-dimensional space, and the probability of the bond length being zero in two dimensions is a lot smaller.

I've isolated one of these instabilities, and at the onset there are indeed two particles close to each other. There are no pair forces between them but the overlap probably sends the angular potential through the roof, since in angle_harmonic the cosine is determined from these distances. In this particular run, the distance between two of the atoms becomes less than 0.003, which, I can imagine, will lead to numerical complications. These two are the ones that then pick up a huge momentum and fly outwards, and due to the bonds, they take the rest with them.

On Mon, Jun 26, 2017 at 4:19 PM, Andrew Jewett wrote:
Hi Stefan

1) also (if you haven't done this yet) try weakening the strength of the 3-body angle potential to see when/if the issue goes away.  If the problem persists with a zero-strength harmonic angle force, that would be worth reporting.

2) another way to approach is to use "fix wall ...harmonic" to squish a 3-D polymer between two plates.  As the distance between the plates decreases, perhaps you will see the unstable behavior you noticed.  If not, then that would be weird.  (The only problem with that idea is that I've noticed that "fix wall" often causes my simulations to become unstable, so doing this might just complicate your problem and not resolve anything.  I have been to busy to dig into why this happens in my simulations.  Its probably my own fault.)

Andrew

On Jun 26, 2017 12:40 PM, "Stefan Paquay" <stefanpaquay@...24...> wrote:
Thanks for thinking with me Steve.

One of my colleagues remarked that potentially it could be due to the short bond length. This makes it possible for particles to overlap due to thermal fluctuations, especially if there's no steric repulsion between neighbours. Of course, if the particles are too close, their angle is not well-defined and they blow up. I'll run some analysis on this to see what the bond lengths are just before blowing up. Increasing the bond length seems to mitigate the issue so I'm thinking this might explain it just fine.

On Mon, Jun 26, 2017 at 10:13 AM, Steve Plimpton wrote:
also, another idea, you can use compute reduce
to print out the max force on any atom at each step.
If you do that every step, you might see a snapshot
or two where things go bad.  Then you can examine
the bad atoms on that step.

Steve

On Mon, Jun 26, 2017 at 6:35 AM, Steve Plimpton wrote:
I don't know any reason 2d would be different than 3d,
unless there is some sort of numerical instability in
the interactions you've specified, e.g. as angles go to 180.

You could print out forces every tilmestep and see if something
is blowing up, leading to bad dynamics.  Once you have a bad
force for a particular config, you should be able to debug
where it comes from.

Steve

On Sat, Jun 24, 2017 at 10:45 PM, Stefan Paquay wrote:
Hi all,

I have noticed some not very intuitive behaviour when using a semi-flexible polymer (by means of harmonic bond and angle potentials) in a 2D simulation. See the attached input script.

I naively assumed that just subtracting out the force in the Z-direction will make this automagically work, because the angles can still relax in the 2D plane (and moreover, if the preferred angle is 180 degrees, everything should be fine).

However, what I see is that in the 2D simulation the polymer dynamics become unstable and the particles are all over the place. The same settings in 3D lead to a perfectly well-behaved situation. Disabling the angular potential removes the instability but, obviously, also doesn't lead to a semi-flexible polymer. Using a planar constraint from USER-MANIFOLD leads to the same issue (not surprisingly).

Therefore my question: Are there some caveats with regards to angle potentials and 2D simulations that I'm not aware of? Maybe it is related to the impossibility of constraining angles at 180 degrees with SHAKE? Is there any way in which this can work?

------------------------------------------------------------------------------
Check out the vibrant tech community on one of the world's most
engaging tech sites, Slashdot.org! http://sdm.link/slashdot
_______________________________________________
lammps-users mailing list
lammps-users@...655....net
https://lists.sourceforge.net/lists/listinfo/lammps-users

------------------------------------------------------------------------------
Check out the vibrant tech community on one of the world's most
engaging tech sites, Slashdot.org! http://sdm.link/slashdot
_______________________________________________
lammps-users mailing list
lammps-users@...655....net
https://lists.sourceforge.net/lists/listinfo/lammps-users

------------------------------------------------------------------------------
Check out the vibrant tech community on one of the world's most
engaging tech sites, Slashdot.org! http://sdm.link/slashdot
_______________________________________________
lammps-users mailing list
lammps-users@...655....net
https://lists.sourceforge.net/lists/listinfo/lammps-users

```units lj
dimension 2
boundary p p p

atom_style molecular

mass * 1.0

bond_style harmonic
bond_coeff 1 50 0.1

angle_style harmonic
angle_coeff 1 500 180

pair_style zero 5.0
pair_coeff * *
pair_modify shift yes

neighbor 6.0 bin
neigh_modify exclude molecule/intra all every 1

min_style hftn
# Super safe time step size:
timestep 0.0001

fix step all nve
fix twod all enforce2d
velocity all create 1.0 412
minimize 1e-8 1e-12 10000 100000

reset_timestep 0

dump traj2 all custom 10000 chains_long_step.dump.bin id mol type x y z fx fy fz

dump traj all custom 50 chains_hi_res.dump.bin id mol type x y z fx fy fz

thermo 500000
run 2000

```

Attachment: chains_flat.data
Description: Binary data