# How to avoid an 'isopeptide bond' to break in GROMACS? **URL:** \ **Category:** User discussions\ **Tags:** simulation-setup\ **Created:** [April 22, 2021, 1:14am UTC](https://gromacs.bioexcel.eu/t/how-to-avoid-an-isopeptide-bond-to-break-in-gromacs/2007 "2021-04-22T01:14:39Z")\ **Posts on this page:** 3\ **Page:** 1
**Author:** ![fidelormz](https://avatars.discourse-cdn.com/v4/letter/f/dc4da7/32.png) [@fidelormz](https://gromacs.bioexcel.eu/u/fidelormz)\ **Post date:** [April 22, 2021, 1:14am UTC](https://gromacs.bioexcel.eu/t/how-to-avoid-an-isopeptide-bond-to-break-in-gromacs/2007/1 "2021-04-22T01:14:39Z")
GROMACS version: 2020.4 GROMACS modification: No Hi, I am trying to simulate a diubiquitin chain, where the NZ of a lysine in one ubiquitin (atom 2464) is covalently linked to the C-ter of the other (atom 1962). In the PDB file this ‘isopeptide bond’ is defined at the end by: CONECT 1962 2464 CONECT 2464 1962 The CONECT lines are eventually ignored when pdb2gmx writes the topology, and apparently when H atoms are added in both molecules somehow the double bonding of C=O at the C-ter breaks down into C-O (x2). So now the “extra oxygen” clashes with NZ, breaking the isopeptide bond I am trying to follow here. Carrying out the simulation this way both of the involved side chains displace both molecules considerably which of course is not what I am looking. Trying to overcome this I’ve followed a few strategies. At first, and ignoring the appearance of the extra oxygen mentioned above, I applied position restraints at the C-ter and NZ atoms in an attempt to mimic the isopeptide bond as much as possible. After failure I tried to restraint the whole residues (GLY\_C-ter and LYS\_NZ) but given that they are in separate chains, numbering issues regarding the atoms result in an error message similar to ‘Atom index is out of bounds’ (even when the position restraints were included in the topology file at the right position). Finally, I managed both molecules as a single chain and again tried it with position restraints and even as ‘frozen groups’; although in both cases the atoms remained well fixed, now an artifact is formed between C-ter atom of the 1st ubiquitin and both N-ter and NZ atoms of the second (giving place to this weird ‘pentavalent carbon’). Is there any other strategy I could try? Or a manner in which I may force this covalent bond to stay? I’m frankly lost. Thank you very much, Fidel. ---
**Author:** ![jalemkul](https://dub1.discourse-cdn.com/flex017/user_avatar/gromacs.bioexcel.eu/jalemkul/32/18_2.png) [@jalemkul](https://gromacs.bioexcel.eu/u/jalemkul)\ **Post date:** [April 22, 2021, 12:35pm UTC](https://gromacs.bioexcel.eu/t/how-to-avoid-an-isopeptide-bond-to-break-in-gromacs/2007/2 "2021-04-22T12:35:43Z")
Branched bonds such as these require the use of `specbond.dat` to define the bond. You will have to create a new lysine `.rtp` entry with a different name and only one HZ proton, and use that residue name in `specbond.dat` as the renamed residue, and then define the bond between NZ and the carbonyl C of the C-terminus. You will also have to choose ‘None’ as the C-terminal patch type, though I don’t know if `pdb2gmx` will complain about that case (missing atoms). But that’s the only way to attempt making such a bond. ---
**Author:** ![fidelormz](https://avatars.discourse-cdn.com/v4/letter/f/dc4da7/32.png) [@fidelormz](https://gromacs.bioexcel.eu/u/fidelormz)\ **Post date:** [April 23, 2021, 1:57am UTC](https://gromacs.bioexcel.eu/t/how-to-avoid-an-isopeptide-bond-to-break-in-gromacs/2007/3 "2021-04-23T01:57:14Z")
Thank you very much for your reply @jalemkul. I tried your suggestion in accordance with the manual. So I created a new lysine (LYS1) which was defined within the `aminoacids.rtp` file (I copied those data from LYS and simply erased the extra H atoms): [LYS1] [atoms] N NH1 -0.47 0 HN H 0.31 1 CA CT1 0.07 2 HA HB 0.09 3 CB CT2 -0.18 4 HB1 HA 0.09 5 HB2 HA 0.09 6 CG CT2 -0.18 7 HG1 HA 0.09 8 HG2 HA 0.09 9 CD CT2 -0.18 10 HD1 HA 0.09 11 HD2 HA 0.09 12 CE CT2 0.21 13 HE1 HA 0.05 14 HE2 HA 0.05 15 NZ NH3 -0.30 16 HZ1 HC 0.33 17 C C 0.51 20 O O -0.51 21 [bonds] CB CA CG CB CD CG CE CD NZ CE N HN N CA C CA C +N CA HA CB HB1 CB HB2 CG HG1 CG HG2 CD HD1 CD HD2 CE HE1 CE HE2 O C NZ HZ1 [impropers] N -C CA HN C CA +N O [cmap] -C N CA C +N Further I declared LYS1 as Protein in the `residuetypes.dat` and updated `specbond.dat` by adding the line: GLY C 1 LYS1 NZ 1 0.2 ISO ISO My guess here is that LYS1 would now appear as an option when I run `pdb2bmx` with `-inter` to define the protonation state of each lysine, but this is not the case. I am using CHARMM27 as forcefield. Is there any other file that I miss to modify? Since I am not adding any new atom or bond types I did not edit the `.itp` files regarding bonded or non-bonded interactions, but perhaps I am supposing it wrong. Besides, by choosing ‘None’ for the C-terminal patch type the error message regarding ‘dangling bond at the terminal end’ actually appeared. So I tried to modify the `c.tdb` file in a similar way that LYS1, tagging this new option as ‘CO-’. But again, no option appeared available in the list after `-inter` was used. I would thank you for any other advice, Fidel.