[gmx-users] Discarded angle between 3 virtual sites

[Jonathan Barnoud]

Hello,

A colleague of mine stumbled upon an issue with angles between 3 virtual 
sites. If the 3 virtual sites are "Virtual Sites N", then the angle is 
removed as being a constant-energy interaction. This can be prevented by 
passing -normvsbds to grompp. Though, using the option is a workaround 
and, from my understanding, the angle should not be removed.

What makes me think that the angle should not be removed is that:
* I do not see why the energy of the interaction would be constant. I 
built, for instance, a test system where the 3 virtual sites forming the 
angle are built from particles that are not connected to each other. 
There, the angle is free to move except for the angles.
* Bonds and dihedrals between the same virtual sites are not removed.
* The angle is removed only if the 3 virtual sites are Virtual Sites N. 
If any of the particles forming the angle is a virtual site of an other 
type, or is not a virtual site, then the angle is not removed.

I experienced the issue with version 2016, 2018 and 2019-beta3. It 
worked as I expected on 4.6.7. I did not try intermediate versions.

 From my understanding, the issue comes from the function 
clean_vsite_angles in src/gromacs/gmxpreprocess/vsite_parm.cpp. When 
looping over the atoms forming a bond, VSITEN are dealt with as 
non-virtual particles:

# line 1209 in 2019rc1
if (vsite_type[atom] != NOTSET && vsite_type[atom] != F_VSITEN)

I guess the problematic angle is meant to be handled latter at line 1255:

/* keep all angles with no virtual sites in them or
            with virtual sites with more than 3 constr. atoms */
if (nvsite == 0 && vsnral > 3)
{
     bKeep = TRUE;
}

At this point nvsite is 0 but vsnral is 0 as well. (I actually do not 
see when the condition can be true, but I am doing a lot of guessing there.)

Did I miss something about virtual sites or is there an issue with the 
andling of the angles between 3 virtual_sitesn?

Thank you,
Jonathan Barnoud