[lammps-users] Simulation of polysaccharides

Dear Lammps users, I’m using Lammps for the first time and I’m facing the problem of simulating polysaccharides. I generated a .pdb file of a monomer and then I calculated the partial charges and optimized geometry using ATB web server. After that, I used moltemplate to create LAMMPS input file. I ran the simulation and all was okay. The problems came when I tried to bind two monomer with moltemplate. I generated the lammps input file of the dimer and, when I ran the simulation, strange and meaningless bond were created.
Have anyone of you ever faced the same problem?
Can you suggest me how to build a correct LAMMPS input file for a polysaccharide?

Thanks in advance for your help.

Nico

What you describe is a complex procedure and sometimes difficult to get right. However, your questions are too vague and at the same time too generic (how to generate a system) and too specific (for polysaccharides).

My suggestion is to build your skills in stages and do some “practice runs” with simple spring-bead polymers and learn how to build connections with those.
run simulations, analyze and check that everything is working correctly. Then move on to something slightly more complex. and so on.

Asking for a “do this, not that” set of instructions for a rather specific kind of system is not likely going to be very successful.
The only case where this is working is simulations of proteins, because they are rather regular, have (mostly) simple single chains and have only a small variation in their internal structure. And even there you can run across rather complex problems if you have modified or non-standard components in there.

Axel.

Dear Lammps users, I’m using Lammps for the first time and I’m facing the problem of simulating polysaccharides. I generated a .pdb file of a monomer and then I calculated the partial charges and optimized geometry using ATB web server. After that, I used moltemplate to create LAMMPS input file. I ran the simulation and all was okay. The problems came when I tried to bind two monomer with moltemplate. I generated the lammps input file of the dimer and, when I ran the simulation, strange and meaningless bond were created.

Dear Nicola

You have found the Achilles’ heel of the ATB server. The ATB molecules are physically realistic, but they cannot be easily modified (except by submitting a new request to the ATB server). However there are strategies you can use. (See below.)

The ATB web server generates ready to use files for molecules (in moltemplate LT) format carefully optimized charge and bonded interactions using DFT. Hopefully the molecules are at least as accurate (if not more accurate) than molecules created with most force-fields (like OPLSAA, GAFF, or CHARMM). However I don’t know of a safe way to modify existing ATB molecules. Of course, you can manually add a bond between atoms using moltemplate (or topotools, or by editing the DATA file with a text editor), but that does not mean that the new molecule will behave accurately. Adding a new bond will only constrain the distance between the two atoms. The angles can rotate and spin freeling. I assume that this is not what you want. You will have to manually add new angles, dihedrals, and impropers as well that contain the pair of atoms you connected, anc constrain these angles. (It may be necessary to modify the atomic charges and angular interactions as well after you add a new bond, assuming you want the molecule to behave accurately.) Again, this is because ATB is not a force field. It is a molecule server. I included some suggestions below for ways to do this. (See part 1)

For comparison, most force-fields (OPLSAA, GAFF, CHARMM, COMPASS, PCFF, DREIDING,…) provide rules for generating new angles, dihedrals, and impropers automatically between bonded atoms. If you add a bond, these rules will automatically create the additional angular interactions that you need (assuming your atom types are selected correctly). But the “GROMOS_54A7_ATB.lt” file (provided by the ATB server) does not include these rules. (I could be wrong, but I am under the impression that GROMOS-based force fields do not include rules for generating angles and dihedrals. If I’m wrong, please let me know. This is something I have been curious about.)

There are several options (that I am aware of):

  1. Use the ATB server to build (or help build) the full-size polymer.

1a) If all of the polymers are the same length, then you can submit a request to ATB for the complete polymer, and use that LT file in your simulation. (If you want to include polymers of different lengths, then do this once for every length of polymer. So if you want to simulate a system containing 100 polymers of length 5, and 100 polymers of length 6, you only have to submit 2 requests to the ATB server: a polymer of length 5, and a polymer of length 6.) There is probably an upper limit to the size of the molecules that ATB can generate, but I’m not sure what that limit is. Try it and find out. Hopefully this works well.

