Internal coordinate interpolation

geometric/nifty.py

@@ -843,10 +843,11 @@ def destroyWorkQueue():
     # Convenience function to destroy the Work Queue objects.
     global WORK_QUEUE, WQIDS
     if hasattr(WORK_QUEUE, '_free'): 
-        print("Freeing Work Queue")
         WORK_QUEUE._free()
     WORK_QUEUE = None
     WQIDS = defaultdict(list)
+    print("Freed Work Queue, sleeping for 5 seconds...")
+    time.sleep(5)
 
 def queue_up(wq, command, input_files, output_files, tag=None, tgt=None, verbose=True, print_time=60):
     """

geometric/optimize.py

@@ -571,6 +571,8 @@ def evaluateStep(self):
             colors['energy'] = "\x1b[92m" if Converged_energy else "\x1b[91m"
             colors['quality'] = "\x1b[91m"
             step_state = StepState.Reject if (Quality < -1.0 or params.transition) else StepState.Poor
+        # 2023-05-10: In the case of a minimization step that increases energy by more than erisemax, reject the step.
+        if not params.transition and not self.farConstraints and self.E-self.Eprev > params.erisemax: step_state = StepState.Reject
         if 'energy' not in colors: colors['energy'] = "\x1b[92m" if Converged_energy else "\x1b[0m"
         if 'quality' not in colors: colors['quality'] = "\x1b[0m"
         colors['grms'] = "\x1b[92m" if Converged_grms else "\x1b[0m"
@@ -674,7 +676,10 @@ def evaluateStep(self):
             elif self.farConstraints:
                 logger.info("\x1b[93mNot rejecting step - far from constraint satisfaction\x1b[0m\n")
             else:
-                logger.info("\x1b[93mRejecting step - quality is lower than %.1f\x1b[0m\n" % (0.0 if params.transition else -1.0))
+                if not params.transition and not self.farConstraints and self.E-self.Eprev > params.erisemax:
+                    logger.info("\x1b[93mRejecting step - energy increases > %.2f during minimization\x1b[0m\n" % params.erisemax)
+                else:
+                    logger.info("\x1b[93mRejecting step - quality is lower than %.1f\x1b[0m\n" % (0.0 if params.transition else -1.0))
                 self.trustprint = "\x1b[1;91mx\x1b[0m"
                 # Store the rejected step.  In case the next step is identical to the rejected one, the next step should be accepted to avoid infinite loops.
                 self.X_rj = self.X.copy()

geometric/params.py

@@ -56,6 +56,8 @@ def __init__(self, **kwargs):
         self.enforce = kwargs.get('enforce', 0.0)
         # Small eigenvalue threshold
         self.epsilon = kwargs.get('epsilon', 1e-5)
+        # Energy increase threshold; 3e-2 -> about 20 kcal/mol
+        self.erisemax = kwargs.get('erisemax', 3e-2)
         # Interval for checking the coordinate system for changes
         self.check = kwargs.get('check', 0)
         # More verbose printout
@@ -222,7 +224,6 @@ def printInfo(self):
         elif self.hessian.startswith('file+last:'):
             logger.info(' Hessian data will be read from file: %s, then computed for the last step.\n' % self.hessian[5:])
 
-
 class NEBParams(object):
     """
     Container for optimization parameters.
@@ -256,6 +257,32 @@ def __init__(self, **kwargs):
         self.trust = max(self.tmin, self.trust)
         self.xyzout = kwargs.get('xyzout', None)
 
+class InterpParams(object):
+    """
+    Container for interpolation parameters.
+    """
+    def __init__(self, **kwargs):
+        # Whether we are optimizing a given trajectory.
+        self.optimize = kwargs.get('optimize', False)
+        if self.optimize:
+            self.nframes = 0
+            self.align_frags = False
+        else:
+            # Number of frames that will be used for the interpolation.
+            self.nframes = kwargs.get('nframes', 50)
+            # Whether we want to align the molecules in reactant and product frames.
+            self.align_frags = kwargs.get('align_frags', False)
+        # Whether we want to align the initial interpolate trajectory before optimization.
+        self.align_system = kwargs.get('align_system', True)
+        # Run in fast mode
+        self.fast = kwargs.get('fast', False)
+        # Verbose printout
+        self.verbose = kwargs.get('verbose', 0)
+        if 'extrapolate' in kwargs:
+            self.extrapolate = tuple(kwargs['extrapolate'])
+        else:
+            self.extrapolate = None
+
 def str2bool(v):
     """ Allows command line options such as "yes" and "True" to be converted into Booleans. """
     if isinstance(v, bool):
