Single polymer chain with 2 values of Lennard Jones 12-6 Sigma

LAMMPS LAMMPS Mailing List Mirror
1

Hello LAMMPS users,

I am trying to simulate a polymer chain using Langevin dynamics (30th July 2016 version of LAMMPS). The monomers interact via bond (harmonic), bond angle (harmonic), dihedral angle (multi-harmonic) and Lennard Jones 12-6 (lj-cut) potentials. I am facing an issue with implementing the Lennard Jones Potential, because it has one value of sigma for monomers closer than 5 repeat units along the chain, and another value for monomers farther than five repeat units along the chain.

I looked at the existing variations of the Lennard Jones 12-6 potentials, as well as special_bonds to add weighting factors, but I could not apply the exactly desired potential. So, I went ahead and modified the pair_lj_cut.cpp with some simple additional lines in the compute function (inside the if (rsq < cutsq[itype][jtype])). I added if-else cases for the force calculation and the energy calculation, and verified that the calculations were correct.

In spite of doing this, I cannot get the simulation to be stable. At some point, the atoms just start flying away (ERROR on proc 0: Bond atoms 662 663 missing on proc 0 at step 4508 (../neigh_bond.cpp:65)). I started my simulation with the lowest possible potential energy structure (bond length, bond angle, dihedral angle all at their most stable values), so the starting structure cannot be an issue. I have also tried playing around with the neighbor lists, but this doesn't help.

Any suggestions regarding this matter will be of great help. Maybe there might be an easy trick that I've overlooked? Thanks in advance.

Thanks,
Kiran Iyer
PhD Student
Department of Chemical Engineering
University of Massachusetts Amherst

2

Are simulation paramerters (partial charges) right ?

Arun

3

Dear Arun,

Thanks for your response. What do you mean by partial charges? If you mean the interaction potentials themselves, then yes. They are correct, and taken from a published article: "Modeling polymer crystallization from solutions" Polymer 41 (2000) 8833-8837.

The chain by itself consists of uncharged monomers.

Thanks,
Kiran Iyer

4

Hello LAMMPS users,

I am trying to simulate a polymer chain using Langevin dynamics (30th
July 2016 version of LAMMPS). The monomers interact via bond (harmonic),
bond angle (harmonic), dihedral angle (multi-harmonic) and Lennard Jones
12-6 (lj-cut) potentials. I am facing an issue with implementing the
Lennard Jones Potential, because it has one value of sigma for monomers
closer than 5 repeat units along the chain, and another value for
monomers farther than five repeat units along the chain.

I looked at the existing variations of the Lennard Jones 12-6
potentials, as well as special_bonds to add weighting factors, but I
could not apply the exactly desired potential. So, I went ahead and
modified the pair_lj_cut.cpp with some simple additional lines in the
compute function (inside the if (rsq < cutsq[itype][jtype])). I added
if-else cases for the force calculation and the energy calculation, and
verified that the calculations were correct.

In spite of doing this, I cannot get the simulation to be stable. At
some point, the atoms just start flying away (ERROR on proc 0: Bond
atoms 662 663 missing on proc 0 at step 4508 (../neigh_bond.cpp:65)). I
started my simulation with the lowest possible potential energy
structure (bond length, bond angle, dihedral angle all at their most
stable values), so the starting structure cannot be an issue. I have
also tried playing around with the neighbor lists, but this doesn't
help.

Any suggestions regarding this matter will be of great help. Maybe there
might be an easy trick that I've overlooked? Thanks in advance.

the most likely reason is, that your modified pair style is incorrect.
what you are trying to do is not as trivial as it may seem and
extrapolating from the typical level of correctness of the code from
people who do hacks like you describe, it is highly probably, that you
are making incorrect assumptions in your code. thus it will break down
at a later point, e.g. after a neighbor list update.

there is a way to implement what you want to do without having to
change source code: you have to give each atom in your polymer chain a
different atom type and then assign suitable parameters to all pairs
of atom types according to your topology condition. this may require
many pair_coeff lines, and due to the repetitive nature, it is
probably best to write a small script/program to generate this part of
the input.

axel.

5

Dear Axel,

Thank you very much for your response. My simulation is stable now.

Thanks and regards,
Kiran Iyer