"""
Contains various constants and parameters for use in Grimme-D3 calculations. The values are extracted from
the original Fortran implementation and tables which can be found at
https://www.chemie.uni-bonn.de/grimme/de/software/dft-d3.
Some of the values have been rescaled to be in units of eV and Å.
References:
- https://doi.org/10.1063/1.3382344
- https://doi.org/10.1002/jcc.21759
"""
from pathlib import Path
import numpy as np
from ..elements import ELEMENT_DICT
# Taken directly from the original Fortran implementation at
# https://www.chemie.uni-bonn.de/grimme/de/software/dft-d3
# Some of the values don't match with the Pyykkö and Atsumi paper, but let's use
# these values to be consistent with the original implementation...
# fmt: off
R_COV = np.array([
0.32, 0.46, 1.20, 0.94, 0.77, 0.75, 0.71, 0.63, 0.64, 0.67,
1.40, 1.25, 1.13, 1.04, 1.10, 1.02, 0.99, 0.96, 1.76, 1.54,
1.33, 1.22, 1.21, 1.10, 1.07, 1.04, 1.00, 0.99, 1.01, 1.09,
1.12, 1.09, 1.15, 1.10, 1.14, 1.17, 1.89, 1.67, 1.47, 1.39,
1.32, 1.24, 1.15, 1.13, 1.13, 1.08, 1.15, 1.23, 1.28, 1.26,
1.26, 1.23, 1.32, 1.31, 2.09, 1.76, 1.62, 1.47, 1.58, 1.57,
1.56, 1.55, 1.51, 1.52, 1.51, 1.50, 1.49, 1.49, 1.48, 1.53,
1.46, 1.37, 1.31, 1.23, 1.18, 1.16, 1.11, 1.12, 1.13, 1.32,
1.30, 1.30, 1.36, 1.31, 1.38, 1.42, 2.01, 1.81, 1.67, 1.58,
1.52, 1.53, 1.54, 1.55
]) # fmt: on
"""
Covalent radii in Ångströms according to `Pyykkö and Atsumi
<https://doi.org/10.1002/chem.200800987>`_.
Values for metals reduced by 10% for use in Grimme-D3. Units are Ångströms.
"""
# Also directly copied from https://www.chemiebn.uni-bonn.de/pctc/mulliken-center/software/dft-d3/dft-d3
# fmt: off
R4R2 = np.array([
1.06224333, 0.82888801, 2.65640046, 2.03933801, 1.92856832,
1.64305726, 1.43499761, 1.3724829 , 1.2638088 , 1.17224667,
3.48508457, 2.89087376, 2.9909978 , 2.58389243, 2.2740204 ,
2.13845222, 1.97347304, 1.82395359, 4.22157869, 3.74458467,
3.49703613, 3.3274206 , 3.2159618 , 2.93504487, 3.07182697,
2.95500811, 2.86483063, 2.7966869 , 2.76544211, 2.69783613,
3.23935665, 2.93208485, 2.68132441, 2.57712005, 2.42939779,
2.28168674, 5.05608556, 4.59006237, 4.21865469, 3.93411462,
3.48575306, 3.27844446, 3.18309242, 3.07781854, 2.99361064,
2.92444844, 2.88011738, 2.95431897, 3.71525741, 3.42280034,
3.16495132, 3.07815323, 2.92805295, 2.78070457, 5.83261529,
5.3747446 , 4.94869507, 4.79924613, 4.74799572, 4.71016829,
4.68842987, 4.66595035, 4.65314893, 4.18033923, 4.65987342,
4.45799636, 4.52070409, 4.48521804, 4.47202839, 4.48392668,
4.1462391 , 4.34297241, 4.07762296, 3.87757775, 3.72481161,
3.54008293, 3.2039037 , 3.11554153, 3.0198106 , 3.06102925,
4.12642191, 3.84417539, 3.58862575, 3.53428553, 3.38157106,
3.22548305, 6.23982524, 5.87914635, 5.03447448, 4.58900926,
4.64162882, 4.57951396, 4.51876857, 4.68334909
]) # fmt: on
"""
sqrt(0.5*sqrt(Z)<r4>/<r2>(Z)) values for Grimme-D3. Units are Ångströms.
"""
# fmt: off
REF_CN = np.array(
[
[+0.9118, +0.0000, -1.0000, -1.0000, -1.0000], # H
[+0.0000, -1.0000, -1.0000, -1.0000, -1.0000], # He
[+0.0000, +0.9865, -1.0000, -1.0000, -1.0000], # Li
[+0.0000, +0.9808, +1.9697, -1.0000, -1.0000], # Be
[+0.0000, +0.9706, +1.9441, +2.9128, +4.5856], # B
[+0.0000, +0.9868, +1.9985, +2.9987, +3.9844], # C
[+0.0000, +0.9944, +2.0143, +2.9903, -1.0000], # N
[+0.0000, +0.9925, +1.9887, -1.0000, -1.0000], # O
[+0.0000, +0.9982, -1.0000, -1.0000, -1.0000], # F
[+0.0000, -1.0000, -1.0000, -1.0000, -1.0000], # Ne
[+0.0000, +0.9684, -1.0000, -1.0000, -1.0000], # Na
[+0.0000, +0.9628, +1.9496, -1.0000, -1.0000], # Mg
[+0.0000, +0.9648, +1.9311, +2.9146, -1.0000], # Al
[+0.0000, +0.9507, +1.9435, +2.9407, +3.8677], # Si
[+0.0000, +0.9947, +2.0102, +2.9859, -1.0000], # P
[+0.0000, +0.9948, +1.9903, -1.0000, -1.0000], # S
[+0.0000, +0.9972, -1.0000, -1.0000, -1.0000], # Cl
[+0.0000, -1.0000, -1.0000, -1.0000, -1.0000], # Ar
[+0.0000, +0.9767, -1.0000, -1.0000, -1.0000], # K
[+0.0000, +0.9831, +1.9349, -1.0000, -1.0000], # Ca
[+0.0000, +1.8627, +2.8999, -1.0000, -1.0000], # Sc
[+0.0000, +1.8299, +3.8675, -1.0000, -1.0000], # Ti
[+0.0000, +1.9138, +2.9110, -1.0000, -1.0000], # V
[+0.0000, +1.8269, 10.6191, -1.0000, -1.0000], # Cr
[+0.0000, +1.6406, +9.8849, -1.0000, -1.0000], # Mn
[+0.0000, +1.6483, +9.1376, -1.0000, -1.0000], # Fe
[+0.0000, +1.7149, +2.9263, +7.7785, -1.0000], # Co
[+0.0000, +1.7937, +6.5458, +6.2918, -1.0000], # Ni
[+0.0000, +0.9576, -1.0000, -1.0000, -1.0000], # Cu
[+0.0000, +1.9419, -1.0000, -1.0000, -1.0000], # Zn
[+0.0000, +0.9601, +1.9315, +2.9233, -1.0000], # Ga
[+0.0000, +0.9434, +1.9447, +2.9186, +3.8972], # Ge
[+0.0000, +0.9889, +1.9793, +2.9709, -1.0000], # As
[+0.0000, +0.9901, +1.9812, -1.0000, -1.0000], # Se
[+0.0000, +0.9974, -1.0000, -1.0000, -1.0000], # Br
[+0.0000, -1.0000, -1.0000, -1.0000, -1.0000], # Kr
[+0.0000, +0.9738, -1.0000, -1.0000, -1.0000], # Rb
[+0.0000, +0.9801, +1.9143, -1.0000, -1.0000], # Sr
[+0.0000, +1.9153, +2.8903, -1.0000, -1.0000], # Y
[+0.0000, +1.9355, +3.9106, -1.0000, -1.0000], # Zr
[+0.0000, +1.9545, +2.9225, -1.0000, -1.0000], # Nb
[+0.0000, +1.9420, 11.0556, -1.0000, -1.0000], # Mo
[+0.0000, +1.6682, +9.5402, -1.0000, -1.0000], # Tc
[+0.0000, +1.8584, +8.8895, -1.0000, -1.0000], # Ru
[+0.0000, +1.9003, +2.9696, -1.0000, -1.0000], # Rh
[+0.0000, +1.8630, +5.7095, -1.0000, -1.0000], # Pd
[+0.0000, +0.9679, -1.0000, -1.0000, -1.0000], # Ag
[+0.0000, +1.9539, -1.0000, -1.0000, -1.0000], # Cd
[+0.0000, +0.9633, +1.9378, +2.9353, -1.0000], # In
[+0.0000, +0.9514, +1.9505, +2.9259, +3.9123], # Sn
[+0.0000, +0.9749, +1.9523, +2.9315, -1.0000], # Sb
[+0.0000, +0.9811, +1.9639, -1.0000, -1.0000], # Te
[+0.0000, +0.9968, -1.0000, -1.0000, -1.0000], # I
[+0.0000, -1.0000, -1.0000, -1.0000, -1.0000], # Xe
[+0.0000, +0.9909, -1.0000, -1.0000, -1.0000], # Cs
[+0.0000, +0.9797, +1.8467, -1.0000, -1.0000], # Ba
[+0.0000, +1.9373, +2.9175, -1.0000, -1.0000], # La
[+2.7991, -1.0000, -1.0000, -1.0000, -1.0000], # Ce
[+0.0000, +2.9425, -1.0000, -1.0000, -1.0000], # Pr
[+0.0000, +2.9455, -1.0000, -1.0000, -1.0000], # Nd
[+0.0000, +2.9413, -1.0000, -1.0000, -1.0000], # Pm
[+0.0000, +2.9300, -1.0000, -1.0000, -1.0000], # Sm
[+0.0000, +1.8286, -1.0000, -1.0000, -1.0000], # Eu
[+0.0000, +2.8732, -1.0000, -1.0000, -1.0000], # Gd
[+0.0000, +2.9086, -1.0000, -1.0000, -1.0000], # Tb
[+0.0000, +2.8965, -1.0000, -1.0000, -1.0000], # Dy
[+0.0000, +2.9242, -1.0000, -1.0000, -1.0000], # Ho
[+0.0000, +2.9282, -1.0000, -1.0000, -1.0000], # Er
[+0.0000, +2.9246, -1.0000, -1.0000, -1.0000], # Tm
[+0.0000, +2.8482, -1.0000, -1.0000, -1.0000], # Yb
[+0.0000, +2.9219, -1.0000, -1.0000, -1.0000], # Lu
[+0.0000, +1.9254, +3.8840, -1.0000, -1.0000], # Hf
[+0.0000, +1.9459, +2.8988, -1.0000, -1.0000], # Ta
[+0.0000, +1.9292, 10.9153, -1.0000, -1.0000], # W
[+0.0000, +1.8104, +9.8054, -1.0000, -1.0000], # Re
[+0.0000, +1.8858, +9.1527, -1.0000, -1.0000], # Os
[+0.0000, +1.8648, +2.9424, -1.0000, -1.0000], # Ir
[+0.0000, +1.9188, +6.6669, -1.0000, -1.0000], # Pt
[+0.0000, +0.9846, -1.0000, -1.0000, -1.0000], # Au
[+0.0000, +1.9896, -1.0000, -1.0000, -1.0000], # Hg
[+0.0000, +0.9267, +1.9302, +2.9420, -1.0000], # Tl
[+0.0000, +0.9383, +1.9356, +2.9081, +3.9098], # Pb
[+0.0000, +0.9820, +1.9655, +2.9500, -1.0000], # Bi
[+0.0000, +0.9815, +1.9639, -1.0000, -1.0000], # Po
[+0.0000, +0.9954, -1.0000, -1.0000, -1.0000], # At
[+0.0000, -1.0000, -1.0000, -1.0000, -1.0000], # Rn
[+0.0000, +0.9705, -1.0000, -1.0000, -1.0000], # Fr
[+0.0000, +0.9662, +1.8075, -1.0000, -1.0000], # Ra
[+0.0000, +2.9070, -1.0000, -1.0000, -1.0000], # Ac
[+0.0000, +2.8844, -1.0000, -1.0000, -1.0000], # Th
[+0.0000, +2.8738, -1.0000, -1.0000, -1.0000], # Pa
[+0.0000, +2.8878, -1.0000, -1.0000, -1.0000], # U
[+0.0000, +2.9095, -1.0000, -1.0000, -1.0000], # Np
[+0.0000, +1.9209, -1.0000, -1.0000, -1.0000], # Pu
]
)
# fmt: on
"""
Reference coordination numbers for Grimme-D3. Values from:
`<https://doi.org/10.1063/1.3382344>`_
"""
K1 = 16
"""
Factor in exponent of coordination number calculation for Grimme-D3.
"""
K2 = 4 / 3
"""
Scale factor for covalent radii in coordination number calculation for Grimme-D3.
"""
K3 = 4
"""
Factor in exponent of Gaussian weight function in reference coordination number
calculation for Grimme-D3.
"""
_REF_C6 = None
_R0_AB = None
[docs]
def load_ref_c6():
"""
Load the reference C6 coefficient values for Grimme-D3. Values from:
`<https://doi.org/10.1063/1.3382344>`_
The returned array has shape (94, 94, 5, 5). The first two indices correspond
to pairs of chemical elements, and the (5, 5) array of values contains the C6
coeffiecients for the corresponding of reference coordination numbers stored in
:data:`REF_CN`. Units are eV*Å^6.
Returns:
ref_c6: numpy.ndarray of shape (94, 94, 5, 5). C6 coefficients.
"""
global _REF_C6
if _REF_C6 is None:
_REF_C6 = np.load(Path(__file__).parent / "d3_c6.npy")
return _REF_C6
[docs]
def load_R0():
"""
Load cut-off radii values for Grimme-D3. Values from:
`<https://doi.org/10.1063/1.3382344>`_
Pairs of chemical elements have their own cut-off radii, corresponding to the
two indices of the returned array. Units are Ångströms.
Returns:
ref_R0: numpy.ndarray of shape (94, 94). Cut-off radii.
"""
global _R0_AB
if _R0_AB is None:
_R0_AB = np.load(Path(__file__).parent / "d3_R0.npy")
return _R0_AB
# fmt: off
DF_DEFAULT_PARAMS = {
'PBE' : {'s6': 1.000, 's8': 0.7875, 'a1': 0.4289, 'a2': 4.4407},
'B1B95' : {'s6': 1.000, 's8': 1.4507, 'a1': 0.2092, 'a2': 5.5545},
'B2GPPLYP': {'s6': 0.560, 's8': 0.2597, 'a1': 0.0000, 'a2': 6.3332},
'B3PW91' : {'s6': 1.000, 's8': 2.8524, 'a1': 0.4312, 'a2': 4.4693},
'BHLYP' : {'s6': 1.000, 's8': 1.0354, 'a1': 0.2793, 'a2': 4.9615},
'BMK' : {'s6': 1.000, 's8': 2.0860, 'a1': 0.1940, 'a2': 5.9197},
'BOP' : {'s6': 1.000, 's8': 3.2950, 'a1': 0.4870, 'a2': 0.5043},
'BPBE' : {'s6': 1.000, 's8': 4.0728, 'a1': 0.4567, 'a2': 4.3908},
'CAMB3LYP': {'s6': 1.000, 's8': 2.0674, 'a1': 0.3708, 'a2': 5.4743},
'LCwPBE' : {'s6': 1.000, 's8': 1.8541, 'a1': 0.3919, 'a2': 5.0897},
'MPW1B95' : {'s6': 1.000, 's8': 1.0508, 'a1': 0.1955, 'a2': 6.4177},
'MPWB1K' : {'s6': 1.000, 's8': 0.9499, 'a1': 0.1474, 'a2': 6.6223},
'mPWLYP' : {'s6': 1.000, 's8': 2.0077, 'a1': 0.4831, 'a2': 4.5323},
'OLYP' : {'s6': 1.000, 's8': 2.6205, 'a1': 0.5299, 'a2': 2.8065},
'OPBE' : {'s6': 1.000, 's8': 3.3816, 'a1': 0.5512, 'a2': 2.9444},
'oTPSS' : {'s6': 1.000, 's8': 2.7495, 'a1': 0.4634, 'a2': 4.3153},
'PBE38' : {'s6': 1.000, 's8': 1.4623, 'a1': 0.3995, 'a2': 5.1405},
'PBEsol' : {'s6': 1.000, 's8': 2.9491, 'a1': 0.4466, 'a2': 6.1742},
'PTPSS' : {'s6': 0.750, 's8': 0.2804, 'a1': 0.0000, 'a2': 6.5745},
'PWB6K' : {'s6': 1.000, 's8': 0.9383, 'a1': 0.1805, 'a2': 7.7627},
'revSSB' : {'s6': 1.000, 's8': 0.4389, 'a1': 0.4720, 'a2': 4.0986},
'SSB' : {'s6': 1.000, 's8': 0.1744, 'a1': -0.0952, 'a2': 5.2170},
'TPSSh' : {'s6': 1.000, 's8': 2.2382, 'a1': 0.4529, 'a2': 4.6550},
'HCTH120' : {'s6': 1.000, 's8': 1.0821, 'a1': 0.3563, 'a2': 4.3359},
'B2PLYP' : {'s6': 0.640, 's8': 0.9147, 'a1': 0.3065, 'a2': 5.0570},
'B3LYP' : {'s6': 1.000, 's8': 1.9889, 'a1': 0.3981, 'a2': 4.4211},
'B97D' : {'s6': 1.000, 's8': 2.2609, 'a1': 0.5545, 'a2': 3.2297},
'BLYP' : {'s6': 1.000, 's8': 2.6996, 'a1': 0.4298, 'a2': 4.2359},
'BP86' : {'s6': 1.000, 's8': 3.2822, 'a1': 0.3946, 'a2': 4.8516},
'DSDBLYP' : {'s6': 0.500, 's8': 0.2130, 'a1': 0.0000, 'a2': 6.0519},
'PBE0' : {'s6': 1.000, 's8': 1.2177, 'a1': 0.4145, 'a2': 4.8593},
'PBE' : {'s6': 1.000, 's8': 0.7875, 'a1': 0.4289, 'a2': 4.4407},
'PW6B95' : {'s6': 1.000, 's8': 0.7257, 'a1': 0.2076, 'a2': 6.3750},
'PWPB95' : {'s6': 0.820, 's8': 0.2904, 'a1': 0.0000, 'a2': 7.3141},
'revPBE0' : {'s6': 1.000, 's8': 1.7588, 'a1': 0.4679, 'a2': 3.7619},
'revPBE38': {'s6': 1.000, 's8': 1.4760, 'a1': 0.4309, 'a2': 3.9446},
'revPBE' : {'s6': 1.000, 's8': 2.3550, 'a1': 0.5238, 'a2': 3.5016},
'rPW86PBE': {'s6': 1.000, 's8': 1.3845, 'a1': 0.4613, 'a2': 4.5062},
'TPSS0' : {'s6': 1.000, 's8': 1.2576, 'a1': 0.3768, 'a2': 4.5865},
'TPSS' : {'s6': 1.000, 's8': 1.9435, 'a1': 0.4535, 'a2': 4.4752}
} # fmt: on
"""
Default Grimme-D3 scaling parameters for a variety of density functionals using Becke-Johnson damping.
Values taken from `<https://www.chemie.uni-bonn.de/grimme/de/software/dft-d3/bj_damping>`_.
"""
[docs]
def get_df_params(params):
if isinstance(params, str):
if params not in DF_DEFAULT_PARAMS:
raise ValueError(f"Default parameters not available for functional `{params}`. See the available default parameters in ppafm.defaults.d3.DF_DEFAULT_PARAMS.")
params = DF_DEFAULT_PARAMS[params]
for key in ["s6", "s8", "a1", "a2"]:
if key not in params:
raise ValueError(f"DFT-D3 parameter dictionary is missing entry for `{key}`")
return params