- /* Exact */
- inv_mass = 0;
- for (j = 0; j < 3*ffparams->iparams[il.iatoms[i]].vsiten.n; j += 3)
- {
- int aj = il.iatoms[i + j + 2];
- coeff = ffparams->iparams[il.iatoms[i+j]].vsiten.a;
- if (moltype->atoms.atom[aj].ptype == eptVSite)
- {
- m_aj = vsite_m[aj];
- }
- else
- {
- m_aj = moltype->atoms.atom[aj].m;
- }
- if (m_aj <= 0)
+ switch (ilist.functionType)
+ {
+ case F_VSITE2:
+ /* Exact */
+ vsite_m[a1] = (cam[1]*cam[2])/(cam[2]*gmx::square(1 - ip.vsite.a) + cam[1]*gmx::square(ip.vsite.a));
+ break;
+ case F_VSITE3:
+ /* Exact */
+ vsite_m[a1] = (cam[1]*cam[2]*cam[3])/(cam[2]*cam[3]*gmx::square(1 - ip.vsite.a - ip.vsite.b) + cam[1]*cam[3]*gmx::square(ip.vsite.a) + cam[1]*cam[2]*gmx::square(ip.vsite.b));
+ break;
+ case F_VSITEN:
+ GMX_RELEASE_ASSERT(false, "VsiteN should not end up in this code path");
+ break;
+ default:
+ /* Use the mass of the lightest constructing atom.
+ * This is an approximation.
+ * If the distance of the virtual site to the
+ * constructing atom is less than all distances
+ * between constructing atoms, this is a safe
+ * over-estimate of the displacement of the vsite.
+ * This condition holds for all H mass replacement
+ * vsite constructions, except for SP2/3 groups.
+ * In SP3 groups one H will have a F_VSITE3
+ * construction, so even there the total drift
+ * estimate shouldn't be far off.
+ */
+ vsite_m[a1] = cam[1];
+ for (int j = 2; j < maxj; j++)