1b) If you have a long polymer, and/or if the ATB server refuses your request, then use the ATB server to create an LT file for a dimer. This LT file contains all of the angles, dihedrals, and impropers necessary to connect your monomers together and build polymers of arbitrary size. Open up the LT file for the dimer created by ATB using a text editor. Go to the “Data Angles” section of the LT file and make note of all of the angles which contain the two atoms in the bond that connects the two monomers. Do the same thing for the “Data Dihedrals” and “Data Impropers” sections. When you create a longer polymer, you will have to include all of these interactions between the atoms in neighboring monomers. You can try doing this manually. There are several ways to do this automatically:

i) For polymers, you can use the “genpoly_lt.py” script to add these interactions. It is documented with examples here:
https://github.com/jewettaij/moltemplate/blob/b9eea0edb051ba3ea3b54cea006a32d63cb346de/doc/doc_genpoly_lt.md#genpoly_ltpy

ii) Alternatively, a more general solution is to create your own (miniature) force-field with rules for generating these kinds of angles, dihedrals, and impropers for the junction between monomers.
The official documentation describing these rules is here:

http://moltemplate.org/doc/moltemplate_manual.pdf#subsection.7.1

…and there are examples how to use these rules here:
http://moltemplate.org/doc/moltemplate_talk_2019-8-15.pdf#page=118
http://moltemplate.org/doc/moltemplate_talk_2013-8-07.pdf#page=54

http://moltemplate.org/examples/force_fields.html

Don’t worry about the fact that the LT files already contain angles, dihedrals, and impropers. The angles created by these rules will not override existing angles, dihedrals, and impropers that already exist in (the “Data Angles”, “Data Dihedrals”, and “Data Impropers” sections of) your LT file. (They have priority.) These rules will only add angles, dihedrals, and impropers if they are missing. However to make sure the polymers you create this way are reasonable, try creating a system with moltemplate that contains only one dimer using these rules, and then open up the “system.data” file generated by moltemplate. The number of angles, dihedrals, and impropers should match the number of angles, dihedrals, and impropers in the dimer LT file that was created directly by the ATB server. (If not, then I can show you how to get around this issue by creating additional atom types.)

Keep in mind that the atomic charges will need to be modified as well to ensure that the final molecule has integer charge. (The procedure(s) described above only explain how to generate angle, dihedral, and improper interactions.)

  1. Try building the molecule using a popular force field.

2a) If the molecules are not too complicated, then you can probably do this with moltemplate, but you must be careful to choose force-field-specific atom types for the atoms in your molecule. (Moltemplate is not yet smart enough to infer atom types from from atomic elements (“atom typing”).) Start with OPLSAA. If you can’t figure out which atom types to use, try GAFF, GAFF2, or DREIDING. The DREIDING force field in moltemplate might not be the most accurate, but it has relatively simple rules for choosing atom types. If you use DREIDING, GAFF or GAFF2, you must use external software (or fix qeq) to calculate the atomic charges. (If you use OPLSAA or COMPASS, the atomic charges are included with the force field, although the charges for that force field are not optimized for your molecule. The COMPASS force field is not available unless you can find an unencrypted “compass.frc” file. If so, I can show you how to convert it to LT format.)

2b) If your molecule is complicated enough that you cannot manually select the correct atom types for all of the atoms in your molecule, then try using other molecule builder software (and convert the files back into LAMMPS or moltemplate format later). First, try using EMC. It can often do atom-typing automatically, and it outputs to LAMMPS data format directly. It may be possible to build the entire system using EMC. Alternatively, if you want to eventually use moltemplate to assemble the final system (from individual molecular parts), then you use EMC to build a system containing only one molecule (for example your polymer) and then convert the LAMMPS DATA file from EMC into moltemplate format using the ltemplify.py tool, and then moltemplate to duplicate that molecule and move it into position. (Ltemplify works with LAMMPS DATA files created by other molecule builders as well.)

If your molecule is complicated, you might consider using the SMIRNOFF force field (which has huge coverage and does not require atom typing), build the system. If you want to use LAMMPS (instead of OpenMM) to run the simulation, then convert the resulting simulation file(s) into LAMMPS format. (You can probably use OpenBabel to do the conversion. I’m sorry to be so vague. I have never tried this. One day, I hope to get the OpenFF force fields like SMIRNOFF working in moltemplate eventually, so no conversion is necessary.)

I hope this gets you started.
I’m sorry that this process is so complicated.

Andrew