@@ -373,6 +400,7 @@ def parse_optimizer_args(*args):
     grp_optparam.add_argument('--reset', type=str2bool, help='Reset approximate Hessian to guess when eigenvalues are under epsilon.\n '
                               'Defaults to True for minimization and False for transition states.\n ')
     grp_optparam.add_argument('--epsilon', type=float, help='Small eigenvalue threshold for resetting Hessian, default 1e-5.\n ')
+    grp_optparam.add_argument('--erisemax', type=float, help='For energy minimization, the maximum value that energy can rise without step being rejected.\n ')
     grp_optparam.add_argument('--check', type=int, help='Check coordinates every <N> steps and rebuild coordinate system, disabled by default.\n')
     grp_optparam.add_argument('--subfrctor', type=int, help='Project out net force and torque components from nuclear gradient.\n'
                               '0 = never project; 1 = auto-detect (default); 2 = always project.')
@@ -516,3 +544,45 @@ def parse_neb_args(*args):
         args_dict['engine'] = 'tera'
 
     return args_dict
+
+def parse_interpolate_args(*args):
+
+    """
+    Read user input from the command line interface.
+    Designed to be called by interpolate.main() passing in sys.argv[1:]
+    """
+
+    parser = ArgumentParserWithFile(add_help=False, formatter_class=argparse.RawTextHelpFormatter, fromfile_prefix_chars='@')
+
+    grp_univ = parser.add_argument_group('universal', 'Relevant to every job')
+    grp_univ.add_argument('input', type=str, help='REQUIRED positional argument: xyz file containing reactant and product structures.\n')
+    grp_univ.add_argument('--optimize', type=str2bool, help='Provide "yes" to optimize an interpolated trajectory.\n'
+                                                            'The following two arguments, nframes and prealign, will be ignored.\n'
+                                                            'The input xyz file must contain 5 or more structures.\n')
+    grp_univ.add_argument('--nframes', type=int, help='Number of frames, needs to be bigger than 5, that will be used to interpolate, default 50\n')
+    grp_univ.add_argument('--fast', type=str2bool, help='Provide "yes" to run in fast mode (no splicing).\n')
+    grp_univ.add_argument('--align_frags', type=str2bool, help='Provide "yes" to align fragments in reactant and product structure prior to interpolation.\n')
+    grp_univ.add_argument('--align_system', type=str2bool, help='Provide "yes" to align the entire system prior to interpolation.\n')
+    grp_univ.add_argument('--extrapolate', type=int, nargs=2, help='Provide two integers to generate extrapolated frames at the head and tail ends.\n'
+                          'The displacement norm between extrapolated frames is the same as for interpolated frames.\n')
+
+    grp_output = parser.add_argument_group('output', 'Control the format and amount of the output')
+    grp_output.add_argument('--prefix', type=str, help='Specify a prefix for log file and temporary directory.\n'
+                            'Defaults to the input file path (incl. file name with extension removed).\n ')
+    grp_output.add_argument('--verbose', type=int, help='Set to positive for more verbose printout.\n'
+                            '0 = Default print level.     1 = Basic info about optimization step.\n'
+                            '2 = Include microiterations. 3 = Lots of printout from low-level functions.\n ')
+
+    grp_help = parser.add_argument_group('help', 'Get help')
+    grp_help.add_argument('-h', '--help', action='help', help='Show this help message and exit')
+
+    # Keep all arguments whose values are not None, so that the setting of default values
+    # in InterpParams() and get_molecule_engine() will work properly.
+    args_dict = {}
+    for k, v in vars(parser.parse_args(*args)).items():
+        if v is not None:
+            args_dict[k] = v
+
+    print("LPW debug:", args_dict)
+
+    return args_dict

geometric/tests/test_openmm.py

