.. include:: autodoc_abbr_options_c.rst .. index:: DFTD3 .. _`sec:dftd3`: Interface to DFTD3 by S. Grimme =============================== .. codeauthor:: Lori A. Burns .. sectionauthor:: Lori A. Burns *Module:* :ref:`Samples ` |PSIfour| contains code to interface to the DFTD3 program of S. Grimme, which is freely downloadable from `http://www.thch.uni-bonn.de/tc/index.php?section=downloads&subsection=getd3&lang=english `_). Installation ~~~~~~~~~~~~ DFTD3 is available as conda package. If using the |PSIfour| binary, the ``dftd3`` executable has already been installed alongside. If using |PSIfour| built from source, and anaconda or miniconda has already been installed (instructions at :ref:`sec:quickconda`), the ``dftd3`` executable can be obtained through ``conda install dftd3``. To be used by |PSIfour|, the program binary (``dftd3``) must be found in your :envvar:`PSIPATH` or :envvar:`PATH` (in that order). If |PSIfour| is unable to execute the binary, an error will be reported. To preferentially use a particular dftd3 compilation, simply adjust its position in the path environment variables. To remove the dftd3 that conda installs alongside |PSIfour|, ``conda remove dftd3``. Alternatively, to build DFTD3 yourself, follow the instructions provided with the source (essentially, unpack the source, edit the Makefile to select a Fortran compiler, and run make). From version 3.1.0 onwards, DFTD3 can be used as-is; for earlier versions, patches are available: :source:`share/scripts/patch_grimme_dftd3.3.0.2`. >>> cd dftd3 >>> ls dftd3.tar patch_grimme_dftd3.3.0.2 >>> tar -xvf dftd3.tar copyc6.f dftd3.f Makefile man.pdf pars.f param >>> patch < patch_grimme_dftd3.3.0.2 patching file dftd3.f >>> make making dftd3.o from dftd3.f ifort -O -c dftd3.f -o dftd3.o making copyc6.o from copyc6.f ifort -O -c copyc6.f -o copyc6.o ifort dftd3.o copyc6.o -o ./dftd3 >>> ls Makefile copyc6.o dftd3.f dftd3.tar param patch_grimme_dftd3.3.0.2 copyc6.f dftd3 dftd3.o man.pdf pars.f Theory ~~~~~~ The local or semilocal character of conventional density functionals necessarily leads to neglect of the long-range correlation interactions which capture attractive van der Waals forces. Initially proposed by Yang [Wu:2002:515]_ and assiduously developed by Grimme, [Grimme:2004:1463]_ [Grimme:2006:1787]_ [Grimme:2010:154104]_ the DFT+Dispersion method appends to the base functional a scaled, damped, and fitted leading term to the well-known dispersion energy series, :math:`E_{disp} = -C_6/R^6 -C_8/R^8 -C_{10}/R^{10}-\cdots`. The DFT-D2 [Grimme:2006:1787]_ variant takes the explicit form below. Here, dispersion coefficients, :math:`C_6^{ij}`, obtained from the geometric mean of tabulated elemental values, are summed over interatomic distances, :math:`R_{ij}`, modulated by a damping function, :math:`f_{damp}(R_{ij})`, that gradually activates the dispersion correction (at a rate characterized by :math:`\alpha_6`) over a distance characterized by the sum of the two atomic vdW radii, :math:`R_{vdW}`, while an overall scaling term, :math:`s_6`, is optimized to be unique to each :math:`E_{xc}` functional. (:math:`\alpha_6` is sometimes allowed to vary as well.) .. math:: E_{disp}^{\text{D2}}=-s_6 \sum_{i,j>i}^{N_{at}} \frac{C_6^{ij}}{(R_{ij})^6} f_{damp}(R_{ij}) :label: DFTD2 .. math:: f_{damp}(R_{ij}) = \frac{1}{1 + e^{- \alpha_6 (R_{ij}/R_{vdW} - 1)}} Grimme recently presented a refined method, DFT-D3, [Grimme:2010:154104]_ which incorporates an additional :math:`R^{-8}` term in the dispersion series and adjusts the :math:`C_{6}^{ij}` combination formula and damping function. The individual atomic :math:`C_6^i` are interpolated from several reference values based upon coordination numbers extracted from the molecular structure, rather than assigned solely by atomic identity as in DFT-D2, and thereby incorporate some awareness of the chemical environment into an otherwise largely heuristic correction. The -D3 dispersion has the following form, where :math:`s_{r,6}` and :math:`s_8` are the customary nonunity parameters fitted for individual functionals. .. math:: E_{disp}^{\text{D3ZERO}}=-\sum_{n=6,8} s_n \sum_{i,j>i}^{N_{at}} \frac{C_n^{ij}}{(R_{ij})^n} f_{damp}(R_{ij}) :label: DFTD3ZERO .. math:: f_{damp}(R_{ij}) = \frac{1}{1 + 6 (R_{ij}/(s_{r,n} R_0^{ij}))^{- \alpha_n}} A modified damping scheme for DFT-D3 using the rational damping form of Becke and Johnson was introduced in [Grimme:2011:1456]_. The parameters fit for individual functionals are now :math:`s_6`, :math:`s_8`, :math:`a_1`, and :math:`a_2`. .. math:: E_{disp}^{\text{D3BJ}}=-\sum_{n=6,8} s_n \sum_{i,j>i}^{N_{at}} \frac{C_n^{ij}}{(R_{ij})^n + (f_{damp})^n} .. math:: f_{damp} = a_1 \sqrt{\frac{C_8^{ij}}{C_6^{ij}}} + a_2 All parameters characterizing the dispersion correction are taken from `http://toc.uni-muenster.de/DFTD3/ `_ or else from the literature. Running DFTD3 ~~~~~~~~~~~~~ A number of *a posteriori* dispersion corrections are available in |PSIfour|. While most are computed within |PSIfours| codebase (-D1, -D2, -CHG, -DAS2009, -DAS2010), the -D3 correction and its variants are available only through the ``DFTD3`` program. Once installed, the ``dftd3``/|PSIfour| interface is transparent, and all corrections are interfaced exactly alike. Dispersion corrections are built into DFT functionals, so appending an *a posteriori* correction to a computation is as simple as ``energy('b2plyp-d')`` *vs.* ``energy('b2plyp')``. For example, the following input file computes (with much redundant work) for water a B3LYP, a B3LYP-D2, and a B3LYP-D3 (zero-damping) energy. :: molecule h2o { O H 1 1.0 H 1 1.0 2 104.5 } set { basis cc-pVDZ } energy('b3lyp') energy('b3lyp-d') energy('b3lyp-d3') Consult the table :ref:`-D Functionals ` to see for each functional what corrections are available and what default parameters define them. The dispersion correction is available after a calculation in the PSI variable :psivar:`DISPERSION CORRECTION ENERGY `. By default, the output from the ``dftd3`` program is suppressed; to see it in the output file, set print > 2. .. _`table:dashd`: .. table:: Variants of S. Grimme's -D correction +-------------------------------------+---------------------------------------------------------------------------------------------------------+---------------------------------------------------------------+ | Extension [#f0]_ | Variant and Computing Program | |scf__dft_dispersion_parameters| | +=====================================+=========================================================================================================+===============================================================+ | -D | alias to -D2P4 | | +-------------------------------------+---------------------------------------------------------------------------------------------------------+---------------------------------------------------------------+ | -D1 | -D1 [#f1]_ within |PSIfour| | | +-------------------------------------+---------------------------------------------------------------------------------------------------------+---------------------------------------------------------------+ | -D2 | alias to -D2P4 | | +-------------------------------------+---------------------------------------------------------------------------------------------------------+---------------------------------------------------------------+ | -D2P4 | -D2 [#f2]_ within |PSIfour| | [:math:`s_6`] | +-------------------------------------+---------------------------------------------------------------------------------------------------------+---------------------------------------------------------------+ | -D2GR | -D2 [#f2]_ through ``dftd3`` | [:math:`s_6`, :math:`\alpha_6`] | +-------------------------------------+---------------------------------------------------------------------------------------------------------+---------------------------------------------------------------+ | -D3 | alias to -D3ZERO | | +-------------------------------------+---------------------------------------------------------------------------------------------------------+---------------------------------------------------------------+ | -D3ZERO | -D3 [#f3]_ w/ original zero-damping through ``dftd3`` | [:math:`s_6`, :math:`s_8`, :math:`s_{r,6}`, :math:`\alpha_6`] | +-------------------------------------+---------------------------------------------------------------------------------------------------------+---------------------------------------------------------------+ | -D3BJ | -D3 [#f4]_ w/ newer Becke-Johnson rational damping through ``dftd3`` | [:math:`s_6`, :math:`s_8`, :math:`a_1`, :math:`a_2`] | +-------------------------------------+---------------------------------------------------------------------------------------------------------+---------------------------------------------------------------+ | -D3M | alias to -D3MZERO | | +-------------------------------------+---------------------------------------------------------------------------------------------------------+---------------------------------------------------------------+ | -D3MZERO | -D3 [#f5]_ w/ reparameterized and more flexible original zero-damping through ``dftd3`` | [:math:`s_6`, :math:`s_8`, :math:`s_{r,6}`, :math:`\beta`] | +-------------------------------------+---------------------------------------------------------------------------------------------------------+---------------------------------------------------------------+ | -D3MBJ | -D3 [#f5]_ w/ reparameterized newer Becke-Johnson rational damping through ``dftd3`` | [:math:`s_6`, :math:`s_8`, :math:`a_1`, :math:`a_2`] | +-------------------------------------+---------------------------------------------------------------------------------------------------------+---------------------------------------------------------------+ .. rubric:: Footnotes .. [#f0] Note that there are functionals with these extensions (*e.g.*, wB97X-D) that, not being Grimme corrections, have nothing to do with this table. .. [#f1] [Grimme:2004:1463]_ .. [#f2] [Grimme:2006:1787]_ .. [#f3] [Grimme:2010:154104]_ .. [#f4] [Grimme:2011:1456]_ .. [#f5] [Smith:2016:2197]_ A few practical examples: * DFT-D2 single point with default parameters (``dftd3`` not called) :: energy('bp86-d') * DFT-D3BJ optimization with default parameters :: optimize('pbe-d3bj') * DFT-D2 optimization with custom s6 parameter :: set dft_dispersion_parameters [1.20] optimize('b3lyp-d2') * DFT-D3ZERO single point (b3lyp) with custom s8 parameter (reset all four values) :: set dft_dispersion_parameters [1.0, 2.0, 1.261, 14.0] energy('b3lyp-d3') If only dispersion corrections (rather than total energies) are of interest, the ``dftd3`` program can be run independently of the scf through the python function :py:func:`~qcdb.interface_dftd3.run_dftd3`. (This function is the same |PSIfour|/``dftd3`` interface that is called during an scf job.) This route is much faster than running a DFT-D energy. * Some set-up:: molecule nene { Ne Ne 1 2.0 } nene.update_geometry() * The same four dispersion corrections/gradients as the section above:: >>> print nene.run_dftd3('bp86', 'd', dertype=0) -7.735e-05 >>> E, G = nene.run_dftd3('pbe', 'd3bj') >>> print G [[0.0, 0.0, -1.1809087569358e-05], [0.0, 0.0, 1.1809087569358e-05]] >>> E, G = nene.run_dftd3('b3lyp', 'd2', {'s6': 1.20}) >>> print E -8.84e-05 >>> E, G = nene.run_dftd3(dashlvl='d3', dashparam={'s8': 2.0, 'alpha6': 14.0, 'sr6': 1.261, 's6': 1.0}) >>> print E -0.00024762 .. autofunction:: qcdb.interface_dftd3.run_dftd3 .. comment print_stdout(' -D correction from Py-side') .. comment eneyne.update_geometry() .. comment E, G = eneyne.run_dftd3('b3lyp', 'd2gr') .. comment compare_values(ref_d2[0], E, 7, 'Ethene-Ethyne -D2') .. comment mA = eneyne.extract_subsets(1) .. comment E, G = mA.run_dftd3('b3lyp', 'd2gr') .. comment compare_values(ref_d2[1], E, 7, 'Ethene -D2') .. comment mB = eneyne.extract_subsets(2) .. comment E, G = mB.run_dftd3('b3lyp', 'd2gr') .. comment compare_values(ref_d2[2], E, 7, 'Ethyne -D2') .. comment #mBcp = eneyne.extract_subsets(2,1) .. comment #E, G = mBcp.run_dftd3('b3lyp', 'd2gr') .. comment #compare_values(ref_d2[2], E, 7, 'Ethyne(CP) -D2') .. comment .. comment E, G = eneyne.run_dftd3('b3lyp', 'd3zero') .. comment compare_values(ref_d3zero[0], E, 7, 'Ethene-Ethyne -D3 (zero)') .. comment mA = eneyne.extract_subsets(1) .. comment E, G = mA.run_dftd3('b3lyp', 'd3zero') .. comment compare_values(ref_d3zero[1], E, 7, 'Ethene -D3 (zero)') .. comment mB = eneyne.extract_subsets(2) .. comment E, G = mB.run_dftd3('b3lyp', 'd3zero') .. comment compare_values(ref_d3zero[2], E, 7, 'Ethyne -D3 (zero)') .. comment .. comment E, G = eneyne.run_dftd3('b3lyp', 'd3bj') .. comment compare_values(ref_d3bj[0], E, 7, 'Ethene-Ethyne -D3 (bj)') .. comment mA = eneyne.extract_subsets(1) .. comment E, G = mA.run_dftd3('b3lyp', 'd3bj') .. comment compare_values(ref_d3bj[1], E, 7, 'Ethene -D3 (bj)') .. comment mB = eneyne.extract_subsets(2) .. comment E, G = mB.run_dftd3('b3lyp', 'd3bj') .. comment compare_values(ref_d3bj[2], E, 7, 'Ethyne -D3 (bj)') .. comment .. comment E, G = eneyne.run_dftd3('b3lyp', 'd3') .. comment compare_values(ref_d3zero[0], E, 7, 'Ethene-Ethyne -D3 (alias)') .. comment E, G = eneyne.run_dftd3('b3lyp', 'd') .. comment compare_values(ref_d2[0], E, 7, 'Ethene-Ethyne -D (alias)') .. comment E, G = eneyne.run_dftd3('b3lyp', 'd2') .. comment compare_values(ref_d2[0], E, 7, 'Ethene-Ethyne -D2 (alias)') .. comment .. comment set basis sto-3g .. comment set scf_type df .. comment set dft_radial_points 50 # use really bad grid for speed since all we want is the -D value .. comment set dft_spherical_points 110 .. comment #set scf print 3 # will print dftd3 program output to psi4 output file .. comment .. comment .. comment print_stdout(' -D correction from C-side') .. comment activate(mA) .. comment energy('b3lyp-d2p4') .. comment compare_values(ref_d2[1], get_variable('DISPERSION CORRECTION ENERGY'), 7, 'Ethene -D2 (calling psi4 Disp class)') .. comment energy('b3lyp-d2gr') .. comment compare_values(ref_d2[1], get_variable('DISPERSION CORRECTION ENERGY'), 7, 'Ethene -D2 (calling dftd3 -old)') .. comment energy('b3lyp-d3zero') .. comment compare_values(ref_d3zero[1], get_variable('DISPERSION CORRECTION ENERGY'), 7, 'Ethene -D3 (calling dftd3 -zero)') .. comment energy('b3lyp-d3bj') .. comment compare_values(ref_d3bj[1], get_variable('DISPERSION CORRECTION ENERGY'), 7, 'Ethene -D3 (calling dftd3 -bj)') .. comment .. comment energy('b3lyp-d2') .. comment compare_values(ref_d2[1], get_variable('DISPERSION CORRECTION ENERGY'), 7, 'Ethene -D2 (alias)') .. comment energy('b3lyp-d3') .. comment compare_values(ref_d3zero[1], get_variable('DISPERSION CORRECTION ENERGY'), 7, 'Ethene -D3 (alias)') .. comment energy('b3lyp-d') .. comment compare_values(ref_d2[1], get_variable('DISPERSION CORRECTION ENERGY'), 7, 'Ethene -D (alias)') .. comment energy('wb97x-d') .. comment compare_values(-0.000834247063, get_variable('DISPERSION CORRECTION ENERGY'), 7, 'Ethene wb97x-d (chg)') .. comment .. comment print_stdout(' non-default -D correction from C-side') .. comment activate(mB) .. comment set dft_dispersion_parameters [0.75] .. comment energy('b3lyp-d2p4') .. comment compare_values(ref_pbe_d2[2], get_variable('DISPERSION CORRECTION ENERGY'), 7, 'Ethene -D2 (calling psi4 Disp class)') .. comment set dft_dispersion_parameters [0.75, 20.0] .. comment energy('b3lyp-d2gr') .. comment compare_values(ref_pbe_d2[2], get_variable('DISPERSION CORRECTION ENERGY'), 7, 'Ethene -D2 (calling dftd3 -old)') .. comment set dft_dispersion_parameters [1.0, 0.722, 1.217, 14.0] .. comment energy('b3lyp-d3zero') .. comment compare_values(ref_pbe_d3zero[2], get_variable('DISPERSION CORRECTION ENERGY'), 7, 'Ethene -D3 (calling dftd3 -zero)') .. comment set dft_dispersion_parameters [1.000, 0.7875, 0.4289, 4.4407] .. comment energy('b3lyp-d3bj') .. comment compare_values(ref_pbe_d3bj[2], get_variable('DISPERSION CORRECTION ENERGY'), 7, 'Ethene -D3 (calling dftd3 -bj)') .. comment .. comment set dft_dispersion_parameters [0.75] .. comment energy('b3lyp-d2') .. comment compare_values(ref_pbe_d2[2], get_variable('DISPERSION CORRECTION ENERGY'), 7, 'Ethene -D2 (alias)') .. comment set dft_dispersion_parameters [1.0, 0.722, 1.217, 14.0] .. comment energy('b3lyp-d3') .. comment compare_values(ref_pbe_d3zero[2], get_variable('DISPERSION CORRECTION ENERGY'), 7, 'Ethene -D3 (alias)') .. comment set dft_dispersion_parameters [0.75] .. comment energy('b3lyp-d') .. comment compare_values(ref_pbe_d2[2], get_variable('DISPERSION CORRECTION ENERGY'), 7, 'Ethene -D (alias)') .. comment activate(mA) .. comment set dft_dispersion_parameters [1.0] .. comment energy('wb97x-d') .. comment compare_values(-0.000834247063, get_variable('DISPERSION CORRECTION ENERGY'), 7, 'Ethene wb97x-d (chg)') .. comment .. comment print_stdout(' non-default -D correction from Py-side') .. comment eneyne.update_geometry() .. comment eneyne.run_dftd3('b3lyp', 'd2gr', {'s6': 0.75}) .. comment compare_values(ref_pbe_d2[0], get_variable('DISPERSION CORRECTION ENERGY'), 7, 'Ethene-Ethyne -D2') .. comment mA = eneyne.extract_subsets(1) .. comment mA.run_dftd3('b3lyp', 'd2gr', {'s6': 0.75}) .. comment compare_values(ref_pbe_d2[1], get_variable('DISPERSION CORRECTION ENERGY'), 7, 'Ethene -D2') .. comment mB = eneyne.extract_subsets(2) .. comment mB.run_dftd3('b3lyp', 'd2gr', {'s6': 0.75}) .. comment compare_values(ref_pbe_d2[2], get_variable('DISPERSION CORRECTION ENERGY'), 7, 'Ethyne -D2') .. comment .. comment eneyne.run_dftd3('b3lyp', 'd3zero', {'s6': 1.0, 's8': 0.722, 'sr6': 1.217, 'alpha6': 14.0}) .. comment compare_values(ref_pbe_d3zero[0], get_variable('DISPERSION CORRECTION ENERGY'), 7, 'Ethene-Ethyne -D3 (zero)') .. comment mA = eneyne.extract_subsets(1) .. comment mA.run_dftd3('b3lyp', 'd3zero', {'s6': 1.0, 's8': 0.722, 'sr6': 1.217, 'alpha6': 14.0}) .. comment compare_values(ref_pbe_d3zero[1], get_variable('DISPERSION CORRECTION ENERGY'), 7, 'Ethene -D3 (zero)') .. comment mB = eneyne.extract_subsets(2) .. comment mB.run_dftd3('b3lyp', 'd3zero', {'s6': 1.0, 's8': 0.722, 'sr6': 1.217, 'alpha6': 14.0}) .. comment compare_values(ref_pbe_d3zero[2], get_variable('DISPERSION CORRECTION ENERGY'), 7, 'Ethyne -D3 (zero)') .. comment .. comment eneyne.run_dftd3('b3lyp', 'd3bj', {'s6': 1.000, 's8': 0.7875, 'a1': 0.4289, 'a2': 4.4407}) .. comment compare_values(ref_pbe_d3bj[0], get_variable('DISPERSION CORRECTION ENERGY'), 7, 'Ethene-Ethyne -D3 (bj)') .. comment mA = eneyne.extract_subsets(1) .. comment mA.run_dftd3('b3lyp', 'd3bj', {'s6': 1.000, 's8': 0.7875, 'a1': 0.4289, 'a2': 4.4407}) .. comment compare_values(ref_pbe_d3bj[1], get_variable('DISPERSION CORRECTION ENERGY'), 7, 'Ethene -D3 (bj)') .. comment mB = eneyne.extract_subsets(2) .. comment mB.run_dftd3('b3lyp', 'd3bj', {'s6': 1.000, 's8': 0.7875, 'a1': 0.4289, 'a2': 4.4407}) .. comment compare_values(ref_pbe_d3bj[2], get_variable('DISPERSION CORRECTION ENERGY'), 7, 'Ethyne -D3 (bj)') .. comment .. comment eneyne.run_dftd3('b3lyp', 'd3', {'s6': 1.0, 's8': 0.722, 'sr6': 1.217, 'alpha6': 14.0}) .. comment compare_values(ref_pbe_d3zero[0], get_variable('DISPERSION CORRECTION ENERGY'), 7, 'Ethene-Ethyne -D3 (alias)') .. comment eneyne.run_dftd3('b3lyp', 'd', {'s6': 0.75}) .. comment compare_values(ref_pbe_d2[0], get_variable('DISPERSION CORRECTION ENERGY'), 7, 'Ethene-Ethyne -D (alias)') .. comment eneyne.run_dftd3('b3lyp', 'd2', {'s6': 0.75}) .. comment compare_values(ref_pbe_d2[0], get_variable('DISPERSION CORRECTION ENERGY'), 7, 'Ethene-Ethyne -D2 (alias)')