How do I model a long soft atom chain in LAMMPS?

Dear all,

I am trying to model a system including a long soft atom chain. About 20 atoms are connected in sequence from one end to another end by the harmonic bond. The chain is “soft” since I will not define any angle/dihedral/improper types. However, I do not expect the chain to tangle itself so a repulsive LJ potential (cutoff < 2^(1/6)sigma) will be imposed among all atoms in the chain. How should I modify this?

I tried “special_bonds lj 1 1 1” command but it looks only impose the LJ potential among the chain atoms separated by no more than 3 bonds.

Any advice will be appreciated. Thanks.

Uke

You need to create a “DATA” file containing the starting conformation of the polymer (list of coordinates) as well as the list of bonded atom pairs. (For most people this also includes a list of 3-body angle, 4-body dihedral, and 4-body improper interactions. I realize you don’t need this yet.) The file format is explained here:

http://lammps.sandia.gov/doc/2001/data_format.html

(there’s a small example here)

You can create the file yourself, or you can use one of these programs to build it for you:
http://lammps.sandia.gov/prepost.html

Some DATA files (as in the above example) also include force-field parameters (“Pair Coeffs”, “Bond Coeffs”, “Angle Coeffs”, “Dihedral Coeffs”,…), but this is optional. Axel and I prefer putting the force-field parameters (the “Coeffs” in the example above) in a separate file. For this reason, both moltemplate and topotools (usually) omit the “Coeffs” from the data file and put them in a LAMMPS “input script” file which you can load independently of your starting conformation. (This way, you can tinker with the force field parameters without having to make a new “DATA” file each time.) This is optional.

Both moltemplate and topotools come with examples for building extremely simple polymers, like the one you are describing:

For topotools:
https://sites.google.com/site/akohlmey/software/topotools/topotools-tutorial—part-1

For moltemplate:
http://moltemplate.org/visual_examples.html#2bead_heteropolymer

The procedure for building these polymers is explained in detail in chapter6 of the moltemplate manual
http://moltemplate.org/doc/index.html
You can make a polymer which is straight and has no sidechains, if you omit the 2nd bead, and eliminate the rule for creating dihedral angle interactions (from “Dihedrals By Type” section of the the “force_field.lt” file for that example),
More examples:

http://moltemplate.org/images/misc/polymers_follow_a_curve.jpg

https://github.com/jewettaij/ndmansfield/blob/master/doc/images/moltemplate_usage/wrap_CG_dsDNA_around_a_curve_from_ndmansfield_LLR.png

If you need to build a polymer that wraps around I curve, I recommend using “genpoly_lt.py” to build the .LT file which moltemplate will read (instead of the procedure shown above):

https://github.com/jewettaij/moltemplate/blob/master/doc/utils/docs_genpoly_lt.txt

Hope this helps.
Cheers
Andrew

I tried “special_bonds lj 1 1 1” command

(This command is orthogonal to what you are trying to do.)

Forgive me for jumping to conclusions.
This is the usual response I give to new LAMMPS users.

Upon reading your email again, I guess the question you were really asking was about non-bonded forces.
You can add them to your simulation by defining a pair_style and assigning repulsive interactions between the particles using the “pair_coeff” command.
Something like this:

pair_style lj/cut 2.5
pair_coeff * * 1.0 1.1225 #1.225~=2^(1/6)

This should create repulsive forces whenever ANY two particles get within a distance of 1.1225 (distance units) from each other. (Unless you turned some of the forces between them off using the “special_bonds” command or the “neigh_modify exclude” command.)

(Alternatively, you can put these numbers in the “Pair Coeffs” section of your data file, if you prefer doing it that way.)

For details, see:
http://lammps.sandia.gov/doc/pair_lj.html
http://lammps.sandia.gov/doc/units.html

Nevertheless, check out the examples from the earlier reply.
Cheers
Andrew

Dear all,

I am trying to model a system including a long soft atom chain. About 20
atoms are connected in sequence from one end to another end by the harmonic
bond. The chain is "soft" since I will not define any
angle/dihedral/improper types. However, I do not expect the chain to tangle
itself so a repulsive LJ potential (cutoff < 2^(1/6)sigma) will be imposed
among all atoms in the chain. How should I modify this?

I tried "special_bonds lj 1 1 1" command but it looks only impose the LJ
potential among the chain atoms separated by no more than 3 bonds.

no, that is not correct. nonbonded interactions are applied to *all*
pairs within the cutoff.
however, 1-2, 1-3, and 1-4 pairs get special treatment, since many
force fields require, that there are no non-bonded interactions
computed, when there are already bonded interactions (bonds, angles,
and/or dihedrals). in those cases, the corresponding pairs are removed
from the neighbor list. in your case, you likely would want to apply:
special_bonds lj/cut 0.0 1.0 1.0
this results in no non-bonded interactions are computed for atoms
connected with a bond.

please have a more careful read of the special_bonds command and also
of a text book describing classical force fields.

axel.