creating long (infinite) polymers with moltemplate

I am looking for a utility to create an infinite polymer chain from the SMILES strings.
I was able to generate the coordinates of a monomer from SMILES using openbabel. Now, got stuck with creating a polymer chain. I can use GaussView/Avagadro for a single chain, but we are looking for creating polymers for a huge database.
If you have any tools (even Moltemplate), can you suggest me the steps.

Hi Karteek

I don’t have any tools that can comprehend SMILES strings.

But if you have used some other tool to obtain the coordinates of a monomer, you can manually create a moltemplate monomer object (similar to “ch2group.lt” for polyethylene. See the “ch2group.lt” file in two examples: 1, 2) If the monomer is simple (as in polyethylene), I typically create the monomer file (eg “ch2group.lt”) by hand using a text editor. If it is complicated, I use a third-party tool to create the monomer, save it in LAMMPS-data format, and convert it to moltemplate format using “ltemplify.py” (as explained in slides 165-185 of this talk.)

Then you can create a moltemplate polymer by using “new” to make many copies of the monomer, and then by adding bonds between the monomers. (See the “butane.lt” file, and the “alkane50.lt” file in those two examples.)

That’s what you would do for a short polymer. I’m not sure what you mean by an “infinite polymer”.

  1. Are you trying to create a polymer with periodic boundary conditions that spans the box? (A polymer with many monomers with the first monomer at one end of the box and the last monomer at the opposite side of the box, with the two monomers bound together.)

a) If so, do you want the initial shape of the polymer to be straight?

(as in this example, see this video. This is a periodic system looking down the length of the long axis. Over time, the length of the box shrinks and the polymer bends and twists. Here’s another example showing the result after the box length has been reduced considerably.)

Either way, place the first polymer at one side of the box, and the last monomer at the other side of the box, and create a bond between them.
If the polymer is initially straight, follow the long alkane chain example, and edit the “Data Boundary” section in the “system.lt” file to insure that the length of the box matches the length of the polymer. And don’t forget to add a bond to the “Data Bonds” section in the polymer file (eg "alkane50.lt) that connects an atom from the last monomer to an atom in the first monomer.

b) Or do you want it to be a bent, space-filling-curve (as in this example)?

There are many ways to create a self-avoiding space-filling curve. Once you have the coordinates of such a curve you can use the “genpoly_lt.py” script to convert this curve into a polymer. (See step iii below.)

i) One of them is to use “ndmansfield” (documentation here) to create an extremely coarse, low-resolution random self-avoiding curve. This will create a curve that lies on a 3D lattice and fills a box. For an idea what this looks like, click on the links above. Since this is unlikely to be a realistic polymer shape, later on you can smooth the polymer and relax it. (see steps ii and iv)

ii) You may want to make this polymer shape more smooth before you use it. One way to do that is to use the the “interpolate_mt.py” script to read the coarse coordinates in the curve generated by “ndmansfield” and then add additional points between the lattice points. This allows you to smooth the curve using interpolation. The “interpolate_mt.py” script also allows you to rescale the x,y,z coordinates (multiply them by a constant. By default, “ndmansfield” generates coordinates which are integers, so you will probably want to use this feature, unless you want all of your peptide bonds to have length 1.) Note: The “interpolate_mt.py” script is included with recent versions of moltemplate and it is documented here.

iii) After that, you can use the “genpoly_lt.py” script to create a polymer that “wraps along the curve”. This script will read the file containing the coordinates and place monomers at each of these locations. (that you generated in steps i and ii) It will also rotate them appropriately, and add bonds connecting them together (between the atoms you specify. (You can also ask it to create angles and dihedrals (etc…) between monomers.) The “genpoly_lt.py” script is included with recent versions of moltemplate and is documented here.

The last step (iii) is explained in slightly more detail here:
http://moltemplate.org/images/misc/polymers_follow_a_curve.png

It is also explained briefly in slides 206-211 from the moltemplate talk:

http://moltemplate.org/doc/moltemplate_talk_2019-8-15.pdf

Note: Slide 211 from that talk outlines the process.

(iv)
You will need to relax the polymer (because initially the polymer is sharply bent). I strongly suggest using a “soft” pair style which will initially allow the monomers to pass through each other (such as the various lj/soft pair_styles, or using pair_style lj/charmm which is documented here, and can be downloaded here.) As the simulation progresses, you can slowly increase the energy penalty for two atoms to occupy the same location until the polymer becomes a self-avoiding chain again. (That’s how this polymer was generated.)

Note: It’s not necessary to use the “genpoly_lt.py” script with curves prepared using “ndmansfield” and “interpolate_mt.py” programs. If you have another way to generate a text file containing a curve with the coordinates of every monomer, you can use that file with “genpoly_lt.py”. “genpoly_lt.py” is a general script and will accept 3d coordinates from any source.

c) Or both? (A long space-filling curve whose end points are connected to opposite sides of the simulation box.)

In that case, follow the procedure (i,ii,iii,iv) explained above in b), but be sure that the first and last monomers are at opposite sides of the box and connect them with a bond. If you are using “ndmansfield” to generate the initial polymer coordinates, then you can use the the “-cyclic yes” and “-cyclic direction 0” arguments (assuming periodicity in the x direction, see documentation) to ensure that the first and last monomers are on opposite ends. Again don’t forget to add a bond connecting them.

Summary:
None of these steps are trivial to use. But at least there are general tools that can do all of the steps in the process, and they are fully documented. (Although perhaps not easy to read.) One day I might create a tutorial to make this easier.

If you have any questions, please email back.

Andrew