Skip to content

Commit

Permalink
Uranyl debug dev.
Browse files Browse the repository at this point in the history
  • Loading branch information
@mgt16-LANL
mgt16-LANL committed May 9, 2024
1 parent 8971f09 commit 1d4bb22
Show file tree
Hide file tree
Showing 6 changed files with 1,059 additions and 422 deletions.
62 changes: 44 additions & 18 deletions architector/io_calc.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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.
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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():
Expand Down Expand Up @@ -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))
Expand All @@ -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()
Expand All @@ -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:
Expand Down Expand Up @@ -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 _:
Expand All @@ -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())
Expand Down
74 changes: 65 additions & 9 deletions architector/io_molecule.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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'],
Expand Down Expand Up @@ -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:
Expand Down Expand Up @@ -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)):
Expand Down Expand Up @@ -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:
Expand All @@ -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())


Expand Down
Loading

0 comments on commit 1d4bb22

Please sign in to comment.