From 98f44d28f09d52cf82a0a2f398bbb55d78647e32 Mon Sep 17 00:00:00 2001 From: Thomas Frederiksen Date: Fri, 3 Apr 2020 02:58:11 +0200 Subject: [PATCH 1/6] enh: initial code for bond completion (bond saturation) --- sisl/geometry.py | 64 +++++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 63 insertions(+), 1 deletion(-) diff --git a/sisl/geometry.py b/sisl/geometry.py index afee4db2e2..22e0e58898 100644 --- a/sisl/geometry.py +++ b/sisl/geometry.py @@ -23,7 +23,7 @@ from .utils.mathematics import fnorm from .quaternion import Quaternion from .supercell import SuperCell, SuperCellChild -from .atom import Atom, Atoms +from .atom import Atom, Atoms, PeriodicTable from .shape import Shape, Sphere, Cube from ._namedindex import NamedIndex @@ -3066,6 +3066,68 @@ def bond_correct(self, ia, atom, method='calc'): raise NotImplementedError( 'Changing bond-length dependent on several lacks implementation.') + def bond_completion(self, nbonds, atom=None, bond=None, idx=None): + """ Return a new geometry with additional atoms added to complete the number of bonds + + This function may be useful to saturate dangling bonds at edges in sp2-carbon structures, + e.g., with H-atoms. + + Parameters + ---------- + nbonds : int + number of bonds requested + atom : `Atom`, optional + the kind of atom to be added, where missing. Defaults to atoms of the same type. + bond, float, optional + bond length to the extra atoms. Defaults to value from `PeriodicTable` + idx : array_like, optional + List of indices for atoms that are to be considered + + Examples + -------- + >>> g = geom.graphene(orthogonal=True).tile(3, 0).tile(4, 1) + >>> g.cell[0] *= 2 + >>> g.bond_completion(3, atom=Atom(1), bond=1.09) + """ + + if idx is not None: + if not isndarray(idx): + selection = _a.asarrayi(idx).ravel() + else: + selection = range(len(self)) + + new_xyz = [] + new_atom = [] + + PT = PeriodicTable() + + for ia in selection: + iaZ = self.atoms[ia].Z + ria = PT.radius(iaZ) + idx = self.close(ia, R=(0.1, 0.1 + 2 * ria))[1] + if len(idx) == nbonds - 1: + # Add one bond + if atom is None: + atom = Atom(iaZ, R=self.atoms[ia].R) + if bond is None: + bond_length = ria + PT.radius(atom.Z) + else: + bond_length = bond + # Compute bond vector + bvec = len(idx) * self.xyz[ia] + for jdx in idx: + bvec -= self.axyz(jdx) + bnorm = bvec.dot(bvec) ** .5 + if bnorm > 0.1: + # only add to geometry if new position is away from ia-atom + bvec *= bond_length / bnorm + new_xyz.append(self.xyz[ia] + bvec) + new_atom.append(atom) + if len(new_xyz) > 0: + return self.add(self.__class__(new_xyz, new_atom)) + else: + return self.copy() + def within(self, shapes, idx=None, idx_xyz=None, ret_xyz=False, ret_rij=False): From 2d0d383ea5ecac8f368ea859a7511010166a0c25 Mon Sep 17 00:00:00 2001 From: Thomas Frederiksen Date: Fri, 3 Apr 2020 10:26:07 +0200 Subject: [PATCH 2/6] maint: bond_completion with .close(..., ret_xyz=True) call --- sisl/geometry.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/sisl/geometry.py b/sisl/geometry.py index 22e0e58898..e05f63cb36 100644 --- a/sisl/geometry.py +++ b/sisl/geometry.py @@ -3104,8 +3104,10 @@ def bond_completion(self, nbonds, atom=None, bond=None, idx=None): for ia in selection: iaZ = self.atoms[ia].Z ria = PT.radius(iaZ) - idx = self.close(ia, R=(0.1, 0.1 + 2 * ria))[1] - if len(idx) == nbonds - 1: + idx = self.close(ia, R=(0.1, 0.1 + 2 * ria), ret_xyz=True) + # We just need second shell coordinates + xyz = idx[1][1] + if len(xyz) == nbonds - 1: # Add one bond if atom is None: atom = Atom(iaZ, R=self.atoms[ia].R) @@ -3114,9 +3116,7 @@ def bond_completion(self, nbonds, atom=None, bond=None, idx=None): else: bond_length = bond # Compute bond vector - bvec = len(idx) * self.xyz[ia] - for jdx in idx: - bvec -= self.axyz(jdx) + bvec = len(xyz) * self.xyz[ia] - np.sum(xyz, axis=0) bnorm = bvec.dot(bvec) ** .5 if bnorm > 0.1: # only add to geometry if new position is away from ia-atom From 01583647bb6be2337056629b71a98cac5ba05c48 Mon Sep 17 00:00:00 2001 From: Thomas Frederiksen Date: Fri, 3 Apr 2020 11:06:20 +0200 Subject: [PATCH 3/6] maint: check bnorm first --- sisl/geometry.py | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/sisl/geometry.py b/sisl/geometry.py index e05f63cb36..24bc405b42 100644 --- a/sisl/geometry.py +++ b/sisl/geometry.py @@ -3108,21 +3108,21 @@ def bond_completion(self, nbonds, atom=None, bond=None, idx=None): # We just need second shell coordinates xyz = idx[1][1] if len(xyz) == nbonds - 1: - # Add one bond - if atom is None: - atom = Atom(iaZ, R=self.atoms[ia].R) - if bond is None: - bond_length = ria + PT.radius(atom.Z) - else: - bond_length = bond # Compute bond vector bvec = len(xyz) * self.xyz[ia] - np.sum(xyz, axis=0) bnorm = bvec.dot(bvec) ** .5 if bnorm > 0.1: # only add to geometry if new position is away from ia-atom + if bond is None: + bond_length = ria + PT.radius(atom.Z) + else: + bond_length = bond bvec *= bond_length / bnorm new_xyz.append(self.xyz[ia] + bvec) - new_atom.append(atom) + if atom is None: + new_atom.append(Atom(iaZ, R=self.atoms[ia].R)) + else: + new_atom.append(atom) if len(new_xyz) > 0: return self.add(self.__class__(new_xyz, new_atom)) else: From bda662de73c0704ab37a861e058a07a0aca0edc3 Mon Sep 17 00:00:00 2001 From: Thomas Frederiksen Date: Fri, 3 Apr 2020 11:13:45 +0200 Subject: [PATCH 4/6] bug: determine which atom to add before looking up radius --- sisl/geometry.py | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/sisl/geometry.py b/sisl/geometry.py index 24bc405b42..acdfd58ae5 100644 --- a/sisl/geometry.py +++ b/sisl/geometry.py @@ -3113,16 +3113,17 @@ def bond_completion(self, nbonds, atom=None, bond=None, idx=None): bnorm = bvec.dot(bvec) ** .5 if bnorm > 0.1: # only add to geometry if new position is away from ia-atom + if atom is None: + this_atom = Atom(iaZ, R=self.atoms[ia].R) + else: + this_atom = atom if bond is None: - bond_length = ria + PT.radius(atom.Z) + bond_length = ria + PT.radius(this_atom.Z) else: bond_length = bond bvec *= bond_length / bnorm new_xyz.append(self.xyz[ia] + bvec) - if atom is None: - new_atom.append(Atom(iaZ, R=self.atoms[ia].R)) - else: - new_atom.append(atom) + new_atom.append(this_atom) if len(new_xyz) > 0: return self.add(self.__class__(new_xyz, new_atom)) else: From 840d75b0e83e24a997dd3479113215ac27d42273 Mon Sep 17 00:00:00 2001 From: Thomas Frederiksen Date: Fri, 3 Apr 2020 11:34:16 +0200 Subject: [PATCH 5/6] enh: let user determine which distance to be considered an existing bond --- sisl/geometry.py | 32 +++++++++++++++++++------------- 1 file changed, 19 insertions(+), 13 deletions(-) diff --git a/sisl/geometry.py b/sisl/geometry.py index acdfd58ae5..bb85307f30 100644 --- a/sisl/geometry.py +++ b/sisl/geometry.py @@ -3066,7 +3066,7 @@ def bond_correct(self, ia, atom, method='calc'): raise NotImplementedError( 'Changing bond-length dependent on several lacks implementation.') - def bond_completion(self, nbonds, atom=None, bond=None, idx=None): + def bond_completion(self, nbonds, bond=None, new_bond=None, atom=None, idx=None): """ Return a new geometry with additional atoms added to complete the number of bonds This function may be useful to saturate dangling bonds at edges in sp2-carbon structures, @@ -3076,10 +3076,12 @@ def bond_completion(self, nbonds, atom=None, bond=None, idx=None): ---------- nbonds : int number of bonds requested - atom : `Atom`, optional - the kind of atom to be added, where missing. Defaults to atoms of the same type. bond, float, optional - bond length to the extra atoms. Defaults to value from `PeriodicTable` + distance between existing atoms to be considered a bond. Defaults to value from `PeriodicTable` + new_bond, float, optional + bond length to added atom. Defaults to value from `PeriodicTable` + atom : `Atom`, optional + kind of atom to be added, where missing. Defaults to atoms of the same type. idx : array_like, optional List of indices for atoms that are to be considered @@ -3102,9 +3104,13 @@ def bond_completion(self, nbonds, atom=None, bond=None, idx=None): PT = PeriodicTable() for ia in selection: - iaZ = self.atoms[ia].Z - ria = PT.radius(iaZ) - idx = self.close(ia, R=(0.1, 0.1 + 2 * ria), ret_xyz=True) + a = self.atoms[ia] + # Determine radius of sphere for bond search + if bond is None: + ria = 0.1 + 2 * PT.radius(a.Z) + else: + ria = 1e-4 + bond + idx = self.close(ia, R=(0.1, ria), ret_xyz=True) # We just need second shell coordinates xyz = idx[1][1] if len(xyz) == nbonds - 1: @@ -3114,16 +3120,16 @@ def bond_completion(self, nbonds, atom=None, bond=None, idx=None): if bnorm > 0.1: # only add to geometry if new position is away from ia-atom if atom is None: - this_atom = Atom(iaZ, R=self.atoms[ia].R) + new_a = Atom(a.Z, R=a.R) else: - this_atom = atom - if bond is None: - bond_length = ria + PT.radius(this_atom.Z) + new_a = atom + if new_bond is None: + bond_length = PT.radius(a.Z) + PT.radius(new_a.Z) else: - bond_length = bond + bond_length = new_bond bvec *= bond_length / bnorm new_xyz.append(self.xyz[ia] + bvec) - new_atom.append(this_atom) + new_atom.append(new_a) if len(new_xyz) > 0: return self.add(self.__class__(new_xyz, new_atom)) else: From a5f5d1183f658854520e6253df20a835bd710978 Mon Sep 17 00:00:00 2001 From: Thomas Frederiksen Date: Fri, 3 Apr 2020 11:44:48 +0200 Subject: [PATCH 6/6] doc: updated example --- sisl/geometry.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/sisl/geometry.py b/sisl/geometry.py index bb85307f30..be2356f86a 100644 --- a/sisl/geometry.py +++ b/sisl/geometry.py @@ -3089,7 +3089,7 @@ def bond_completion(self, nbonds, bond=None, new_bond=None, atom=None, idx=None) -------- >>> g = geom.graphene(orthogonal=True).tile(3, 0).tile(4, 1) >>> g.cell[0] *= 2 - >>> g.bond_completion(3, atom=Atom(1), bond=1.09) + >>> g.bond_completion(3, bond=1.42, new_bond=1.09, atom=Atom(1)) """ if idx is not None: @@ -3111,7 +3111,7 @@ def bond_completion(self, nbonds, bond=None, new_bond=None, atom=None, idx=None) else: ria = 1e-4 + bond idx = self.close(ia, R=(0.1, ria), ret_xyz=True) - # We just need second shell coordinates + # We just need second-shell coordinates xyz = idx[1][1] if len(xyz) == nbonds - 1: # Compute bond vector