Merge pull request #216 from hjnpark/master

minor changes in neb.py, params.py, README.md. Documentation is updated.

README.md

@@ -7,9 +7,9 @@
 
 This is a geometry optimization code for molecular structures.
 The code works by calling external software for the energy and 
-gradient through wrapper functions.  Q-Chem, TeraChem, Psi4, 
-Molpro, Gaussian 09/16, and CFOUR are supported quantum chemistry
-codes through the command line interface.  The PySCF and
+gradient through wrapper functions. Q-Chem, TeraChem, Psi4, 
+Molpro, Gaussian 09/16, CFOUR, and QUICK are supported quantum chemistry
+codes through the command line interface. The PySCF and
 QCArchive packages also provide interfaces to geomeTRIC for
 optimization. MM optimizations using OpenMM and Gromacs are
 also supported through the command line interface.

docs/source/how-it-works.rst

@@ -196,7 +196,7 @@ This is designed to minimize the number of times the IC step is converted to Car
 The above image shows the relationship between :math:`r_{IC}` (x-axis) and the RMSD (left y-axis; blue) for the structure shown (the taxol example at the 10th optimization cycle).
 The plot is generated by scanning the value of :math:`r_{IC}`, optimizing :math:`\lambda` for each value, the calculating :math:`\boldsymbol{\delta}^{(x)}` and the RMSD.
 
-The optimization of :math:`\lambda` for a target step size (i.e. :math:`|\boldsymbol{\delta}^{(q)}| \rightarrow r_{IC}`) follows `Hebden et al. <https://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.716.2997>`_
+The optimization of :math:`\lambda` for a target step size (i.e. :math:`|\boldsymbol{\delta}^{(q)}| \rightarrow r_{IC}`) follows `Hebden et al. <https://www.numerical.rl.ac.uk/media/reports/AERE_TP515.pdf>`_
 The derivative of the step size with respect to :math:`\lambda_n` at iteration :math:`n` is given by:
 
 .. math::

docs/source/options.rst

@@ -344,11 +344,18 @@ These options control the NEB method.
 
 ....
 
+``chain_coords``
+
+Name of the coordinate file containing multiple frames for NEB. This is a **required** positional argument along with the ``input`` file.
+It will override the molecular structure in the ``input`` file, which should contain the structure of the first image.
+
+....
+
 ``--maxg [0.05]``
 
 ``--avgg [0.025]``
 
-The convergence occurs when the maximum and average RMS-gradient of all images fall below these thresholds.
+The convergence occurs when the maximum and average RMS-gradient (in ev/Ang) of all images fall below these thresholds.
 
 ....
 
@@ -358,12 +365,6 @@ This value will be used to build the guessed Hessian for the input chain.
 
 ....
 
-``--neb_history [1]``
-
-Number of chain histories to save in memory. Note that each chain is memory intensive.
-
-....
-
 ``--neb_maxcyc [100]``
 
 This sets the maximum number of NEB iterations.
@@ -383,7 +384,7 @@ The climbing images will be selected in descending order from the highest energy
 
 This option determines force components that will be projected. The default value is ``0`` which is NEB.
 The other two available bands are hybrid (``1``) and plain (``2``) band.
-The details of the force components of each band can be found in the :ref:` NEB page <neb>`.
+The details of the force components of each band can be found in the :ref:`NEB page <neb>`.
 
 ....
 
@@ -408,14 +409,11 @@ The trust radius is the maximum allowed Cartesian displacement of an optimizatio
 
 Depending on the quality of individual optimization steps, the trust radius can be increased from its current value up to the ``--tmax`` value, or it can be decreased down to a minimum value.
 
-The minimum trust radius cannot be user-set; its value is ``0.0`` for transition state and MECI jobs, and the smaller of the ``drms`` convergence criteria and ``1.2e-3`` for energy minimizations.
-The purpose of the minimum trust radius is to prevent unpredictable behavior when the trust radius becomes extremely small (e.g. if the step is so small that the energy change is smaller than the SCF convergence criteria).
-
 ....
 
 ``--skip [yes/no]``
 
-Setting it to ``yes`` will skip Hessian updates that will introduce negative eigenvalues.
+Setting it to ``yes`` will skip Hessian updates that would introduce negative eigenvalues instead of resetting it. By default, it will reset the Hessian when negative Hessian eigenvalues are detected.
 
 ....
 

docs/source/transition.rst

@@ -18,7 +18,7 @@ The estimated TS structures generated by these approximate methods are sometimes
 
 At the start of a TS optimization in geomeTRIC, the full Hessian matrix is calculated by central finite difference of the gradient.
 The local calculation may take a long time as the number of gradients needed is :math:`6 \times N_{\mathrm{atoms}}`.
