Dear all,
I’m trying to set up a DNA strand. I have named the 5’ and 3’ residues DG5 and DC3 in the pdb. I’ve used -ter option in pdb2gmx and specified the 5TER and 3TER endings, respectively. The process ends successfully, as seen in the enclosed text, however, the total charge of the chain is -29.140 e, which doesn’t look like a rounding error. Do you know what could be wrong with my input? The force field files are from charmm36-jul2021.ff.tgz from MacKerell’s web page, and I did not modify them, so I assume they are correct.
Thanks for your time.
Best,
Ramon
Identified residue DG51 as a starting terminus.
Identified residue DC331 as a ending terminus.
8 out of 8 lines of specbond.dat converted successfully
Select start terminus type for DG5-1
0: NH3+
1: NH2
2: 5MET
3: 5PHO
4: 5POM
5: 5TER
6: None
5
Start terminus DG5-1: 5TER
Select end terminus type for DC3-31
0: COO-
1: COOH
2: CT1
3: CT2
4: 3TER
5: None
4
End terminus DC3-31: 3TER
Opening force field file ./charmm36-jul2021.ff/aminoacids.arn
Checking for duplicate atoms....
Generating any missing hydrogen atoms and/or adding termini.
Now there are 31 residues with 981 atoms
Making bonds...
Number of bonds was 1059, now 1057
Generating angles, dihedrals and pairs...
Before cleaning: 2550 pairs
Before cleaning: 2790 dihedrals
Keeping all generated dihedrals
Making cmap torsions...
There are 2790 dihedrals, 84 impropers, 1913 angles
2457 pairs, 1057 bonds and 0 virtual sites
Total mass 9471.066 a.m.u.
Total charge -29.140 e
Writing topology
Writing coordinate file...
--------- PLEASE NOTE ------------
You have successfully generated a topology from: chainJ.pdb.
The Charmm36-jul2021 force field and the tip3p water model are used.
--------- ETON ESAELP ------------
GROMACS reminds you: "Take Dehydrated Water On Your Desert Trips" (Space Quest III)
most probably I’m responsible for the error (I wrote uhe script for converting the CHARMM toppar files to GROMACS format), and also most probably the culprits are the termini. Could you please attach the topology generated by pdb2gmx?
Kind regards, thank you for the bug reports and sorry for the error,
Hi Andras,
Thanks for your quick answer. Here you have the topol.top. Tell me if you need any other file or information.
Best regards,
Ramon topol.top (267.3 KB)
Hi again,
I have just noticed that if I choose None for both termini types I get a more reasonable total charge of -30.000 e. Is that how I should proceed?
I attach the new topol.top in case you want to compare them. topol.top (267.1 KB)
We need to remove any terminus-specific nucleotides. There should be none for CHARMM; this is not like the AMBER force field where such things are needed. Residues should be named ADE, CYT, GUA, THY for DNA and one should use 5TER and 3TER. I did not come across this in my tests because I never use e.g. DG5, etc. because such residues are noncompliant with the CHARMM force field.
If @awacha will please confirm, these need to be removed from your scripts so they are not written to .rtp files.
in my defence (:-)): residues DX3 and DX5 (X in ACGTU) existed in earlier GROMACS ports of the CHARMM FF (at least in the feb2021 version), and they are exactly the same as in the most recent (jul2021) one, created with my scripts.
CHARMM toppar files (namely top_all36_na.rtf) only contains five base residues: GUA, ADE, CYT, THY and URA. These are ribose nucleotides, for constructing RNA. Deoxy residues are made from these in CHARMM by applying the DEOX patch. As GROMACS lacks such advanced patching functionality, deoxy residues (DC, DA, DG, DT and DU) were supplied in former ports of the FF. Technically, these are simply CYT, ADE, GUA, THY and URA, with the DEOX patch already applied. Although these do not exist in the upstream CHARMM, I wouldn’t recommend removing these. 5’ and 3’-terminal deoxy residues are a different case: they are made by applying the 5TER and 3TER patches to the deoxy residues. As termini enties already exist in the na.c.tdb and na.n.tdb files, these ready-made residues should definitely go.
The cause of the error you, @rcrehuet encountered is that applying the 5TER terminus with pdb2gmx to the DG5 residue means that the 5TER patch is applied a second time. This results in an additional HN5 atom (see the first topol.top file you sent earlier), increasing the total charge of your molecule by +0.43.
For the time being, selecting None for the 5’-terminus is probably the easiest choice. The other one would be to name your first residue simply DG and remove the HN5 atom before handing the file to pdb2gmx.
The deoxy residues (DX) should stay indeed, but I still favor removing the 5’- and 3’-specific versions. I don’t ever recall adding these so they probably predated me working with the (old) port conversion script. As we move to greater interoperability between CHARMM and GROMACS (more on that tomorrow at the webinar!) I would like to minimize differences and ad hoc changes, and therefore remove such customized residues, especially when “default” patching leads to errors such as these.
Dear Andras and Justin,
Thanks for your extensive explanations and for delving into the issue in such detail. I am glad that my question will help improve a bit the force field terms. This is one of the less user-friendly aspects of Gromacs (and probably other engines) so any effort to make it more systematic is really welcome.
All the best,
Ramon