The selection of transition state optimization method?

8 posts / 0 new
Last post
tyrion
The selection of transition state optimization method?

I use DL-Find optimiser to find the transition state. When the "initial_hessian" and "update_method" are set to "onepoint" and "bofill" respectively, the informations are as below:

MM Energies (kcal/mol) total : -1.00089967E+05
short range: 1.28339379E+04 coulombic : -1.18583958E+05 3body : 0.00000000E+00
bond : 7.65787712E+02 angle : 1.79850102E+03 str-bend : 0.00000000E+00
per dihed : 3.09576402E+03 harm dihed : 0.00000000E+00 angle-angle: 0.00000000E+00
ang-ang-tor: 0.00000000E+00 urey-brdly : 0.00000000E+00 inversion : 0.00000000E+00
Contribution to energy from turbomole: -3774.171508 (a.u.)
Contribution to energy from dl_poly: -159.499669 (a.u.)
Contribution to energy from additional MM energy terms: 0.000000 (a.u.)
-----------------------------------------------------------------------------------
QM/MM Energy: -3933.671177 (a.u.)
-----------------------------------------------------------------------------------
Energy calculation finished, energy: -3.933671177E+03
Finite-difference Hessian: variable 2370/2943 direction= 1
Energy difference to midpoint: 1.73E-05 H
Abs. Gradient difference to midpoint: 7.01E-03

The process of "Finite-difference Hessian: variable 2370/2943 direction= 1" consumed 8 days, which was been too long...
What the selection of "initial_hessian" and "update_method" are the best choice for the transition state optimization?
Can anyone help me ? Thanks for replying.

tyrion
Sorry to forget to say that

Sorry to forget to say that my server is 24 cores with 48 threads and the version of chemshell is 3.6.0 on linux.

tomkeal
tomkeal's picture
Hi,

Hi,

For large active regions where a Hessian calculation is prohibitively expensive we recommend using the Dimer method in DL-FIND instead, which avoids the calculation of a Hessian entirely. For an example see the tutorial page at: http://www.cse.scitech.ac.uk/ccg/software/chemshell/tutorial/dl-find_ts....

Tom

tyrion
Dear Tom

Dear Tom

Thanks for your suggestion. I am sorry to express unclearly. My systems is performed with the QM/MM calculation, including 87 atoms in the QM region and 987 atoms in the MM region. should I use the Dimer method instead P-RFO for my current system? In my script, the "active_atoms" argument refers to MM region (987 atoms). Is this setting correct? Thank you again!

tomkeal
tomkeal's picture
Yes, for a QM/MM system you

Yes, for a QM/MM system you should as a general rule use the Dimer method, in which case it is OK to have a large active region. If you did use P-RFO you would have to select a much smaller active region as suggested below.

Tom

tyrion
Dear Tom

Dear Tom

OK! I will try the Dimer method for searching the transition state, Thank you very much!!

uccaagh
Dear Tyrion,

Dear Tyrion,

You can greatly reduce the time of Hessian calculations by limiting the active size of your calculation with the keyword 'active_atoms= { }' with a list of atoms you wish to calculate the hessian for (perhaps try just the QM atoms first?) e.g. active_atoms= {15-25 40-45}

If you really want to calculate the Hessian for your whole system perhaps try benchmarking a few steps with initial_hessian=twopoint

tyrion
Dear uccaagh,

Dear uccaagh,

Thank you so much for your reply, I did not express it clearly. I major in QM/MM calculation. I use Chemshll to combine DL-poly for MM region and Turbomole for QM region, and there are 87 atoms in my QM region. In other words, I set the active_atoms=qmatoms, and the atoms among the qmatoms are as below:

