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.