Completing inter-monomer bonds

User discussions
,
1

GROMACS version: 2020.5
GROMACS modification: No
I am currently working on constructing a polymer topology in GROMACS, and I am having difficulty correctly defining the inter-monomer bonds. My aim is to connect the carbon atom C2 in one monomer with the oxygen atom O4 in the next monomer.

Below is the current [bonds] section from my topology file:
[bonds]
; ai aj funct r k
1 2 1 1.5380e-01 2.5179e+05 ; C2 - C3
1 9 1 1.4320e-01 2.5824e+05 ; C2 - O5
1 15 1 1.0960e-01 2.7757e+05 ; C2 - H9
2 3 1 1.4230e-01 2.6501e+05 ; C3 - O2
2 4 1 1.5380e-01 2.5179e+05 ; C3 - C4
2 16 1 1.0970e-01 2.7665e+05 ; C3 - H10
3 13 1 9.7300e-02 3.1079e+05 ; O2 - H4
4 5 1 1.4230e-01 2.6501e+05 ; C4 - O3
4 6 1 1.5380e-01 2.5179e+05 ; C4 - C5
4 17 1 1.0970e-01 2.7665e+05 ; C4 - H11
5 14 1 9.7300e-02 3.1079e+05 ; O3 - H5
6 7 1 1.4320e-01 2.5824e+05 ; C5 - O4
6 8 1 1.5380e-01 2.5179e+05 ; C5 - C6
6 18 1 1.0970e-01 2.7665e+05 ; C5 - H12
8 9 1 1.4320e-01 2.5824e+05 ; C6 - O5
8 10 1 1.5240e-01 2.6192e+05 ; C6 - C7
8 19 1 1.0970e-01 2.7665e+05 ; C6 - H13
10 11 1 1.2180e-01 5.3363e+05 ; C7 - O6
10 12 1 1.2180e-01 5.3363e+05 ; C7 - O7

Based on my understanding, I need to add two lines to define the inter-monomer bond: one for C2 to O4 in the next monomer and another for O4 in the next monomer back to C2 (maybe a line expressing O4 - -C2 and other C2 - +O4?). Also, I am unsure about the exact syntax and parameters (r and k) required . Should I take the r and k values after running ACPYPE for my capping group? (I deleted them before it, just using them for RESP charges calculation but not for parmchk2, tleap, and acpype steps) If so, should I delete the atoms, bonds, pairs, angles, etc from my capping groups or just delete the atoms but keep some bonds and pairs, etc related with the inter-monomer bond?

Could someone please guide me on how to properly define these inter-monomer bonds? Please take into consideration that the [bonds] section I give you here is just for the monomeric repetitive unit (not the head or tail monomeric units). Any advice or reference to relevant documentation would be greatly appreciated.

Thank you for your assistance.

2

Hi,

Bonds are not directional, so no need to have the “back” bond.

This should be solvable in 3 ways:

  • adding lines in specbond.dat (note that the actual distances have to be within +/- 10% of the values declared there)
  • modifying the .rtp and pretending it’s a protein (might produce orthogonal issues with generic amino acid processing)
  • using the add_bond() function from Gromologist to do the linking after the topologies of individual residues are generated as separate [ moleculetype ]s
3

If I modify the .rtp file, when should I do it? These files are usually in the folder gromacs/top/, but in my case I already used GLYCAM_06j-1 in tleap (so there is no .ff folder in which to modify any aminoacids.rtp file). Alternatively, could you please help me to understand how to use add_bond() function from Gromologist? I have a master top file calling all the itp files (each one a building block GLC, GLA, GLZ, GLX, and GLT). The GLC is repetitive subunit, GLA is head, GLZ is tail, always connecting C2 to the O4 in the next subunit. GLX is intermediate like GLC but with a linker in COO group and GLT in GLC is also intermediate but truncated (also GLX and GLT are connected to each other)

Could I put in my master top file an extra [bonds] section (extra because each subunit has one) with lines like this?

; Include monomer topologies
#include “./GLC.itp”
#include “./GLT.itp”
#include “./GLX.itp”
#include “./GLA.itp”
#include “./GLZ.itp”

[ bonds ]
; Define inter-subunit bonds
; ai aj funct r k
1 7 1 0.1538 251790 ; GLA-C2 - GLC-O4
1 9 1 0.1538 251790 ; GLC-C2 - GLZ-O4
1 7 1 0.1538 251790 ; GLC-C2 - GLC-O4
1 7 1 0.1538 251790 ; GLC-C2 - GLX-O4
1 7 1 0.1538 251790 ; GLX-C2 - GLC-O4
1 7 1 0.1538 251790 ; GLC-C2 - GLT-O4
1 7 1 0.1538 251790 ; GLT-C2 - GLC-O4
10 39 1 0.1538 251790 ; GLT-C7 - GLX-N4

4

Sorry, no idea if you have solved it already, but add_bond needs just atom indices (1-based) from the selected molecule, from the other molecule, and the other molecule entry as a parameter. These won’t be inter-chain bonds, but rather the tool will merge two chains into one and then add the bonds intra-chain.

If you read your topology as t, then t.molecules will list the molecules defined in the topology, and you can refer to them as t.molecules[0], t.molecules[1] etc.

And you can do, say, t.molecules[0].add_bond(20, 42, t.molecules[1]) to add a bond between atom 20 of molecule 0 and atom 42 of molecule 1. After that, you just save the topology with t.save_top('newtopo.top').