GROMACS version:gromacs-2022.3
GROMACS modification: Yes
I am attempting to implement a Hunter-Sanders model for pi-pi interactions. To that end I have added DUM atoms into a structure on the aromatic rings in my protein and ligand. I then used a script to insert reference to them into the protein.itp and ligand.itp files (I also added intra-residue/ligand [ exclusions ] and [ virtual_sites ] or the DUM atoms).
My first simulations were a disaster so I reran the trajectory with energygroups. I discovered that the virtual particles were forming 1,4 interactions with my lipids, molecules of which they are not (supposed to be) members.
I am not sure what to make of this. I was thinking that there are no [ bonds ] or [ constraints ] defined between the DUM atoms on the virtual sites and the corresponding proteins and ligand. Is it possible that [ pairs ] are being generated automatically between the DUM and my lipids?
If so, how do I stop this? The topology is assembled on a per molecule basis so [ exclusions ] are only within a molecule (right?). And I don’t think that I want to add [ constraints ] or [ bonds ].
An update:
I ran gmx_fda and it does not show 1,4 interactions between the dum atoms on the virtual sites and the other molecules like POPC. One note, gmx_fda requires tpr file made with version 2020.4 which doesn’t recognize the [ virtual_sites1 ] directive I had in my .itp file. I needed to change that to [ virtual_sitesn ]. The documentation makes it seem that the latter can replicate the former.
I am not willing to discount the dummy-lipid 1,4 interactions that I saw in the interaction energy dump, given how badly my lipid bilayer behaves (note it behaves fine without the dummy atoms on the virtual sites.)
What led you to this conclusion? 1-4 interactions are exclusively intramolecular. If you have some kind of virtual site in a protein or small molecule, it’s cannot engage in 1-4 interactions with lipids. That would only be possible if the virtual sites were part of the lipid molecule definition.
Hi @jalemkul,
Before I dive in let me just state my thanks for your response and that no further response is requested in the following, but would still be appreciated.
What led me to that conclusion is that I used gmx make_ndx to make a group of the Dummy atoms on the virtual sites. I then added that group and the existing POPC index group to the energygroups line in my .mdp file and reran the trajectory with gmx mdrun...-rerun. When I look into the resulting .edr file, I see dummy-POPC 1,4 interactions listed, which surprised me.
Your response confirms for me that that is indeed suspicious.
At first I thought that this was the consequence of whatever was causing the aforementioned disastrous simulation, maybe some index error that was considering the virtual atoms to be part of the wrong molecule or something.
Taking a hint from my previous use of virtual_sitesn in gmx_fda, I changed to using that instead of virtual_sites1 in a new simulation. The interesting thing is that disaster with the shrinking height (z direction) of the box and the squashing of the membrane is averted, but the rerun .edr still shows 1,4 interactions.
I’ve rechecked my .itp files and they are all indexed correctly and consistently with the input .pdb file. So the bad 1,4 interactions (if they even exist or are erroneously being claimed to by the rerun edr) are, if not independent of, not as strongly correlated with the squish I was seeing.
So I have two different issues:
- squish when using
virtual_sites1 - 1,4 interactions showing in the rerun edr.
The first seems to be fixed by switching to virtual_sitesn. I diffed the gmx dump (only learned later about gmx check) of the two .tprs, one with virtual_sites1 the other with n. The only oddity I see is that on the line of the dumped tpr in which VSITE1 is defined the line isn’t broken before defining the next functype
functype[3137]=CONSTR, dA= 1.33000001e-01, dB= 1.33000001e-01
functype[3138]=VSITE1, functype[3139]=VSITE3OUT, a= 0.00000000e+00, b= 0.00000000e+00, c= 2.32668877e+00
functype[3140]=VSITE3OUT, a= 0.00000000e+00, b= 0.00000000e+00, c= 2.50276303e+00
But the format looks fine for VSITEN in the other .tpr dump. If I have time, I’ll investigate further. But for now I have my fix.
In regard to second issue, now that I think that it isn’t directly related to the disastrous results of the first simulation and I’ve seen by gmx_fda that it might not be so in the force calculation, I am open to the possibility that it is a reporting error and not actually affecting the trajectory. Again, if I have time, I’ll investigate further. But for now, if virtual_siten trajectory is reasonable and given gmx_fda shows no intermolecular 1,4 forces and I only checked the rerun edr because of the first issue, I’m willing to table this.