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

Are simulation paramerters (partial charges) right ?

Arun

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

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.

Dear Axel,

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

Thanks and regards,

Kiran Iyer