.. tip:: All input files can be downloaded: :download:`Files `. opt ====== .. contents:: :local: This keyword controls the search of geometry minimum, transistion state, and reaction path. Options ------------ .. option:: type .. list-table:: :stub-columns: 1 :widths: 5 20 * - Value - ``Min`` Search the geometry minimum. * - - ``NEB`` Use the nudged elastic band (NEB) algorithm to search the reaction path and transition state. * - - ``Dimer`` Use the dimer algorithm to search the transition state given the reactant and product geometry. * - - ``QST2`` Use the QST2 algorithm to search the transition state given the reactant and product geometry. * - - ``TS`` Search the transition state from a single structure. * - Default - ``Min`` Determine what type of geometry optimization you want to do. For ``Min`` and ``TS``, the initial structure should be given in ``mol``. The optimization process is output to ``x-opt-traj.xyz`` and the final minimum or transition state is output to ``x-opt.xyz``. For ``NEB``, ``Dimer``, and ``QST2``, 2 structures have to be given in ``mol`` and ``mol2``, which represent a reactant and product pose, respectively. Usually, one can first use ``NEB`` to rapidly find a reasonable path and transition state. If ``NEB`` is hard to converge, then use the structures from ``NEB`` result to do a ``Dimer`` or ``QST2`` search. ``Dimer`` and ``QST2`` require reactant and product structures of high quality. For ``NEB``, the reaction path is output to ``x-opt-traj.xyz`` and the final transition state is output to ``x-opt.xyz``. For ``Dimer`` and ``QST2``, the transition state is output to ``x-opt.xyz``. The ``x-opt-traj.xyz`` is **NOT** reaction path! It is just the optimization process as ``Min`` or ``TS``. .. option:: max_step .. list-table:: :stub-columns: 1 :widths: 5 20 * - Value - An integer * - Default - ``1000`` The maximum number of geometry optimization steps. .. option:: energy_cov .. list-table:: :stub-columns: 1 :widths: 5 20 * - Value - A real number * - Default - ``1.E-4`` The energy convergence threshold. When the energy change is smaller than this value, this energy condition is satisfied. .. option:: grad_cov .. list-table:: :stub-columns: 1 :widths: 5 20 * - Value - A real number * - Default - ``1.E-3`` The gradient convergence threshold. This actually determines 4 convergence thresholds: .. list-table:: * - Maximum gradient component - ``grad_cov`` * - RMS gradient: - ``grad_cov`` * 2/3 * - Maximum atomic displacement - ``grad_cov`` * 4 * - RMS atomic displacement - ``grad_cov`` * 8/3 When all these 4 conditions are met, this gradient condition is satisfied. .. option:: max_dr .. list-table:: :stub-columns: 1 :widths: 5 20 * - Value - A real number * - Default - ``0.5`` The maximum atomic displacement in an optimization step. If the molecule is highly flexible (Mathematically, the potential energy surface is very flat), or the structure (especially transition state) is very close to the stationary point but not converge, setting a smaller ``max_dr`` like ``0.1`` is very useful. .. option:: num_images .. list-table:: :stub-columns: 1 :widths: 5 20 * - Value - An integer * - Default - ``10`` The number of images for NEB transition state search. This number canNOT be set **too small**, say, ``5``. .. option:: neb_k .. list-table:: :stub-columns: 1 :widths: 5 20 * - Value - A real number * - Default - ``0.01`` The force constant for NEB transition state search. For a specific system, the optimal number of ``neb_k`` should be chose by trail-and-error. .. option:: fix_atoms .. list-table:: :stub-columns: 1 :widths: 5 20 * - Value - Atom range * - Default - None Assign the atoms that are fixed during geometry optimization. For example: .. code-block:: bash :linenos: opt fix_atom 2 5-9 23 26 end The atoms 2,5,6,7,8,9,23,26 will be fixed during geometry optimization. .. option:: fix_bond .. list-table:: :stub-columns: 1 :widths: 5 20 * - Value - 2 integers * - Default - None Assign the bond that are fixed during geometry optimization. For example: .. code-block:: bash :linenos: opt fix_bond 1 4 fix_bond 2 6 end The bonds (1,4) and (2,6) will be fixed during geometry optimization. .. option:: fix_angle .. list-table:: :stub-columns: 1 :widths: 5 20 * - Value - 3 integers * - Default - None Assign the angle that are fixed during geometry optimization. For example: .. code-block:: bash :linenos: opt fix_angle 1 4 5 fix_angle 2 6 7 end The angles (1, 4, 5) and (2, 6, 7) will be fixed during geometry optimization. .. option:: fix_torsion .. list-table:: :stub-columns: 1 :widths: 5 20 * - Value - 4 integers * - Default - None Assign the torsion that are fixed during geometry optimization. For example: .. code-block:: bash :linenos: opt fix_torsion 1 4 5 9 fix_torsion 2 6 7 12 end The torsions (1, 4, 5, 9) and (2, 6, 7, 12) will be fixed during geometry optimization. .. attention:: Currently, the transition state search algorithm ``QST2`` and ``TS`` do **NOT** support constraints. If you want to search a transition state with constraints, please use ``NEB`` or ``Dimer``. Theoretical Background -------------------------------- Minimum ^^^^^^^^^^^^^^^^ The minimum is defined as a stable isomer on its potential energy surface (PES) of a molecule. The gradients on all atoms are zero. The optimization of minimum depends strongly on the initial structure. **For different starters, one can get different isomers.** Transition State ^^^^^^^^^^^^^^^^^^^^^^ Transition state is a short-lived configuration of atoms that in maximum on one direction but minimum on other directions. The gradients on all atoms are also zero. In Qbics, one can use NEB or dimer method to search the transition state. Only (unoptimized) reactant and product structures are needed. No exact Hessian needs to be computed. A good strategy is: - Use cheap method, like xTB, to find a reasonable path and transition state with NEB (``type neb``). - Then, use standard DFT method to refine the transition state with dimer (``type dimer``) or QST2 (``type QST2``), even the previous NEB result is not converged. This strategy is shown below: .. figure:: figs/optk.jpg Input Examples -------------------- Example: Minimum Structure of Aspirin ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ Search the minimum structure of aspirin at B3LYP/def2-SVP level of theory: .. code-block:: bash :linenos: :caption: opt-1.inp basis def2-svp end scf charge 0 spin2p1 1 end mol O 1.23330 0.55400 0.77920 O -0.69520 -2.71480 -0.75020 O 0.79580 -2.18430 0.86850 O 1.78130 0.81050 -1.48210 C -0.08570 0.60880 0.44030 C -0.79270 -0.55150 0.12440 C -0.72880 1.84640 0.41330 C -2.14260 -0.47410 -0.21840 C -2.07870 1.92380 0.07060 C -2.78550 0.76360 -0.24530 C -0.14090 -1.85360 0.14770 C 2.10940 0.67150 -0.31130 C 3.53050 0.59960 0.16350 H -0.18510 2.75450 0.65930 H -2.72470 -1.36050 -0.45640 H -2.57970 2.88720 0.05060 H -3.83740 0.82380 -0.50900 H 3.72900 1.41840 0.85930 H 4.20450 0.69690 -0.69240 H 3.71050 -0.36590 0.64260 H -0.25550 -3.59160 -0.73370 end task opt b3lyp end Example: Minimum Structure with Constraints ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ Search the minimum structure of a molecule "1UML" with a bond and a torsion fixed at xTB level of theory: .. tabs:: .. tab:: opt-2.inp .. code-block:: bash :linenos: :caption: opt-2.inp opt fix_bond 7 51 fix_torsion 2 9 12 13 end mol 1uml.xyz end task opt xtb end .. tab:: 1uml.xyz .. code-block:: bash :linenos: :caption: 1uml.xyz 60 1UMl(0) C 0.38066735 2.55219474 1.46267636 N -0.33789766 3.82795956 1.53541195 C -1.09121603 3.87593593 2.64608144 N -0.94245652 2.73117700 3.32171417 C -0.03158842 1.89319979 2.55453525 C 0.37212443 0.52262484 2.93349994 O 0.95733416 -0.14814113 2.10415165 N 0.07166899 0.06707847 4.15265909 C -0.24927084 4.91745505 0.56953160 C 0.93256260 5.84926075 0.89138538 O 0.85638839 6.40073825 2.18023567 C -0.21178237 4.41479892 -0.88585063 C -1.30451123 3.40488529 -1.25181254 N -0.92687172 2.67905416 -2.46663252 C -1.02012131 3.09677073 -3.75965819 C -0.51157942 2.12142921 -4.62384570 C -0.09088312 1.03387236 -3.77904056 C -0.36650560 1.41729423 -2.43324566 C 0.47064744 -0.18584082 -4.05073192 C 0.77451739 -1.04742178 -2.99644766 C 0.51501197 -0.67068355 -1.66606604 C -0.06269010 0.57406126 -1.40267529 N 0.81935074 -1.45803325 -0.60027087 C 1.28805430 -2.72903213 -0.64795569 O 1.47475540 -3.39721237 -1.66461437 C 1.57993549 -3.25675705 0.74748896 C 1.64225207 -4.78720972 0.77417719 C -0.79428279 -4.92640828 1.45233103 C 0.27286077 -5.43783743 0.71500759 C 0.09621230 -6.58712911 -0.05216185 C -1.13149600 -7.23287171 -0.06650889 C -2.18992415 -6.73793686 0.69577410 C -2.02072844 -5.58468177 1.44785223 H 1.08446564 2.21827120 0.70020968 H -1.72176781 4.71224102 2.94804163 H -0.42535733 0.67320222 4.82086735 H 0.33691344 -0.89030993 4.42449590 H -1.16218006 5.52526169 0.70631806 H 0.98220194 6.64417179 0.12373872 H 1.88374201 5.29787879 0.83772723 H 0.29674389 7.26495989 2.15187259 H 0.77353707 3.94775022 -1.06475242 H -0.26304168 5.28737023 -1.55851059 H -2.24791634 3.94699892 -1.42947786 H -1.48837270 2.67328764 -0.45339564 H -1.43235349 4.05491839 -4.07609186 H -0.44710423 2.17375337 -5.71065901 H 0.67778198 -0.47993189 -5.07960408 H 1.21706689 -2.02142178 -3.20534792 H -0.27079576 0.87203324 -0.37507411 H 0.68015960 -1.04879454 0.33459425 H 0.94277827 -2.86688017 1.55772521 H 2.59161906 -2.87935288 0.99367456 H 2.12787894 -5.07775312 1.72353104 H 2.29189074 -5.17230270 -0.02720370 H -0.66818790 -4.01186770 2.03180634 H 0.92422969 -6.97980706 -0.64230710 H -1.26922150 -8.12751976 -0.67371391 H -3.14853209 -7.25672910 0.70037925 H -2.85015656 -5.19299484 2.03675744 Check the constraints during optimization: .. figure:: figs/opt-2.gif Example: Transition State of S\ :sub:`N`\ 2 Reaction ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ Search the transion state of an S\ :sub:`N`\ 2 reaction with NEB algorithm at B3LYP/def2-SVP level of theory: .. code-block:: bash :linenos: :caption: ts-1.inp basis # Define basis set. def2-svp end opt type NEB # Type: Min, NEB, Dimer, QST2, TS num_images 10 # The number of images for NEB calculations. neb_k 0.01 # The force constant for NEB calculations. end scf charge -1 # The net charge. spin2p1 1 # 2S+1 end xtb chrg -1 end mol C -2.25147439 4.89406277 -1.00469981 H -1.89481996 3.88525277 -1.00469981 H -1.89480154 5.39846096 -0.13104831 H -3.32147439 4.89407596 -1.00469981 Cl -1.66479756 5.72372709 -2.44173406 Cl -2.67350651 4.09697871 0.73250622 end mol2 C -2.36845504 4.69197207 -0.60149770 H -1.76657311 4.00286639 -1.15626927 H -1.80200132 5.57659799 -0.39786281 H -3.23625780 4.94775799 -1.17280492 Cl -1.66479756 5.72372709 -2.44173406 Cl -2.86278952 3.94963672 0.91579319 end task opt b3lyp # opt xtb # You can also try this. end The reaction path is given in ``ts-1-opt-traj.xyz``: .. figure:: figs/ts-1.gif The energies can be found in the output file ``ts-1.out``: .. code-block:: bash :linenos: :caption: ts-1.out NEB path updated in "ts-1-opt-traj.xyz": ---------------------------------------------------- # Energy Dist Gtang Gperp ---------------------------------------------------- 0 -960.06873748 0.13968 0.00000 0.00000 1 -960.06864683 0.11091 0.00029 0.00035 2 -960.06738628 0.07888 0.00032 0.00028 3 -960.06508634 0.07078 0.00008 0.00030 4 -960.06240481 0.09854 -0.00028 0.00036 5 -960.05912811 0.18516 -0.00087 0.00032 6 -960.05781339 0.21042 -0.00025 0.00051 7 -960.06424460 0.26756 -0.00057 0.00032 8 -960.06880165 0.00000 0.00000 0.00000 9 -960.05738746 0.06067 0.00014 0.00026 ---------------------------------------------------- Geometry optimization step 34: Current energy: -960.05738746 Delta Energy: 8.34686E-08; Target: 1.00000E-04; Converged? Yes Max displacement: 2.30167E-04; Target: 4.00000E-03; Converged? Yes RMS displacement: 1.07545E-04; Target: 2.66667E-03; Converged? Yes Max gradient: 5.45965E-04; Target: 1.00000E-03; Converged? Yes RMS gradient: 2.56874E-04; Target: 6.66667E-04; Converged? Yes Stationary point has reached. In the table, structure ``0`` and ``1`` are the reactant and product, respectively, and structure ``6`` is the transition state, which is also given in ``ts-1-opt.xyz``. You can change ``b3lyp`` to ``xtb`` to perform the calculation in a much more rapid way. You can also try ``type dimer``, which is computationally cheaper than ``NEB`` but requires reactant and product structures must be of high quality. Example: Transion State of Decarboxylation Reaction ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ Search the transion state of the following decarboxylation reaction with NEB algorithm at B3LYP/def2-SVP level of theory: .. tabs:: .. tab:: ts-2.inp .. code-block:: bash :linenos: :caption: ts-2.inp mol a1.xyz end mol2 a2.xyz end scf charge -1 spin2p1 1 end opt type neb # You can also try "dimer". num_images 15 neb_k 0.01 end basis def2-svp end task opt b3lyp end .. tab:: a1.xyz .. code-block:: bash :linenos: :caption: a1.xyz 16 E = -587.26844619 Hartree; coordinates in Angstrom. C 0.04106534 -2.18782379 0.01873379 C -1.35798345 -2.26546988 0.01269198 C -2.17792255 -1.12924287 -0.02333896 C -1.53026482 0.10833393 -0.03797859 C -0.12071906 0.20403755 -0.02314596 C 0.67915283 -0.93908160 -0.00278636 C 0.47143440 -3.57161137 0.04441853 H -3.26722478 -1.21160460 -0.04107296 H -2.12778938 1.02510522 -0.06470926 H 0.34489770 1.19455049 -0.03304010 H 1.77014421 -0.92296985 -0.00764617 O -1.74251880 -3.55390987 0.03659370 N -0.56119829 -4.37040148 0.05822977 C 1.96372187 -4.05774618 0.02851820 O 2.15749685 -5.25776354 0.25204089 O 2.75531543 -3.11946977 -0.21744012 .. tab:: a2.xyz .. code-block:: bash :linenos: :caption: a2.xyz 16 E = -34.49382101 Hartree; coordinates in Angstrom. C -0.31300710 -2.30106195 -0.01428127 C -1.70090129 -2.13170314 -0.01142600 C -2.23889881 -0.83944478 -0.02217974 C -1.41495826 0.23848161 -0.03526503 C -0.02069520 0.06834558 -0.03813291 C 0.51990509 -1.17609015 -0.02785645 C 0.29527044 -3.71579140 -0.00257527 H -3.29999458 -0.70170793 -0.02007771 H -1.83005138 1.22465134 -0.04345765 H 0.61932516 0.92576433 -0.04849442 H 1.58298567 -1.29755818 -0.03013007 O -2.48723821 -3.16847622 0.00114316 N 0.74811153 -4.76914425 0.00614076 C 2.27925143 -4.13711065 0.09555385 O 1.94633923 -5.30045357 0.44105892 O 2.61216378 -2.97376825 -0.24995180 Structures in ``a1.xyz`` and ``a2.xyz`` are shown below. They are put arbitrarily without optimization: .. tabs:: .. tab:: a1.xyz .. figure:: figs/a1.jpg .. tab:: a2.xyz .. figure:: figs/a2.jpg The optimized transtion state ``ts-2-opt.xyz`` and path ``ts-2-opt-traj.xyz`` are shown below: .. tabs:: .. tab:: ts-2-opt.xyz .. figure:: figs/a.jpg .. tab:: ts-2-opt-traj.xyz .. figure:: figs/ts-2.gif Example: Transion State of Pd(OAc)-Catalyzed Nucleopalladation ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ Search the transion state of the following Pd(OAc)-catalyzed nucleopalladation with DIMER algorithm at B3LYP/def2-SVP level of theory: .. tabs:: .. tab:: ts-3.inp .. code-block:: bash :linenos: :caption: ts-3.inp mol a1.xyz end mol2 a2.xyz end scf charge 0 spin2p1 1 end opt type dimer # You can also use qst2 end basis def2-svp end pseudopotential def2-ecp # Since Pd is used, you need to use ECP. end task opt b3lyp end .. tab:: a1.xyz .. code-block:: bash :linenos: :caption: a1.xyz 43 element x y z Pd 0.24654 -1.25332 0.40518 N 2.21910 -1.24512 -0.02293 O 4.02373 -2.55737 -0.59458 N 0.71331 0.97790 0.25853 C 2.90583 -2.40807 -0.13902 C 4.14710 0.14624 -0.66749 H 4.76219 -0.71420 -0.85891 C 1.86165 3.54336 0.19530 H 2.32306 4.52041 0.16979 C 2.82401 -0.02669 -0.28507 C 0.56127 3.36146 0.56593 H -0.05640 4.19293 0.87529 C 0.02466 2.06497 0.55234 H -1.01171 1.95668 0.81328 C 2.10089 -3.60617 0.39103 H 2.29784 -3.68424 1.46574 H 2.50312 -4.49514 -0.09449 C 3.99404 2.53651 -0.49534 H 4.44575 3.51746 -0.52470 C 0.63638 -3.41684 0.12343 C 2.02781 1.13223 -0.06025 C 2.63885 2.42073 -0.14140 C 4.72289 1.41237 -0.76883 H 5.76763 1.48006 -1.03642 C -0.35167 -3.32488 1.08900 H -1.30757 -3.58022 0.86748 H -0.08593 -3.27603 2.15625 H 0.32333 -3.63444 -0.90750 N -1.81017 -1.32804 0.64184 O -2.91585 -2.79869 -0.74520 C -2.83467 -1.72049 -0.17802 C -2.07940 -0.18984 1.24478 O -1.43547 0.22090 2.19254 C -3.75316 -0.59348 -0.30558 C -3.63015 1.69613 0.35364 H -3.18847 2.48978 0.92320 C -4.76410 -0.32523 -1.20052 H -5.18473 -1.11145 -1.79977 C -3.23750 0.40158 0.51050 C -4.56763 1.99386 -0.63755 H -4.81215 3.01413 -0.85747 C -5.15401 0.98013 -1.36177 H -5.89024 1.23011 -2.11222 .. tab:: a2.xyz .. code-block:: bash :linenos: :caption: a2.xyz 43 element x y z Pd 0.26512 -0.64711 0.61399 N 2.06999 -0.92815 -0.05458 O 3.65856 -2.55368 -0.43497 N 0.99509 1.45229 0.64166 C 2.54673 -2.20957 -0.08507 C 4.14972 0.20649 -0.71873 H 4.62578 -0.72848 -0.96208 C 2.30494 3.86864 0.23224 H 2.82051 4.80502 0.07331 C 2.84319 0.19559 -0.25035 C 1.03545 3.83626 0.74375 H 0.50815 4.74282 0.99632 C 0.42408 2.59308 0.95757 H -0.55466 2.53340 1.41954 C 1.49284 -3.20295 0.42350 H 1.66001 -3.32238 1.50139 H 1.64757 -4.17628 -0.05294 C 4.28051 2.61430 -0.57548 H 4.82459 3.53688 -0.70941 C 0.10935 -2.62297 0.15412 C 2.24272 1.45031 0.10134 C 2.96211 2.66317 -0.08830 C 4.85325 1.40272 -0.86406 H 5.87106 1.36062 -1.22258 C -1.07037 -3.17687 0.95906 H -1.41934 -4.15247 0.60366 H -0.83933 -3.23046 2.02995 H -0.13100 -2.69707 -0.92522 N -2.14266 -2.22101 0.73663 O -3.47014 -3.49326 -0.64534 C -3.14336 -2.41981 -0.19876 C -1.88727 -0.84977 0.97034 O -1.25091 -0.41692 1.95149 C -3.68247 -1.07281 -0.49566 C -3.19644 1.22329 0.10017 H -2.60860 1.95263 0.63571 C -4.76405 -0.68826 -1.26341 H -5.36566 -1.42688 -1.76954 C -2.89679 -0.12547 0.17462 C -4.28294 1.60635 -0.67387 H -4.53744 2.65330 -0.74569 C -5.05313 0.66777 -1.34699 H -5.89479 0.99766 -1.93734 Structures in ``a1.xyz`` and ``a2.xyz`` are shown below. They are put arbitrarily without optimization: .. tabs:: .. tab:: a1.xyz .. figure:: figs/a12.jpg .. tab:: a2.xyz .. figure:: figs/a22.jpg The optimized transtion state is ``ts-3-opt.xyz``, shown below: .. figure:: figs/ats2.jpg