---------------------------------------------------------------------------------------------------------------------------------------------------------------------
set qm_atoms { 3905 3906 3907 3908 3909 3910 3911 3912 3913 3914 3915 3916 3917 3918 3919 3920 3921 3922 3923 3924 3925 3926 3927 3928 3929 3930 3931 3932 3933 3934 3935 3936 3937 3938 3939 3940 3941 3942 3943 3944 946 947 948 949 950 951 952 953 954 955 956 2050 2051 2052 2053 2054 2055 2056 2057 2058 2059 2060 2693 2694 2695 2696 2697 2698 2699 2700 1814 1815 1816 1817 1818 1819 1820 1821 1822 31249 31250 31251 }
set active_atom_numbers { 73 77 78 116 117 118 119 120 121 122 123 153 154 155 157 158 340 341 343 679 734 736 737 925 926 928 929 930 931 932 933 934 935 936 937 938 939 940 941 942 943 944 945 946 947 948 949 950 951 952 953 954 955 956 957 958 959 960 961 962 963 964 965 966 967 968 969 970 971 972 973 974 975 976 977 978 979 980 981 982 983 984 985 986 987 988 989 990 991 992 993 994 995 996 997 998 999 1000 1001 1002 1003 1004 1005 1007 1008 1013 1014 1015 1016 1368 1369 1370 1371 1375 1377 1378 1379 1380 1381 1382 1383 1384 1385 1386 1387 1388 1389 1390 1391 1392 1393 1394 1395 1396 1397 1398 1399 1400 1401 1402 1404 1405 1406 1407 1408 1409 1410 1411 1412 1413 1414 1419 1430 1431 1432 1433 1455 1456 1457 1546 1547 1548 1549 1550 1570 1572 1573 1574 1575 1576 1577 1578 1579 1580 1581 1582 1583 1584 1585 1586 1641 1642 1643 1644 1646 1647 1648 1649 1653 1660 1661 1662 1663 1664 1665 1666 1667 1668 1669 1670 1671 1672 1673 1674 1675 1676 1677 1678 1679 1680 1681 1682 1683 1684 1685 1688 1700 1701 1702 1706 1709 1762 1764 1765 1766 1767 1768 1769 1770 1777 1778 1779 1780 1781 1782 1783 1784 1785 1786 1787 1788 1789 1790 1791 1792 1793 1794 1795 1796 1797 1798 1799 1800 1801 1802 1803 1804 1805 1806 1807 1808 1809 1810 1811 1812 1813 1814 1815 1816 1817 1818 1819 1820 1821 1822 1823 1824 1825 1826 1827 1828 1829 1830 1831 1833 1834 1835 1836 1837 1838 1839 1841 1842 1843 1844 1845 1846 1847 1848 1849 1850 1851 1852 1853 1854 1855 1856 1857 1858 1859 1860 1897 1898 1899 1906 1907 1908 1909 1910 1911 1912 1913 1914 1915 1916 1917 1918 1919 1920 1921 1922 1923 1924 1925 1926 1927 1928 1929 1930 1931 1932 1933 1934 1935 1936 1937 1938 1939 1940 1941 1942 1944 1948 1949 2011 2012 2013 2015 2016 2017 2018 2019 2020 2021 2022 2023 2024 2025 2026 2027 2028 2029 2030 2031 2032 2033 2034 2035 2036 2037 2038 2039 2040 2041 2042 2043 2044 2045 2046 2047 2048 2049 2050 2051 2052 2053 2054 2055 2056 2057 2058 2059 2060 2061 2062 2063 2064 2065 2066 2067 2068 2069 2070 2071 2072 2073 2074 2075 2076 2077 2078 2079 2080 2081 2082 2083 2084 2085 2086 2087 2088 2089 2090 2091 2092 2093 2094 2095 2096 2097 2098 2099 2100 2101 2102 2106 2107 2108 2109 2110 2111 2112 2113 2114 2115 2119 2120 2122 2125 2126 2127 2128 2144 2186 2187 2191 2192 2193 2194 2195 2198 2428 2434 2438 2453 2454 2457 2458 2468 2469 2470 2471 2472 2473 2474 2475 2476 2477 2478 2479 2480 2481 2482 2483 2484 2485 2486 2487 2488 2489 2490 2491 2492 2493 2494 2495 2496 2497 2498 2499 2500 2501 2502 2503 2504 2505 2506 2507 2508 2509 2510 2511 2512 2513 2514 2515 2516 2517 2518 2519 2520 2521 2522 2523 2524 2525 2526 2527 2528 2529 2530 2531 2532 2533 2534 2535 2536 2537 2538 2539 2540 2541 2542 2543 2544 2545 2546 2547 2548 2549 2550 2551 2552 2553 2554 2555 2556 2557 2558 2559 2560 2571 2572 2573 2574 2575 2576 2577 2578 2587 2643 2644 2645 2647 2648 2649 2650 2651 2653 2654 2655 2656 2657 2658 2659 2660 2661 2662 2663 2664 2665 2666 2667 2668 2669 2670 2671 2672 2673 2674 2675 2676 2677 2678 2679 2680 2681 2682 2683 2684 2685 2686 2687 2688 2689 2690 2691 2692 2693 2694 2695 2696 2697 2698 2699 2700 2701 2702 2703 2704 2705 2706 2707 2708 2709 2710 2711 2712 2717 2718 2719 2720 2721 2722 2723 2724 2725 2726 2727 2728 2729 2730 2731 2732 2733 2734 2735 2736 2737 2738 2739 2740 2741 2753 2756 2757 2758 2767 2768 2769 3905 3906 3907 3908 3909 3910 3911 3912 3913 3914 3915 3916 3917 3918 3919 3920 3921 3922 3923 3924 3925 3926 3927 3928 3929 3930 3931 3932 3933 3934 3935 3936 3937 3938 3939 3940 3941 3942 3943 3944 4033 4034 4035 4342 4343 4344 5161 5162 5434 5435 5436 6379 6380 6381 6452 6793 6794 6795 6974 7105 7106 7107 8164 8166 8386 8387 8388 9433 9434 9435 9610 9611 9612 9694 9695 9696 9985 9986 9987 10312 10313 10314 10690 10691 10692 10891 10892 10893 10912 1091210914 11290 11398 11399 11400 11428 11429 11430 11506 11507 11508 11824 11825 11826 11845 11846 11847 12961 12962 12963 12994 12995 12996 13105 13106 13107 13171 13172 13173 14062 14063 14064 14320 14321 14322 14347 14348 14349 14389 1438914391 14626 14627 14628 14641 14642 14643 14836 14837 14838 15139 15140 15141 15241 15242 15243 15616 15618 16021 16022 16023 16537 16538 16539 16576 16577 16578 16666 16667 16668 17221 17222 17223 17434 17435 17436 17692 17693 17694 1786517867 18064 18065 18066 18175 18176 18322 18323 18324 19162 19163 19164 19444 19445 19540 19541 19542 19624 19625 19626 20194 20195 20196 20842 20844 20893 20894 20895 21163 21164 21165 21178 21179 21180 22099 22100 22101 22238 22507 2250722509 22702 22703 22704 23555 23556 23731 23732 23733 24243 25189 25190 25191 25339 25340 25341 25846 25848 25987 25988 25989 26572 26573 26704 26705 26706 27250 27251 27252 27391 27392 27393 28000 28001 28002 28015 28016 28017 28693 2869328695 28882 28883 28884 29236 29237 29238 29749 29751 30238 30239 30240 30529 30530 30531 30646 30647 30928 30929 30930 31249 31250 31251 31399 31400 31401 31594 31595 31596 31780 31781 31782 31924 31925 31926 32050 32051 32052 32729 3279632798 32799 33427 33428 33429 33763 33764 33765 33832 33833 33834 33936 33976 33977 33978 34273 34274 34275 34336 34338 34387 34388 34989 }
---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
Among the qmatoms, the "set qm_atoms" argument represents for my QM region atoms, which contains 87 atoms. And the "set active_atom_numbers" argument refers to the atoms in my MM region. And my setting script are as below:
----------------------------------------------------------------------------------------------------------------------------
set systemname 49.81
set ambertop 6A4Z_solv.prmtop
source ${systemname}.charges
source qmatoms
set qm_atoms [ expand_range $qm_atoms ]
dl-find \
coords=${systemname}.c \
result=${systemname}_result.c \
list_option= full \
tolerance=0.00135 delta=0.01 \
update_method=bofill maxupdate=10000 \
optimiser=prfo \
initial_hessian=onepoint \
maxcycle=8888 maxene=8000 maxstep = 0.1 \
active_atoms= $active_atom_numbers\
theory= hybrid : [ list \
coupling= shift \
qm_region= $qm_atoms \
cutoff= 0 \
list_option= medium \
atom_charges= $atom_charges \
qm_theory= turbomole : [ list \
hamiltonian= b3-lyp \
basis= 6-31g* \
charge= 0 \
accuracy= medium \
scftype=hf \
mult= 6 \
scratchdir=tmp \
read_control = yes \
harmonic= yes ] \
mm_theory= dl_poly : [ list \
amber_prmtop_file=${ambertop} \
scale14 = [ list [ expr 1 / 1.2 ] 0.5 ] \
conn= ${systemname}.c \
use_pairlist = no \
save_dl_poly_files = yes \
list_option=none ] ]
-------------------------------------------------------------------------------------------------------------------------
The suggestion you provided is that I should set "active_atoms= $qm_atoms"? Or is there anything else needing to revise in my script? Thank you!
Thanks for your helping!

Log in or register to post comments