-Using the :ref:`Work Queue <installcctools>` library, you could parallelize the gradient calculations by distributing the calculations on remote machines.
+Using the :ref:`Work Queue <installcctools>` library or :ref:`BigChem <installbigchem>`, you could parallelize the gradient calculations by distributing the calculations on remote machines.
 Alternatively, if you have an exact or approximate Hessian matrix computed externally, you may provide it as a text file and skip the finite difference Hessian calculation.
 At the conclusion of the optimization, the user may optionally request an additional Hessian calculation and vibrational analysis to check the number of imaginary modes in the final optimized structure; in most cases, the desired number of imaginary modes is 1.
 
@@ -267,6 +267,22 @@ The worker node must have the QC software installed with the environment variabl
 If successful, the worker will establish a connection to the master and begin to accept gradient jobs.
 Parallelization is achieved by running multiple workers on one or more nodes (you can run workers locally too).
 
+Similarly, `BigChem <https://github.com/mtzgroup/bigchem>`_ can be used to parallelize the gradient calculations for the Hessian. First, the Redis server can be started on the head node::
+
+    redis-server --bind 0.0.0.0 --daemonize yes --logfile redis.log
+
+Before we start the workers, the environment variables need to be set::
+
+    export BIGCHEM_BROKER_URL="redis://your_head_node_hostname_or_ip/0"
+    export BIGCHEM_BACKEND_URL="redis://your_head_node_hostname_or_ip/0"
+
+Now the workers can be started on any computing nodes that can reach the head node::
+
+    celery -A bigchem.tasks worker --without-heartbeat --without-mingle --without-gossip --loglevel=INFO
+
+While both the server and workers are active, use the ``--bigchem yes`` flag to enable BigChem for parallel Hessian calculations.
+
+
 The Hessian calculation and vibrational analyses should give the same results as if you had requested them directly from the quantum chemistry code.
 After the vibrational analysis, the Gibbs free energy corrections are computed using an ideal gas / rigid rotor / harmonic oscillator approximation (imaginary frequency modes are ignored).
 The free energy calculation may be customized by passing ``--thermo <temp> <pres>`` and providing the temperature and pressure.

examples/constraints.txt

@@ -66,7 +66,7 @@
 #  
 
 $freeze
-bond 5 6
+distance 5 6
 xyz 5 xyz
 xy 5-11,13,35
 $set

geometric/neb.py

@@ -1323,7 +1323,7 @@ def evaluate(self, trial):
             return cnorm - self.target
 
 
-def recover(chain_hist, forceCart, result=None):
+def recover(chain_hist, result=None):
     """
     Recover from a failed optimization.
 
@@ -1332,9 +1332,6 @@ def recover(chain_hist, forceCart, result=None):
     chain_hist : list
         List of previous Chain objects;
         the last element is the current chain
-    forceCart : bool
-        Whether to use Cartesian coordinates or
-        adopt the IC system of the current chain
     result : dict
         Dictionary with energies and gradients
 
@@ -1380,7 +1377,6 @@ def BFGSUpdate(Y, old_Y, G, old_G, H, params):
     ndg = np.array(Dg).flatten() / np.linalg.norm(np.array(Dg))
     nhdy = np.dot(H, Dy).flatten() / np.linalg.norm(np.dot(H, Dy))
     if verbose:
-        # HP: 2023-2-15: Not sure what is nhdy is for. I changed np.array(H*dy) to np.dot(H, dy)
         logger.info("Denoms: %.3e %.3e \n" % ((Dg.T * Dy)[0, 0], (Dy.T * H * Dy)[0, 0]))
         logger.info("Dots: %.3e %.3e \n" % (np.dot(ndg, ndy), np.dot(ndy, nhdy)))
     H += Mat1 - Mat2
@@ -1399,7 +1395,6 @@ def updatehessian(old_chain, chain, HP, HW, Y, old_Y, GW, old_GW, GP, old_GP, La
     """
     This function updates the Hessians and check their eigenvalues.
     """
-    H_reset = False
     HP_bak = HP.copy()
     HW_bak = HW.copy()
     BFGSUpdate(Y, old_Y, GP, old_GP, HP, params)
@@ -1413,7 +1408,6 @@ def updatehessian(old_chain, chain, HP, HW, Y, old_Y, GW, old_GW, GP, old_GP, La
             HP = HP_bak.copy()
             HW = HW_bak.copy()
         else:
-            H_reset = True
             logger.info(
                 "Eigenvalues below %.4e (%.4e) - will reset the Hessian \n"
                 % (params.epsilon, np.min(Eig1))
@@ -1422,7 +1416,7 @@ def updatehessian(old_chain, chain, HP, HW, Y, old_Y, GW, old_GW, GP, old_GP, La
 
     del HP_bak
     del HW_bak
-    return chain, Y, GW, GP, HP, HW, old_Y, old_GP, old_GW, H_reset
+    return chain, Y, GW, GP, HP, HW, old_Y, old_GP, old_GW
 
 
 def qualitycheck(old_chain, new_chain, trust, Quality, ThreLQ, ThreRJ, ThreHQ, Y, GW, GP, old_Y, old_GW, old_GP, params_tmax):
@@ -1711,7 +1705,7 @@ def OptimizeChain(chain, engine, params):
         # |      Update the Hessian Matrix      |#
         # =======================================#
 
-        chain, Y, GW, GP, HP, HW, Y_prev, GP_prev, GW_prev, _ = updatehessian(chain_prev, chain, HP, HW, Y, Y_prev, GW,
+        chain, Y, GW, GP, HP, HW, Y_prev, GP_prev, GW_prev = updatehessian(chain_prev, chain, HP, HW, Y, Y_prev, GW,
                                                                  GW_prev, GP, GP_prev, LastForce, params, None)
     return chain, optCycle
 

geometric/params.py

@@ -482,11 +482,12 @@ def parse_neb_args(*args):
     grp_nebparam = parser.add_argument_group('nebparam', 'Control the NEB calculation')
     grp_nebparam.add_argument('--maxg', type=float, help='Converge when the maximum RMS-gradient of all images falls below this threshold (default 0.05 ev/Ang).\n ')
     grp_nebparam.add_argument('--avgg', type=float, help='Converge when the average RMS-gradient of all images falls below this threshold (default 0.025 ev/Ang).\n ')
-    grp_nebparam.add_argument('--guessk', type=float, help='Guess Hessian eigenvalue for displacements (default 0.05).\n ')
+    grp_nebparam.add_argument('--guessk', type=float, help='Guess the initial Hessian eigenvalue for displacements (default 0.05).\n ')
     #HP 5/10/2024: guessw will be enabled once IC NEB is implemented.
     #grp_nebparam.add_argument('--guessw', type=float, help='Guess weight for chain coordinates (default 0.1).\n ')
     grp_nebparam.add_argument('--nebk', type=float, help='NEB spring constant in units of kcal/mol/Ang^2 (default 1.0).\n ')
-    grp_nebparam.add_argument('--neb_history', type=int, help='Chain history to keep in memory; note chains are very memory intensive, >1 GB each (default 1).\n ')
+    #HP 1/16/2025: neb_history is commented out because rebuilding the Hessian based on changes in IC isn't available.
+    #grp_nebparam.add_argument('--neb_history', type=int, help='Chain history to keep in memory; note chains are very memory intensive, >1 GB each (default 1).\n ')
     grp_nebparam.add_argument('--neb_maxcyc', type=int, help='Maximum number of chain optimization cycles to perform (default 100).\n ')
     grp_nebparam.add_argument('--climb', type=float, help='Activate climbing image for max-energy points when max RMS-gradient falls below this threshold (default 0.5).\n ')
     grp_nebparam.add_argument('--ncimg', type=int, help='Number of climbing images to expect (default 1).\n ')
@@ -497,7 +498,7 @@ def parse_neb_args(*args):
     grp_nebparam.add_argument('--trust', type=float, help='Starting trust radius (default 0.1). \n ')
     grp_nebparam.add_argument('--tmax', type=float, help='Maximum trust radius (default 0.3).\n ')
     grp_nebparam.add_argument('--tmin', type=float, help='Minimum trust radius, do not reject steps trust radius is below this threshold.\n ')
-    grp_nebparam.add_argument('--skip', type=str2bool, help='Skip Hessian updates that would introduce negative eigenvalues.\n ')
+    grp_nebparam.add_argument('--skip', type=str2bool, help='Setting it to ``yes`` will skip Hessian updates that would introduce negative eigenvalues instead of resetting it. By default, it will reset the Hessian when negative Hessian eigenvalues are detected.\n ')
     grp_nebparam.add_argument('--epsilon', type=float, help='Small eigenvalue threshold for resetting Hessian, default 1e-5.\n ')
     grp_nebparam.add_argument('--wqport', type=int, help='Work Queue port used to distribute singlepoint calculations. Workers must be started separately.\n ')
     grp_nebparam.add_argument('--bigchem', type=str2bool, help='Provide "Yes" to use BigChem for performing the NEB calculation in parallel. \n'

geometric/qcf_neb.py

@@ -350,7 +350,7 @@ def nextchain(info_dict):
                 % (params.epsilon, np.min(Eig))
             )
 
-            chain, Y, GW, GP, HW, HP = recover([chain_prev], LastForce, result_prev)
+            chain, Y, GW, GP, HW, HP = recover([chain_prev], result_prev)
             Ys = []
             GWs = []
             GPs = []