Merge pull request #22 from lanl/Oxo_handling

Oxo handling

README.md

@@ -160,7 +160,9 @@ inputDict = {
     # in case the fill_ligand and specified ligands list cannot fully map to the coordination environment.
     "secondary_fill_ligand": "water",
     # or integer index in reference to the ligand list!!
-    "force_trans_oxos":False, # Force trans configurations for oxos (Useful for actinyls)
+    "force_trans_oxos":True, # Force trans configurations for oxos (Useful for actinyls)
+    # Will only be activated when actinides are present - otherwise will not force trans oxos.
+    "override_oxo_opt":True, # Override no relaxation of oxo groups (not generally suggested)
     "lig_assignment":'bruteforce', # or "similarity" - How to automatically assign ligand types.
 
     ######### Sanity check parameters ########

architector/complex_construction.py

@@ -239,22 +239,25 @@ def final_eval(self,single_point=False):
         self.initMol.dist_sanity_checks(params=self.parameters,assembly=single_point)
         self.initMol.graph_sanity_checks(params=self.parameters,assembly=single_point)
         if self.assembled:
-            if self.parameters['ff_preopt'] and (not single_point):
+            if self.parameters['ff_preopt']:
                 if self.parameters['debug']:
                     print('Doing UFF - pre-optimization before final evaluation.')
                 calculator = CalcExecutor(self.complexMol, parameters=self.parameters,
                                 final_sanity_check=False,
-                                relax=True, assembly=single_point, 
-                                ff_preopt_run=self.parameters['ff_preopt'])
+                                relax=True, assembly=False, 
+                                ff_preopt_run=self.parameters['ff_preopt'],
+                                fix_m_neighbors=True)
                 if calculator:
                     self.complexMol.ase_atoms.set_positions(calculator.mol.ase_atoms.get_positions())
+                elif self.parameters['debug']:
+                    print('Pre-opt failed!')
             if self.parameters['debug']:
-                print("Final Evaluation - Opt Molecule/Single point")
+                print("Final Evaluation - Opt Molecule/Single point, Single Point? {}".format(single_point))
             self.calculator = CalcExecutor(self.complexMol,
                                            parameters=self.parameters,
                                            final_sanity_check=self.parameters['full_sanity_checks'],
-                                           relax=single_point,
-                                           assembly=single_point)
+                                           relax=(not single_point),
+                                           assembly=False)
             if (self.parameters['debug']) and (self.calculator):
                 print('Finished method: {} {}.'.format(self.calculator.method,
                                                     self.calculator.relax))

architector/data/avg_m_oxo_dists_csd.csv

@@ -0,0 +1,28 @@
+metal,avg,min,max,std,count,
+Am,0.76788761,0.760014514,0.773881513,0.007070752,4,
+Ce,0.814035047,0.814035047,0.814035047,0,1,
+Co,0.959395961,0.952031159,0.966760764,0.010415403,2,
+Cr,0.87392861,0.799180314,1.018453462,0.023193201,332,
+Fe,0.945075794,0.907370912,1.031026613,0.038593856,24,
+Hf,0.849062547,0.849062547,0.849062547,0,1,
+Ir,0.969059148,0.932336913,0.987423129,0.031802388,3,
+La,1.06431122,1.052084396,1.076538045,0.017291341,2,
+Mn,0.877267486,0.767598761,1.028880997,0.065610529,43,
+Mo,0.846951502,0.759104478,1.198433111,0.022023247,3690,
+Nb,0.82760043,0.783656402,1.108295549,0.03621571,119,
+Ni,1.186873115,1.083468208,1.290278022,0.146236622,2,
+Np,0.764864742,0.722526454,1.051528513,0.028816006,188,
+Os,0.899358307,0.843020833,1.026595988,0.016334881,327,
+Pa,0.750843127,0.750843127,0.750843127,0,1,
+Pu,0.744027505,0.719301711,0.776125597,0.010125398,84,
+Re,0.876609345,0.569944295,1.201640147,0.028364752,3123,
+Ru,0.917002467,0.856749575,0.971583731,0.023489772,79,
+Ta,0.842807854,0.821873648,0.971339713,0.036170194,15,
+Tc,0.880331045,0.819714799,0.99984932,0.022365732,314,
+Th,0.810283252,0.810283252,0.810283252,0,1,
+Ti,0.828192864,0.809741878,0.899145729,0.013356636,56,
+U,0.764211927,0.583212372,1.112612324,0.040440426,3345,
+V,0.81503691,0.773877719,1.210517607,0.022300605,2577,
+W,0.862678567,0.758994226,1.226126834,0.026110934,835,
+Zr,0.831452232,0.831452232,0.831452232,0,1,
+Others,0.868044177,0.812749222,0.999335112,0.025628015,583.4230769,