@@ -19,8 +19,10 @@ def test_dlc_openmm_water3(localizer):
     Optimize the geometry of three water molecules using standard delocalized internal coordinates.
     The coordinate system will break down and have to be rebuilt.
     """
-    progress = geometric.optimize.run_optimizer(engine='openmm', pdb=os.path.join(datad,'water3.pdb'), coordsys='dlc', input='tip3p.xml',
-                                                converge=['gmax', '1.0e-5'], check=50)
+    # LPW 2023-09-21 Increased check from 50 to 100; this was "perturbed" by the improvement to InternalCoordinates.newCartesian() .
+    progress = geometric.optimize.run_optimizer(engine='openmm', pdb=os.path.join(datad,'water3.pdb'), coordsys='dlc', input='tip3p.xml', 
+                                                converge=['gmax', '1.0e-5'], check=100) 
+
     # The results here are in Angstrom
     #
     ref = np.array([[ 1.19172917, -1.71174316,  0.79961878],
@@ -46,7 +48,8 @@ def test_dlc_openmm_water3(localizer):
     xdiff = (progress.xyzs[-1] - ref).flatten()
     rmsd, maxd = geometric.optimize.calc_drms_dmax(progress.xyzs[-1], ref, align=True)
     rmsd2, maxd2 = geometric.optimize.calc_drms_dmax(progress.xyzs[-1], ref2, align=True)    
-    print("RMS / Max displacement from reference:", rmsd, maxd)
+    print("RMS / Max displacement from reference 1:", rmsd, maxd)
+    print("RMS / Max displacement from reference 2:", rmsd2, maxd2)
     # This test is a bit stochastic and doesn't converge to the same minimized geometry every time.
     # Check that the energy is 0.01 a.u. above reference. Not really the qm_energy, this is a misnomer
     assert progress.qm_energies[-1] < (e_ref + 0.01)
@@ -62,14 +65,14 @@ def test_dlc_openmm_water12(localizer):
     The coordinate system is expected to break down and the optimizer should skip the optimization step
     after rebuilding the coordinate system.
     """
-    progress = geometric.optimize.run_optimizer(engine='openmm', pdb=os.path.join(datad,'water12.pdb'), 
+    progress = geometric.optimize.run_optimizer(engine='openmm', pdb=os.path.join(datad,'water12.pdb'),
                                                 coordsys='dlc', input='tip3p.xml', maxiter=20, converge=['maxiter'])
-
-    have_skip_step = False
-    for line in open('tip3p.log').readlines():
-        if 'Skipping optimization step' in line:
-            have_skip_step = True
-    assert have_skip_step
+    # LPW 2023-09-21: The coordinate system no longer breaks down after the improvement to newCartesian().
+    # have_skip_step = False
+    # for line in open('tip3p.log').readlines():
+    #     if 'Skipping optimization step' in line:
+    #         have_skip_step = True
+    # assert have_skip_step
 
 
 @addons.using_openmm

geometric/xml_helper.py

@@ -18,9 +18,16 @@ def read_coors_from_xml(M, state_xml):
                 M.xyzs[0][aid, 2] = z*10.0
     return root;
 
+def write_coors_to_xml(M, state_xml, out_xml):
+    if type(state_xml) is str:
+        with open(state_xml) as fobj:
+            xml = fobj.read()
+        root = ET.fromstring(xml)
+    elif type(state_xml) is ET.Element:
+        root = copy.deepcopy(state_xml)
+    else:
+        raise IOError("Expected second argument to write_coors_to_xml to be file name or xml.etree.ElementTree.Element type")
 
-def write_coors_to_xml(M, root_template, filename):
-    root = copy.deepcopy(root_template)   
     for child in root:
         if child.tag == "Positions":
             for aid in range(M.na):
@@ -30,8 +37,8 @@ def write_coors_to_xml(M, root_template, filename):
                 child[aid].attrib['x'] = "%.16e" % x
                 child[aid].attrib['y'] = "%.16e" % y
                 child[aid].attrib['z'] = "%.16e" % z
-    # write
-    fout = open(filename, "w")
+    # write to the output XML file.
+    fout = open(out_xml, "w")
     print(ET.tostring(root).decode("utf-8"), file=fout)
     fout.close()
 

setup.py

@@ -18,8 +18,9 @@
     long_description_content_type="text/markdown",
     entry_points={'console_scripts': [
         'geometric-optimize = geometric.optimize:main',
-        'run-ase = geometric.ase_engine:main',
-	    'geometric-neb = geometric.neb:main',
+        'geometric-neb = geometric.neb:main',
+        'geometric-interpolate = geometric.interpolate:main',
+        'run-ase = geometric.ase_engine:main'
     ]},
     install_requires=[
         'numpy>=1.11',