diff --git a/architector/io_calc.py b/architector/io_calc.py index 5d46c0b..a4d3463 100644 --- a/architector/io_calc.py +++ b/architector/io_calc.py @@ -18,7 +18,7 @@ import ase from ase.io import Trajectory from ase.optimize import LBFGSLineSearch -from ase.constraints import (FixAtoms, FixBondLengths) +from ase.constraints import (FixAtoms, FixBondLengths, FixInternals) has_sella = True try: from sella import Internals @@ -91,7 +91,7 @@ def __init__(self, structure, parameters={}, init_sanity_check=False, detect_spin_charge=False, fix_m_neighbors=False, fix_indices=None, default_params=params, ase_opt_method=None, ase_opt_kwargs={}, species_run=False, intermediate=False, skip_spin_assign=False, save_trajectories=False, - store_results=False, use_constraints=False, + store_results=False, use_constraints=False, trans_oxo_triples=[], calculator=None, debug=False): """CalcExecutor is the class handling all calculations of full metal-ligand complexes. @@ -165,6 +165,8 @@ def __init__(self, structure, parameters={}, init_sanity_check=False, elif isinstance(self.in_struct, architector.io_molecule.Molecule): if len(self.in_struct.ase_atoms.constraints) > 0: self.mol.ase_atoms.set_constraint(self.in_struct.ase_atoms.constraints) + else: + self.mol.ase_atoms.constraints = [] self.method = method default_params = params.copy() default_params['debug'] = debug @@ -196,8 +198,8 @@ def __init__(self, structure, parameters={}, init_sanity_check=False, self.force_generation = False self.force_oxo_relax = False self.replace_trics = False + self.trans_oxo_triples = trans_oxo_triples self.results = None - self.detect_spin_charge = detect_spin_charge if len(parameters) > 0: for key,val in parameters.items(): @@ -304,8 +306,10 @@ def calculate(self): # verbosity=0) # Difference of more than 1. Still perform a ff_preoptimization if requested. if (np.abs(self.mol.xtb_charge - self.mol.charge) > 1): - if ((not self.override_oxo_opt) or (self.assembly)) and (not self.force_oxo_relax): - self.relax = False # E.g - don't relax if way off in oxdiation states (III) vs (V or VI) + if len(self.trans_oxo_triples) > 0: + pass + elif ((not self.override_oxo_opt) or (self.assembly)) and (not self.force_oxo_relax): + self.relax = False # E.g - don't relax if way off in oxidation states (III) vs (V or VI) elif self.assembly: # FF more stable for highly charged assembly complexes. self.method = 'GFN-FF' uhf_vect = np.zeros(len(self.mol.ase_atoms)) @@ -323,25 +327,44 @@ def calculate(self): if not obabel_ff_requested: self.mol.ase_atoms.calc = calc if self.relax: - if self.parameters.get("freeze_molecule_add_species",False) and self.species_run: + cs = self.mol.ase_atoms.constraints + if self.parameters.get("freeze_molecule_add_species", + False) and self.species_run: if self.parameters['debug']: print('Fixing first component!') fix_inds = self.mol.find_component_indices(component=0) c = FixAtoms(indices=fix_inds.tolist()) - self.mol.ase_atoms.set_constraint(c) - elif isinstance(self.fix_indices,list): + cs.append(c) + elif isinstance(self.fix_indices, list): c = FixAtoms(indices=self.fix_indices) - self.mol.ase_atoms.set_constraint(c) + cs.append(c) if len(self.mol.ase_constraints) > 0: - fix_list = [[x[0],x[1]] for x in self.mol.ase_constraints] + fix_list = [[x[0], x[1]] for x in self.mol.ase_constraints] c = FixBondLengths(fix_list) - self.mol.ase_atoms.set_constraint(c) + cs.append(c) + if len(self.trans_oxo_triples) > 0: + bonds = [] + angles_degs = [] + for triple in self.trans_oxo_triples: + if not isinstance(self.fix_indices,list): + tlist = [] + else: + tlist = self.fix_indices + if (triple[0] not in tlist) or (triple[1] not in tlist): + bonds.append([self.mol.ase_atoms.get_distance(triple[0],triple[1]), + [triple[0],triple[1]]]) + bonds.append([self.mol.ase_atoms.get_distance(triple[2],triple[1]), + [triple[2],triple[1]]]) + angles_degs.append([180.0,list(triple)]) + c = FixInternals(bonds=bonds, angles_deg=angles_degs) + cs.append(c) + self.mol.ase_atoms.set_constraint(cs) with arch_context_manage.make_temp_directory( prefix=self.parameters['temp_prefix']) as _: if has_sella and self.ase_opt_kwargs.get('sella_internal_trics',False): if self.debug: print('Adding Internals....') - ints = Internals(self.mol.ase_atoms,allow_fragments=True) + ints = Internals(self.mol.ase_atoms, allow_fragments=True) del self.ase_opt_kwargs['sella_internal_trics'] self.replace_trics = True ints.find_all_bonds() @@ -352,12 +375,12 @@ def calculate(self): self.init_energy = copy.deepcopy(self.mol.ase_atoms.get_total_energy()) if self.parameters['save_trajectories']: if self.logfile is not None: - dyn = self.opt_method(self.mol.ase_atoms, + dyn = self.opt_method(self.mol.ase_atoms, trajectory='temp.traj', logfile=self.logfile, **self.ase_opt_kwargs) else: - dyn = self.opt_method(self.mol.ase_atoms, + dyn = self.opt_method(self.mol.ase_atoms, trajectory='temp.traj', **self.ase_opt_kwargs) else: @@ -389,10 +412,11 @@ def calculate(self): self.init_energy = 10000 self.calc_time = time.time() - self.calc_time # Remove constraint - if self.parameters.get("freeze_molecule_add_species",False) and (self.species_run): + if self.parameters.get("freeze_molecule_add_species", + False) and (self.species_run): if self.parameters['debug']: print('Removing fixing first component!') - self.mol.ase_atoms.set_constraint() + self.mol.ase_atoms.set_constraint() else: with arch_context_manage.make_temp_directory( prefix=self.parameters['temp_prefix']) as _: @@ -419,9 +443,11 @@ def calculate(self): try: self.init_energy = io_obabel.obmol_energy(self.mol) out_atoms, energy = io_obabel.obmol_opt(self.mol, center_metal=True, - fix_m_neighbors=self.fix_m_neighbors,fix_indices=self.fix_indices, # Note - fixing metal neighbors in UFF + fix_m_neighbors=self.fix_m_neighbors, + fix_indices=self.fix_indices, # Note - fixing metal neighbors in UFF # Done to maintain metal center symmetry - return_energy=True) + trans_oxo_triples=self.trans_oxo_triples, + return_energy=True) self.successful = True self.energy = energy self.mol.ase_atoms.set_positions(out_atoms.get_positions()) diff --git a/architector/io_molecule.py b/architector/io_molecule.py index f998c6f..75f75ff 100644 --- a/architector/io_molecule.py +++ b/architector/io_molecule.py @@ -1585,6 +1585,36 @@ def get_disjoint_molecules(self): mols.append(newmol) return mols + def detect_trans_oxos(self): + """Detect and flag trans_oxos. + + Returns: + trans_oxo_triples: list(tuple (3)) + List of tuples of indices of oxo1-metal-oxo2. + """ + if len(self.graph) < 1: + self.create_mol_graph() + info_dict = self.split_ligs() + oxo_inds = [] + for i,lig in enumerate(info_dict['lig_smiles']): + if lig == '[O-2]': + oxo_inds.append(i) + trans_oxo_triples = [] + if isinstance(info_dict['metal_ind'],(int,float)): + info_dict['metal_ind'] = [info_dict['metal_ind']] + for m_ind in info_dict['metal_ind']: + toxinds = [] + for oxind in oxo_inds: + if m_ind in info_dict['bound_metal_inds'][oxind]: + toxinds.append(oxind) + if len(toxinds) > 1: + for pair in itertools.combinations(toxinds,2): + ind0 = info_dict['original_lig_inds'][pair[0]][0] + ind1 = info_dict['original_lig_inds'][pair[1]][0] + angle = self.ase_atoms.get_angle(ind0,m_ind,ind1) + if angle > 179: + trans_oxo_triples.append((ind0,m_ind,ind1)) + return trans_oxo_triples def functionalize_3D(self, functional_groups=['carboxyllic_acid'], @@ -1643,7 +1673,6 @@ def functionalize_3D(self, Keep the existing structure frozen, by default True """ from architector.io_calc import CalcExecutor - new_functional_groups = [] for i,fg in enumerate(functional_groups): if fg in io_ptable.functional_groups_dict: @@ -1754,10 +1783,22 @@ def functionalize_3D(self, self.atom_types += funct_mol.atom_types if funct_mol.charge is not None: self.charge += funct_mol.charge # Add charge for charged functional groups. - if isinstance(self.xtb_charge,int): - self.xtb_charge += funct_mol.charge + if isinstance(self.xtb_charge, (int,float)): + if isinstance(funct_mol.xtb_charge, (int,float)): + self.xtb_charge += funct_mol.xtb_charge + else: + self.xtb_charge += funct_mol.charge + else: + if isinstance(funct_mol.xtb_charge, (int,float)): + self.xtb_charge = funct_mol.xtb_charge + if funct_mol.uhf is not None: + self.uhf += funct_mol.uhf + if isinstance(self.xtb_uhf, (int,float)): + if isinstance(funct_mol.xtb_uhf, (int,float)): + self.xtb_uhf += funct_mol.xtb_uhf else: - self.xtb_charge = funct_mol.charge + if isinstance(funct_mol.xtb_uhf, (int,float)): + self.xtb_uhf = funct_mol.xtb_uhf self.BO_dict.update(newligbodict) self.create_graph_from_bo_dict() elif isinstance(idx,(np.ndarray,list)): @@ -1876,10 +1917,22 @@ def functionalize_3D(self, self.atom_types += funct_mol.atom_types if funct_mol.charge is not None: self.charge += funct_mol.charge # Add charge for charged functional groups. - if isinstance(self.xtb_charge,int): - self.xtb_charge += funct_mol.charge + if isinstance(self.xtb_charge, (int,float)): + if isinstance(funct_mol.xtb_charge, (int,float)): + self.xtb_charge += funct_mol.xtb_charge + else: + self.xtb_charge += funct_mol.charge + else: + if isinstance(funct_mol.xtb_charge, (int,float)): + self.xtb_charge = funct_mol.xtb_charge + if funct_mol.uhf is not None: + self.uhf += funct_mol.uhf + if isinstance(self.xtb_uhf, (int,float)): + if isinstance(funct_mol.xtb_uhf, (int,float)): + self.xtb_uhf += funct_mol.xtb_uhf else: - self.xtb_charge = funct_mol.charge + if isinstance(funct_mol.xtb_uhf, (int,float)): + self.xtb_uhf = funct_mol.xtb_uhf self.BO_dict.update(newligbodict) self.create_graph_from_bo_dict() else: @@ -1890,13 +1943,16 @@ def functionalize_3D(self, freeze_indices = None if uff_opt: tmpmol = self.write_mol2('temp.mol2',writestring=True) - obmol_opt = CalcExecutor(tmpmol, relax=True, method='UFF', fix_indices=freeze_indices) + obmol_opt = CalcExecutor(tmpmol, relax=True, method='UFF', fix_indices=freeze_indices, + trans_oxo_triples=self.detect_trans_oxos()) self.ase_atoms.set_positions(obmol_opt.mol.ase_atoms.get_positions()) if xtb_opt: tmpmol = self.write_mol2('temp.mol2',writestring=True) obmol_opt = CalcExecutor(tmpmol, relax=True, method='GFN2-xTB', - fix_indices=freeze_indices) + fix_indices=freeze_indices, + trans_oxo_triples=self.detect_trans_oxos() + ) self.ase_atoms.set_positions(obmol_opt.mol.ase_atoms.get_positions()) diff --git a/architector/io_obabel.py b/architector/io_obabel.py index c964225..60aba8d 100644 --- a/architector/io_obabel.py +++ b/architector/io_obabel.py @@ -617,7 +617,7 @@ def convert_ase_obmol(ase_atoms): def obmol_opt(structure, center_metal=False, fix_m_neighbors=True, - return_energy=False, fix_indices=None): + return_energy=False, fix_indices=None, trans_oxo_triples=[]): """obmol_opt take in a structure and optimize with openbabel return as ase atoms as default Will default to MMFF94 if it is applicable - otherwise it is UFF. @@ -634,6 +634,8 @@ def obmol_opt(structure, center_metal=False, fix_m_neighbors=True, Return the energy of the optimized structure?, default False fix_indices : list, optional Fix the atoms at these indices (0-indexed), defualt None + trans_oxo_triples : list (tuples), optional, + Fix the angle of these indices (0-indexed) and bond lengths, defualt None Returns ---------- @@ -646,9 +648,9 @@ def obmol_opt(structure, center_metal=False, fix_m_neighbors=True, OBMol = convert_ase_obmol(structure) elif isinstance(structure,str): if 'TRIPOS' in structure: - OBMol = convert_mol2_obmol(structure,readstring=True) + OBMol = convert_mol2_obmol(structure, readstring=True) elif structure[-5:] == '.mol2': - OBMol = convert_mol2_obmol(structure,readstring=False) + OBMol = convert_mol2_obmol(structure, readstring=False) elif isinstance(structure,architector.io_molecule.Molecule): mol2str = structure.write_mol2('cool.mol2', writestring=True) OBMol = convert_mol2_obmol(mol2str, readstring=True) @@ -661,11 +663,20 @@ def obmol_opt(structure, center_metal=False, fix_m_neighbors=True, else: FF = ob.OBForceField.FindForceField('UFF') - if isinstance(fix_indices,list): + if isinstance(fix_indices, list): constr = ob.OBFFConstraints() for j in fix_indices: constr.AddAtomConstraint(int(j+1)) - FF.Setup(OBMol,constr) + if len(trans_oxo_triples) > 0: + for triple in trans_oxo_triples: + if (triple[0] not in fix_indices) or (triple[2] not in fix_indices): + constr.AddAngleConstraint(int(triple[0])+1,int(triple[1])+1,int(triple[2])+1,180.0) + at1 = OBMol.GetAtom(int(triple[1])+1) + dist0 = at1.GetDistance(int(triple[0])+1) + constr.AddDistanceConstraint(int(triple[0])+1,int(triple[1])+1,dist0) + dist2 = at1.GetDistance(int(triple[2])+1) + constr.AddDistanceConstraint(int(triple[2])+1,int(triple[1])+1,dist2) + FF.Setup(OBMol, constr) elif fix_m_neighbors: _,anums,graph = get_OBMol_coords_anums_graph(OBMol, return_coords=False, get_types=False) syms = [io_ptable.elements[x] for x in anums] @@ -681,7 +692,17 @@ def obmol_opt(structure, center_metal=False, fix_m_neighbors=True, elif len(mets) == 0: constr = ob.OBFFConstraints() # print('No Metals present for FF optimization.') - FF.Setup(OBMol,constr) + FF.Setup(OBMol, constr) + elif len(trans_oxo_triples) > 0: + constr = ob.OBFFConstraints() + for triple in trans_oxo_triples: + constr.AddAngleConstraint(int(triple[0])+1,int(triple[1])+1,int(triple[2])+1,180.0) + at1 = OBMol.GetAtom(int(triple[1])+1) + dist0 = at1.GetDistance(int(triple[0])+1) + constr.AddDistanceConstraint(int(triple[0])+1,int(triple[1])+1,dist0) + dist2 = at1.GetDistance(int(triple[2])+1) + constr.AddDistanceConstraint(int(triple[2])+1,int(triple[1])+1,dist2) + FF.Setup(OBMol, constr) else: FF.Setup(OBMol) @@ -1019,7 +1040,7 @@ def obmol_lig_split(mol2string, # Key block for catching where coordinating atoms were deprotonated ### WORKING -> Does not work great for nitrogen compounds. for l, atom in enumerate(ob.OBMolAtomIter(ligobmol)): - if l in coord_atom_list: + if (l in coord_atom_list) and (len(lig) > 1): total_val = (io_ptable.valence_electrons[atom.GetAtomicNum()] + atom.GetTotalValence()) close = np.argmin(np.abs(np.array(io_ptable.filled_valence_electrons)-total_val)) if atom.GetFormalCharge() > 0: # Coordinating atoms rarely +ve @@ -1058,6 +1079,15 @@ def obmol_lig_split(mol2string, elif np.abs(total_val - 16) < 2: # octet breaker near 16 newcharge = int(atom.GetFormalCharge()-(16-total_val)) atom.SetFormalCharge(newcharge) + elif (len(lig) == 1): # Single atom ligand. + atom.SetImplicitHCount(0) # No implicit hydrogens + total_val = (io_ptable.valence_electrons[atom.GetAtomicNum()] + atom.GetTotalValence()) + close = np.argmin(np.abs(np.array(io_ptable.filled_valence_electrons)-total_val)) + if atom.GetFormalCharge() > 0: # Coordinating atoms rarely +ve + atom.SetFormalCharge(0) + else: + newcharge = int(atom.GetFormalCharge()-(io_ptable.filled_valence_electrons[close]-total_val)) + atom.SetFormalCharge(newcharge) if (not allow_radicals) and (ligobmol.GetTotalSpinMultiplicity() > 1): for l, atom in enumerate(ob.OBMolAtomIter(ligobmol)): if l not in coord_atom_list: @@ -1108,9 +1138,12 @@ def obmol_lig_split(mol2string, return ligand_smiles,coord_atom_lists else: info_dict = dict() - if len(met_inds) > 0: + if len(met_inds) == 1: info_dict['metal'] = io_ptable.elements[anums[met_inds[0]]] info_dict['metal_ind'] = met_inds[0] + elif len(met_inds) > 1: + info_dict['metal'] = [io_ptable.elements[anums[x]] for x in met_inds] + info_dict['metal_ind'] = met_inds else: info_dict['metal'] = None info_dict['metal_ind'] = None @@ -1127,14 +1160,23 @@ def obmol_lig_split(mol2string, info_dict['lig_coord_ats'] = lig_coord_ats info_dict['original_lig_inds'] = ligs_inds info_dict['mapped_smiles_inds'] = mapped_smiles_inds + bound_metal_inds = [] + for lig in ligs_inds: + inds = np.array(lig) + m_inds = [] + for m_ind in met_inds: + if np.any(graph[inds, m_ind] == 1): + m_inds.append(m_ind) + bound_metal_inds.append(m_inds) + info_dict['bound_metal_inds'] = bound_metal_inds elif calc_coord_atoms: - for i,lig_obmol in enumerate(lig_obmols): - _,anums,_ = get_OBMol_coords_anums_graph(lig_obmol,get_types=False) + for i, lig_obmol in enumerate(lig_obmols): + _, anums, _ = get_OBMol_coords_anums_graph(lig_obmol, get_types=False) lig_coord_ats.append(','.join([io_ptable.elements[x] for x in np.array(anums)[np.array(coord_atom_lists[i])]])) info_dict['lig_coord_ats'] = lig_coord_ats else: info_dict['lig_coord_ats'] = None - return ligand_smiles,coord_atom_lists, info_dict + return ligand_smiles, coord_atom_lists, info_dict def get_canonical_label(obmol): diff --git a/development/Secondary_Solvation_Shell.ipynb b/development/Secondary_Solvation_Shell.ipynb deleted file mode 100644 index 36a9a7a..0000000 --- a/development/Secondary_Solvation_Shell.ipynb +++ /dev/null @@ -1,252 +0,0 @@ -{ - "cells": [ - { - "cell_type": "markdown", - "id": "4302c967", - "metadata": {}, - "source": [ - "# 7 - Secondary Solvation Shell\n", - "\n", - "Beyond simple ligands searched from SMILES, we want to dig into functionalizations of ligands!\n", - "\n", - "In this tutorial we will look at the ligand functionalization routines in Architector, covering:\n", - "\n", - "**(A)** Viewing default functional groups present in architector by name!\n", - "\n", - "**(B)** Identifying functionalization sites on ligands.\n", - "\n", - "**(C)** Adding both single and multiple functionalizations to a single ligand." - ] - }, - { - "cell_type": "code", - "execution_count": 1, - "id": "2a424827", - "metadata": {}, - "outputs": [], - "source": [ - "from architector import (build_complex,\n", - " view_structures)" - ] - }, - { - "cell_type": "code", - "execution_count": 2, - "id": "d3df2a8a", - "metadata": {}, - "outputs": [], - "source": [ - "inputDict = {'core':{'metal':'Fe','coreType':'octahedral'},\n", - " 'ligands':[{'smiles':'O','coordList':[0]}],\n", - " 'parameters':dict()}" - ] - }, - { - "cell_type": "code", - "execution_count": 3, - "id": "693d1541", - "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Total valid symmetries for core octahedral: 1\n", - "GENERATING CONFORMATIONS for O\n", - "CONFORMERS GENERATED for O\n", - "ASSEMBLING COMPLEX\n", - "LIGAND: O\n", - "FINDING CORRECT CONFORMER\n", - "LIGAND: O\n", - "FINDING CORRECT CONFORMER\n", - "LIGAND: O\n", - "FINDING CORRECT CONFORMER\n", - "LIGAND: O\n", - "FINDING CORRECT CONFORMER\n", - "LIGAND: O\n", - "FINDING CORRECT CONFORMER\n", - "LIGAND: O\n", - "FINDING CORRECT CONFORMER\n", - "Initial Sanity: True\n", - "Complex sanity after adding ligands: True\n", - "Complex class generated: True\n", - "OPTIMIZING MOLECULE\n", - " Step[ FC] Time Energy fmax\n", - "*Force-consistent energies used in optimization.\n", - "BFGSLineSearch: 0[ 0] 16:37:24 -867.170038* 4.4804\n", - "BFGSLineSearch: 1[ 2] 16:37:24 -868.184000* 1.3493\n", - "BFGSLineSearch: 2[ 3] 16:37:24 -869.721016* 1.7982\n", - "BFGSLineSearch: 3[ 4] 16:37:24 -870.628923* 3.0902\n", - "BFGSLineSearch: 4[ 6] 16:37:24 -871.176430* 0.6051\n", - "BFGSLineSearch: 5[ 8] 16:37:24 -871.201601* 0.5484\n", - "BFGSLineSearch: 6[ 9] 16:37:25 -871.269653* 0.3177\n", - "BFGSLineSearch: 7[ 11] 16:37:25 -871.272075* 0.1762\n", - "BFGSLineSearch: 8[ 12] 16:37:25 -871.274558* 0.0569\n", - "ComplexSanity: True\n" - ] - } - ], - "source": [ - "out = build_complex(inputDict)" - ] - }, - { - "cell_type": "code", - "execution_count": 4, - "id": "05e712e2", - "metadata": {}, - "outputs": [], - "source": [ - "start = out[list(out.keys())[0]]['mol2string']" - ] - }, - { - "cell_type": "code", - "execution_count": 5, - "id": "c107170d", - "metadata": {}, - "outputs": [ - { - "data": { - "application/3dmoljs_load.v0": "
\n

You appear to be running in JupyterLab (or JavaScript failed to load for some other reason). You need to install the 3dmol extension:
\n jupyter labextension install jupyterlab_3dmol

\n
\n", - "text/html": [ - "
\n", - "

You appear to be running in JupyterLab (or JavaScript failed to load for some other reason). You need to install the 3dmol extension:
\n", - " jupyter labextension install jupyterlab_3dmol

\n", - "
\n", - "" - ] - }, - "metadata": {}, - "output_type": "display_data" - } - ], - "source": [ - "view_structures(start)" - ] - }, - { - "cell_type": "code", - "execution_count": 6, - "id": "b4a73928", - "metadata": {}, - "outputs": [], - "source": [ - "from architector.io_crest import add_explicit_solvents" - ] - }, - { - "cell_type": "code", - "execution_count": 7, - "id": "1ce1b935", - "metadata": {}, - "outputs": [ - { - "name": "stderr", - "output_type": "stream", - "text": [ - "Note: The following floating-point exceptions are signalling: IEEE_DENORMAL\n", - "ERROR STOP \n", - "\n", - "Error termination. Backtrace:\n", - "#0 0x106a114dc\n", - "#1 0x106a12395\n", - "#2 0x106a13b04\n", - "#3 0x10601e7e9\n", - "#4 0x105fccac0\n", - "#5 0x105fce65e\n", - "#6 0x106125c9f\n", - "#7 0x1061260fe\n" - ] - }, - { - "ename": "CalledProcessError", - "evalue": "Command 'crest solute.xyz --qcg solvent.xyz --chrg 3 --uhf 5 --nsolv 6 --T 12 --ensemble --alpb water --mdtime 50 --mddump 200 > crest.out' returned non-zero exit status 1.", - "output_type": "error", - "traceback": [ - "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", - "\u001b[0;31mCalledProcessError\u001b[0m Traceback (most recent call last)", - "\u001b[0;32m/var/folders/hm/58cd4wn126gg0_gbqn15c21800144l/T/ipykernel_40018/2338525459.py\u001b[0m in \u001b[0;36m\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mstructs\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0menergies\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0madd_explicit_solvents\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mstart\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m", - "\u001b[0;32m~/software/architector/architector/io_crest.py\u001b[0m in \u001b[0;36madd_explicit_solvents\u001b[0;34m(complex_mol2string, n_solvents, solvent)\u001b[0m\n\u001b[1;32m 239\u001b[0m \u001b[0mfile1\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mwrite\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmolcontrol_str\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 240\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 241\u001b[0;31m \u001b[0msub\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mrun\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mexecStr\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mshell\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mTrue\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mcheck\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mTrue\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 242\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 243\u001b[0m \u001b[0mconformerList\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mxtb_energies\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mread_conformers\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"final_ensemble.xyz\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", - "\u001b[0;32m/opt/anaconda3/envs/architector/lib/python3.9/subprocess.py\u001b[0m in \u001b[0;36mrun\u001b[0;34m(input, capture_output, timeout, check, *popenargs, **kwargs)\u001b[0m\n\u001b[1;32m 526\u001b[0m \u001b[0mretcode\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mprocess\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mpoll\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 527\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0mcheck\u001b[0m \u001b[0;32mand\u001b[0m \u001b[0mretcode\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 528\u001b[0;31m raise CalledProcessError(retcode, process.args,\n\u001b[0m\u001b[1;32m 529\u001b[0m output=stdout, stderr=stderr)\n\u001b[1;32m 530\u001b[0m \u001b[0;32mreturn\u001b[0m \u001b[0mCompletedProcess\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mprocess\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mretcode\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mstdout\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mstderr\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", - "\u001b[0;31mCalledProcessError\u001b[0m: Command 'crest solute.xyz --qcg solvent.xyz --chrg 3 --uhf 5 --nsolv 6 --T 12 --ensemble --alpb water --mdtime 50 --mddump 200 > crest.out' returned non-zero exit status 1." - ] - } - ], - "source": [ - "structs,energies = add_explicit_solvents(start)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "4942d607", - "metadata": {}, - "outputs": [], - "source": [ - "ls" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "e4904017", - "metadata": {}, - "outputs": [], - "source": [] - } - ], - "metadata": { - "kernelspec": { - "display_name": "Python 3 (ipykernel)", - "language": "python", - "name": "python3" - }, - "language_info": { - "codemirror_mode": { - "name": "ipython", - "version": 3 - }, - "file_extension": ".py", - "mimetype": "text/x-python", - "name": "python", - "nbconvert_exporter": "python", - "pygments_lexer": "ipython3", - "version": "3.9.7" - } - }, - "nbformat": 4, - "nbformat_minor": 5 -} diff --git a/development/n_lig_combos.py b/development/n_lig_combos.py deleted file mode 100644 index 30c1f15..0000000 --- a/development/n_lig_combos.py +++ /dev/null @@ -1,132 +0,0 @@ - # def assemble_random_complex(self): - # """Assemble a complex with random combinations of ligand conformers.""" - - # # Variables - # constraintList = [] - - # mol = io_molecule.Molecule() # Initialize molecule. - # mol.load_ase(io_obabel.smiles2Atoms(self.coreSmiles),atom_types=[self.complexMol[0].symbol]) - # init_sane = True - # # Iterate over ligands - # for i,ligand in enumerate(self.ligandList): - # if init_sane: - # # Find correct conformer to use - # conf_list = ligand.conformerList - # cc = np.random.choice(np.arange(len(conf_list)),1)[0] - # conformer = conf_list[cc] - # mol.append_ligand({'ase_atoms':conformer,'bo_dict':ligand.BO_dict, - # 'atom_types':ligand.atom_types}) - # sane = mol.graph_sanity_check(params=self.parameters,assembly=True) - # if not sane: - # init_sane = False - # else: - # if self.parameters['debug']: - # print('Random conformer failed on ligand {}!'.format(self.ligandList[i-1].smiles)) - # if self.parameters['debug']: - # print('Initial Sanity: ', init_sane) - # init_mol2str = mol.write_mol2('Charge: {} Unpaired_Electrons: {} XTB_Unpaired_Electrons: {} .mol2'.format( - # int(sum(self.complexMol.get_initial_charges())), - # int(io_xtb_calc.calc_suggested_spin(self.complexMol,parameters=self.parameters)), - # int(sum(self.complexMol.get_initial_magnetic_moments()))), writestring=True) - # energy = np.inf - # init_energy = np.inf - # mol2str = None - # if init_sane: - # # Add contraints for core - Removed b/c unnecessary - # fixCore = ase_con.FixAtoms(indices=[0]) - # constraintList.append(fixCore) - # mol.ase_atoms.set_constraint(constraintList) # Fix core metal only. - # sane = io_molecule.sanity_check(mol.ase_atoms,params=self.parameters,assembly=True) - # if self.parameters['debug']: - # print('Complex sanity after adding ligands: ',sane) - # if sane: - # _ = mol.ase_atoms.copy() - # # Add the calculator. - # io_xtb_calc.set_XTB_calc(mol.ase_atoms, self.parameters, assembly=False) - # try: - # with arch_context_manage.make_temp_directory(prefix=self.parameters['temp_prefix']) as _: - # init_energy = mol.ase_atoms.get_total_energy() - # except Exception as e: - # if self.parameters['debug']: - # print(e) - # sane = False - # if sane: - # try: - # with arch_context_manage.make_temp_directory(prefix=self.parameters['temp_prefix']) as _: - # if self.parameters['save_trajectories']: - # dyn = BFGSLineSearch(mol.ase_atoms, maxstep=3000, trajectory='temp.traj') - # else: - # dyn = BFGSLineSearch(mol.ase_atoms, maxstep=3000) - # dyn.run(fmax=0.1) - # if self.parameters['save_trajectories']: - # traj = Trajectory('temp.traj') - # newposits = mol.ase_atoms.get_positions() - mol.ase_atoms.get_positions()[0] # Recenter - # mol.ase_atoms.set_positions(newposits) - # sane1 = io_molecule.sanity_check(mol.ase_atoms, params=self.parameters) # Check close atoms or separated atoms - # mol.ase_atoms.set_positions(mol.ase_atoms.get_positions()) # Set to final configuration - # sane2 = mol.graph_sanity_check(params=self.parameters, assembly=False) # Check for graph-based elongated bonds - # if sane1 and sane2: - # sane = True - # else: - # sane = False - # energy = mol.ase_atoms.get_total_energy() - # except Exception as e: # Catch final relaxation failures. - # if self.parameters['debug']: - # print(e) - # print('Failed final relaxation.') - # sane = False - # energy = 100000000 - # if self.parameters['save_trajectories']: - # for i,ats in enumerate(traj): - # self.parameters['ase_db'].write(ats,geo_step=i) - # mol2str = mol.write_mol2('Charge: {} Unpaired_Electrons: {} XTB_Unpaired_Electrons: {} .mol2'.format( - # int(sum(self.complexMol.get_initial_charges())), - # int(io_xtb_calc.calc_suggested_spin(self.complexMol,parameters=self.parameters)), - # int(sum(self.complexMol.get_initial_magnetic_moments()))), writestring=True) - # return mol,mol2str,energy,init_mol2str,init_energy,sane - # else: - # return mol,mol2str,energy,init_mol2str,init_energy,sane - - #### Random ligand conformer generation routine. Potentially useful but not in production. - # if inputDict['parameters']['n_lig_combos']>1: - # mol2strs = [] - # mols = [] - # xtb_energies1 = [] - # init_mol2strs = [] - # init_xtb_energies = [] - # for k in range(2*inputDict['parameters']['n_lig_combos']): - # mol,mol2str,energy,init_mol2str,init_energy,sane = structs[i].assemble_random_complex() - # if sane: - # mols.append(mol) - # mol2strs.append(mol2str) - # xtb_energies1.append(energy) - # init_mol2strs.append(init_mol2str) - # init_xtb_energies.append(init_energy) - # if len(mol2strs) == inputDict['parameters']['n_lig_combos']: - # break - # for k in range(len(mol2strs)): - # outkey = keys[i] + '_random_complex_' + str(k) - # ordered_conf_dict[outkey] = {'ase_atoms':mol.ase_atoms, - # # 'complex':structs[i], # Debug - # 'total_charge':np.sum(mols[k].ase_atoms.get_initial_charges()), - # 'n_unpaired_electrons': np.sum(mols[k].ase_atoms.get_initial_magnetic_moments()), - # 'calc_n_unpaired_electrons': io_xtb_calc.calc_suggested_spin(structs[i].complexMol, - # parameters=inputDict['parameters']), - # 'metal_ox':inputDict['parameters']['metal_ox'], - # 'init_energy':init_xtb_energies[k], - # 'energy':xtb_energies1[k], - # 'mol2string':mol2strs[k], 'init_mol2string':init_mol2strs[k], - # 'energy_sorted_index': energy_sorted_inds[i], - # 'inputDict':inputDict - # } - # if inputDict['parameters']['return_full_complex_class']: # Return whole complex class (all ligand geometries!) - # ordered_conf_dict[outkey].update({'full_complex_class':structs[i] - # }) - # if inputDict['parameters']['return_timings']: - # tdict.update({'core_preprocess_time':core_preprocess_time, - # 'symmetry_preprocess_time':symmetry_preprocess_time, - # 'total_liggen_time':np.sum([x.total_liggen_time for x in structs[i].ligandList]), - # 'total_complex_assembly_time':structs[i].assemble_total_time, - # 'final_relaxation_time':structs[i].final_eval_total_time, - # 'sum_total_conformer_time_spent':fin_time2 - int_time1 - # }) \ No newline at end of file diff --git a/development/uranium_clusters.ipynb b/development/uranium_clusters.ipynb new file mode 100644 index 0000000..81054c2 --- /dev/null +++ b/development/uranium_clusters.ipynb @@ -0,0 +1,897 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "from architector import build_complex,view_structures,convert_io_molecule\n", + "out = build_complex(\n", + " {'core': {'metal': 'U', 'coreCN': 4}, 'ligands': \n", + " ['oxo','oxo','oxygen_edge'],\n", + " 'parameters': {'metal_ox': 6, 'metal_spin': 0, \n", + " 'assemble_method': 'UFF', 'full_method': 'UFF'}})" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [ + { + "data": { + "application/3dmoljs_load.v0": "
\n

3Dmol.js failed to load for some reason. Please check your browser console for error messages.

\n
\n", + "text/html": [ + "
\n", + "

3Dmol.js failed to load for some reason. Please check your browser console for error messages.

\n", + "
\n", + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "view_structures(out)" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [], + "source": [ + "import copy\n", + "mol = convert_io_molecule(out[list(out.keys())[0]]['mol2string'])\n", + "mol1 = convert_io_molecule(out[list(out.keys())[0]]['mol2string'])\n", + "mol2 = copy.deepcopy(mol)" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "data": { + "application/3dmoljs_load.v0": "
\n

3Dmol.js failed to load for some reason. Please check your browser console for error messages.

\n
\n", + "text/html": [ + "
\n", + "

3Dmol.js failed to load for some reason. Please check your browser console for error messages.

\n", + "
\n", + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "view_structures(mol,labelinds=True)" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [], + "source": [ + "mol2.remove_atoms([1,2])" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [], + "source": [ + "mol.functionalize_3D(functional_groups=[mol2],\n", + " functionalization_inds=[[1, 2]],\n", + " bond_orders=[[1, 1]],\n", + " functional_group_mol_inds=[[0, 0]],\n", + " remove_hydrogens_when_adding=[[False,False]],\n", + " core_frozen=False,\n", + " xtb_opt=True)" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [ + { + "data": { + "application/3dmoljs_load.v0": "
\n

3Dmol.js failed to load for some reason. Please check your browser console for error messages.

\n
\n", + "text/html": [ + "
\n", + "

3Dmol.js failed to load for some reason. Please check your browser console for error messages.

\n", + "
\n", + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "view_structures(mol,labelinds=True)" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [ + { + "data": { + "application/3dmoljs_load.v0": "
\n

3Dmol.js failed to load for some reason. Please check your browser console for error messages.

\n
\n", + "text/html": [ + "
\n", + "

3Dmol.js failed to load for some reason. Please check your browser console for error messages.

\n", + "
\n", + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "view_structures(mol1,labelinds=True)" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [], + "source": [ + "mol.functionalize_3D(functional_groups=[mol1],\n", + " functionalization_inds=[[0, 0]],\n", + " bond_orders=[[1, 1]],\n", + " functional_group_mol_inds=[[1, 2]],\n", + " remove_hydrogens_when_adding=[[False,False]],\n", + " core_frozen=False,\n", + " xtb_opt=True)" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [ + { + "data": { + "application/3dmoljs_load.v0": "
\n

3Dmol.js failed to load for some reason. Please check your browser console for error messages.

\n
\n", + "text/html": [ + "
\n", + "

3Dmol.js failed to load for some reason. Please check your browser console for error messages.

\n", + "
\n", + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "view_structures(mol,labelinds=True)" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": {}, + "outputs": [], + "source": [ + "mol.functionalize_3D(functional_groups=[mol1],\n", + " functionalization_inds=[[5, 5]],\n", + " bond_orders=[[1, 1]],\n", + " functional_group_mol_inds=[[1, 2]],\n", + " remove_hydrogens_when_adding=[[False,False]],\n", + " core_frozen=False,\n", + " xtb_opt=True)" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": {}, + "outputs": [ + { + "data": { + "application/3dmoljs_load.v0": "
\n

3Dmol.js failed to load for some reason. Please check your browser console for error messages.

\n
\n", + "text/html": [ + "
\n", + "

3Dmol.js failed to load for some reason. Please check your browser console for error messages.

\n", + "
\n", + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "view_structures(mol,labelinds=True)" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": {}, + "outputs": [ + { + "data": { + "application/3dmoljs_load.v0": "
\n

3Dmol.js failed to load for some reason. Please check your browser console for error messages.

\n
\n", + "text/html": [ + "
\n", + "

3Dmol.js failed to load for some reason. Please check your browser console for error messages.

\n", + "
\n", + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "view_structures(mol,w=400,h=400)" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": {}, + "outputs": [], + "source": [ + "from architector.io_arch_dock import add_non_covbound_species" + ] + }, + { + "cell_type": "code", + "execution_count": 22, + "metadata": {}, + "outputs": [], + "source": [ + "out = add_non_covbound_species(mol,\n", + " parameters={'species_list':[mol]*2,\n", + " \"species_add_copies\":100, # Number \"randomized\" generations to do.\n", + " \"species_location_method\":'random', # Use random for actiniyls \n", + " \"species_relax\":False, # Don't relax (none of the methods are super accurate)\n", + " \"species_intermediate_relax\":False, # Don't relax intermediate\n", + " \"species_intermediate_method\":'UFF', # Use UFF for speed\n", + " \"species_method\":'UFF', # Use UFF for speed\n", + " \"debug\":False}) # Turn on if you want to see what it's doing in the background })" + ] + }, + { + "cell_type": "code", + "execution_count": 26, + "metadata": {}, + "outputs": [ + { + "data": { + "application/3dmoljs_load.v0": "
\n

3Dmol.js failed to load for some reason. Please check your browser console for error messages.

\n
\n", + "text/html": [ + "
\n", + "

3Dmol.js failed to load for some reason. Please check your browser console for error messages.

\n", + "
\n", + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "view_structures(out[1][:12])" + ] + }, + { + "cell_type": "code", + "execution_count": 24, + "metadata": {}, + "outputs": [ + { + "data": { + "application/3dmoljs_load.v0": "
\n

3Dmol.js failed to load for some reason. Please check your browser console for error messages.

\n
\n", + "text/html": [ + "
\n", + "

3Dmol.js failed to load for some reason. Please check your browser console for error messages.

\n", + "
\n", + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "view_structures(out[0])" + ] + }, + { + "cell_type": "code", + "execution_count": 28, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "0" + ] + }, + "execution_count": 28, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "out[0].uhf" + ] + }, + { + "cell_type": "code", + "execution_count": 29, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "24.0" + ] + }, + "execution_count": 29, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "out[0].charge" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "base", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.10.12" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +}