Merge pull request #14 from lanl/FF_preopt

Add FF pre-optimization option for final-generated architector structures.

Author: mgt16-LANL

README.md

@@ -30,8 +30,9 @@ pip install -e .
 1. See tutorials for basic introduction to capabilties and code examples: `documentation/tutorials/`
 2. Reference for core and ligand geometry labels see: `documentation/view_default_core_ligand_types.ipynb`
 3. Utility for aiding in determining ligand coordination sites see: `utils/ligand_viewing_coordinating_atom_selecting.ipynb`
+4. Reference for debugging cases where no structures are generated see: `documentation/Debugging_Guide.ipynb`
 
-* Note that ligands used in (3) can even be drawn in [Avogadro](https://avogadro.cc/) and copied as SMILES strings into this analysis.
+* Note that ligands used in (3) can even be drawn in [Avogadro](https://avogadro.cc/) and or most other 2D molecular editors (e.g. ChemDraw) and copied as SMILES strings into this analysis.
 * If other analyses are used to determine the coordinating atom indices we can't guarantee the generated structure will match what was input. If generating complexes with new ligands we HIGHLY recommend using the utility in (3)
 
 ## XTB (backend) Potentially Useful References:
@@ -115,6 +116,11 @@ inputDict = {
     "xtb_max_iterations": 250, # Max iterations for xtb SCF.
     "force_generation":False, # Whether to force the construction to proceed without xtb energies - defaults to UFF evaluation
     # in cases of XTB outright failure. Will still enforce sanity checks on output structures.
+    "ff_preopt":False, # Whether to force forcefield (FF) pre-optimization of fully-built complexes before final xtb evalulation.
+    # FF preoptimization helps especially in cases where longer carbon chains are present that tend to overlap.
+    # This option will also decrease the maximum force tolerance for xtb relaxation to 0.2 and set assemble_method='GFN-FF'
+    # By default for acceleration.
+
 
     # Covalent radii and vdw radii of the metal if nonstandard radii requested.
     "vdwrad_metal": vdwrad_metal,

architector/complex_construction.py

@@ -19,7 +19,6 @@
 import architector.io_align_mol as io_align_mol
 import architector.io_crest as io_crest
 import architector.io_symmetry as io_symmetry
-
 from architector.io_calc import CalcExecutor
 
 class Ligand:
@@ -240,11 +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['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'])
+                if calculator:
+                    self.complexMol.ase_atoms.set_positions(calculator.mol.ase_atoms.get_positions())
             if self.parameters['debug']:
                 print("Final Evaluation - Opt Molecule/Single point")
-            self.calculator = CalcExecutor(self.complexMol,parameters=self.parameters,
-                                            final_sanity_check=self.parameters['full_sanity_checks'],
-                                            relax=single_point,assembly=single_point)
+            self.calculator = CalcExecutor(self.complexMol,
+                                           parameters=self.parameters,
+                                           final_sanity_check=self.parameters['full_sanity_checks'],
+                                           relax=single_point,
+                                           assembly=single_point)
+            if (self.parameters['debug']) and (self.calculator):
+                print('Finished method: {} {}.'.format(self.calculator.method,
+                                                    self.calculator.relax))
             if self.parameters['debug'] and (not self.calculator.successful):
                 print('Failed final relaxation. - Retrying with UFF/XTB')
                 print(self.initMol.write_mol2('cool.mol2', writestring=True))

architector/io_calc.py

@@ -70,7 +70,7 @@ def __init__(self, structure, parameters={}, init_sanity_check=False,
                 final_sanity_check=False, relax=False, assembly=False,
                 method='GFN2-xTB', xtb_solvent='none', xtb_accuracy=1.0,
                 xtb_electronic_temperature=300, xtb_max_iterations=250,
-                fmax=0.1, maxsteps=1000,
+                fmax=0.1, maxsteps=1000, ff_preopt_run=False,
                 detect_spin_charge=False, fix_m_neighbors=False,
                 default_params=params, ase_opt_method=None,
                 calculator=None):
@@ -107,6 +107,8 @@ def __init__(self, structure, parameters={}, init_sanity_check=False,
             Max force in eV/Angstrom, by default 0.1
         maxsteps : int, optional
             Total number of optimization steps to take, by default 1000
+        ff_preopt_run : bool, optional
+            Perform forcefield pre-optimization?, by default False
         detect_spin_charge : bool, optional
             Detect the charge and spin with openbabel routines?, by default False
         fix_m_neighbors : bool, optional
@@ -130,6 +132,7 @@ def __init__(self, structure, parameters={}, init_sanity_check=False,
         self.calculator = calculator
         self.relax = relax
         self.assembly = assembly
+        self.ff_preopt_run = ff_preopt_run
         self.xtb_solvent = xtb_solvent
         self.xtb_accuracy = xtb_accuracy
         self.xtb_electronic_temperature = xtb_electronic_temperature
@@ -148,7 +151,10 @@ def __init__(self, structure, parameters={}, init_sanity_check=False,
                 self.relax = False
                 self.method = self.assemble_method
             else:
-                self.method = self.full_method
+                if self.ff_preopt_run:
+                    self.method = 'UFF'
+                else:
+                    self.method = self.full_method
         if ase_opt_method is None: # Default to LBFGSLineSearch
             self.opt_method = LBFGSLineSearch
         else:

architector/io_obabel.py

@@ -379,7 +379,8 @@ def build_3D(OBmol,addHydrogens=True):
 
 
 def generate_obmol_conformers(smiles, rmsd_cutoff=0.4, conf_cutoff=3000, energy_cutoff=50.0, 
-        confab_verbose = False, output_format='mol2', neutralize=False, functionalizations=None):
+        confab_verbose = False, output_format='mol2', neutralize=False, functionalizations=None,
+        return_energies = False):
     """generate_obmol_conformers 
     generate conformers with openbabel for given smiles
     using confab conformer generation routine
@@ -405,11 +406,15 @@ def generate_obmol_conformers(smiles, rmsd_cutoff=0.4, conf_cutoff=3000, energy_
         neutralize smiles?, by default False
     functionalizations : dict, optional
         add functionalizations?, by default None
+    return_energies : bool, optional
+        return the FF energies in addition to the conformers generated
 
     Returns
     -------
     output_strings : list (str)
         list of conformers generated as whatever format desired
+    output_energies : list (float)
+        forcefield energies
     """
     obmol = get_obmol_smiles(smiles,
                              neutralize=neutralize,
@@ -419,18 +424,36 @@ def generate_obmol_conformers(smiles, rmsd_cutoff=0.4, conf_cutoff=3000, energy_
         FF = ob.OBForceField.FindForceField("MMFF94")
         FF.Setup(obmol) # Make sure setup works OK
     else:
-        FF = ob.OBForceField.FindForceField("MMFF94")
+        FF = ob.OBForceField.FindForceField("UFF")
         FF.Setup(obmol) # Make sure setup works OK
+    # Run Diverse conformer generation
     FF.DiverseConfGen(rmsd_cutoff, conf_cutoff, energy_cutoff, confab_verbose)
     FF.GetConformers(obmol)
     confs_to_write = obmol.NumConformers()
     obconversion = ob.OBConversion()
     obconversion.SetOutFormat(output_format)
     output_strings = []
+    output_energies = []
     for conf_num in range(confs_to_write):
         obmol.SetConformer(conf_num)
+        if return_energies: # Calculate FF energies
+            if mmff94_ok:
+                FF = ob.OBForceField.FindForceField("MMFF94")
+                FF.Setup(obmol) # Make sure setup works OK
+            else:
+                FF = ob.OBForceField.FindForceField("UFF")
+                FF.Setup(obmol) # Make sure setup works OK
+            energy = FF.Energy()
+            if mmff94_ok:
+                energy = energy * units.kcal / units.mol
+            else:
+                energy = energy * units.kJ / units.mol
+            output_energies.append(energy)
         output_strings.append(obconversion.WriteString(obmol))
-    return output_strings
+    if return_energies:
+        return output_strings,output_energies
+    else:
+        return output_strings
 
 
 def functionalize(OBmol,functional_group='C',smiles_inds=[0]):

architector/io_process_input.py

@@ -744,6 +744,7 @@ def inparse(inputDict):
             "fmax":0.1, # eV/Angstrom - max force for relaxation.
             "maxsteps":1000, # Steps involved in relaxation
             "force_generation":False, # Whether to force the construction to proceed without xtb energies - defaults to UFF
+            "ff_preopt":False, # Perform forcefield pre-optimization of full complex? - default False.
             # In cases of XTB outright failure.
             # For very large speedup - use GFN-FF, though this is much less stable (especially for Lanthanides)
             # Or UFF
@@ -853,6 +854,12 @@ def inparse(inputDict):
             coreTypes = [x for x in newinpDict['coreTypes'] if x in core_geo_class.trans_geos]
             newinpDict['coreTypes'] = coreTypes
 
+        # FF-preoptimization shifts
+        if outparams['ff_preopt']:
+            outparams['assemble_method'] = 'GFN-FF'
+            outparams['assemble_sanity_checks'] = False # Allow for close/overlapping atoms 
+            outparams['fmax'] = 0.2 # Slightly decrease accuracy to encourage difficult convergence
+ 
         # # Load logger
         # if outparams['logg']
         newinpDict['parameters'] = outparams

documentation/Debugging_Guide.ipynb

@@ -49,6 +49,7 @@
     "Solutions:\n",
     "1. Play with sanity check cutoffs. You can turn off final an intermediate sanity checks with \"full_sanity_checks\":False, and \"assemble_sanity_checks\":False, respectively. You can also tune the actual cutoff values as well : \"full_graph_sanity_cutoff\":1.7 - can be tuned higher to allow for looser constraints, \"full_smallest_dist_cutoff\":0.55 - can be tuned lower for closer interatomic distances, and \"full_min_dist_cutoff\":3.5 - can be tuned higher to allow for more isolated atoms.\n",
     "2. Metal van-der-waals radii (\"vdwrad_metal\") and covalent radii (\"covrad_metal\") can be manually tuned in cases where a very distinct oxiation state is used. \n",
+    "3. You might have ligands with trailing carbon chains that tend to overlap when being tile to the metal center surface. For example, the TODGA ligand can have this issue. Try setting \"ff_preopt\":True - this will generate the full complex, relax it with UFF - and only then perform the final relaxation method (usually GFN2-xTB). The UFF pre-optimization dramatically increases the chance for successful XTB relaxation in tests, but does introduce more bias into the generated structures.\n",
     "3. A last-dich effort can be made to simply force Architector to return basically anything with \"force_generation\":True. This can be useful to determine if something is going wrong internally in Architector. It is not recommended to use this option for production (multiple-structures runs)."
    ]
   },