Problem recognising Iodine atom in GROMACS

User discussions
1

GROMACS version: 2022
GROMACS modification: No
CGenFF Version :4.6
Python script for CHARMM to gmx conversion : for python 2.7
Force-field : CHARMM36-jul 2022.

Hey everyone. I am trying to simulate the dynamics of an iodine-containing ligand with a well studied protein.

While following Dr. Lemkul’s tutorials for simulating Protein-Ligand interactions, I am unable to generate the “ions.tpr” file. GROMACS throws up a fatal error :

“Duplicate atom index in virtual_sites 3”.

The [virtual_sites3] subsection is present in lig.itp file. Further, the lig.gro file has 2 (dummy?) atomic positions named LP1 and LP2, which I presume is suspicious. Upon checking the said subsection in lig.itp, I find the specific mention of the atomic index corresponding to the Iodine atom. I further presume that the lone pairs of Iodine atoms are causing this issue.

Upon searching on the web, I find people reporting this issue being caused in older versions of GROMACS or some versions of the python conversion script. Given that I use very recent versions of all involved programs and scripts, I am unable to find the cause and hence, solution of the error.

Any help will be greatly appreciated.

2

There are no virtual sites in the ligand present in the tutorial, so please provide us more information about what you’re doing, including the relevant topology snippet that defines the virtual sites.

The conversion script uses a hack for an old version of GROMACS (arguably exploiting a bug) and it likely needs to be updated, but it is rather simple to construct the virtual site. If you tell us more about what you’re working with, it should be easy to solve.

3

Respected sir, thank you for replying.

I am working on simulating the dynamics of a well-studied ligand containing 2 Iodine heteroatoms, with a well-studied protein. The ligand under consideration is different from that in your tutorial.

As I try to add ions to the protein-ligand complex, by executing the following :

Screenshot (23)

GROMACS throws up a fatal error :

Screenshot (28)

I am unable to solve this error. On this portal, I find posts suggesting that this error might arise due to usage of old versions of the scripts involved. Since I am using fairly recent versions of all scripts involved, I am unable to find the root cause and hence solution of this error.

Upon some analysis, I find that the .itp file has a sub-section named [virtual_sites3] :

Screenshot (27)
Screenshot (27)765×128 2.86 KB

where 15 and 16 are the indices of the two Iodine atoms.

Also, the .gro file also has two (dummy?) atomic positions named LP1 and LP2 respectively :

Screenshot (29)

In my limited experience with using GROMACS, I have not yet encountered such observations. Are these observations a cause for concern?

On a sidenote, I also tried to simulate the protein-ligand dynamics without the Iodine atoms, which worked without generating any error.

4

You have one virtual site for each I atom, to model the σ hole associated with halogens. See https://www.sciencedirect.com/science/article/pii/S0968089616304576?via%3Dihub

The script you used is indeed outdated, designed for GROMACS versions prior to 2020. We fixed this long ago. Please make sure you’ve got the latest version of the conversion script, available from either Alex MacKerell’s site (MacKerell Lab) or my GitHub (which is where all future development will occur (GitHub - Lemkul-Lab/cgenff_charmm2gmx: Python scripts to convert CGenFF stream files to GROMACS format)

5

Respected sir,
I now used the conversion script you just provided. However, the same error and observations occur.

I believe that the two scripts I used (the script I used originally and the script you just shared) are the same. Hence the similarity in the observations.

6

I don’t see how that’s possible. You have a [virtual_sites3] directive, which we explicitly do not produce. You should get a [virtual_sites2] directive. Please check what code your conversion script has. It should be:

        if (self.nvsites > 0):
            func=2
            f.write("[ virtual_sites2 ]\n")
            f.write("; Site   from              funct a\n")
            for atomi in range (0,self.nvsites):
                vsite = 0
                at1 = 0
                at2 = 0
                # find atom name matches
                for ai in range (0, self.natoms):
                    if (self.G.node[ai]['name'] == self.G.node[atomi]['vsite']):
                        vsite = ai
                    if (self.G.node[ai]['name'] == self.G.node[atomi]['at1']):
                        at1 = ai
                    if (self.G.node[ai]['name'] == self.G.node[atomi]['at2']):
                        at2 = ai
                dist=self.G.node[atomi]['dist']*-1  # invert sign for GROMACS convention
                f.write("%5d %5d %5d %5d %8.3f\n" % (vsite+1, at1+1, at2+1, func, dist))
            f.write("\n")

Make sure you are using a Python3 version. We discontinued development of the Python2 version due to compatibility issues.

7

Using Python3 and the relevant scripts solved the issue. Thank you so much for your help.

Warm Regards