architector/io_calc.py

@@ -56,12 +56,15 @@
 "full_charge": None, # Assign charge to the complex (overrides ligand charges and metal_ox)!
 "full_method":"GFN2-xTB", # Which xtb method to use for final cleaning/evaulating conformers.
 "assemble_method":"GFN2-xTB",
+"ff_preopt":False,
+"override_oxo_opt":False,
 "fmax":0.1,
 "maxsteps":1000,
 "vdwrad_metal":None,
 "covrad_metal":None,
 "scaled_radii_factor":None,
 "force_generation":False,
+"species_run":False,
 "debug":False,
 }
 
@@ -72,7 +75,7 @@ def __init__(self, structure, parameters={}, init_sanity_check=False,
                 xtb_electronic_temperature=300, xtb_max_iterations=250,
                 fmax=0.1, maxsteps=1000, ff_preopt_run=False,
                 detect_spin_charge=False, fix_m_neighbors=False,
-                default_params=params, ase_opt_method=None,
+                default_params=params, ase_opt_method=None, species_run=False,
                 calculator=None):
         """CalcExecutor is the class handling all calculations of full metal-ligand complexes.
 
@@ -117,6 +120,8 @@ def __init__(self, structure, parameters={}, init_sanity_check=False,
             default parameters to evaluate to, by default params
         ase_opt_method : ase.optimize Optimizer, optional
             Valid ase style optimizer if desired, by default None
+        species_run : bool, optional
+            Flag if this run is performed for a "species" to be added.
         calculator : ase.calculators Calculator, optional
             Valid ase style calculator if desired, by default None
         """
@@ -140,21 +145,32 @@ def __init__(self, structure, parameters={}, init_sanity_check=False,
         self.fmax = fmax
         self.fix_m_neighbors = fix_m_neighbors
         self.maxsteps = maxsteps
+        self.species_run = species_run
         self.force_generation = False
 
         self.detect_spin_charge = detect_spin_charge
         if len(parameters) > 0:
             for key,val in parameters.items():
                 setattr(self,key,val)
-            if assembly:
-                self.init_sanity_check = True
-                self.relax = False
-                self.method = self.assemble_method
+
+        if assembly:
+            self.init_sanity_check = True
+            self.relax = False
+            self.method = self.assemble_method
+        elif species_run:
+            self.method = self.species_xtb_method
+            self.relax = self.species_relax
+        elif len(parameters) > 0:
+            if self.ff_preopt_run:
+                self.method = 'UFF'
+                self.relax=True
             else:
-                if self.ff_preopt_run:
-                    self.method = 'UFF'
-                else:
-                    self.method = self.full_method
+                self.method = self.full_method
+        else:
+            if self.ff_preopt_run:
+                self.method = 'UFF'
+                self.relax=True
+
         if ase_opt_method is None: # Default to LBFGSLineSearch
             self.opt_method = LBFGSLineSearch
         else:
@@ -198,8 +214,10 @@ def calculate(self):
                            max_iterations=self.xtb_max_iterations,
                            electronic_temperature=self.xtb_electronic_temperature,
                            accuracy=self.xtb_accuracy)
-                if np.abs(self.mol.xtb_charge - self.mol.charge) > 1: # Difference of more than 1.
-                    self.relax = False # E.g - don't relax if way off in oxdiation states (III) vs (V or VI)
+                # 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):
+                        self.relax = False # E.g - don't relax if way off in oxdiation states (III) vs (V or VI)
                     if self.assembly: # FF more stable for highly charged assembly complexes.
                         self.method = 'GFN-FF'
                 uhf_vect = np.zeros(len(self.mol.ase_atoms))

architector/io_lig.py

@@ -19,6 +19,8 @@
 """
 
 import numpy as np
+import os
+import pandas as pd
 import itertools
 import scipy
 import architector.io_ptable as io_ptable
@@ -42,6 +44,20 @@
 
 warnings.filterwarnings('ignore') # Supress numpy warnings.
 
+def get_oxo_refdict():
+    """get_oxo_refdict
+    Pull the metal-oxo distance reference dictionary.
+
+    Returns
+    -------
+    ref_df : dict
+        dictionary with {'type':angles(List)}
+    """
+    filepath = os.path.abspath(os.path.join(__file__, "..", "data", "avg_m_oxo_dists_csd.csv"))
+    ref_df = pd.read_csv(filepath)
+    refdict = dict(zip(ref_df.metal, ref_df.avg))
+    return refdict
+
 def set_XTB_calc(ase_atoms):
     """set_XTB_calc 
     assign xtb calculator to atoms instance!
