# DUM atoms on virtual sites forming 1,4 interactions with other molecules
**URL:** https://gromacs.bioexcel.eu/t/dum-atoms-on-virtual-sites-forming-1-4-interactions-with-other-molecules/6121
**Category:** User discussions
**Tags:** simulation-setup
**Created:** [April 5, 2023, 9:04pm UTC](https://gromacs.bioexcel.eu/t/dum-atoms-on-virtual-sites-forming-1-4-interactions-with-other-molecules/6121 "2023-04-05T21:04:11Z")
**Posts on this page:** 4
**Page:** 1
### Author:  [@aheim](https://gromacs.bioexcel.eu/u/aheim)
#### Post date: [April 5, 2023, 9:04pm UTC](https://gromacs.bioexcel.eu/t/dum-atoms-on-virtual-sites-forming-1-4-interactions-with-other-molecules/6121/1 "2023-04-05T21:04:11Z")
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].
---
### Author:  [@aheim](https://gromacs.bioexcel.eu/u/aheim)
#### Post date: [April 6, 2023, 10:00pm UTC](https://gromacs.bioexcel.eu/t/dum-atoms-on-virtual-sites-forming-1-4-interactions-with-other-molecules/6121/2 "2023-04-06T22:00:23Z")
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.)
---
### Author:  [@jalemkul](https://gromacs.bioexcel.eu/u/jalemkul)
#### Post date: [April 8, 2023, 5:58pm UTC](https://gromacs.bioexcel.eu/t/dum-atoms-on-virtual-sites-forming-1-4-interactions-with-other-molecules/6121/3 "2023-04-08T17:58:07Z")
> [@aheim](#):
>
> 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.
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.
---
### Author:  [@aheim](https://gromacs.bioexcel.eu/u/aheim)
#### Post date: [April 12, 2023, 7:29pm UTC](https://gromacs.bioexcel.eu/t/dum-atoms-on-virtual-sites-forming-1-4-interactions-with-other-molecules/6121/4 "2023-04-12T19:29:01Z")
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:
1. squish when using `virtual_sites1`
2. 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 `.tpr`s, 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`
```auto
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.