Single ion virtual sites using [ virtual_sites1 ] in Gromacs 2021

GROMACS version: 2021
GROMACS modification: No


Gromacs 2021 supports virtual_sites1! Would someone like to explain to me the possible benefit of this, could it be used on a single site, for instance on a single ions?

I have tried it with solvated Na+ ions (see my modifications below), trying to divide the sodium ions LJ potential 50%:50% between the real Na+ sites and their virtual site (by using the same LJ-sigma and half the LJ-epsilon value for each site). I does not work… have I misunderstood the purpose of virtual_sites1 ?

From my ions.itp

[ moleculetype ]
; molname nrexcl
Na 1

[ atoms ]
; id at type res nr residu name at name cg nr charge
1 Na 1 Na Na 1 1.00000 22.98977
2 DNa 1 Na DNa 1 0.00000 0.00000

[ virtual_sites1 ]
; Site from funct
2 1 1 ; is the last 1 needed? Unclear from the docs

[ exclusions ]
1 2
2 1

From my ffnonbonded.itp

[ atomtypes ]
Na 11 22.98977 1.00000 A 0.215954 0.737725 ; Aqueous sodium ion with 0.5epsilon value. Originally 1.47545
DNa 0 0.00000 0.00000 V 0.215954 0.737725 ; Virtual sodium site with 0.5
epsilon value



The purpose is to be able to have multiple interaction sites on the same atom location. This can be useful for free-energy calculations in particular.

In your example you simply double the LJ interactions, as you now have two interactions instead of 1. You should also exclude LJ interactions between the atom and the site, but the GROMACS kernels mask out LJ and Coulomb interactions at zero distance, so this actually works.

That is in line with my general understanding, but it seems that I do need to add an [ exclusions ] section for Na+ ion and its virtual site (see my first post) since if not I get extreme Total potential energies (10^30 kJ/mol).

It also seems like I cannot transfer all of the charge (+1.0) and/or the LJ-sigma and epsilon values from the real Na+ site to the virtual site without the simulation crashing, which I initially thought would be possible. However, using exclusions and transferring half the charge and ~half the LJ-epsilon value (but it obviously depends on the mixing rules… ) but using the orig LJ-sigma value for both the real and virtual Na+, seems to produce the similar simulation results as in a normal simulation without using any virtual single ion sites,



I now realize that the LJ combination rules can affectt how you should redistribute the LJ parameters.

Yes, that is what I meant with ‘it depends on the mixing rules’, although I did not think of it at first either.

Anyway, I think I understand how it works now.

Hi mholmboe i have been struggling with this virtual sites for several months now. I am trying to use virtual_sites but i first worry is how to include in the pdb or gro file? Could you please hint me on how this actually works?


my suggestion to know how gro files and topology files looks like is to run the following simple test.

  1. take a pdb file of a simple molecule (e.i lysozime)
  2. run gmx pdb2gmx -f …pdb -o -p with and without the option -vsite
  3. look at the difference in the outputs


1 Like

It is hard to answer without knowing exactly what you want to do, but you could start by trying what alevilla is suggesting below.
Michael Holmboe