Skip to content

Commit

Permalink
Merge pull request #4 from duartegroup/dev
Browse files Browse the repository at this point in the history 
Faster metallocage building
  • Loading branch information
@t-young31
t-young31 authored Jul 5, 2020
2 parents 15b3aad + 5e2ad20 commit d5c33fc
Show file tree
Hide file tree
Showing 27 changed files with 639 additions and 588 deletions.
4 changes: 2 additions & 2 deletions README.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -31,10 +31,10 @@ conda install --file requirements.txt

## Installation

For installation instructions see the [installation instructions](https://duartegroup.github.io/cgbind/install.html).
For detailed instructions see the [installation instructions](https://duartegroup.github.io/cgbind/install.html).
If the requirements are already satisfied
```
python setup.py install
pip install cgbind
```

***
Expand Down
108 changes: 36 additions & 72 deletions cgbind/atoms.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,20 +4,27 @@

class Atom:

def __init__(self, atomic_symbol, x, y, z):
def __init__(self, atomic_symbol, x=0.0, y=0.0, z=0.0, coord=None):

self.label = atomic_symbol
self.coord = np.array([float(x), float(y), float(z)])

# Coordinate takes priority over x, y, z values
if coord is not None:
self.coord = coord


heteroatoms = ['O', 'N', 'S', 'P', 'F', 'Cl']

metals = ['Li', 'Be', 'Na', 'Mg', 'Al', 'K', 'Ca', 'Sc', 'Ti', 'V', 'Cr', 'Mn', 'Fe', 'Co', 'Ni', 'Cu', 'Zn', 'Ga',
'Rb', 'Sr', 'Y', 'Zr', 'Nb', 'Mo', 'Tc', 'Ru', 'Rh', 'Pd', 'Ag', 'Cd', 'In', 'Sn', 'Cs', 'Ba', 'La', 'Ce',
'Pr', 'Nd', 'Pm', 'Sm', 'Eu', 'Gd', 'Tb', 'Dy', 'Ho', 'Er', 'Tm', 'Yb', 'Lu', 'Hf', 'Ta', 'W', 'Re', 'Os',
'Ir', 'Pt', 'Au', 'Hg', 'Tl', 'Pb', 'Bi', 'Po', 'Fr', 'Ra', 'Ac', 'Th', 'Pa', 'U', 'Np', 'Pu', 'Am', 'Cm',
'Bk', 'Cf', 'Es', 'Fm', 'Md', 'No', 'Lr', 'Rf', 'Db', 'Sg', 'Bh', 'Hs', 'Mt', 'Ds', 'Rg', 'Cn', 'Nh', 'Fl',
'Mc', 'Lv']
metals = ['Li', 'Be', 'Na', 'Mg', 'Al', 'K', 'Ca', 'Sc', 'Ti', 'V', 'Cr',
'Mn', 'Fe', 'Co', 'Ni', 'Cu', 'Zn', 'Ga', 'Rb', 'Sr', 'Y', 'Zr',
'Nb', 'Mo', 'Tc', 'Ru', 'Rh', 'Pd', 'Ag', 'Cd', 'In', 'Sn', 'Cs',
'Ba', 'La', 'Ce', 'Pr', 'Nd', 'Pm', 'Sm', 'Eu', 'Gd', 'Tb', 'Dy',
'Ho', 'Er', 'Tm', 'Yb', 'Lu', 'Hf', 'Ta', 'W', 'Re', 'Os', 'Ir',
'Pt', 'Au', 'Hg', 'Tl', 'Pb', 'Bi', 'Po', 'Fr', 'Ra', 'Ac', 'Th',
'Pa', 'U', 'Np', 'Pu', 'Am', 'Cm', 'Bk', 'Cf', 'Es', 'Fm', 'Md',
'No', 'Lr', 'Rf', 'Db', 'Sg', 'Bh', 'Hs', 'Mt', 'Ds', 'Rg', 'Cn',
'Nh', 'Fl', 'Mc', 'Lv']

avg_bond_lengths = {
'HH': 0.74,
Expand All @@ -42,69 +49,23 @@ def __init__(self, atomic_symbol, x, y, z):
'HI': 1.61
}

atomic_masses = {
'H': 1.01,
'B': 10.81,
'C': 12.01,
'N': 14.01,
'O': 16.0,
'F': 19.0,
'Si': 28.09,
'P': 32.07,
'S': 32.07,
'Cl': 35.45,
'Br': 79.90,
'I': 126.90
}
atomic_masses = {'H': 1.01, 'B': 10.81, 'C': 12.01, 'N': 14.01, 'O': 16.0,
'F': 19.0, 'Si': 28.09, 'P': 32.07, 'S': 32.07, 'Cl': 35.45,
'Br': 79.90, 'I': 126.90}


# Van der Waals radii in Å taken from
# http://www.webelements.com/periodicity/van_der_waals_radius/

# Van der Waals radii in Å taken from http://www.webelements.com/periodicity/van_der_waals_radius/

vdw_radii = {
'H': 1.20,
'He': 1.40,
'Li': 1.82,
'Be': 1.53,
'B': 1.92,
'C': 1.70,
'N': 1.55,
'O': 1.52,
'F': 1.47,
'Ne': 1.54,
'Na': 2.27,
'Mg': 1.73,
'Al': 1.84,
'Si': 2.10,
'P': 1.80,
'S': 1.80,
'Cl': 1.75,
'Ar': 1.88,
'K': 2.75,
'Ca': 2.31,
'Ni': 1.63,
'Cu': 1.40,
'Zn': 1.39,
'Ga': 1.87,
'Ge': 2.11,
'As': 1.85,
'Se': 1.90,
'Br': 1.85,
'Kr': 2.02,
'Rb': 3.03,
'Sr': 2.49,
'Pd': 1.63,
'Ag': 1.72,
'Cd': 1.58,
'In': 1.93,
'Sn': 2.17,
'Sb': 2.06,
'Te': 2.06,
'I': 1.98,
'Xe': 2.16,
'Cs': 3.43,
'Ba': 2.49,
'Pt': 1.75,
'Au': 1.66}
vdw_radii = {'H': 1.20, 'He': 1.40, 'Li': 1.82, 'Be': 1.53, 'B': 1.92,
'C': 1.70, 'N': 1.55, 'O': 1.52, 'F': 1.47, 'Ne': 1.54,
'Na': 2.27, 'Mg': 1.73, 'Al': 1.84, 'Si': 2.10, 'P': 1.80,
'S': 1.80, 'Cl': 1.75, 'Ar': 1.88, 'K': 2.75, 'Ca': 2.31,
'Ni': 1.63, 'Cu': 1.40, 'Zn': 1.39, 'Ga': 1.87, 'Ge': 2.11,
'As': 1.85, 'Se': 1.90, 'Br': 1.85, 'Kr': 2.02, 'Rb': 3.03,
'Sr': 2.49, 'Pd': 1.63, 'Ag': 1.72, 'Cd': 1.58, 'In': 1.93,
'Sn': 2.17, 'Sb': 2.06, 'Te': 2.06, 'I': 1.98, 'Xe': 2.16,
'Cs': 3.43, 'Ba': 2.49, 'Pt': 1.75, 'Au': 1.66}

atoms = ['H', 'He', 'Li', 'Be', 'B', 'C', 'N', 'O', 'F', 'Ne', 'Na', 'Mg', 'Al', 'Si', 'P', 'S', 'Cl', 'Ar',
'K', 'Ca', 'Sc', 'Ti', 'V', 'Cr', 'Mn', 'Fe', 'Co', 'Ni', 'Cu', 'Zn', 'Ga', 'Ge', 'As', 'Se', 'Br',
Expand Down Expand Up @@ -188,15 +149,17 @@ def get_atomic_mass(atom):
if atom.label in atomic_masses.keys():
atom_mass = atomic_masses[atom.label]
else:
logger.warning(f"Couldn't find the atomic mass for {atom.label}. Guessing at 10")
logger.warning(f"Couldn't find the atomic mass for {atom.label}. "
f"Guessing at 10")
atom_mass = 10

return atom_mass


def get_vdw_radii(atom):
"""
Get the van der Wall radii of an atom givne its atomic symbol e.g. Kr --> ~2 Å
Get the van der Wall radii of an atom givne its atomic symbol
e.g. Kr --> ~2 Å
:param atom: (cgbind.atoms.Atom)
:return: (float) VdW radius in Å
Expand All @@ -205,7 +168,8 @@ def get_vdw_radii(atom):
if atom.label in vdw_radii.keys():
vdv_radii = vdw_radii[atom.label]
else:
logger.error(f"Couldn't find the VdV radii for {atom.label}. Guessing at 1.5")
vdv_radii = 1.5
logger.warning(f"Couldn't find the VdV radii for {atom.label}. "
f"Guessing at 2.0")
vdv_radii = 2.0

return vdv_radii
120 changes: 50 additions & 70 deletions cgbind/build.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -47,6 +47,7 @@ def get_fitted_linker_coords_and_cost(linker, template_x_coords, coords_to_fit,

if best_linker_coords is None:
logger.warning('Fitted linker coords could not be found')
best_linker_coords = linker.get_coords()

return best_linker_coords, min_cost

Expand Down Expand Up @@ -92,7 +93,7 @@ def get_kfitted_coords_and_cost(linker, template_x_coords, coords_to_fit, return
return new_linker_coords, cost


def get_linker_atoms_to_add_and_cost(linker, template_linker, curr_coords):
def get_linker_atoms_and_cost(linker, template_linker, current_atoms, x_coords=None):
"""
Get the xyzs of a linker that is fitted to a template_linker object and the associated cost function – i.e. the
repulsion to the current cage structure
Expand All @@ -103,28 +104,28 @@ def get_linker_atoms_to_add_and_cost(linker, template_linker, curr_coords):
:return: list(list)), float
"""

if x_coords is None:
x_coords = linker.get_xmotif_coordinates()

# Ensure the shift amount dr is set
if linker.dr is None:
logger.error('Cannot build a cage dr was None')
return None, 9999999.9

# Copy the linker so no attributes are modified in place..
new_linker = deepcopy(linker)
return linker.atoms, 9999999.9

# Expand the template by an about dr
shifted_coords = get_shifted_template_x_motif_coords(linker_template=template_linker, dr=new_linker.dr)
x_coords = [new_linker.atoms[atom_id].coord for motif in new_linker.x_motifs for atom_id in motif.atom_ids]
shifted_coords = get_shifted_template_x_motif_coords(linker_template=template_linker,
dr=linker.dr)

linker_coords, cost = get_fitted_linker_coords_and_cost(linker=new_linker, template_x_coords=shifted_coords,
coords_to_fit=x_coords, curr_coords=curr_coords)
linker_coords, cost = get_fitted_linker_coords_and_cost(linker=linker,
template_x_coords=shifted_coords,
coords_to_fit=x_coords,
curr_coords=[atom.coord for atom in current_atoms])
logger.info(f'Repulsive + fitting cost for adding the linker is {cost:.5f}')

if linker_coords is None:
logger.warning('Linker coords were None')
return None, cost
atoms = [Atom(linker.atoms[i].label, coord=linker_coords[i])
for i in range(linker.n_atoms)]

new_linker.set_atoms(coords=linker_coords)
return new_linker.atoms, cost
return atoms, cost


def cost_fitted_x_motifs(dr, linker, linker_template, x_coords):
Expand All @@ -145,36 +146,6 @@ def cost_fitted_x_motifs(dr, linker, linker_template, x_coords):
return cost


def calc_linker_cost(init_conformer, linker, x_motifs, template_linker):
"""
For a set of linker xyzs e.g. one conformer and a linker object together with a list of x motifs in the
structure return cost function associated with fitting the X motifs to the template
:param conformer: (cgbind.molecule.BaseStruct)
:param linker: (Linker)
:param x_motifs: (list(Xmotif))
:param template_linker: (Template.Linker)
:return: (tuple(Linker, float)) New linker and cost
"""
coords = init_conformer.get_coords()
x_coords = [coords[atom_id] for motif in x_motifs for atom_id in motif.atom_ids]

# Minimise the cost function as a function of dr in Å
min_result = minimize(cost_fitted_x_motifs, x0=np.array([1.0]), args=(linker, template_linker, x_coords),
method='L-BFGS-B', tol=1e-3, bounds=Bounds(-100.0, 10.0))

conformer = BaseStruct()
conformer.dr = min_result.x[0]
conformer.x_motifs = x_motifs
conformer.set_atoms(init_conformer.atoms)

x_motif_atoms = [atom_id for motif in x_motifs for atom_id in motif.atom_ids]
conformer.x_atoms = [atom_id for atom_id in linker.x_atoms if atom_id in x_motif_atoms]
conformer.cost = min_result.fun

return conformer


def build_homoleptic_cage(cage, max_cost):
"""
Construct the geometry (atoms) of a homoleptic cage
Expand All @@ -189,74 +160,83 @@ def build_homoleptic_cage(cage, max_cost):
logger.info(f'Have {len(cage.linkers[0].possibilities)} linkers to fit')

min_cost, best_linker = 99999999.9, None
atoms = []

# For all the possible linker conformer / Xmotif set possibilities
for linker in cage.linkers[0].possibilities:
cage_cost = 0.0

for i, template_linker in enumerate(cage.cage_template.linkers):
linker_atoms, cost = get_linker_atoms_to_add_and_cost(linker, template_linker,
curr_coords=[atom.coord for atom in atoms])
cage_cost += cost
# Atoms for and cost in building this cage
atoms, cage_cost = [], 0.0

if linker_atoms is None:
logger.error('Failed to get linker')
break
# Coordinates of the X motif atoms in this linker - used to rotate
x_coords = linker.get_xmotif_coordinates()

for i, template_linker in enumerate(cage.cage_template.linkers):
linker_atoms, cost = get_linker_atoms_and_cost(linker,
template_linker,
atoms,
x_coords)
cage_cost += cost
atoms += linker_atoms

if cage_cost < min_cost:
min_cost = cage_cost
best_linker = deepcopy(linker)

if cage_cost < max_cost:
logger.info(f'Total L-L repulsion + fit to template in building cage is {cage_cost:.2f}')
logger.info(f'Total L-L repulsion + fit to template in building '
f'cage is {cage_cost:.2f}')
break

else:
atoms = []

if len(atoms) == 0:
# If there is no break due to a small repulsion then build the best
# possible cage
if cage_cost > max_cost:
if best_linker is None:
logger.error('Could not achieve the required cost threshold for building the cage')
logger.error('Could not achieve the required cost threshold for '
'building the cage')
return None
else:
logger.warning('Failed to reach the threshold. Returning the cage that minimises the L-L repulsion')
logger.warning('Failed to reach the threshold. Returning the cage '
'that minimises the L-L repulsion')
atoms = []
for i, template_linker in enumerate(cage.cage_template.linkers):
linker_atoms, _ = get_linker_atoms_and_cost(best_linker,
template_linker,
atoms)
atoms += linker_atoms

# Set the delta r for the whole cage
cage.dr = best_linker.dr
atoms = []
for i, template_linker in enumerate(cage.cage_template.linkers):
linker_atoms, _ = get_linker_atoms_to_add_and_cost(best_linker, template_linker,
curr_coords=[atom.coord for atom in atoms])
atoms += linker_atoms

# Add the metals from the template shifted by dr
for metal in cage.cage_template.metals:
metal_coord = cage.dr * metal.shift_vec / np.linalg.norm(metal.shift_vec) + metal.coord
atoms.append(Atom(cage.metal, x=metal_coord[0], y=metal_coord[1], z=metal_coord[2]))
atoms.append(Atom(cage.metal, coord=metal_coord))

cage.set_atoms(atoms)
return None


def build_heteroleptic_cage(cage, max_cost):
logger.info('Building a heteroleptic cage')
logger.warning('Due to the very large space that needs to be minimised only the *best* linker conformer is used')
logger.warning('Due to the very large space that needs to be minimised '
'only the *best* linker conformer is used')

added_linkers, atoms = [], []

for i, linker in enumerate(cage.linkers):

linker.set_ranked_linker_possibilities(metal=cage.metal)

linker_atoms, cost = get_linker_atoms_to_add_and_cost(linker.possibilities[0], cage.cage_template.linkers[i],
curr_coords=[atom.coord for atom in atoms])
linker_atoms, cost = get_linker_atoms_and_cost(linker.possibilities[0],
cage.cage_template.linkers[i],
atoms)

logger.info(f'L-L repulsion + fit to template in building cage is {cost:.2f}')
atoms += linker_atoms
linker.dr = linker.possibilities[0].dr

logger.warning('Heteroleptic cages will have the average dr of all linkers - using the average')
logger.warning('Heteroleptic cages will have the average dr of all linkers '
'- using the average')
cage.dr = np.average(np.array([linker.dr for linker in cage.linkers]))

# Add the metals from the template shifted by dr
Expand Down
Loading

0 comments on commit d5c33fc

Please sign in to comment.