@@ -817,6 +833,15 @@ def clean_conformation_ff(X, OBMol, catoms, shape, graph,
     fail = False
     # Set metal to zero
     metal_coords = (X[last_atom_index-1,0],X[last_atom_index-1,1],X[last_atom_index-1,2])
+    # Check if oxo
+    if (len(catoms) == 1) and (len(atomic_numbers) == 2) and (atomic_numbers[0] == 8):
+        msym = io_ptable.elements[atomic_numbers[1]]
+        refdict = get_oxo_refdict()
+        bondls = float(refdict.get(msym,refdict.get('Others',None)))
+        bondls = bondls * (io_ptable.rcov1[io_ptable.elements.index(msym)] + io_ptable.rcov1[8])
+        oxo = True
+    else:
+        oxo = False
     # set coordinates using OBMol to keep bonding info
     for i, atom in enumerate(openbabel.OBMolAtomIter(OBMol)):
         atom.SetVector(X[i, 0]-metal_coords[0], X[i, 1]-metal_coords[1], X[i, 2]-metal_coords[2])
@@ -825,6 +850,9 @@ def clean_conformation_ff(X, OBMol, catoms, shape, graph,
         constr = openbabel.OBFFConstraints()
         ff = openbabel.OBForceField.FindForceField('UFF')
         constr.AddAtomConstraint(int(last_atom_index))
+        # Oxo set distance based on CSD
+        if oxo:
+            constr.AddDistanceConstraint(1,2,bondls)
         s = ff.Setup(OBMol,constr)
         if not s:
             if debug:
@@ -840,6 +868,9 @@ def clean_conformation_ff(X, OBMol, catoms, shape, graph,
                 constr = openbabel.OBFFConstraints()
                 ff = openbabel.OBForceField.FindForceField('UFF')
                 constr.AddAtomConstraint(int(last_atom_index))
+                # Oxo set distance based on CSD
+                if oxo:
+                    constr.AddDistanceConstraint(1,2,bondls)
                 s = ff.Setup(OBMol,constr)
                 for i in range(100):
                     ff.SteepestDescent(10)
@@ -849,7 +880,7 @@ def clean_conformation_ff(X, OBMol, catoms, shape, graph,
             # Reorder the coordinating atom assignment to more closely match unconstrained
             # UFF relaxation to the desired geometry.
             tligcoordList, catoms = reorder_ligcoordList(coords, catoms, shape, ligcoordList, isCp=isCp)
-            if add_angle_constraints:
+            if add_angle_constraints and (len(atomic_numbers) > 2):
                 if debug:
                     print('Finished initial UFF relaxation without angle constraints.')
                 ff = openbabel.OBForceField.FindForceField('UFF')
@@ -1521,6 +1552,8 @@ def get_aligned_conformer(ligsmiles, ligcoordList, corecoordList, metal='Fe',
         crowding_penalty = 1000
         if debug:
             print('Failed Relaxation!')
+    if ligsmiles == '[O-2]':
+        sane=True
     minval += crowding_penalty
     return outatoms, minval, sane, final_relax, bo_dict, atypes, tligcoordList
 

architector/io_process_input.py

@@ -426,7 +426,7 @@ def inparse(inputDict):
                             newdict.update({'smiles':ligDict['smiles']})
                             if (ligDict['smiles'] == '[H-]')or (ligDict['smiles'] == '[O-2]'): # Handle hydrides
                                 newinpDict['parameters']['relax'] = False
-                                newinpDict['parameters']['metal_spin'] = 0 # Set to low_spin
+                                # newinpDict['parameters']['metal_spin'] = 0 # Set to low_spin
                             if (ligDict['ligType'] in core_geo_class.liglist_geo_map_dict.keys()) or (ligDict['ligType'] == 'mono'):
                                 newdict.update({'ligType':ligDict['ligType']})
                             else:
@@ -444,7 +444,7 @@ def inparse(inputDict):
                             newdict.update({'smiles':ligDict['smiles']})
                             if (ligDict['smiles'] == '[H-]') or (ligDict['smiles'] == '[O-2]'): # Handle hydrides
                                 newinpDict['parameters']['relax'] = False
-                                newinpDict['parameters']['metal_spin'] = 0 # Set to low-spin for oxos
+                                # newinpDict['parameters']['metal_spin'] = 0 # Set to low-spin for oxos
                             if isinstance(ligDict['coordList'][0],list):
                                 newdict.update({'coordList':ligDict['coordList']})
                                 # Assign ligType arbitrarily -> will be ignored by io_symmetry if list of lists
@@ -598,7 +598,7 @@ def inparse(inputDict):
                     tligdict = {'smiles':lig_smiles,'coordList':coord_atoms[i]}
                     if (tligdict['smiles'] == '[H-]') or (tligdict['smiles'] == '[O-2]'): # Handle hydrides
                         newinpDict['parameters']['relax'] = False
-                        newinpDict['parameters']['metal_spin'] = 0 # Set to low-spin for oxos
+                        # newinpDict['parameters']['metal_spin'] = 0 # Set to low-spin for oxos
                         newinpDict['parameters']['full_spin'] = None
                         tligdict.update({'ligType':'mono'})
                     elif len(tligdict['coordList']) == 1: # Monodentate
@@ -709,7 +709,8 @@ def inparse(inputDict):
             # in case the fill_ligand and specified ligands list cannot fully map to the coordination environment.
             "secondary_fill_ligand": "water",
             # or integer index in reference to the ligand list!!
-            "force_trans_oxos":False, # Force trans configurations for oxos (Useful for actinyls)
+            "force_trans_oxos":True, # Force trans configurations for oxos (Useful for actinyls)
+            "override_oxo_opt":True, # Override no relaxation of oxo groups (not generally suggested)
             "lig_assignment":'bruteforce', # or "similarity" How to automatically assign ligand types.
 
             # Cutoff parameters
@@ -770,13 +771,18 @@ def inparse(inputDict):
         outparams.update(default_parameters) 
         # If acinide default to forcing trans oxos if 2 present (actinyls)
         # Can still be turned off if desired.
-        if newinpDict['parameters']['is_actinide']:
+        if (newinpDict['parameters']['is_actinide']) and \
+            (newinpDict['parameters'].get('force_trans_oxos',False)):
             count = 0 
             for lig in newinpDict['ligands']:
                 if lig['smiles'] == '[O-2]':
                     count += 1
-            if count == 2:
+            if (count == 2):
                 newinpDict['parameters']['force_trans_oxos'] = True
+            else:
+                newinpDict['parameters']['force_trans_oxos'] = False
+        else:
+            newinpDict['parameters']['force_trans_oxos'] = False
 
         if ('parameters' in newinpDict): # Load input parameters as changes to the defaults.
             if isinstance(newinpDict['parameters'],dict):
@@ -860,7 +866,7 @@ def inparse(inputDict):
 
         # FF-preoptimization shifts
         if outparams['ff_preopt']:
-            outparams['assemble_method'] = 'GFN-FF'
+            outparams['assemble_method'] = 'GFN-FF' # Less likely to crash with short interatomic distances.
             outparams['assemble_sanity_checks'] = False # Allow for close/overlapping atoms 
             outparams['fmax'] = 0.2 # Slightly decrease accuracy to encourage difficult convergence
 
@@ -871,6 +877,9 @@ def inparse(inputDict):
         # Initialize seed.
         if isinstance(newinpDict['parameters']['seed'],(int,float,np.float64,np.int64)):
             np.random.seed(int(newinpDict['parameters']['seed']))
+
+        if newinpDict['parameters']['debug']:
+            print(newinpDict)
 
         return newinpDict
     else:
@@ -1027,7 +1036,7 @@ def inparse_2D(inputDict):
                         newdict.update({'smiles':ligDict['smiles']})
                         if (ligDict['smiles'] == '[H-]')or (ligDict['smiles'] == '[O-2]'): # Handle hydrides
                             newinpDict['parameters']['relax'] = False
-                            newinpDict['parameters']['metal_spin'] = 0 # Set to low_spin
+                            # newinpDict['parameters']['metal_spin'] = 0 # Set to low_spin
                         if (ligDict['ligType'] == 'mono') or (isinstance(ligDict['ligType'],str)):
                             newdict.update({'ligType':ligDict['ligType']})
                         else:
@@ -1045,7 +1054,7 @@ def inparse_2D(inputDict):
                         newdict.update({'smiles':ligDict['smiles']})
                         if (ligDict['smiles'] == '[H-]') or (ligDict['smiles'] == '[O-2]'): # Handle hydrides
                             newinpDict['parameters']['relax'] = False
-                            newinpDict['parameters']['metal_spin'] = 0 # Set to low-spin for oxos
+                            # newinpDict['parameters']['metal_spin'] = 0 # Set to low-spin for oxos
                         if isinstance(ligDict['coordList'][0],list):
                             newdict.update({'coordList':ligDict['coordList']})
                             # Assign ligType arbitrarily -> will be ignored by io_symmetry if list of lists