AtomGroup.guess_bonds misleading documentation, could pass in fudge_factor

[lilyminium]

Expected behavior

mdanalysis/package/MDAnalysis/core/groups.py

Lines 2918 to 2926 in 971248a

	def guess_bonds(self, vdwradii=None):
	"""Guess bonds that exist within this :class:`AtomGroup` and add them to
	the underlying :attr:`~AtomGroup.universe`.
	Parameters
	----------
	vdwradii : dict, optional
	Dict relating atom types: vdw radii

AtomGroup.guess_bonds guesses bonds, only.

Actual behavior

It guesses bonds, angles, and dihedrals.

If the method is going to guess angles, it could guess improper dihedrals too from the angles.

Code to reproduce the behavior

Show us how to reproduce the failiure. If you can, use trajectory files from the test data.

mdanalysis/package/MDAnalysis/core/groups.py

Lines 2937 to 2961 in 971248a

	from ..topology.core import guess_bonds, guess_angles, guess_dihedrals
	from .topologyattrs import Bonds, Angles, Dihedrals
	def get_TopAttr(u, name, cls):
	"""either get *name* or create one from *cls*"""
	try:
	return getattr(u._topology, name)
	except AttributeError:
	attr = cls([])
	u.add_TopologyAttr(attr)
	return attr
	# indices of bonds
	box = self.dimensions if self.dimensions.all() else None
	b = guess_bonds(self.atoms, self.atoms.positions, vdwradii=vdwradii, box=box)
	bondattr = get_TopAttr(self.universe, 'bonds', Bonds)
	bondattr.add_bonds(b, guessed=True)
	a = guess_angles(self.bonds)
	angleattr = get_TopAttr(self.universe, 'angles', Angles)
	angleattr.add_bonds(a, guessed=True)
	d = guess_dihedrals(self.angles)
	diheattr = get_TopAttr(self.universe, 'dihedrals', Dihedrals)
	diheattr.add_bonds(d)

Currently version of MDAnalysis

• Which version are you using? current develop branch

Also

It would be useful to be able to pass in fudge_factor and lower_bound kwargs to this method (to pass onto guessers.guess_bonds)

[richardjgowers]

So impropers are an odd case in that most impropers aren't used (as in have a potential associated with them) - they're mostly used to keep things planar. So guessing them all would produce a lot of useless entries.

Passing in the other kwargs to guess_bonds makes sense

[soma2000-lang]

@lilyminium Is this issue still open for working .Then I would like to take up this issue.

[xhgchen]

I am interested in looking into this further but unsure what changes need to be made exactly. Would someone be able to explain a little more clearly? I have looked around for the terms and found fudge_factor described in the guessing page (https://userguide.mdanalysis.org/stable/formats/guessing.html) but I am unsure what lower_bound means. Thanks!

[lilyminium]

@xhgchen the function guess_bonds (below) takes in the keyword arguments fudge_factor and lower_bound. AtomGroup.guess_bonds calls this function but doesn't have these keyword arguments itself. A fix to this issue would make those keyword arguments available to AtomGroup.guess_bonds.

mdanalysis/package/MDAnalysis/topology/guessers.py

Lines 157 to 262 in 971248a

	def guess_bonds(atoms, coords, box=None, **kwargs):
	r"""Guess if bonds exist between two atoms based on their distance.
	Bond between two atoms is created, if the two atoms are within
	.. math::
	d < f \cdot (R_1 + R_2)
	of each other, where :math:`R_1` and :math:`R_2` are the VdW radii
	of the atoms and :math:`f` is an ad-hoc *fudge_factor*. This is
	the `same algorithm that VMD uses`_.
	Parameters
	----------
	atoms : AtomGroup
	atoms for which bonds should be guessed
	coords : array
	coordinates of the atoms (i.e., `AtomGroup.positions)`)
	fudge_factor : float, optional
	The factor by which atoms must overlap eachother to be considered a
	bond. Larger values will increase the number of bonds found. [0.55]
	vdwradii : dict, optional
	To supply custom vdwradii for atoms in the algorithm. Must be a dict
	of format {type:radii}. The default table of van der Waals radii is
	hard-coded as :data:`MDAnalysis.topology.tables.vdwradii`. Any user
	defined vdwradii passed as an argument will supercede the table
	values. [``None``]
	lower_bound : float, optional
	The minimum bond length. All bonds found shorter than this length will
	be ignored. This is useful for parsing PDB with altloc records where
	atoms with altloc A and B maybe very close together and there should be
	no chemical bond between them. [0.1]
	box : array_like, optional
	Bonds are found using a distance search, if unit cell information is
	given, periodic boundary conditions will be considered in the distance
	search. [``None``]
	Returns
	-------
	list
	List of tuples suitable for use in Universe topology building.
	Warnings
	--------
	No check is done after the bonds are guessed to see if Lewis
	structure is correct. This is wrong and will burn somebody.
	Raises
	------
	:exc:`ValueError` if inputs are malformed or `vdwradii` data is missing.
	.. _`same algorithm that VMD uses`:
	http://www.ks.uiuc.edu/Research/vmd/vmd-1.9.1/ug/node26.html
	.. versionadded:: 0.7.7
	.. versionchanged:: 0.9.0
	Updated method internally to use more :mod:`numpy`, should work
	faster. Should also use less memory, previously scaled as
	:math:`O(n^2)`. *vdwradii* argument now augments table list
	rather than replacing entirely.
	"""
	# why not just use atom.positions?
	if len(atoms) != len(coords):
	raise ValueError("'atoms' and 'coord' must be the same length")
	fudge_factor = kwargs.get('fudge_factor', 0.55)
	vdwradii = tables.vdwradii.copy() # so I don't permanently change it
	user_vdwradii = kwargs.get('vdwradii', None)
	if user_vdwradii: # this should make algo use their values over defaults
	vdwradii.update(user_vdwradii)
	# Try using types, then elements
	atomtypes = atoms.types
	# check that all types have a defined vdw
	if not all(val in vdwradii for val in set(atomtypes)):
	raise ValueError(("vdw radii for types: " +
	", ".join([t for t in set(atomtypes) if
	not t in vdwradii]) +
	". These can be defined manually using the" +
	" keyword 'vdwradii'"))
	lower_bound = kwargs.get('lower_bound', 0.1)
	if box is not None:
	box = np.asarray(box)
	# to speed up checking, calculate what the largest possible bond
	# atom that would warrant attention.
	# then use this to quickly mask distance results later
	max_vdw = max([vdwradii[t] for t in atomtypes])
	bonds = []
	pairs, dist = distances.self_capped_distance(coords,
	max_cutoff=2.0*max_vdw,
	min_cutoff=lower_bound,
	box=box)
	for idx, (i, j) in enumerate(pairs):
	d = (vdwradii[atomtypes[i]] + vdwradii[atomtypes[j]])*fudge_factor
	if (dist[idx] < d):
	bonds.append((atoms[i].index, atoms[j].index))
	return tuple(bonds)

[xhgchen]

With all the new contributors coming in, I just want to say that I think this is a great starter issue to work on, especially since Lily has taken the time to explain things so clearly :)