.. tip:: All input files can be downloaded: :download:`Files `. .. tip:: For more information of this section, please refer to these pages: - :doc:`../keywords/opt` - :doc:`../keywords/xtb` Geometry Optimization ========================================= .. contents:: :local: This tutorial will lead you step by step to do geometry optimization calculations using Qbics. Assume you have read :doc:`dft` or :doc:`charmm`. For optimization of transition states, please refer to :doc:`ts`. Geometry optimization is to find the minimum energy structure of a molecule. It is a common task in computational chemistry and materials. As long as you can do energy calculation, simply chaning ``energy`` to ``opt`` in the input file, you can do geometry optimization. Example: DFT for CH\ :sub:`3`\ CH\ :sub:`2`\ SO\ :sub:`3`\ H ----------------------------------------------------------------------- For DFT, we can do geometry optimization at B3LYP-D3BJ/def2-svp level of theory for the molecule CH\ :sub:`3`\ CH\ :sub:`2`\ SO\ :sub:`3`\ H. The input file is as follows: .. code-block:: :caption: ch3ch2so3h.inp :linenos: basis def2-svp end scf charge 0 spin2p1 1 end grimmedisp type bj end mol S 0.00000 0.00000 0.00000 O 1.37599 0.00000 0.00000 O -0.71199 1.18363 0.00000 O -0.45996 -0.69162 -1.43483 H -0.82001 -0.11945 -2.10909 C -0.62121 -1.12008 1.07039 H 0.14029 -1.49066 1.80382 H -1.47518 -0.71991 1.67461 C -1.14768 -2.34547 0.38391 H -0.34399 -2.85090 -0.10958 H -1.58848 -2.99835 1.10803 H -1.88611 -2.06128 -0.33641 end task opt b3lyp end During optimization, Qbics will write the latest optimized structure to ``ch3ch2so3h-opt.xyz``, and the optimization process will be written to ``ch3ch2so3h-opt-traj.xyz``. You can visualize the optimization process using `Qbics-MolStar `_. After optimization, we can see the energies by this command: .. code-block:: bash $ grep "E = " ch3ch2so3h-opt-traj.xyz | cat -n 1 E = -703.24635624 Hartree; coordinates in Angstrom. 2 E = -703.25408101 Hartree; coordinates in Angstrom. 3 E = -703.27695102 Hartree; coordinates in Angstrom. 4 E = -703.28018008 Hartree; coordinates in Angstrom. 5 E = -703.28163730 Hartree; coordinates in Angstrom. 6 E = -703.28278993 Hartree; coordinates in Angstrom. 7 E = -703.28366187 Hartree; coordinates in Angstrom. 8 E = -703.28479281 Hartree; coordinates in Angstrom. 9 E = -703.28607847 Hartree; coordinates in Angstrom. 10 E = -703.28606122 Hartree; coordinates in Angstrom. 11 E = -703.28666766 Hartree; coordinates in Angstrom. 12 E = -703.28677457 Hartree; coordinates in Angstrom. 13 E = -703.28709275 Hartree; coordinates in Angstrom. 14 E = -703.28735283 Hartree; coordinates in Angstrom. 15 E = -703.28758537 Hartree; coordinates in Angstrom. 16 E = -703.28770103 Hartree; coordinates in Angstrom. 17 E = -703.28775039 Hartree; coordinates in Angstrom. 18 E = -703.28777717 Hartree; coordinates in Angstrom. 19 E = -703.28779907 Hartree; coordinates in Angstrom. 20 E = -703.28782446 Hartree; coordinates in Angstrom. 21 E = -703.28783816 Hartree; coordinates in Angstrom. 22 E = -703.28786067 Hartree; coordinates in Angstrom. 23 E = -703.28786999 Hartree; coordinates in Angstrom. 24 E = -703.28787815 Hartree; coordinates in Angstrom. 25 E = -703.28788727 Hartree; coordinates in Angstrom. 26 E = -703.28789946 Hartree; coordinates in Angstrom. 27 E = -703.28792996 Hartree; coordinates in Angstrom. 28 E = -703.28798086 Hartree; coordinates in Angstrom. 29 E = -703.28805214 Hartree; coordinates in Angstrom. 30 E = -703.28821936 Hartree; coordinates in Angstrom. 31 E = -703.28858252 Hartree; coordinates in Angstrom. 32 E = -703.28826689 Hartree; coordinates in Angstrom. 33 E = -703.28902797 Hartree; coordinates in Angstrom. 34 E = -703.28916009 Hartree; coordinates in Angstrom. 35 E = -703.28936878 Hartree; coordinates in Angstrom. 36 E = -703.28949234 Hartree; coordinates in Angstrom. 37 E = -703.28966812 Hartree; coordinates in Angstrom. 38 E = -703.28983363 Hartree; coordinates in Angstrom. 39 E = -703.28987269 Hartree; coordinates in Angstrom. 40 E = -703.29037808 Hartree; coordinates in Angstrom. 41 E = -703.29064316 Hartree; coordinates in Angstrom. 42 E = -703.29099480 Hartree; coordinates in Angstrom. 43 E = -703.29143705 Hartree; coordinates in Angstrom. 44 E = -703.29170664 Hartree; coordinates in Angstrom. 45 E = -703.29220829 Hartree; coordinates in Angstrom. 46 E = -703.29260499 Hartree; coordinates in Angstrom. 47 E = -703.29312866 Hartree; coordinates in Angstrom. 48 E = -703.29323150 Hartree; coordinates in Angstrom. 49 E = -703.29367916 Hartree; coordinates in Angstrom. 50 E = -703.29379655 Hartree; coordinates in Angstrom. 51 E = -703.29399614 Hartree; coordinates in Angstrom. 52 E = -703.29413105 Hartree; coordinates in Angstrom. 53 E = -703.29421491 Hartree; coordinates in Angstrom. 54 E = -703.29430457 Hartree; coordinates in Angstrom. 55 E = -703.29434373 Hartree; coordinates in Angstrom. 56 E = -703.29437098 Hartree; coordinates in Angstrom. 57 E = -703.29437738 Hartree; coordinates in Angstrom. 58 E = -703.29437298 Hartree; coordinates in Angstrom. 59 E = -703.29436155 Hartree; coordinates in Angstrom. 60 E = -703.29435907 Hartree; coordinates in Angstrom. 61 E = -703.29435586 Hartree; coordinates in Angstrom. 62 E = -703.29435027 Hartree; coordinates in Angstrom. 63 E = -703.29435261 Hartree; coordinates in Angstrom. The energy has converged to ``-703.29435261 Hartree``, which is the final energy of the optimized structure. It has taken 63 steps to converge. Example: xTB Pre-optimization befor DFT Geomtry Optimization -------------------------------------------------------------- In the above example, the starting geometry is very bad, so it took 63 steps. Therefore, we can use :doc:`../keywords/xtb` to do a pre-optimization before DFT optimization, which is very fast: .. code-block:: :caption: ch3ch2so3h-xtb.inp :linenos: basis def2-svp end scf charge 0 spin2p1 1 end grimmedisp type bj end mol S 0.00000 0.00000 0.00000 O 1.37599 0.00000 0.00000 O -0.71199 1.18363 0.00000 O -0.45996 -0.69162 -1.43483 H -0.82001 -0.11945 -2.10909 C -0.62121 -1.12008 1.07039 H 0.14029 -1.49066 1.80382 H -1.47518 -0.71991 1.67461 C -1.14768 -2.34547 0.38391 H -0.34399 -2.85090 -0.10958 H -1.58848 -2.99835 1.10803 H -1.88611 -2.06128 -0.33641 end task opt xtb end This optimization only take several seconds. The final optimized structure is written to ``ch3ch2so3h-xtb-opt.xyz``. We can use it as the starting geometry for DFT optimization: .. code-block:: :caption: ch3ch2so3h-dft.inp :linenos: basis def2-svp end scf charge 0 spin2p1 1 end grimmedisp type bj end mol ch3ch2so3h-xtb-opt.xyz end task opt xtb end This task takes only 17 steps to converge to ``-703.29434128 Hartree``, which is very close to the one from direct DFT optmization, but is much faster! Example: DFT Geomtry Optimization with Constraints -------------------------------------------------------------- We can apply some constraints for geometry optimization. First, let us visualize this molecule with `Qbics-MolStar `_. Just drag ``ch3ch2so3h.inp`` into `Qbics-MolStar `_, then right-click :guilabel:`All`, then click :guilabel:`Add Component`, then click :guilabel:`Label` in :guilabel:`Representation`, then click :guilabel:`+ Create Component`, you can see clearly the atomic indices: .. figure:: figs/a7.png Now, we want to fix the bond 4-5 and torsion 9-6-1-4, we can use the following input: .. code-block:: :caption: ch3ch2so3h-cons.inp :linenos: basis def2-svp end scf charge 0 spin2p1 1 end grimmedisp type bj end opt fix_bond 4 5 fix_torsion 9 6 1 4 end mol S 0.00000 0.00000 0.00000 O 1.37599 0.00000 0.00000 O -0.71199 1.18363 0.00000 O -0.45996 -0.69162 -1.43483 H -0.82001 -0.11945 -2.10909 C -0.62121 -1.12008 1.07039 H 0.14029 -1.49066 1.80382 H -1.47518 -0.71991 1.67461 C -1.14768 -2.34547 0.38391 H -0.34399 -2.85090 -0.10958 H -1.58848 -2.99835 1.10803 H -1.88611 -2.06128 -0.33641 end task opt b3lyp end This is very easy: just add constraints to ``opt...end``. Some exmaples of constraints are: - ``fix_atoms 1-4 8 10-11``: Fix atoms: ``1 2 3 4 8 10 11``; - ``fix_bond 4 5``: Fix the bond 4-5; - ``fix_angle 5 7 9``: Fix the angle 5-7-9; - ``fix_torsion 1 3 4 5``: Fix the torsion 1-3-4-5. There can be an arbitrary number of constraints. After calculation, it converge to ``-703.28739992 Hartree``, with fixed bond length and torsion. This can also be seen with `Qbics-MolStar `_. Drag ``ch3ch2so3h-cons-opt-traj.xyz`` into `Qbics-MolStar `_. Then, click atom 4 and 5, and then right-click one atom, click :guilabel:`Measure Distance` to shown the bond length; for the torsion, do the similar. Then, click the icon at the left-up corner, and click :guilabel:`Start` to play the trajectory. You can see that the bond length and torsion are successfully constrained. For comparison, we also shoa a bond length that is not constrained. .. figure:: figs/a8.gif Example: MM for Solvated Protein 5VBM ----------------------------------------- For a MM system, a geometry optimization is also necessary. For the solvated protein 5VBM system given in :doc:`charmm`, we can optimize the geometry by simply changing ``energy`` to ``opt``: .. code-block:: :caption: 5vbm-sol.inp :linenos: charmm # Here: # toppar_water_ions.str: the parameter file for water and ions. # par_all36m_prot.prm: the parameter file for proteins. # Both are obtained from website: https://mackerell.umaryland.edu/charmm_ff.shtml. parameters toppar_water_ions.str par_all36m_prot.prm topology 5vbm-sol.psf scaling14 1.0 rcutoff 12.0 rswitch 10.0 use_PBC # Add this line to turn on PBC. Lbox 70 70 70 # Define the box size. electrostatic pme # Use PME for electrostatics. PMEk 72 72 72 end opt max_step 1000 end mol 5vbm-sol.pdb end task opt charmm end However, since this system contains 40353 atoms, it will take very long to reach the strict stationary point, which is often unnecessary if it is used as an initial structure for molecular dynamic simulations. Thus, we use ``max_step 1000`` in ``opt...end`` to limit the maximum step of optimization to 1000. In its output ``5vbm-sol.out``, we can see that just in 100 steps, the enregy is optimized from the ridiculous value ``11507712447667.06 kcal/mol`` to ``-146029.06 kcal/mol``, which is much more reasonable. Finally, the optimized energy is ``-170315.90 kcal/mol``.