Simulation of novel terminal DNA residues with Charmm36

User discussions
1

GROMACS version: Any, but I am using 2020.1
GROMACS modification: No

Hi guys, hopefully you will be able to help me out with this. I am trying to use gromacs to simulate a DNA nanostructure, in which some of the 3’ terminal residues are modified with a cholesterol, attached via a tetra ethylene glycol linker (TEG-Cl), using the charmm36 ffs.
I can generate the charmm36 parameters for the TEG-Cl using cgenff, or swissparam, and these have been pretty sufficient for some previous, similar work I have done, which was performed in NAMD.

I have been able to generate NAMD style .psf topologies using psfgen by ammending the swissparam generated .rtf to my nucleic acids topology and applying it as a 3’ terminal patch. However I would really like to be able to achieve this in gromacs (The dream would be to build the system in charmn-gui).

So I have generated the parameters and topologies for my TEG-Cl, oriented it correctly in VMD and concatenated the pbd files. I have no problems using the swissparam topologies to treat the TEG-Cl as a ligand, but am struggling to register it as part of the DNA chain.

It is my understanding that I will need to edit the merged.rtp in my charmm ff folder to add an entry for the TEG-Cl, using the swissparam outputted .itp as a base. I am running into a little trouble here as the swissparam .itp has a few more columns than my merged.rtp.
Will I also need to modify ffbonded.itp if I want to have the TEG-Cl bonded the rest of the DNA? and are there any other parts of the forcefields I would need to adapt?
let me know if there is any more info that would be useful, I can ofc provide the coordinate & parameter files if you think that would be useful.
All the best,
Jonah

2

Hi Jonah,

One way is to follow your idea, adding an .rtp entry to the FF and going through pdb2gmx. The extra columns in .itp are mostly obsolete for .rtp, since the latter can get e.g. masses for each atom type from ffnonbonded.itp (atom and residue numbering is context-dependent anyway), but in .rtp you need to provide all bonds and impropers (can be listed in a friendly format with list_bonds() in Gromologist, see below) - angles, dihedrals and 1-4 pairs can be easily inferred from bonds. For pdb2gmx to recognize things as part of the chain, you might need a local copy of residuetypes.dat that registers your custom residue as a DNA residue: that’s a frequently overlooked step.

Another way, if you have both topologies already generated (e.g. as two .itps within a single .top), is to remove the excess atoms (like terminal hydrogens) and introduce a bond between the two residues with Gromologist (look up Editing topologies for examples). This is assuming one of the residues has a terminal phosphate that you can use for linking.

Once done, run your system through gmx grompp (or use find_missing_ff_params() from Gromologist) to see which parameters are not included in your FF files. Both will complain about each bonded term that is not defined, so that you can include them by analogy and modify [ bondtypes ], [ angletypes ] etc accordingly.

Let me know if any of that helps.

Best regards,
Miłosz

3

Hi Milosz,

I’d not heard of the gromologist tool, this looks very useful, and I will spend some time getting to grips with it and get back to you.
Thank you for you time,
Jonah

4

Indeed I made it to automatize common non-standard topology editing routines - it assumes your force field is nothing out of ordinary, but that should be the case here. If you find any issues or have questions, please let me know.

Best,
Miłosz

5

Hi Milosz,

the gromolgist tool is amazing, thank you for this it is exactly what I need.
So I have followed your suggestions, combined my top & pdb files, opened in gromologist, and made the required bond. Then I renumbered the atoms & residues and outputted fresh topologys and coordinates, great.
I put this together in grompp and it was throwing up atomtype not found errors for everything in the novel residue, but I solved this by using clear_ff_params() in gromologist, and re-adding an include statement in the outputted .top for the ligand.itp from swissparam (caled tegchol.itp).

Then, as you said, using grompp or find_missing_ff_params() I can see there is a a bond type, along with a few U-B and dihedral types missing. I have been able to infer the parameters relatively easily, but I’m a little stuck at what (i hope) may be the last hurdle - getting these parameters recognised by grompp.

After combining in gromologist, my .top looks like this;
#define _FF_CHARMM

; Include ff parameters
#include “ffparams.itp”
#include “tegchol.itp” ; swissparam generated parameters for novel residue

; Include DNA_chain_A topology
#include “DNA_chain_A_ed.itp” ; pdb2gmx generated parameters for DNA chain

[ system ]
; Name
ring_chainA

[ molecules ]
; Compound #mols
DNA_chain_A 1

grompp has shown given me line numbers in the DNA_chain_A.itp that have the missing parameters, eg.

[ bonds ]
; ai aj funct c0 c1 c2 c3
571 572 1
571 573 1
573 620 1 ← location of error

so I have been scrolling to the next section within [ bonds ], and adding in parameters for the corresponding bond

; ai aj fu b0 kb, b0 kb
620 690 1 0.1093 287014.9 0.1093 287014.9
621 691 1 0.0972 469365.3 0.0972 469365.3
573 620 1 0.1418 303937.5 0.1418 303937.5 ← added line.

but this does not seem to change the error message I am receiving from grompp. I have also tried adding in these new parameters to the ffparams.itp and tegchol.itp files that are called in my topology but without much joy.
I feel there maybe something I am missing here?
Also, when editing parameter files to add new bonds, does the order matter? The only reference I can find in the gromacs manual is that this is a “free” format, which makes me think no?

Best regards,
Jonah

6

I think I’ve cracked it, my error was adding in a line on the .itp file with the new bond definitions, as oppose to simply adding the parameters to the line specified by grompp. I will leave this up as a testament to my own impatience, and possibly help to anyone else working on something similar.

7

Hi Jonah,

Glad you liked it! I believe there’s a warning message when running, but be sure to double-check for stray #ifdef/#endif commands from the merger of two #ifdef’ed [ position_restraint ] sections in the merged molecule - this is the only thing that remains to be made fully error-proof in one way or another.

Regarding your question, are you sure you need to duplicate the bond parameters? This is usually used for alchemical transformations to specify A/B state values; I don’t think it would somehow affect the recognition of the parameter as incorrect, but just to be sure.

Otherwise, if you’re #including the correct files, either editing the [ bonds ] (in the molecule params) or [ bondtypes ] (in the FF files) section should work, so without looking directly at your files I can’t really tell what might have gone wrong.

Order used to matter for dihedral wildcards only - in versions up to 201something, wildcards (X) would be used preferentially if they were listed before the more specific dihedral. If you have two identical entries, there will be a warning saying which one is used and which one ignored.

Best,
Miłosz

8

You’ve helped me solve in a day what I’ve been toiling on for most of the week, many many thanks.

9

Great, happy to have helped! Remember same struggles back in the day, exactly what prompted me to automatize this.